Program Listing for File mesh_intersection.hpp

↰ Return to documentation for file (src/pcms/transfer/mesh_intersection.hpp)

#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