Mortar Transfer Function

core/specfem/assembly/nonconforming_interfaces/dim2/impl/compute_intersection.cpp

@@ -165,3 +165,161 @@ specfem::assembly::nonconforming_interfaces_impl::compute_intersection(
 
   return intersections;
 }
+
+void specfem::assembly::nonconforming_interfaces_impl::set_transfer_functions(
+    const Kokkos::View<
+        specfem::point::global_coordinates<specfem::dimension::type::dim2> *,
+        Kokkos::HostSpace> &element1,
+    const Kokkos::View<
+        specfem::point::global_coordinates<specfem::dimension::type::dim2> *,
+        Kokkos::HostSpace> &element2,
+    const specfem::mesh_entity::type &edge1,
+    const specfem::mesh_entity::type &edge2,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &mortar_quadrature,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &element_quadrature,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function1,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function1_prime,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function2,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function2_prime) {
+
+  // ======= ensure array shapes are correct =======
+  const auto nquad_mortar = mortar_quadrature.extent(0);
+  const auto ngll = element_quadrature.extent(0);
+  if ((transfer_function1.extent(0) != nquad_mortar ||
+       transfer_function1_prime.extent(0) != nquad_mortar ||
+       transfer_function1.extent(1) != ngll ||
+       transfer_function1_prime.extent(1) != ngll) ||
+      (transfer_function2.extent(0) != nquad_mortar ||
+       transfer_function2_prime.extent(0) != nquad_mortar ||
+       transfer_function2.extent(1) != ngll ||
+       transfer_function2_prime.extent(1) != ngll)) {
+    std::ostringstream oss;
+    oss << "Incompatible dimensions of `Kokkos::View`s in "
+           "specfem::assembly::nonconforming_interfaces_impl:"
+           ":set_transfer_functions:\n";
+    if (transfer_function1.extent(0) != nquad_mortar ||
+        transfer_function1_prime.extent(0) != nquad_mortar ||
+        transfer_function2.extent(0) != nquad_mortar ||
+        transfer_function2_prime.extent(0) != nquad_mortar) {
+      oss << "Mortar quadrature has " << nquad_mortar
+          << " quadrature points, which should match the first axis of the "
+             "transfer function tensor.";
+    }
+    if (transfer_function1.extent(1) != ngll ||
+        transfer_function1_prime.extent(1) != ngll ||
+        transfer_function2.extent(1) != ngll ||
+        transfer_function2_prime.extent(1) != ngll) {
+      oss << "Edge quadrature (element ngll) has " << ngll
+          << " quadrature points, which should match the second axis of the "
+             "transfer function tensor.";
+    }
+    oss << "\n Shape of\n"
+        << "        transfer_function1: (" << transfer_function1.extent(0)
+        << ", " << transfer_function1.extent(1) << ")\n"
+        << "  transfer_function1_prime: (" << transfer_function1_prime.extent(0)
+        << ", " << transfer_function1_prime.extent(1) << ")\n"
+        << "        transfer_function2: (" << transfer_function2.extent(0)
+        << ", " << transfer_function2.extent(1) << ")\n"
+        << "  transfer_function2_prime: (" << transfer_function2_prime.extent(0)
+        << ", " << transfer_function2_prime.extent(1) << ")\n";
+    throw std::runtime_error(oss.str());
+  }
+
+  // ======= populate transfer function and deriv =======
+
+  const auto intersections =
+      specfem::assembly::nonconforming_interfaces_impl::compute_intersection(
+          element1, element2, edge1, edge2, mortar_quadrature);
+
+  for (int iquad = 0; iquad < nquad_mortar; iquad++) {
+    {
+      const auto [hxi, hpxi] =
+          specfem::quadrature::gll::Lagrange::compute_lagrange_interpolants(
+              intersections[iquad].first, ngll, element_quadrature);
+      for (int igll = 0; igll < ngll; igll++) {
+        transfer_function1(iquad, igll) = hxi(igll);
+        transfer_function1_prime(iquad, igll) = hpxi(igll);
+      }
+    }
+    {
+      const auto [hxi, hpxi] =
+          specfem::quadrature::gll::Lagrange::compute_lagrange_interpolants(
+              intersections[iquad].second, ngll, element_quadrature);
+      for (int igll = 0; igll < ngll; igll++) {
+        transfer_function2(iquad, igll) = hxi(igll);
+        transfer_function2_prime(iquad, igll) = hpxi(igll);
+      }
+    }
+  }
+}
+
+void specfem::assembly::nonconforming_interfaces_impl::set_transfer_functions(
+    const Kokkos::View<
+        specfem::point::global_coordinates<specfem::dimension::type::dim2> *,
+        Kokkos::HostSpace> &element1,
+    const Kokkos::View<
+        specfem::point::global_coordinates<specfem::dimension::type::dim2> *,
+        Kokkos::HostSpace> &element2,
+    const specfem::mesh_entity::type &edge1,
+    const specfem::mesh_entity::type &edge2,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &mortar_quadrature,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &element_quadrature,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function1,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function2) {
+
+  // ======= ensure array shapes are correct =======
+  const auto nquad_mortar = mortar_quadrature.extent(0);
+  const auto ngll = element_quadrature.extent(0);
+  if ((transfer_function1.extent(0) != nquad_mortar ||
+       transfer_function1.extent(1) != ngll) ||
+      (transfer_function2.extent(0) != nquad_mortar ||
+       transfer_function2.extent(1) != ngll)) {
+    std::ostringstream oss;
+    oss << "Incompatible dimensions of `Kokkos::View`s in "
+           "specfem::assembly::nonconforming_interfaces_impl:"
+           ":set_transfer_functions:\n";
+    if (transfer_function1.extent(0) != nquad_mortar ||
+        transfer_function2.extent(0) != nquad_mortar) {
+      oss << "Mortar quadrature has " << nquad_mortar
+          << " quadrature points, which should match the first axis of the "
+             "transfer function tensor.";
+    }
+    if (transfer_function1.extent(1) != ngll ||
+        transfer_function2.extent(1) != ngll) {
+      oss << "Edge quadrature (element ngll) has " << ngll
+          << " quadrature points, which should match the second axis of the "
+             "transfer function tensor.";
+    }
+    oss << "\n Shape of\n"
+        << "  transfer_function1: (" << transfer_function1.extent(0) << ", "
+        << transfer_function1.extent(1) << ")\n"
+        << "  transfer_function2: (" << transfer_function2.extent(0) << ", "
+        << transfer_function2.extent(1) << ")";
+    throw std::runtime_error(oss.str());
+  }
+
+  // ======= populate transfer function =======
+
+  const auto intersections =
+      specfem::assembly::nonconforming_interfaces_impl::compute_intersection(
+          element1, element2, edge1, edge2, mortar_quadrature);
+
+  for (int iquad = 0; iquad < nquad_mortar; iquad++) {
+    {
+      const auto hxi = std::get<0>(
+          specfem::quadrature::gll::Lagrange::compute_lagrange_interpolants(
+              intersections[iquad].first, ngll, element_quadrature));
+      for (int igll = 0; igll < ngll; igll++) {
+        transfer_function1(iquad, igll) = hxi(igll);
+      }
+    }
+    {
+      const auto hxi = std::get<0>(
+          specfem::quadrature::gll::Lagrange::compute_lagrange_interpolants(
+              intersections[iquad].second, ngll, element_quadrature));
+      for (int igll = 0; igll < ngll; igll++) {
+        transfer_function2(iquad, igll) = hxi(igll);
+      }
+    }
+  }
+}

core/specfem/assembly/nonconforming_interfaces/dim2/impl/compute_intersection.hpp

@@ -40,4 +40,70 @@ std::vector<std::pair<type_real, type_real> > compute_intersection(
     const specfem::mesh_entity::type &edge2,
     const Kokkos::View<type_real *, Kokkos::HostSpace> &mortar_quadrature);
 
+/**
+ * @brief Populates the transfer function for a given intersection.
+ * The transfer function is a linear map from the element edge function basis to
+ * the mortar basis.
+ *
+ * @param element1 - global coordinate representation of the first element.
+ * @param element2 - global coordinate representation of the second element.
+ * @param edge1 - edge on the first element
+ * @param edge2 - edge on the second element
+ * @param mortar_quadrature - a list of local (mortar space) knots describing
+ *             　　　　　　　　　the quadrature rule on the interval [-1,1]
+ * @param transfer_function1 - the nquad_mortar x ngll transfer tensor to
+ * populate, mapping fields on edge1 to the mortar.
+ * @param transfer_function2 - the nquad_mortar x ngll transfer tensor to
+ * populate, mapping fields on edge2 to the mortar.
+ */
+void set_transfer_functions(
+    const Kokkos::View<
+        specfem::point::global_coordinates<specfem::dimension::type::dim2> *,
+        Kokkos::HostSpace> &element1,
+    const Kokkos::View<
+        specfem::point::global_coordinates<specfem::dimension::type::dim2> *,
+        Kokkos::HostSpace> &element2,
+    const specfem::mesh_entity::type &edge1,
+    const specfem::mesh_entity::type &edge2,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &mortar_quadrature,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &element_quadrature,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function1,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function2);
+
+/**
+ * @brief Populates the transfer function for a given intersection.
+ * The transfer function is a linear map from the element edge function basis to
+ * the mortar basis.
+ *
+ * @param element1 - global coordinate representation of the first element.
+ * @param element2 - global coordinate representation of the second element.
+ * @param edge1 - edge on the first element
+ * @param edge2 - edge on the second element
+ * @param mortar_quadrature - a list of local (mortar space) knots describing
+ *             　　　　　　　　　the quadrature rule on the interval [-1,1]
+ * @param transfer_function1 - the nquad_mortar x ngll transfer tensor to
+ * populate, mapping fields on edge1 to the mortar.
+ * @param transfer_function1_prime - derivative of transfer_function1 w.r.t.
+ * edge coordinate.
+ * @param transfer_function2 - the nquad_mortar x ngll transfer tensor to
+ * populate, mapping fields on edge2 to the mortar.
+ * @param transfer_function2_prime - derivative of transfer_function2 w.r.t.
+ * edge coordinate.
+ */
+void set_transfer_functions(
+    const Kokkos::View<
+        specfem::point::global_coordinates<specfem::dimension::type::dim2> *,
+        Kokkos::HostSpace> &element1,
+    const Kokkos::View<
+        specfem::point::global_coordinates<specfem::dimension::type::dim2> *,
+        Kokkos::HostSpace> &element2,
+    const specfem::mesh_entity::type &edge1,
+    const specfem::mesh_entity::type &edge2,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &mortar_quadrature,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &element_quadrature,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function1,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function1_prime,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function2,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function2_prime);
+
 } // namespace specfem::assembly::nonconforming_interfaces_impl

core/specfem/assembly/nonconforming_interfaces/dim2/impl/compute_intersection_in_assembly.hpp

@@ -8,6 +8,39 @@
 
 namespace specfem::assembly::nonconforming_interfaces_impl {
 
+/**
+ * @brief A helper function that retrieves the element global coordinates and
+ * edge orientation of a given adjacency in the adjacency graph.
+ *
+ * @param mesh - the assembly::mesh object
+ * @param edge - the adjacency graph edge of the intersection.
+ *               mesh.graph()[edge] should reference the edge we want to
+ *               construct the intersection between.
+ * @return std::tuple<
+ * Kokkos::View<
+ * specfem::point::global_coordinates<specfem::dimension::type::dim2> *,
+ * Kokkos::HostSpace>,
+ * Kokkos::View<
+ * specfem::point::global_coordinates<specfem::dimension::type::dim2> *,
+ * Kokkos::HostSpace>,
+ * specfem::mesh_entity::type, specfem::mesh_entity::type> - in order: the
+ * cooordinates of the source element, the coordinates of the target element,
+ * the orientation on the source element in the intersection, and the
+ * orientation on the target element in the intersection.
+ */
+template <typename EdgeType>
+std::tuple<
+    Kokkos::View<
+        specfem::point::global_coordinates<specfem::dimension::type::dim2> *,
+        Kokkos::HostSpace>,
+    Kokkos::View<
+        specfem::point::global_coordinates<specfem::dimension::type::dim2> *,
+        Kokkos::HostSpace>,
+    specfem::mesh_entity::type, specfem::mesh_entity::type>
+_expand_edge_index(
+    const specfem::assembly::mesh<specfem::dimension::type::dim2> &mesh,
+    const EdgeType &edge);
+
 /**
  * @brief Computes the intersection between two elements, returning the knots of
  * the "mortar", the codimension 1 element joining them at their intersection,
@@ -30,5 +63,62 @@ namespace specfem::assembly::nonconforming_interfaces_impl {
 template <typename EdgeType>
 std::vector<std::pair<type_real, type_real> > compute_intersection(
     const specfem::assembly::mesh<specfem::dimension::type::dim2> &mesh,
-    const EdgeType &edge, const Kokkos::View<type_real *, Kokkos::HostSpace> &mortar_quadrature);
+    const EdgeType &edge,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &mortar_quadrature);
 } // namespace specfem::assembly::nonconforming_interfaces_impl
+
+/**
+ * @brief Populates the transfer function for a given intersection.
+ * The transfer function is a linear map from the element edge function basis to
+ * the mortar basis.
+ *
+ * @param mesh - the assembly::mesh object
+ * @param edge - the adjacency graph edge of the intersection.
+ *               mesh.graph()[edge] should reference the edge we want to
+ *               construct the intersection between.
+ * @param mortar_quadrature - a list of local (mortar space) knots describing
+ *             　　　　　　　　　the quadrature rule on the interval [-1,1]
+ * @param transfer_function1 - the nquad_mortar x ngll transfer tensor to
+ * populate, mapping fields on edge1 to the mortar.
+ * @param transfer_function2 - the nquad_mortar x ngll transfer tensor to
+ * populate, mapping fields on edge2 to the mortar.
+ */
+template <typename EdgeType>
+void set_transfer_functions(
+    const specfem::assembly::mesh<specfem::dimension::type::dim2> &mesh,
+    const EdgeType &edge,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &mortar_quadrature,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &element_quadrature,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function1,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function2);
+
+/**
+ * @brief Populates the transfer function for a given intersection.
+ * The transfer function is a linear map from the element edge function basis to
+ * the mortar basis.
+ *
+ * @param mesh - the assembly::mesh object
+ * @param edge - the adjacency graph edge of the intersection.
+ *               mesh.graph()[edge] should reference the edge we want to
+ *               construct the intersection between.
+ * @param mortar_quadrature - a list of local (mortar space) knots describing
+ *             　　　　　　　　　the quadrature rule on the interval [-1,1]
+ * @param transfer_function1 - the nquad_mortar x ngll transfer tensor to
+ * populate, mapping fields on edge1 to the mortar.
+ * @param transfer_function1_prime - derivative of transfer_function1 w.r.t.
+ * edge coordinate.
+ * @param transfer_function2 - the nquad_mortar x ngll transfer tensor to
+ * populate, mapping fields on edge2 to the mortar.
+ * @param transfer_function2_prime - derivative of transfer_function2 w.r.t.
+ * edge coordinate.
+ */
+template <typename EdgeType>
+void set_transfer_functions(
+    const specfem::assembly::mesh<specfem::dimension::type::dim2> &mesh,
+    const EdgeType &edge,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &mortar_quadrature,
+    const Kokkos::View<type_real *, Kokkos::HostSpace> &element_quadrature,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function1,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function1_prime,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function2,
+    Kokkos::View<type_real **, Kokkos::HostSpace> &transfer_function2_prime);