.. _program_listing_file_src_pcms_transfer_mesh_intersection.hpp:

Program Listing for File mesh_intersection.hpp
==============================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_pcms_transfer_mesh_intersection.hpp>` (``src/pcms/transfer/mesh_intersection.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #ifndef PCMS_TRANSFER_MESH_INTERSECTION_HPP
   #define PCMS_TRANSFER_MESH_INTERSECTION_HPP
   
   #include <pcms/configuration.h>
   #include <pcms/localization/point_search.h>
   #include <pcms/localization/queue_visited.hpp>
   #include <Omega_h_fail.hpp>
   #include <Omega_h_int_scan.hpp>
   #include <r3d.hpp>
   #include <Omega_h_mesh.hpp>
   #include <Omega_h_for.hpp>
   #include <Kokkos_MathematicalFunctions.hpp>
   
   namespace pcms
   {
   // Gather the vertex coordinates of a simplex element (triangle for Dim==2,
   // tetrahedron for Dim==3) into an r3d simplex, ready for
   // r3d::intersect_simplices.
   template <int Dim>
   [[nodiscard]] OMEGA_H_INLINE r3d::Few<r3d::Vector<Dim>, Dim + 1>
   get_vert_coords_of_elem(const Omega_h::Reals& coords,
                           const Omega_h::LOs& elems2nodes, const int id)
   {
     const auto elm_verts = Omega_h::gather_verts<Dim + 1>(elems2nodes, id);
   
     const Omega_h::Matrix<Dim, Dim + 1> elm_vert_coords =
       Omega_h::gather_vectors<Dim + 1, Dim>(coords, elm_verts);
   
     r3d::Few<r3d::Vector<Dim>, Dim + 1> r3d_vector;
     for (int i = 0; i < Dim + 1; ++i) {
       for (int d = 0; d < Dim; ++d) {
         r3d_vector[i][d] = elm_vert_coords[i][d];
       }
     }
   
     return r3d_vector;
   }
   
   template <int Dim>
   [[nodiscard]] OMEGA_H_INLINE bool simplices_have_degenerate_overlap(
     const r3d::Few<r3d::Vector<Dim>, Dim + 1>& a,
     const r3d::Few<r3d::Vector<Dim>, Dim + 1>& b)
   {
     // Each iteration treats one simplex as the half-space owner. Both are
     // needed: a face of `a` may separate the pair while no face of `b` does.
     for (int side = 0; side < 2; ++side) {
       const auto& owner = (side == 0) ? a : b;
       const auto& other = (side == 0) ? b : a;
   
       // Face `f` is the face opposite vertex `f`.
       for (int f = 0; f < Dim + 1; ++f) {
         // Vertices spanning the face, in index order with `f` skipped.
         int face_verts[Dim];
         int n_face_verts = 0;
         for (int v = 0; v < Dim + 1; ++v) {
           if (v != f) {
             face_verts[n_face_verts++] = v;
           }
         }
         const auto& origin = owner[face_verts[0]];
   
         // Face normal, oriented to point away from the opposite vertex so the
         // owner's interior lies on the negative side.
         r3d::Vector<Dim> normal;
         if constexpr (Dim == 2) {
           const auto edge = owner[face_verts[1]] - origin;
           normal[0] = edge[1];
           normal[1] = -edge[0];
         } else {
           const auto e1 = owner[face_verts[1]] - origin;
           const auto e2 = owner[face_verts[2]] - origin;
           normal[0] = e1[1] * e2[2] - e1[2] * e2[1];
           normal[1] = e1[2] * e2[0] - e1[0] * e2[2];
           normal[2] = e1[0] * e2[1] - e1[1] * e2[0];
         }
   
         Omega_h::Real opposite_side = 0.0;
         Omega_h::Real normal_sq = 0.0;
         for (int d = 0; d < Dim; ++d) {
           opposite_side += normal[d] * (owner[f][d] - origin[d]);
           normal_sq += normal[d] * normal[d];
         }
         // A degenerate face gives no usable axis.
         if (normal_sq == 0.0) {
           continue;
         }
         const Omega_h::Real orientation = (opposite_side > 0.0) ? -1.0 : 1.0;
   
         bool all_outside = true;
         for (int v = 0; v < Dim + 1 && all_outside; ++v) {
           Omega_h::Real signed_distance = 0.0;
           for (int d = 0; d < Dim; ++d) {
             signed_distance += normal[d] * (other[v][d] - origin[d]);
           }
           all_outside = (orientation * signed_distance) >= 0.0;
         }
         if (all_outside) {
           return true;
         }
       }
     }
     return false;
   }
   
   struct IntersectionResults
   {
     Omega_h::LOs tgt2src_offsets;
     Omega_h::LOs tgt2src_indices;
   };
   
   class FindIntersections
   {
   private:
     Omega_h::Mesh& source_mesh_;
     Omega_h::Mesh& target_mesh_;
   
   public:
     FindIntersections(Omega_h::Mesh& source_mesh, Omega_h::Mesh& target_mesh)
       : source_mesh_(source_mesh), target_mesh_(target_mesh)
     {
     }
   
     template <int Dim>
     void adjBasedIntersectSearch(const Kokkos::View<const LO*>& start_elements,
                                  const Omega_h::LOs& tgt2src_offsets,
                                  Omega_h::Write<Omega_h::LO>& nIntersections,
                                  Omega_h::Write<Omega_h::LO>& tgt2src_indices,
                                  bool is_count_only, bool use_prefilter = true);
   };
   
   IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh,
                                        Omega_h::Mesh& target_mesh,
                                        bool use_prefilter = true);
   
   IntersectionResults intersectTargets(
     Omega_h::Mesh& source_mesh, Omega_h::Mesh& target_mesh,
     const GridPointSearchVariant& source_search, bool use_prefilter = true);
   } // namespace pcms
   #endif // PCMS_TRANSFER_MESH_INTERSECTION_HPP