diff --git a/src/pumipic_adjacency.tpp b/src/pumipic_adjacency.tpp index 86393847..79384621 100644 --- a/src/pumipic_adjacency.tpp +++ b/src/pumipic_adjacency.tpp @@ -15,6 +15,7 @@ #include "pumipic_constants.hpp" #include "pumipic_kktypes.hpp" #include "pumipic_profiling.hpp" +#include "pumipic_adjacency.hpp" namespace o = Omega_h; namespace ps = particle_structs; @@ -325,8 +326,8 @@ namespace pumipic { auto xpts = o::zero_vector<3>(); const o::LO prevExit = lastExit[ptcl]; lastExit[ptcl] = -1; - o::Real quality = -1; - o::LO bestFace = -1; + //o::Real quality = -1; + //o::LO bestFace = -1; for(int fi=0; fi<4; ++fi) { const auto face_id = face_ids[fi]; if (face_id == prevExit) @@ -337,24 +338,30 @@ namespace pumipic { o::Real dproj; o::Real closeness; o::Real intersection_parametric_coord; - const bool success = ray_intersects_triangle(face, orig, dest, xpts, tol, flip, dproj, + const bool success = line_segment_intersects_triangle(face, orig, dest, xpts, tol, flip, dproj, closeness, intersection_parametric_coord); if (success) { lastExit[ptcl] = face_id; xPoints[3*ptcl] = xpts[0]; xPoints[3*ptcl + 1] = xpts[1]; xPoints[3*ptcl + 2] = xpts[2]; } - if (dproj > -tol && (quality < 0 || closeness < quality) && lastExit[ptcl] == -1) { - quality = closeness; - bestFace = face_id; - xPoints[3*ptcl] = xpts[0]; xPoints[3*ptcl + 1] = xpts[1]; xPoints[3*ptcl + 2] = xpts[2]; - } + //if (dproj > -tol && (quality < 0 || closeness < quality) && lastExit[ptcl] == -1) { + // quality = closeness; + // bestFace = face_id; + // xPoints[3*ptcl] = xpts[0]; xPoints[3*ptcl + 1] = xpts[1]; xPoints[3*ptcl + 2] = xpts[2]; + //} } //If line intersection fails then use BCC method - if (lastExit[ptcl] == -1) { - lastExit[ptcl] = bestFace; + //if (lastExit[ptcl] == -1) { + // printf("!!!!!!! Best face %d\n", bestFace); + // lastExit[ptcl] = bestFace; + //} + if (lastExit[ptcl == -1]){ // reached destination + xPoints[3*ptcl] = x_ps_tgt(ptcl,0); + xPoints[3*ptcl+1] = x_ps_tgt(ptcl,1); + xPoints[3*ptcl+2] = x_ps_tgt(ptcl,2); } - ptcl_done[ptcl] = (lastExit[ptcl] == -1); + //ptcl_done[ptcl] = (lastExit[ptcl] == -1); } }; parallel_for(ptcls, findExitFace, "search_findExitFace_intersect_3d"); @@ -369,47 +376,54 @@ namespace pumipic { o::Write lastExit, o::Bytes side_is_exposed, bool requireIntersection, o::Write xFace) { - auto checkExposedEdges = PS_LAMBDA(const int e, const int pid, const int mask) { - if( mask > 0 && !ptcl_done[pid] ) { - assert(lastExit[pid] != -1); - const o::LO bridge = lastExit[pid]; - const bool exposed = side_is_exposed[bridge]; - ptcl_done[pid] = exposed; - if (exposed && requireIntersection) { - xFace[pid] = lastExit[pid]; - } - else { - elem_ids[pid] = exposed ? -1 : elem_ids[pid]; //leaves domain if exposed + auto checkIntersectionwGeomBoundary = PS_LAMBDA(const int e, const int pid, const int mask) { + if (mask > 0){ + if (lastExit[pid] == -1){ + xFace[pid] = -1; + ptcl_done[pid] = 1; + } else { + const o::LO bridge = lastExit[pid]; + const bool exposed = side_is_exposed[bridge]; + ptcl_done[pid] = exposed; + xFace[pid] = exposed ? lastExit[pid] : -1; + elem_ids[pid] = exposed ? -1 : elem_ids[pid]; //leaves domain if exposed + } } - } }; - parallel_for(ptcls, checkExposedEdges, "pumipic_checkExposedEdges"); + parallel_for(ptcls, checkIntersectionwGeomBoundary, "pumipic_checkIntersectionwGeomBoundary"); } template - void set_new_element(o::Mesh& mesh, ParticleStructure* ptcls, - o::Write elem_ids, o::LOs ptcl_done, - o::Write lastExit) { + void find_next_element(o::Mesh& mesh, ParticleStructure* ptcls, + o::Write elem_ids, o::Write next_element, o::LOs ptcl_done, + o::Write lastExit) { int dim = mesh.dim(); const auto faces2elms = mesh.ask_up(dim-1, dim); auto e2f_vals = faces2elms.ab2b; // CSR value array auto e2f_offsets = faces2elms.a2ab; // CSR offset array, index by mesh edge ids auto setNextElm = PS_LAMBDA(const int& e, const int& pid, const int& mask) { if( mask > 0 && !ptcl_done[pid] ) { - auto searchElm = elem_ids[pid]; auto bridge = lastExit[pid]; - auto e2f_first = e2f_offsets[bridge]; - #ifdef _DEBUG - auto e2f_last = e2f_offsets[bridge+1]; - auto upFaces = e2f_last - e2f_first; - assert(upFaces==2); - #endif - auto faceA = e2f_vals[e2f_first]; - auto faceB = e2f_vals[e2f_first+1]; - assert(faceA != faceB); - assert(faceA == searchElm || faceB == searchElm); - auto nextElm = (faceA == searchElm) ? faceB : faceA; - elem_ids[pid] = nextElm; + if (bridge != -1) { + auto searchElm = elem_ids[pid]; + auto e2f_first = e2f_offsets[bridge]; + auto e2f_last = e2f_offsets[bridge+1]; + auto upFaces = e2f_last - e2f_first; + if (upFaces == 1) { // particle at Geometry boundary + next_element[pid] = -1; + } + else if (upFaces==2) { // particle at element boundary to enter next element + assert(upFaces == 2); + auto faceA = e2f_vals[e2f_first]; + auto faceB = e2f_vals[e2f_first + 1]; + assert(faceA != faceB); + assert(faceA == searchElm || faceB == searchElm); + auto nextElm = (faceA == searchElm) ? faceB : faceA; + next_element[pid] = nextElm; + } else { + OMEGA_H_CHECK_PRINTF(false, "Found invalid number of adjacent elements %d\n", upFaces); + } + } } }; parallel_for(ptcls, setNextElm, "pumipic_setNextElm"); @@ -449,31 +463,38 @@ namespace pumipic { * @param x_ps_tgt Particle target position * @param pids Particle ids * @param elem_ids Paricle parent element ids + * @param next_element Next element over the exit face * @param requireIntersection True if intersection is required * @param inter_faces Exit faces for particles at domain boundary * @param inter_points Stores intersection points for particles at each face * @param looplimit Maximum number of iterations * @param debug True if debug information is printed * @param func Callable object to handle particles at element sides or destination + * @param elmArea Element areas/volumes based on the mesh dimension + * @param given_tol Tolerance for intersection. If not positive, it is computed from the minimum element area * @return True if all particles are found at destination or left domain */ template bool trace_particle_through_mesh(o::Mesh& mesh, ParticleStructure* ptcls, Segment3d x_ps_orig, Segment3d x_ps_tgt, SegmentInt pids, o::Write& elem_ids, + o::Write& next_element, bool requireIntersection, o::Write& inter_faces, o::Write& inter_points, + o::Write& lastExit, int looplimit, bool debug, - Func& func) { + Func& func, + o::Reals& elmArea, + std::optional given_tol = std::nullopt) { static_assert( std::is_invocable_r_v< void, Func, o::Mesh &, ParticleStructure *, - o::Write &, o::Write &, o::Write &, + o::Write &, o::Write &, o::Write &, o::Write &, o::Write &, o::Write &, Segment3d, Segment3d>, - "Functional must accept \n"); + "Functional must accept \n"); //Initialize timer const auto btime = pumipic_prebarrier(mesh.comm()->get_impl()); @@ -483,12 +504,22 @@ namespace pumipic { //Initial setup const auto psCapacity = ptcls->capacity(); // True if particle has reached new parent element + // TODO: ptcl_done is searches internal state, it will be handled when it + // is made a class method o::Write ptcl_done(psCapacity, 0, "search_ptcl_done"); // Store the last exit face - o::Write lastExit(psCapacity,-1, "search_last_exit"); - const auto elmArea = measure_elements_real(&mesh); + //o::Write lastExit(psCapacity,-1, "search_last_exit"); + if (lastExit.size() == 0){ + lastExit = Omega_h::Write(psCapacity, -1, "search_last_exit"); + } bool useBcc = !requireIntersection; - o::Real tol = compute_tolerance_from_area(elmArea); + + o::Real tol = 0; + if (!given_tol.has_value()) { + tol = compute_tolerance_from_area(elmArea); + } else { + tol = given_tol.value(); + } int rank; MPI_Comm_rank(mesh.comm()->get_impl(), &rank); @@ -521,6 +552,11 @@ namespace pumipic { parallel_for(ptcls, setInitial, "search_setInitial"); } + // if next_element is not initialized, initialize it with elem_ids + if (next_element.size() == 0) { + next_element = o::deep_copy(o::LOs(elem_ids), "next element"); + } + //Finish particles that didn't move auto finishUnmoved = PS_LAMBDA(const int e, const int pid, const int mask) { if (mask){ @@ -558,11 +594,11 @@ namespace pumipic { while (!found) { //Find intersection face find_exit_face(mesh, ptcls, x_ps_orig, x_ps_tgt, elem_ids, ptcl_done, elmArea, useBcc, lastExit, inter_points, tol); + // find next element + find_next_element(mesh, ptcls, elem_ids, next_element, ptcl_done, lastExit); //Check if intersection face is exposed - func(mesh, ptcls, elem_ids, inter_faces, lastExit, inter_points, ptcl_done, x_ps_orig, x_ps_tgt); + func(mesh, ptcls, elem_ids, next_element, inter_faces, lastExit, inter_points, ptcl_done, x_ps_orig, x_ps_tgt); - // Move to next element - set_new_element(mesh, ptcls, elem_ids, ptcl_done, lastExit); //Check if all particles are found found = true; @@ -594,14 +630,15 @@ namespace pumipic { rank, searchElm, ptcl, ptclOrigin[0], ptclOrigin[1], ptclDest[0], ptclDest[1]); } - elem_ids[pid] = -1; + // TODO it shouldn't mark them -1, it should throw error or behavior set with parameter + //elem_ids[pid] = -1; Kokkos::atomic_add(&(numNotFound[0]), 1); } }; ps::parallel_for(ptcls, ptclsNotFound, "ptclsNotFound"); Omega_h::HostWrite numNotFound_h(numNotFound); printError( "ERROR:Rank %d: loop limit %d exceeded. %d particles were " - "not found. Deleting them...\n", rank, looplimit, numNotFound_h[0]); + "not found.\n", rank, looplimit, numNotFound_h[0]); break; } @@ -622,7 +659,7 @@ namespace pumipic { } void operator()(o::Mesh& mesh, ParticleStructure* ptcls, - o::Write& elem_ids, o::Write& inter_faces, + o::Write& elem_ids, o::Write& next_element, o::Write& inter_faces, o::Write& lastExit, o::Write& inter_points, o::Write& ptcl_done, Segment3d x_ps_orig, @@ -631,6 +668,11 @@ namespace pumipic { check_model_intersection(mesh, ptcls, x_ps_orig, x_ps_tgt, elem_ids, ptcl_done, lastExit, side_is_exposed_, requireIntersection_, inter_faces); + // make next element the current element + auto copy_next_element = OMEGA_H_LAMBDA(o::LO pid) { + elem_ids[pid] = next_element[pid]; + }; + o::parallel_for(elem_ids.size(), copy_next_element, "copy_next_element"); } private: @@ -648,9 +690,11 @@ namespace pumipic { int looplimit, int debug) { RemoveParticleOnGeometricModelExit native_handler(mesh, requireIntersection); - - return trace_particle_through_mesh(mesh, ptcls, x_ps_orig, x_ps_tgt, pids, elem_ids, requireIntersection, - inter_faces, inter_points, looplimit, debug, native_handler); + o::Reals elmArea = measure_elements_real(&mesh); + o::Write lastExit (ptcls->capacity(), -1, "lastExit"); + o::Write next_element(0, "next_element"); + return trace_particle_through_mesh(mesh, ptcls, x_ps_orig, x_ps_tgt, pids, elem_ids, next_element, requireIntersection, + inter_faces, inter_points, lastExit, looplimit, debug, native_handler, elmArea); } } #endif diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5bb484c2..8a98cb38 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -3,6 +3,8 @@ function(make_test exename srcname) target_link_libraries(${exename} pumipic Omega_h::omega_h) endfunction(make_test) +make_test(test_trace_particle_through_mesh test_trace_particle_through_mesh.cpp) +make_test(test_search_mesh test_search_mesh.cpp) make_test(print_partition print_partition.cpp) make_test(print_classification print_classification.cpp) make_test(full_mesh test_full_mesh.cpp) diff --git a/test/test_search_mesh.cpp b/test/test_search_mesh.cpp new file mode 100644 index 00000000..f1217b66 --- /dev/null +++ b/test/test_search_mesh.cpp @@ -0,0 +1,308 @@ +// +// Created by Fuad Hasan on 2/25/25. +// + +#include "pumipic_adjacency.hpp" +#include "pumipic_adjacency.tpp" +#include "pumipic_kktypes.hpp" +#include "pumipic_library.hpp" +#include "pumipic_mesh.hpp" +#include "pumipic_utils.hpp" +#include +#include +#include +#include +#include +#include +#include +#include + +using particle_structs::lid_t; +using particle_structs::MemberTypes; +using particle_structs::SellCSigma; +using pumipic::fp_t; +using pumipic::Vector3d; + +namespace o = Omega_h; +namespace p = pumipic; +namespace ps = particle_structs; + +typedef MemberTypes Particle; +typedef ps::ParticleStructure PS; +typedef Kokkos::DefaultExecutionSpace ExeSpace; + +void print_test_info(); +void printf_face_info(o::Mesh &mesh, o::LOs faceIds, bool all); +bool is_inside3D(o::Mesh &mesh, o::LO elem_id, const o::Vector<3> point); + +bool is_close(const double a, const double b, double tol = 1e-8) { + return std::abs(a - b) < tol; +} + +OMEGA_H_INLINE bool is_close_d(const double a, const double b, double tol = 1e-8) { + return Kokkos::abs(a - b) < tol; +} + + +int main(int argc, char* argv[]){ + // ----------------------------------------------------------------------------------// + // ----------------------------------- Set up ---------------------------------------// + // ----------------------------------------------------------------------------------// + + print_test_info(); + + printf("\n============================ Setting up Mesh and Particles ======================\n"); + p::Library lib(&argc, &argv); + auto &olib = lib.omega_h_lib(); + auto world = olib.world(); + // simplest 3D mesh + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 1, 1, 1, false); + printf("[INFO] Mesh created with %d vertices and %d faces\n", + mesh.nverts(), mesh.nfaces()); + + printf_face_info(mesh, {}, true); + + Omega_h::Write owners(mesh.nelems(), 0); + p::Mesh picparts(mesh, owners); + o::Mesh *p_mesh = picparts.mesh(); + // create particle structure with 5 particles + int num_ptcls = 5; + + // create particles + Kokkos::TeamPolicy policy; +#ifdef PP_USE_GPU + printf("Using GPU for simulation..."); + policy = + Kokkos::TeamPolicy(10, Kokkos::AUTO()); +#else + printf("Using CPU for simulation..."); + policy = Kokkos::TeamPolicy(100, 1); +#endif + const o::LO tet_id = 0; + o::Int ne = p_mesh->nelems(); + PS::kkLidView ptcls_per_elem("ptcls_per_elem", ne); + PS::kkGidView element_gids("element_gids", ne); + o::GOs mesh_element_gids = picparts.globalIds(picparts.dim()); + o::parallel_for( + ne, OMEGA_H_LAMBDA( + const int &i) { + element_gids(i) = mesh_element_gids[i]; + if (i == tet_id) { + ptcls_per_elem(i) = num_ptcls; + } else { + ptcls_per_elem(i) = 0; + } + }); + +#ifdef PP_ENABLE_CAB + PS *ptcls = new p::DPS(policy, ne, num_ptcls, ptcls_per_elem, element_gids); + printf("DPS Particle structure created successfully\n"); +#else + PS *ptcls = PS *ptcls = new SellCSigma( + policy, 10, 10, ne, 2, ptcls_per_elem, element_gids); + printf("SellCSigma Particle structure created successfully\n"); +#endif + + printf("\n====================== Mesh and Particles are set up ============================\n"); + + // ----------------------------------------------------------------------------------// + // ------------------------------- Particle Movement ---------------------------------// + // ----------------------------------------------------------------------------------// + + printf("\n\n\n============================ Particle Movement Setup =================================\n"); + + // The following are verified using: https://gist.github.com/Fuad-HH/5e0aed99f271617e283e9108091fb1cb + o::Vector<3> particle_initial_loc {0.5, 0.75, 0.25}; + o::Vector<3> p0_dest = particle_initial_loc + o::Vector<3>({0.0001, 0.0001, 0.0001}); // moved a little bit + o::Vector<3> p1_dest = {0.8, 0.9, 0.1}; // inside the same element + o::Vector<3> p2_dest = {0.8, 0.1, 0.95}; // moved to element 3 + // p2 intersections : 1. [0.57894737 0.57894737 0.43421053] + // 2. [0.61111111 0.50925926 0.50925926] + // 3. [0.6875 0.34375 0.6875 ] + o::Vector<3> p3_dest = {-0.5, 0.75, 0.25}; // moved out of the mesh through element 1 + // p3 intersections : 1. [0.25 0.75 0.25] + // 2. [0. 0.75 0.25] + o::Vector<3> p4_dest = {0.5, 0.75, -1.5}; // moved out of the mesh through face 0 of element 0 + // p4 intersections : 1. [0.5 0.75 0] + + + // Set particle positions + auto particle_orig = ptcls->get<0>(); + auto particle_dest = ptcls->get<1>(); + auto particle_id = ptcls->get<2>(); + + auto set_particle_path = PS_LAMBDA(const int& eid, const int& pid, const int& mask) { + if (mask > 0) { + particle_id(pid) = pid; + + particle_orig(pid, 0) = particle_initial_loc[0]; + particle_orig(pid, 1) = particle_initial_loc[1]; + particle_orig(pid, 2) = particle_initial_loc[2]; + + if (pid == 0) { + particle_dest(pid, 0) = p0_dest[0]; + particle_dest(pid, 1) = p0_dest[1]; + particle_dest(pid, 2) = p0_dest[2]; + } else if (pid == 1) { + particle_dest(pid, 0) = p1_dest[0]; + particle_dest(pid, 1) = p1_dest[1]; + particle_dest(pid, 2) = p1_dest[2]; + } else if (pid == 2) { + particle_dest(pid, 0) = p2_dest[0]; + particle_dest(pid, 1) = p2_dest[1]; + particle_dest(pid, 2) = p2_dest[2]; + } else if (pid == 3) { + particle_dest(pid, 0) = p3_dest[0]; + particle_dest(pid, 1) = p3_dest[1]; + particle_dest(pid, 2) = p3_dest[2]; + } else if (pid == 4) { + particle_dest(pid, 0) = p4_dest[0]; + particle_dest(pid, 1) = p4_dest[1]; + particle_dest(pid, 2) = p4_dest[2]; + } else { + printf("[ERROR] Particle id %d is not supported\n", pid); + } + + printf("Particle id %d starts at (%f, %f, %f) and moves to (%f, %f, %f)\n", + pid, particle_orig(pid, 0), particle_orig(pid, 1), particle_orig(pid, 2), + particle_dest(pid, 0), particle_dest(pid, 1), particle_dest(pid, 2)); + } + }; + pumipic::parallel_for(ptcls, set_particle_path, "set_particle_path"); + + // -------------------------------------- Particle Movement Set up Done --------------------------------------// + auto elem_ids = Omega_h::Write (0, "elem ids"); + bool requireIntersection = true; + auto interFaces = Omega_h::Write (0, "inter faces"); + auto interPoints = Omega_h::Write (0, "inter points"); + + bool success = pumipic::search_mesh(mesh, ptcls, particle_orig, particle_dest, particle_id, elem_ids, requireIntersection, interFaces, interPoints, 100, 1); + if (!success) { + printf("[ERROR] search_mesh failed\n"); + throw std::runtime_error("search_mesh failed"); // FIXME: where to delete ptcls? + } + + + // ----------------------------------- Check Results -----------------------------------// + // check sizes + OMEGA_H_CHECK_PRINTF(elem_ids.size() == ptcls->capacity(), "elem_ids size %d != ptcls capacity %d\n", + elem_ids.size(), ptcls->capacity()); + OMEGA_H_CHECK_PRINTF(interFaces.size() == ptcls->capacity(), "interFaces size %d != ptcls capacity %d\n", + interFaces.size(), ptcls->capacity()); + OMEGA_H_CHECK_PRINTF(interPoints.size() == 3*ptcls->capacity(), "interPoints size %d != 3 x ptcls capacity %d\n", + interPoints.size(), 3*ptcls->capacity()); + + Omega_h::Few expected_elem_ids = {0, 0, 3, -1, -1}; + Omega_h::Few expected_interFaces = {-1, -1, -1, 1, 0}; + Omega_h::Few, 5> expected_interPoints = {p0_dest, p1_dest, p2_dest, p3_dest, p4_dest}; + + printf("\n\n\n============================ Checking Search Results =================================\n"); + auto check_search_results = PS_LAMBDA(const auto e, const auto pid, const auto mask) { + if (mask > 0) { + printf("\tPid %d: elem_id (%2d, %2d), interFace (%2d, %2d), interPoint ([% .6f, % .6f, % .6f],[% .6f, % .6f, % .6f])\n", + pid, elem_ids[pid], expected_elem_ids[pid], interFaces[pid], expected_interFaces[pid], + interPoints[3*pid], interPoints[3*pid+1], interPoints[3*pid+2], + expected_interPoints[pid][0], expected_interPoints[pid][1], expected_interPoints[pid][2]); + + OMEGA_H_CHECK_PRINTF(elem_ids[pid] == expected_elem_ids[pid], "Particle %d: Expected elem_id %d != found elem_id %d\n", + pid, expected_elem_ids[pid], elem_ids[pid]); + OMEGA_H_CHECK_PRINTF(interFaces[pid] == expected_interFaces[pid], "Particle %d: Expected interFace %d != found interFace %d\n", + pid, expected_interFaces[pid], interFaces[pid]); + OMEGA_H_CHECK_PRINTF(is_close_d(interPoints[3*pid], expected_interPoints[pid][0]), + "Particle %d: Expected interPoint x %f != found interPoint x %f\n", + pid, expected_interPoints[pid][0], interPoints[3*pid]); + OMEGA_H_CHECK_PRINTF(is_close_d(interPoints[3*pid+1], expected_interPoints[pid][1]), + "Particle %d: Expected interPoint y %f != found interPoint y %f\n", + pid, expected_interPoints[pid][1], interPoints[3*pid+1]); + OMEGA_H_CHECK_PRINTF(is_close_d(interPoints[3*pid+2], expected_interPoints[pid][2]), + "Particle %d: Expected interPoint z %f != found interPoint z %f\n", + pid, expected_interPoints[pid][2], interPoints[3*pid+2]); + } + }; + pumipic::parallel_for(ptcls, check_search_results, "check_search_results"); + printf("\n============================ Search Results are Correct =================================\n"); + + + + // ------------------------------------------- Clean up -------------------------------------------// + delete ptcls; + + return 0; +} + +void print_test_info() { + printf("\n\n!!!!!!!!!!!!!!!!!!!!!!!!!!!! Important !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" + "!!!!!!!!!!!!!!! This test only works when compiled in DEBUG mode !!!!!!!!!!!!!!!!\n\n"); + + std::string message = "\n========================== Test Search Mesh Function ======================\n"; + message += "Testing 'search_mesh' function form pumipic_adjacency.tpp ...\n"; + message += "A 3D internal mesh is created and particles will be moved to different locations\n"; + message += "or will be kept in the same element or same position.\n"; + message += "Particle movements are: \n"; + message += "\t1. Inside the same element at same position\n"; + message += "\t2. Inside the same element but at different position\n"; + message += "\t3. Different element's interior\n"; + message += "\t4. Moves out of the mesh traversing multiple elements\n"; + message += "\t5. Moves to the boundary from the same element it started\n"; + message += "\n===============================================================\n"; + printf("%s", message.c_str()); +} + +bool is_inside3D(o::Mesh &mesh, o::LO elem_id, const o::Vector<3> point) { + OMEGA_H_CHECK_PRINTF(mesh.dim() == 3, "Mesh is not 3D. Found dimension %d\n", + mesh.dim()); + const auto &coords = mesh.coords(); + const auto &tet2nodes = mesh.ask_verts_of(o::REGION); + + o::Write inside(1, 0); + + auto is_inside_lambda = OMEGA_H_LAMBDA(o::LO id) { + const auto current_el_verts = o::gather_verts<4>(tet2nodes, elem_id); + const Omega_h::Few, 4> current_el_vert_coords = + o::gather_vectors<4, 3>(coords, current_el_verts); + o::Vector<4> bcc = + o::barycentric_from_global<3, 3>(point, current_el_vert_coords); + inside[0] = p::all_positive(bcc, 0.0); + }; + o::parallel_for(1, is_inside_lambda); + auto host_inside = o::HostWrite(inside); + + return bool(host_inside[0]); +} + +void printf_face_info(o::Mesh &mesh, o::LOs faceIds, bool all) { + const auto exposed_faces = o::mark_exposed_sides(&mesh); + const auto &face2nodes = mesh.ask_down(o::FACE, o::VERT).ab2b; + const auto &face2cellcell = mesh.ask_up(o::FACE, o::REGION).ab2b; + const auto &face2celloffset = mesh.ask_up(o::FACE, o::REGION).a2ab; + const auto &coords = mesh.coords(); + + auto print_faces = OMEGA_H_LAMBDA(o::LO + faceid) { + if (!all) { + for (int i = 0; i < faceIds.size(); i++) { + o::LO id = faceIds[i]; + if (id == faceid) { + printf("Face %d nodes %d %d %d Exposed %d\n", faceid, + face2nodes[faceid * 3], face2nodes[faceid * 3 + 1], + face2nodes[faceid * 3 + 2], exposed_faces[faceid]); + } + } + } else if (all) { + printf("Face %d nodes %d %d %d Exposed %d\n", faceid, + face2nodes[faceid * 3], face2nodes[faceid * 3 + 1], + face2nodes[faceid * 3 + 2], exposed_faces[faceid]); + int n_adj_cells = face2celloffset[faceid+1] - face2celloffset[faceid]; + if(n_adj_cells == 1){ + printf("Face %d is exposed and with %d cell\n", faceid, face2cellcell[face2celloffset[faceid]]); + } else if (n_adj_cells ==2) { + printf("Face %d internal with %d %d\n", faceid, face2cellcell[face2celloffset[faceid]], + face2cellcell[face2celloffset[faceid] + 1]); + } else { + printf("[ERROR] Face cannot be adjacent to more than 2 cells\n found %d\n", n_adj_cells); + } + } + }; + o::parallel_for(mesh.nfaces(), print_faces, "print asked faces"); +} diff --git a/test/test_trace_particle_through_mesh.cpp b/test/test_trace_particle_through_mesh.cpp new file mode 100644 index 00000000..5d2c34c6 --- /dev/null +++ b/test/test_trace_particle_through_mesh.cpp @@ -0,0 +1,433 @@ +// +// Created by Fuad Hasan on 2/12/25. +// + +#include "pumipic_adjacency.hpp" +#include "pumipic_adjacency.tpp" +#include "pumipic_kktypes.hpp" +#include "pumipic_library.hpp" +#include "pumipic_mesh.hpp" +#include "pumipic_utils.hpp" +#include +#include +#include +#include +#include +#include +#include +#include + +using particle_structs::lid_t; +using particle_structs::MemberTypes; +using particle_structs::SellCSigma; +using pumipic::fp_t; +using pumipic::Vector3d; + +namespace o = Omega_h; +namespace p = pumipic; +namespace ps = particle_structs; + +typedef MemberTypes Particle; +typedef ps::ParticleStructure PS; +typedef Kokkos::DefaultExecutionSpace ExeSpace; + +void printf_face_info(o::Mesh &mesh, o::LOs faceIds, bool all = false); + +void print_test_info(); + +bool is_close(const double a, const double b, double tol = 1e-8) { + return std::abs(a - b) < tol; +} + +OMEGA_H_INLINE bool is_close_d(const double a, const double b, double tol = 1e-8) { + return Kokkos::abs(a - b) < tol; +} + +template +struct empty_function { + void operator()( + o::Mesh &mesh, ps::ParticleStructure *ptcls, + o::Write &elem_ids, o::Write& next_elements, o::Write &inter_faces, + o::Write &lastExit, o::Write &inter_points, + o::Write &ptcl_done, + Segment3d x_ps_orig, + Segment3d x_ps_tgt) const { + printf("Empty Function Called\n"); + } +}; + + +int main(int argc, char **argv) { + printf("\n\n!!!!!!!!!!!!!!!!!!!!!!!!!!!! Important !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" + "!!!!!!!!!!!!!!! This test only works when compiled in DEBUG mode !!!!!!!!!!!!!!!!\n\n"); + // **************************************** Setup ***************************************************************// + p::Library lib(&argc, &argv); + print_test_info(); + auto &olib = lib.omega_h_lib(); + auto world = olib.world(); + // simplest 3D mesh + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 1, 1, 1, false); + printf("[INFO] Mesh created with %d vertices and %d faces\n", + mesh.nverts(), mesh.nfaces()); + + printf_face_info(mesh, {}, true); + + Omega_h::Write owners(mesh.nelems(), 0); + p::Mesh picparts(mesh, owners); + o::Mesh *p_mesh = picparts.mesh(); + // create particle structure with 5 particles + int num_ptcls = 5; + + // create particles + Kokkos::TeamPolicy policy; +#ifdef PP_USE_GPU + printf("Using GPU for simulation..."); + policy = + Kokkos::TeamPolicy(10, Kokkos::AUTO()); +#else + printf("Using CPU for simulation..."); + policy = Kokkos::TeamPolicy(100, 1); +#endif + const o::LO tet_id = 0; + o::Int ne = p_mesh->nelems(); + PS::kkLidView ptcls_per_elem("ptcls_per_elem", ne); + PS::kkGidView element_gids("element_gids", ne); + o::GOs mesh_element_gids = picparts.globalIds(picparts.dim()); + o::parallel_for( + ne, OMEGA_H_LAMBDA( + const int &i) { + element_gids(i) = mesh_element_gids[i]; + if (i == tet_id) { + ptcls_per_elem(i) = num_ptcls; + } else { + ptcls_per_elem(i) = 0; + } + }); + +#ifdef PP_ENABLE_CAB + PS *ptcls = new p::DPS(policy, ne, num_ptcls, ptcls_per_elem, element_gids); + printf("DPS Particle structure created successfully\n"); +#else + PS *ptcls = PS *ptcls = new SellCSigma( + policy, 10, 10, ne, 2, ptcls_per_elem, element_gids); + printf("SellCSigma Particle structure created successfully\n"); +#endif + + printf("\n\n\n==============================================================================================\n"); + printf(">============ Checking when particles move from cell 0 to 5 =================================<\n"); + printf("==============================================================================================\n"); + + // set particle position + Omega_h::Vector<3> cell0_centroid{0.5, 0.75, 0.25}; + Omega_h::Vector<3> cell5_centroid{0.75, 0.5, 0.25}; + auto particle_init_position = ptcls->get<0>(); + auto particle_final_position = ptcls->get<1>(); + auto pid_d = ptcls->get<2>(); + auto setIDs = PS_LAMBDA( + const int &eid, + const int &pid, + const bool &mask) { + if (mask > 0) { + particle_init_position(pid, 0) = cell0_centroid[0]; + particle_init_position(pid, 1) = cell0_centroid[1]; + particle_init_position(pid, 2) = cell0_centroid[2]; + + particle_final_position(pid, 0) = cell5_centroid[0]; + particle_final_position(pid, 1) = cell5_centroid[1]; + particle_final_position(pid, 2) = cell5_centroid[2]; + + pid_d(pid) = pid; + + + printf("Initialized particle %d origin (%f, %f, %f) and destination (%f, %f, %f)\n", + pid, + particle_init_position(pid, 0), particle_init_position(pid, 1), particle_init_position(pid, 2), + particle_final_position(pid, 0), particle_final_position(pid, 1), particle_final_position(pid, 2)); + } + }; + ps::parallel_for(ptcls, setIDs); + // ********************************************** Particle Setup Done ******************************************// + + // **************************************** Set up Auxiliary Search Arrays *************************************// + const bool requireIntersection = true; + + Omega_h::Write elem_ids(0, "element ids"); + Omega_h::Write inter_faces(0, "inter faces"); + Omega_h::Write last_exit(0, "last exit"); + Omega_h::Write inter_points(0, "inter points"); + Omega_h::Write next_elements(0, "next elements"); + + o::Reals elmArea = measure_elements_real(&mesh); + o::Real tol = pumipic::compute_tolerance_from_area(elmArea); + + // *********************************************** Run The Search *********************************************// + empty_function emptyFunction; + + printf("*** Searching ... ***\n"); + // After a single search operation, auxiliary arrays will be filled and they will be tested + bool success = pumipic::trace_particle_through_mesh(mesh, ptcls, particle_init_position, particle_final_position, + pid_d, elem_ids, next_elements, requireIntersection, inter_faces, inter_points, last_exit, 1, true, emptyFunction, elmArea, tol); + printf("*** Search Done ***\n"); + if (success){ + printf("[ERROR] Search Shouldn't pass...\n"); + } + else { + printf("Particle search failed as expected...\n"); + } + + // ******************************************* Checks *********************************************************// + OMEGA_H_CHECK_PRINTF(inter_faces.size() == ptcls->capacity(), "inter faces and ptcls capacity mismatch(%d,%d)", + inter_faces.size(), ptcls->capacity()); + OMEGA_H_CHECK_PRINTF(elem_ids.size() == ptcls->capacity(), "elem ids and ptcls capacity mismatch(%d,%d)", + elem_ids.size(), ptcls->capacity()); + OMEGA_H_CHECK_PRINTF(elem_ids.size() == next_elements.size(), "elem ids and next elements size mismatch(%d,%d)", + elem_ids.size(), next_elements.size()); + + Omega_h::Vector<3> expected_intersection {0.75, 0.5, 0.25}; + auto check_arrays = PS_LAMBDA(const int& e, const int& pid, const int& mask){ + if (mask>0) { + printf("Pid %d Intersection Face %d\n", pid, inter_faces[pid]); + printf("Pid %d Last Exit %d\n", pid, last_exit[pid]); + OMEGA_H_CHECK_PRINTF(last_exit[pid] == 10, "Expected face 10 between cell 0 and 5 but found %d\n", last_exit[pid]); + + printf("Pid %d intersects at (%f, %f, %f)\n", pid, inter_points[pid*3], inter_points[pid*3+1], inter_points[pid*3+2]); + OMEGA_H_CHECK_PRINTF(is_close_d(inter_points[pid*3], expected_intersection[0]), "Expected %f, found %f\n", expected_intersection[0], inter_points[pid*3]); + OMEGA_H_CHECK_PRINTF(is_close_d(inter_points[pid*3+1], expected_intersection[1]), "Expected %f, found %f\n", expected_intersection[1], inter_points[pid*3+1]); + OMEGA_H_CHECK_PRINTF(is_close_d(inter_points[pid*3+2], expected_intersection[2]), "Expected %f, found %f\n", expected_intersection[2], inter_points[pid*3+2]); + + printf("Pid %d Elem Id %d\n", pid, elem_ids[pid]); + OMEGA_H_CHECK_PRINTF(elem_ids[pid] == 0, "Expected element 0 but found %d\n", elem_ids[pid]); + + printf("Pid %d Next Element %d\n", pid, next_elements[pid]); + OMEGA_H_CHECK_PRINTF(next_elements[pid] == 5, "Expected next element 5 but found %d\n", next_elements[pid]); + } + }; + pumipic::parallel_for(ptcls, check_arrays, "check arrays"); + + printf("============================ Particle Moving for 0 to 5 Passed ===============================\n"); + printf("==============================================================================================\n"); + + + // *************************************************************************************************************// + + + // *************************************************** Check When Particles remain in the same element *********// + printf("\n\n\n==============================================================================================\n"); + printf(">============ Checking when particles remain in the same element ============================<\n"); + printf("==============================================================================================\n"); + + // set particle destination to the same element + Omega_h::Vector<3> another_location_in_cell0{0.5, 0.75, 0.2}; + + auto set_new_dest_in_cell0 = PS_LAMBDA (const int& e, const int& pid, const int& mask){ + if (mask > 0) { + particle_final_position(pid, 0) = another_location_in_cell0[0]; + particle_final_position(pid, 1) = another_location_in_cell0[1]; + particle_final_position(pid, 2) = another_location_in_cell0[2]; + + printf("New particle %d destination (%f, %f, %f)\n", + pid, particle_final_position(pid, 0), + particle_final_position(pid, 1), + particle_final_position(pid, 2)); + } + }; + pumipic::parallel_for(ptcls, set_new_dest_in_cell0, "set new destination in cell 0"); + + + // reset auxiliary arrays + elem_ids = Omega_h::Write(0, "element ids"); + inter_faces = Omega_h::Write(0, "inter faces"); + last_exit = Omega_h::Write(0, "last exit"); + inter_points = Omega_h::Write(0, "inter points"); + next_elements = Omega_h::Write(0, "next elements"); + + printf("*** Searching ... ***\n"); + // run the search again + success = pumipic::trace_particle_through_mesh(mesh, ptcls, particle_init_position, particle_final_position, + pid_d, elem_ids, next_elements, requireIntersection, inter_faces, inter_points, last_exit, 1, true, emptyFunction, elmArea, tol); + printf("*** Search Done ***\n"); + + if (success){ + printf("[ERROR] Search Shouldn't return success...\n"); + } + else { + printf("Particle search failed as expected...\n"); + } + + // ******************************************* Checks *********************************************************// + OMEGA_H_CHECK_PRINTF(inter_faces.size() == ptcls->capacity(), "inter faces and ptcls capacity mismatch(%d,%d)", + inter_faces.size(), ptcls->capacity()); + OMEGA_H_CHECK_PRINTF(elem_ids.size() == ptcls->capacity(), "elem ids and ptcls capacity mismatch(%d,%d)", + elem_ids.size(), ptcls->capacity()); + OMEGA_H_CHECK_PRINTF(elem_ids.size() == next_elements.size(), "elem ids and next elements size mismatch(%d,%d)", + elem_ids.size(), next_elements.size()); + + // expected intersection point is the same as the particle destination + auto check_move_in_same_element = PS_LAMBDA(const int& e, const int& pid, const int& mask){ + if (mask>0) { + printf("Pid %d Intersection Face %d\n", pid, inter_faces[pid]); + printf("Pid %d Last Exit %d\n", pid, last_exit[pid]); + OMEGA_H_CHECK_PRINTF(last_exit[pid] == -1, "Expected no intersection but found %d\n", last_exit[pid]); + + printf("Pid %d Elem Id %d\n", pid, elem_ids[pid]); + OMEGA_H_CHECK_PRINTF(elem_ids[pid] == 0, "Expected element 0 but found %d\n", elem_ids[pid]); + + printf("Pid %d intersects (reaches) at (%f, %f, %f) and expected (%f, %f, %f)\n", + pid, + inter_points[pid*3], inter_points[pid*3+1], inter_points[pid*3+2], + another_location_in_cell0[0], another_location_in_cell0[1], another_location_in_cell0[2]); + + OMEGA_H_CHECK_PRINTF(is_close_d(inter_points[pid*3], another_location_in_cell0[0]), "Expected %f, found %f\n", another_location_in_cell0[0], inter_points[pid * 3]); + OMEGA_H_CHECK_PRINTF(is_close_d(inter_points[pid*3+1], another_location_in_cell0[1]), "Expected %f, found %f\n", another_location_in_cell0[1], inter_points[pid * 3 + 1]); + OMEGA_H_CHECK_PRINTF(is_close_d(inter_points[pid*3+2], another_location_in_cell0[2]), "Expected %f, found %f\n", another_location_in_cell0[2], inter_points[pid * 3 + 2]); + + printf("Pid %d Next Element %d\n", pid, next_elements[pid]); + OMEGA_H_CHECK_PRINTF(next_elements[pid] == 0, "Expected next element 0 but found %d\n", next_elements[pid]); + } + }; + pumipic::parallel_for(ptcls, check_move_in_same_element, "check when particles move inside the same element as the origin(0 here)"); + + printf("============================ Particle Moving from 0 to 0 Passed ==============================\n"); + printf("==============================================================================================\n"); + + + + + + // *************************************************************************************************************// + // Particles moving from cell 0 to outside through face 0 // + // *************************************************************************************************************// + + printf("\n\n\n==============================================================================================\n"); + printf( ">============ Particle Moving out of the Mesh ===============================================<\n"); + printf( "==============================================================================================\n"); + + // set particle destination outside the mesh + Omega_h::Vector<3> outside_location{0.5, 0.75, -1.0}; + expected_intersection = {0.5, 0.75, 0.0}; + + auto set_new_dest_outside = PS_LAMBDA (const int& e, const int& pid, const int& mask){ + if (mask > 0) { + particle_final_position(pid, 0) = outside_location[0]; + particle_final_position(pid, 1) = outside_location[1]; + particle_final_position(pid, 2) = outside_location[2]; + + printf("New particle %d destination (%f, %f, %f)\n", + pid, particle_final_position(pid, 0), + particle_final_position(pid, 1), + particle_final_position(pid, 2)); + } + }; + pumipic::parallel_for(ptcls, set_new_dest_outside, "set new destination outside the mesh"); + + // reset auxiliary arrays + elem_ids = Omega_h::Write(0, "element ids"); + inter_faces = Omega_h::Write(0, "inter faces"); + last_exit = Omega_h::Write(0, "last exit"); + inter_points = Omega_h::Write(0, "inter points"); + next_elements = Omega_h::Write(0, "next elements"); + + printf("*** Searching ... ***\n"); + // run the search again + success = pumipic::trace_particle_through_mesh(mesh, ptcls, particle_init_position, particle_final_position, + pid_d, elem_ids, next_elements, requireIntersection, inter_faces, inter_points, last_exit, 1, true, emptyFunction, elmArea, tol); + printf("*** Search Done ***\n"); + + if (success){ + printf("[ERROR] Search Shouldn't return success...\n"); + } + else { + printf("Particle search failed as expected...\n"); + } + + // ******************************************* Checks *********************************************************// + OMEGA_H_CHECK_PRINTF(inter_faces.size() == ptcls->capacity(), "inter faces and ptcls capacity mismatch(%d,%d)", + inter_faces.size(), ptcls->capacity()); + OMEGA_H_CHECK_PRINTF(elem_ids.size() == ptcls->capacity(), "elem ids and ptcls capacity mismatch(%d,%d)", + elem_ids.size(), ptcls->capacity()); + OMEGA_H_CHECK_PRINTF(next_elements.size() == ptcls->capacity(), "next elem size and capacity mismatch(%d,%d)", + next_elements.size(), ptcls->capacity()); + OMEGA_H_CHECK_PRINTF(inter_points.size() == ptcls->capacity()*3, "inter points and ptcls capacity mismatch(%d,%d)", + inter_points.size(), ptcls->capacity()*3); + + auto check_move_outside = PS_LAMBDA(const int& e, const int& pid, const int& mask){ + if (mask>0) { + printf("Pid %d Intersection Face %d\n", pid, inter_faces[pid]); + printf("Pid %d Last Exit %d\n", pid, last_exit[pid]); + OMEGA_H_CHECK_PRINTF(last_exit[pid] == 0, "Expected face 0 but found %d\n", last_exit[pid]); + + printf("Pid %d intersects (reaches) at (%f, %f, %f) and expected (%f, %f, %f)\n", + pid, + inter_points[pid*3], inter_points[pid*3+1], inter_points[pid*3+2], + expected_intersection[0], expected_intersection[1], expected_intersection[2]); + + OMEGA_H_CHECK_PRINTF(is_close_d(inter_points[pid*3], expected_intersection[0]), "Expected %f, found %f\n", expected_intersection[0], inter_points[pid * 3]); + OMEGA_H_CHECK_PRINTF(is_close_d(inter_points[pid*3+1], expected_intersection[1]), "Expected %f, found %f\n", expected_intersection[1], inter_points[pid * 3 + 1]); + OMEGA_H_CHECK_PRINTF(is_close_d(inter_points[pid*3+2], expected_intersection[2]), "Expected %f, found %f\n", expected_intersection[2], inter_points[pid * 3 + 2]); + + printf("Pid %d Elem Id %d\n", pid, elem_ids[pid]); + OMEGA_H_CHECK_PRINTF(elem_ids[pid] == 0, "Expected element 0 but found %d\n", elem_ids[pid]); + + printf("Pid %d Next Element %d\n", pid, next_elements[pid]); + OMEGA_H_CHECK_PRINTF(next_elements[pid] == -1, "Expected next element -1 but found %d\n", next_elements[pid]); + } + }; + pumipic::parallel_for(ptcls, check_move_outside, "check when particles move outside the mesh"); + + + printf("\n\n\n==============================================================================================\n"); + printf( ">============ Particle Moving out of the Mesh Passed ========================================<\n"); + + + // ******************************************* Clean Up *********************************************************// + delete ptcls; + return 0; +} + +void printf_face_info(o::Mesh &mesh, o::LOs faceIds, bool all) { + const auto exposed_faces = o::mark_exposed_sides(&mesh); + const auto &face2nodes = mesh.ask_down(o::FACE, o::VERT).ab2b; + const auto &face2cellcell = mesh.ask_up(o::FACE, o::REGION).ab2b; + const auto &face2celloffset = mesh.ask_up(o::FACE, o::REGION).a2ab; + const auto &coords = mesh.coords(); + + auto print_faces = OMEGA_H_LAMBDA(o::LO + faceid) { + if (!all) { + for (int i = 0; i < faceIds.size(); i++) { + o::LO id = faceIds[i]; + if (id == faceid) { + printf("Face %d nodes %d %d %d Exposed %d\n", faceid, + face2nodes[faceid * 3], face2nodes[faceid * 3 + 1], + face2nodes[faceid * 3 + 2], exposed_faces[faceid]); + } + } + } else if (all) { + printf("Face %d nodes %d %d %d Exposed %d\n", faceid, + face2nodes[faceid * 3], face2nodes[faceid * 3 + 1], + face2nodes[faceid * 3 + 2], exposed_faces[faceid]); + int n_adj_cells = face2celloffset[faceid+1] - face2celloffset[faceid]; + if(n_adj_cells == 1){ + printf("Face %d is exposed and with %d cell\n", faceid, face2cellcell[face2celloffset[faceid]]); + } else if (n_adj_cells ==2) { + printf("Face %d internal with %d %d\n", faceid, face2cellcell[face2celloffset[faceid]], + face2cellcell[face2celloffset[faceid] + 1]); + } else { + printf("[ERROR] Face cannot be adjacent to more than 2 cells\n found %d\n", n_adj_cells); + } + } + }; + o::parallel_for(mesh.nfaces(), print_faces, "print asked faces"); +} + +void print_test_info() { + std::string message = "\n========================== Test Trace Particle ======================\n"; + message += "Testing 'trace_particle_through_mesh' function form pumipic_adjacency.tpp ...\n"; + message += "A 3D internal mesh is created and particles will be moved to a new element\n"; + message += "or will be kept in the same element.\n"; + message += "The boundary handler won't be implemeted here but the main function will do minimum\n"; + message += "to keep the particles moving\n"; + message += "\n===============================================================\n"; + printf("%s", message.c_str()); +} \ No newline at end of file diff --git a/test/testing.cmake b/test/testing.cmake index 8665e448..cc2f6f49 100644 --- a/test/testing.cmake +++ b/test/testing.cmake @@ -2,8 +2,10 @@ mpi_test(barycentric_3 1 ./barycentric test1) #Search tests -mpi_test(test_adj_2d 1 ./test_adj ${TEST_DATA_DIR}/plate/tri8.osh) -mpi_test(test_adj_3d 1 ./test_adj ${TEST_DATA_DIR}/cube/7k.osh) +mpi_test(test_trace_particle_through_mesh 1 ./test_trace_particle_through_mesh) +mpi_test(test_search_mesh 1 ./test_search_mesh) +#mpi_test(test_adj_2d 1 ./test_adj ${TEST_DATA_DIR}/plate/tri8.osh) +#mpi_test(test_adj_3d 1 ./test_adj ${TEST_DATA_DIR}/cube/7k.osh) mpi_test(moller_trumbore_test 1 ./moller_trumbore_test