Skip to content

Commit

Permalink
test find_exit_face without bcc based search
Browse files Browse the repository at this point in the history
  • Loading branch information
@Fuad-HH
Fuad-HH committed Oct 24, 2024
1 parent 4cb57af commit afc8c74
Showing 1 changed file with 74 additions and 24 deletions.
98 changes: 74 additions & 24 deletions test/test_find_exit_face.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,37 @@ typedef MemberTypes<Vector3d, Vector3d, int> Particle;
typedef ps::ParticleStructure<Particle> PS;
typedef Kokkos::DefaultExecutionSpace ExeSpace;

OMEGA_H_DEVICE bool all_positive(const o::Vector<3>& a, Omega_h::Real tol = EPSILON) {
auto isPos = 1;
for (o::LO i = 0; i < a.size(); ++i) {
const auto gtez = o::are_close(a[i], 0.0, tol, tol) || a[i] > 0;
isPos = isPos && gtez;
}
return isPos;
}

bool is_inside(o::Mesh &mesh, o::LO elem_id, const o::Vector<2> point) {
const auto &coords = mesh.coords();
const auto &faces2nodes = mesh.ask_verts_of(o::FACE);

o::Write<o::LO> inside(2, 0);

auto is_inside_lambda = OMEGA_H_LAMBDA(o::LO id){
const auto current_el_verts = o::gather_verts<3>(faces2nodes, elem_id);
const Omega_h::Few<Omega_h::Vector<2>, 3> current_el_vert_coords =
o::gather_vectors<3, 2>(coords, current_el_verts);
// check if the particle is in the element
o::Vector<3> bcc = o::barycentric_from_global<2, 2>(
point, current_el_vert_coords);
// if all bcc are positive, the point is inside the element
inside[0] = all_positive(bcc,0.0);
};
o::parallel_for(1, is_inside_lambda);
auto host_inside = o::HostWrite(inside);

return bool(host_inside[0]);
}

void print_2D_mesh_edges(o::Mesh &mesh) {
printf("------------------ Mesh Info ------------------\n");
const auto edge2node = mesh.ask_down(o::EDGE, o::VERT).ab2b;
Expand All @@ -81,7 +112,7 @@ void print_2D_mesh_edges(o::Mesh &mesh) {
printf("------------------ Mesh Info ------------------\n");
}

bool test_2D_intersection(const std::string mesh_fname, p::Library *lib) {
bool test_2D_intersection(const std::string mesh_fname, p::Library *lib, bool useBcc) {
printf("Test: 2D intersection...\n");
o::Library &olib = lib->omega_h_lib();
// load mesh
Expand Down Expand Up @@ -141,36 +172,48 @@ bool test_2D_intersection(const std::string mesh_fname, p::Library *lib) {

// const o::Vector<3> start_point {0.0, 2.0, -3.0};
// const o::Vector<3> end_point {0.0, -4.0, -3.0};
const o::Vector<3> start_point{2.52054, -3.19384, 0};
const o::Vector<3> start_point2{-3.03604, -3.193840, 0};
const o::Vector<3> end_point{-0.294799, 3.07299, 0};
const o::Vector<2> start_point{2.52054, -3.19384};
const o::Vector<2> start_point2{-3.03604, -3.193840};
const o::Vector<2> end_point{-0.294799, 3.07299};
if (!is_inside(mesh, 1, start_point)) {
printf("Error: Start point is not inside the expected element.\n");
exit(1);
}
if (!is_inside(mesh, 1, start_point2)) {
printf("Error: Start point 1 is not inside the expected element.\n");
exit(1);
}
if (is_inside(mesh, 1, end_point)) {
printf("Error: End point is incorrectly inside the expected element.\n");
exit(1);
}

printf("Start point y: %f %f %f\n", start_point[0], start_point[1],
start_point[2]);
printf("Start point x: %f %f %f\n", start_point2[0], start_point2[1],
start_point2[2]);
printf("End point: %f %f %f\n", end_point[0], end_point[1], end_point[2]);
printf("Start point y: %f %f\n", start_point[0],
start_point[1]);
printf("Start point x: %f %f\n", start_point2[0],
start_point2[1]);
printf("End point: %f %f\n", end_point[0], end_point[1]);

auto set_initial_and_final_position =
PS_LAMBDA(const int &e, const int &pid, const int &mask) {
// see test calculation notebook
if (pid == 0) {
particle_init_position(pid, 0) = start_point[0];
particle_init_position(pid, 1) = start_point[1];
particle_init_position(pid, 2) = start_point[2];
particle_init_position(pid, 2) = 0.0;

particle_final_position(pid, 0) = end_point[0];
particle_final_position(pid, 1) = end_point[1];
particle_final_position(pid, 2) = end_point[2];
particle_final_position(pid, 2) = 0.0;
}
if (pid == 1) {
particle_init_position(pid, 0) = start_point2[0];
particle_init_position(pid, 1) = start_point2[1];
particle_init_position(pid, 2) = start_point2[2];
particle_init_position(pid, 2) = 0.0;

particle_final_position(pid, 0) = end_point[0];
particle_final_position(pid, 1) = end_point[1];
particle_final_position(pid, 2) = end_point[2];
particle_final_position(pid, 2) = 0.0;
}
};

Expand All @@ -179,12 +222,10 @@ bool test_2D_intersection(const std::string mesh_fname, p::Library *lib) {

// search for intersection
const auto psCapacity = ptcls->capacity();
o::Write<o::LO> elem_ids;
o::Write<o::Real> xpoints_d(3 * psCapacity, 0.0, "intersection points");
o::Write<o::LO> xface_id(psCapacity, 0.0, "intersection faces");
o::Write<o::LO> xface_id(psCapacity, -1, "intersection faces");

const auto elmArea = measure_elements_real(&mesh);
bool useBcc = true;
o::Real tol = p::compute_tolerance_from_area(elmArea);
auto inter_points =
o::Write<o::Real>(2 * psCapacity, 0.0, "intersection points");
Expand All @@ -197,17 +238,15 @@ bool test_2D_intersection(const std::string mesh_fname, p::Library *lib) {
particle_final_position, elem_ids_2d, ptcl_done, elmArea,
useBcc, lastExit, inter_points, tol);

printf("Last Exit Face: ");

auto lastExit_h = o::HostRead<o::LO>(lastExit);
bool passed_right = (lastExit_h[0] == 7);
bool passed_left = (lastExit_h[1] == 2);

if (!passed_right) {
printf("Expected exit face: 7, got %d\n", lastExit_h[0]);
printf("ERROR: Expected exit face: 7, got %d\n", lastExit_h[0]);
}
if (!passed_left) {
printf("Expected exit face: 2, got %d\n", lastExit_h[1]);
printf("ERROR: Expected exit face: 2, got %d\n", lastExit_h[1]);
}

// clean up
Expand All @@ -223,11 +262,22 @@ int main(int argc, char **argv) {
exit(1);
}
std::string mesh_fname = argv[1];
bool passed = test_2D_intersection(mesh_fname, &lib);
if (!passed) {
printf("Test failed. Run in verbose mode to view details.\n");

bool passed_bcc = test_2D_intersection(mesh_fname, &lib, true);
if (!passed_bcc) {
printf("ERROR: Test failed with BCC. Run in verbose mode to view details.\n");
}

bool passed_woBcc = test_2D_intersection(mesh_fname, &lib, false);

This comment has been minimized.

Sign in to view

Copy link
@Fuad-HH

Fuad-HH Oct 24, 2024

Author Contributor

@jacobmerson I have included the test to check if the exit face-finding algorithm works without the bcc based method. It also fails.
if (!passed_woBcc) {
printf("ERROR: Test failed without BCC. Run in verbose mode to view details.\n");
}

if (passed_bcc && passed_woBcc) {
printf("All tests passed!\n");
return 0;
} else {
return 1;
}

return 0;
}

0 comments on commit afc8c74

Please sign in to comment.