SCOREC / pumi-pic

support libraries for unstructured mesh particle in cell simulations on GPUs and CPUs
BSD 3-Clause "New" or "Revised" License
36 stars 15 forks source link

Added line_segment flag and Unit Test for Moller Trumbore Line Triangle Intersection Function #131

Closed Fuad-HH closed 2 weeks ago

Fuad-HH commented 2 weeks ago

The moller_trumbore_line_triangle function returns true if a ray (given the origin and destination coordinates) intersects a triangle. But, for some applications, we need it to be intersected only if the intersection point lies within the line segment. In the new version, it takes a flag which if set to true checks if the the line segment intersects within origin to destination, not extrapolated.

In the test, the line enters the tetrahedral (id = 12) and stops inside it. It has faces 28, 29, 33, and 39. Face 28 (face nodes are 5,11,4) is the top face where the given line segment doesn't intersect but when treated as a ray it intersects.

Ray origin: 0.000000 -0.200000 -0.500000 Ray Destination: 0.000000 -0.200000 0.900000

Face 28 Nodes: 5 (1.000000 -1.000000 1.000000), 11 (0.000000 0.000000 1.000000), 4 (-1.000000 -1.000000 1.000000)

An interactive plot for the tested line and faces is attached here. The HTML file can be downloaded and opened in a browser to view it. https://drive.google.com/file/d/1H1K41VDQc-BldMaXLUYcWaWQxLxGBkzi/view?usp=sharing

jacobmerson commented 2 weeks ago

We may want to wrap the current version in a function that does that final check so it will not be called if not needed.

Fuad-HH commented 2 weeks ago

If I do it then the line segment length has to be calculated twice effectively since the original one doesn't return the length of the line segment. Moreover, the original function doesn't return t which is needed to decide if it intersected. One solution is having a lower-level function that gets the direction and triangle and returns the needed things that we can use in two different wrappers. Like the following:

  1. Function ray_passes_through_triangle gets the direction and triangle to find the parameter t and intersects flag and return them.
  2. Two wrappers get the origin and destination, find the direction to give it to ray_passes_through_triangle and save the length.
  3. The wrappers decide if intersect based on their own scheme.
jacobmerson commented 2 weeks ago

I think if you do the following it's the least amount of work and will have the effect we want.

// new function that gives back the parametric coordinate
// put in the comment that we are using moller_trumbore method here
 OMEGA_H_DEVICE bool ray_intersects_triangle(const o::Few<o::Vector<3>, 3>& faceVerts,
                                                    const o::Vector<3>& orig, const o::Vector<3>& dest,
                                                    o::Vector<3>& xpoint, const o::Real tol, const o::LO flip,
                                                    o::Real& dproj, o::Real& closeness, double& intersection_parametric_coord);

// current function that really should be ray/line intersection
// consider marking as deprecated https://en.cppreference.com/w/cpp/language/attributes/deprecated
 OMEGA_H_DEVICE bool moller_trumbore_line_triangle(const o::Few<o::Vector<3>, 3>& faceVerts,
                                                    const o::Vector<3>& orig, const o::Vector<3>& dest,
                                                    o::Vector<3>& xpoint, const o::Real tol, const o::LO flip,
                                                    o::Real& dproj, o::Real& closeness)
{
 double parametric_coord;
 return ray_triangle_intersection(faceVerts,
                                                     orig, dest,
                                                    xpoint, tol, const o::LO flip,
                                                    dproj, closeness, parametric_coord);
}

// new function that just checks line segment triangle intersection
// might as well pass back the parametric coord since someone may want to use it
 OMEGA_H_DEVICE bool line_segment_intersects_triangle(const o::Few<o::Vector<3>, 3>& faceVerts,
                                                    const o::Vector<3>& orig, const o::Vector<3>& dest,
                                                    o::Vector<3>& xpoint, const o::Real tol, const o::LO flip,
                                                    o::Real& dproj, o::Real& closeness, double& intersection_parametric_coord)
{
          line_triangle_intersection(faceVerts,
                                                     orig, dest,
                                                    xpoint, tol, const o::LO flip,
                                                    dproj, closeness, intersection_parametric_coord);
          return intersection_parametric_coord > 0 && intersection_parametric_coord <1;
}

New functions could be called ray_triangle_intersection and line_segment_triangle_intersection instead, but the suggested name makes it more clear what the bool return value actually means.

Fuad-HH commented 2 weeks ago

Thank you so much for the suggestion. I have added the changes in the above commit.

Fuad-HH commented 2 weeks ago

These tests need the newly added meshes in the pumipic-data pull request 2

jacobmerson commented 2 weeks ago

I merged the changes in pumi-pic data

jacobmerson commented 2 weeks ago

@Fuad-HH can you also update the usage of the moller_trombole function inside of search_mesh? We will most likely want to change that to line_intersects version after, but let's start by using the version with the more obvious name.

jacobmerson commented 2 weeks ago

Otherwise this looks great.

jacobmerson commented 2 weeks ago

@cwsmith this is ready to review. This does not change the behavior of the search method.