[DO NOT MERGE] 3D weakly_conforming -> nonconforming Override Test

core/specfem/algorithms/locate_point/dim3/locate_point_impl.cpp

@@ -591,9 +591,9 @@ std::pair<std::pair<type_real, type_real>, bool> locate_point(
                      Kokkos::HostSpace>
       coorg("coorg", mesh.ngnod);
   for (int i = 0; i < mesh.ngnod; i++) {
-    coorg(i).x = mesh.h_control_node_coordinates(0, ispec, i);
-    coorg(i).y = mesh.h_control_node_coordinates(1, ispec, i);
-    coorg(i).z = mesh.h_control_node_coordinates(2, ispec, i);
+    coorg(i).x = mesh.h_control_node_coordinates(ispec, i, 0);
+    coorg(i).y = mesh.h_control_node_coordinates(ispec, i, 1);
+    coorg(i).z = mesh.h_control_node_coordinates(ispec, i, 2);
   }
 
   // initial guess of 0 (center of edge)

core/specfem/assembly/assembly/dim3/assembly.cpp

@@ -76,6 +76,12 @@ specfem::assembly::assembly<specfem::element::dimension_tag::dim3>::assembly(
     ngllz,     nglly, ngllx, this->element_intersections, this->jacobian_matrix,
     this->mesh
   };
+
+  this->nonconforming_interfaces = {
+    ngllz,     nglly, ngllx, this->element_intersections, this->jacobian_matrix,
+    this->mesh
+  };
+
   this->fields = { this->mesh, this->element_types, simulation };
 
   this->mpi_interfaces = { mesh.adjacency_graph,

core/specfem/assembly/assembly/dim3/assembly.hpp

@@ -20,6 +20,8 @@
 #include "specfem/receivers.hpp"
 #include "specfem/source.hpp"
 
+#include "specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp"
+
 namespace specfem::io {
 class reader;
 }
@@ -125,6 +127,9 @@ template <> struct assembly<specfem::element::dimension_tag::dim3> {
   specfem::assembly::conforming_interfaces<dimension_tag>
       conforming_interfaces; ///< Conforming interfaces between different media
                              ///< (e.g., fluid-solid interface)
+  specfem::assembly::nonconforming_interfaces<dimension_tag>
+      nonconforming_interfaces; ///< Nonconforming interfaces between different
+                                ///< media (e.g., fluid-solid interface)
 
   specfem::assembly::boundary_values<dimension_tag>
       boundary_values; ///< Field
@@ -166,23 +171,22 @@ template <> struct assembly<specfem::element::dimension_tag::dim3> {
    * assignment if exists)
    * @param flux_scheme_config Flux scheme rules for nonconforming interfaces
    */
-  assembly(
-      const specfem::mesh::mesh<dimension_tag> &mesh,
-      const specfem::quadrature::quadratures &quadratures,
-      std::vector<std::shared_ptr<specfem::sources::source<dimension_tag> > >
-          &sources,
-      const std::vector<
-          std::shared_ptr<specfem::receivers::receiver<dimension_tag> > >
-          &receivers,
-      const std::vector<specfem::enums::wavefield> &stypes, const type_real t0,
-      const type_real dt, const int max_timesteps, const int max_sig_step,
-      const int nsteps_between_samples,
-      const specfem::simulation::type simulation,
-      const bool allocate_boundary_values,
-      const std::shared_ptr<specfem::io::reader> &property_reader,
-      const specfem::element_coupling::flux_scheme_configuration
-          &flux_scheme_config =
-              specfem::element_coupling::flux_scheme_configuration());
+  assembly(const specfem::mesh::mesh<dimension_tag> &mesh,
+           const specfem::quadrature::quadratures &quadratures,
+           std::vector<std::shared_ptr<specfem::sources::source<dimension_tag>>>
+               &sources,
+           const std::vector<
+               std::shared_ptr<specfem::receivers::receiver<dimension_tag>>>
+               &receivers,
+           const std::vector<specfem::enums::wavefield> &stypes,
+           const type_real t0, const type_real dt, const int max_timesteps,
+           const int max_sig_step, const int nsteps_between_samples,
+           const specfem::simulation::type simulation,
+           const bool allocate_boundary_values,
+           const std::shared_ptr<specfem::io::reader> &property_reader,
+           const specfem::element_coupling::flux_scheme_configuration
+               &flux_scheme_config =
+                   specfem::element_coupling::flux_scheme_configuration());
 
   assembly() = default;
 

core/specfem/assembly/element_intersections/dim3/element_intersections.cpp

@@ -1,5 +1,6 @@
 #include "specfem/assembly/element_intersections.hpp"
 #include "specfem/assembly/element_types.hpp"
+#include "specfem/element_connections/tags.hpp"
 #include "specfem/enums.hpp"
 #include "specfem/macros.hpp"
 #include "specfem/mesh.hpp"
@@ -34,8 +35,8 @@ specfem::assembly::element_intersections<
 
   // Collect self/coupled face vectors for each tag combination
   struct CollectedFaces {
-    std::vector<specfem::mesh_entity::face<dimension_tag> > self_collect;
-    std::vector<specfem::mesh_entity::face<dimension_tag> > coupled_collect;
+    std::vector<specfem::mesh_entity::face<dimension_tag>> self_collect;
+    std::vector<specfem::mesh_entity::face<dimension_tag>> coupled_collect;
   };
 
   specfem::tag_dispatch::Storage<CollectedFaces, IntersectionCombinations>
@@ -45,15 +46,25 @@ specfem::assembly::element_intersections<
         constexpr auto coupled_medium = specfem::element_coupling::attributes<
             TagsType::dimension_tag, TagsType::interface_tag>::coupled_medium();
 
-        std::vector<specfem::mesh_entity::face<dimension_tag> > self_collect;
-        std::vector<specfem::mesh_entity::face<dimension_tag> > coupled_collect;
+        std::vector<specfem::mesh_entity::face<dimension_tag>> self_collect;
+        std::vector<specfem::mesh_entity::face<dimension_tag>> coupled_collect;
 
         const auto &graph = mesh.graph();
 
         // Filter out corresponding connections
         auto filter = [&graph](const auto &edge) {
+          if constexpr (TagsType::connection_tag ==
+                        specfem::element_connections::type::weakly_conforming) {
+            return false;
+          } else if constexpr (TagsType::connection_tag ==
+                               specfem::element_connections::type::
+                                   nonconforming) {
+            return true;
+          }
           return graph[edge].connection == TagsType::connection_tag;
-        };
+        }; // TODO yell at Hanson if this is still in when the PR is reviewed.
+           // This is to test the nonconforming kernel on a conforming mesh
+           // prior to a NC mesher.
 
         // Create a filtered graph view
         const auto &nc_graph = boost::make_filtered_graph(graph, filter);
@@ -164,7 +175,7 @@ specfem::assembly::element_intersections<
 specfem::assembly::face_view_from_collected_faces(
     const std::string &label,
     const std::vector<
-        specfem::mesh_entity::face<specfem::element::dimension_tag::dim3> >
+        specfem::mesh_entity::face<specfem::element::dimension_tag::dim3>>
         &collected_faces,
     const specfem::mesh_entity::element<specfem::element::dimension_tag::dim3>
         &element) {

core/specfem/assembly/fields/data_access/impl/load_access_functions.hpp

@@ -256,4 +256,48 @@ load_after_simd_dispatch(const std::false_type, const IndexType &index,
   return;
 }
 
+template <bool on_device, typename IndexType, typename ContainerType,
+          typename... AccessorTypes,
+          typename std::enable_if_t<
+              (specfem::data_access::is_index_type<IndexType>::value &&
+               specfem::data_access::is_chunk_face<IndexType>::value &&
+               (specfem::data_access::is_field<AccessorTypes>::value && ...)),
+              int> = 0>
+KOKKOS_FORCEINLINE_FUNCTION void
+load_after_simd_dispatch(const std::false_type, const IndexType &index,
+                         const ContainerType &field,
+                         AccessorTypes &...accessors) {
+  static_assert(!IndexType::using_simd,
+                "IndexType must not use SIMD in this overload");
+  constexpr static auto MediumTag =
+      std::tuple_element_t<0, std::tuple<AccessorTypes...>>::medium_tag;
+  // Check that all accessors have the same medium tag
+  specfem::assembly::fields_impl::check_accessor_compatibility<
+      AccessorTypes...>();
+  const auto &current_field = field.template get_field<MediumTag>();
+
+  constexpr static int ncomponents = specfem::element::attributes<
+      std::tuple_element_t<0, std::tuple<AccessorTypes...>>::dimension_tag,
+      std::tuple_element_t<0, std::tuple<AccessorTypes...>>::medium_tag>::
+      components;
+
+  specfem::execution::for_each_level(
+      index.get_iterator(),
+      [&](const typename IndexType::iterator_type::index_type &iterator_index) {
+        const auto point_index = iterator_index.get_local_index();
+        const int iglob =
+            field.template get_iglob<on_device, MediumTag>(point_index);
+        for (int icomp = 0; icomp < ncomponents; ++icomp) {
+          (specfem::assembly::fields_impl::base_load_accessor<
+               on_device, AccessorTypes::data_class>(
+               iglob, icomp, current_field,
+               accessors(point_index.iface, point_index.ipoint_i,
+                         point_index.ipoint_j, icomp)),
+           ...);
+        }
+      });
+
+  return;
+}
+
 } // namespace specfem::assembly::simulation_field_impl

core/specfem/assembly/fields/data_access/load_on_device.hpp

@@ -22,7 +22,7 @@ KOKKOS_FORCEINLINE_FUNCTION void load_on_device(const IndexType &index,
                                                 AccessorTypes &...accessors) {
 
   constexpr static auto MediumTag =
-      std::tuple_element_t<0, std::tuple<AccessorTypes...> >::medium_tag;
+      std::tuple_element_t<0, std::tuple<AccessorTypes...>>::medium_tag;
 
   // Check that all accessors have the same medium tag
   fields_impl::check_accessor_compatibility<AccessorTypes...>();
@@ -40,15 +40,16 @@ template <typename IndexType, typename ContainerType, typename... AccessorTypes,
               ((specfem::data_access::is_index_type<IndexType>::value) &&
                (specfem::data_access::is_point<IndexType>::value ||
                 specfem::data_access::is_chunk_element<IndexType>::value ||
-                specfem::data_access::is_chunk_edge<IndexType>::value) &&
+                specfem::data_access::is_chunk_edge<IndexType>::value ||
+                specfem::data_access::is_chunk_face<IndexType>::value) &&
                (specfem::data_access::is_field<AccessorTypes>::value && ...)),
               int> = 0>
 KOKKOS_FORCEINLINE_FUNCTION void load_on_device(const IndexType &index,
                                                 const ContainerType &field,
                                                 AccessorTypes &...accessors) {
 
   constexpr static auto MediumTag =
-      std::tuple_element_t<0, std::tuple<AccessorTypes...> >::medium_tag;
+      std::tuple_element_t<0, std::tuple<AccessorTypes...>>::medium_tag;
 
   // Check that all accessors have the same medium tag
   fields_impl::check_accessor_compatibility<AccessorTypes...>();

core/specfem/assembly/nonconforming_interfaces.hpp

@@ -1,30 +1,4 @@
 #pragma once
 
-#include "specfem/element_connections.hpp"
-#include "specfem/element_coupling.hpp"
-#include "specfem/element_coupling/tags.hpp"
-#include "specfem/enums.hpp"
-
-namespace specfem::assembly {
-
-namespace nonconforming_interfaces_impl {
-
-template <specfem::element::dimension_tag DimensionTag,
-          specfem::element_coupling::interface_tag InterfaceTag,
-          specfem::element::boundary_tag BoundaryTag,
-          specfem::element_connections::type ConnectionTag,
-          specfem::element_coupling::flux_scheme_tag FluxSchemeTag>
-struct interface_container;
-
-} // namespace nonconforming_interfaces_impl
-
-/**
- * @brief Information on coupled interfaces between two mediums
- * @tparam DimensionTag Dimension of spectral elements
- */
-template <specfem::element::dimension_tag DimensionTag>
-struct nonconforming_interfaces;
-
-} // namespace specfem::assembly
-
 #include "nonconforming_interfaces/dim2/nonconforming_interfaces.hpp"
+#include "nonconforming_interfaces/dim3/nonconforming_interfaces.hpp"

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

@@ -3,7 +3,7 @@
 
 #include "specfem/assembly/element_intersections.hpp"
 #include "specfem/assembly/mesh.hpp"
-#include "specfem/data_access.hpp"
+#include "specfem/assembly/nonconforming_interfaces/fwd.hpp"
 #include "specfem/element_coupling/flux_scheme_configuration.hpp"
 #include "specfem/element_coupling/tags.hpp"
 #include "specfem/enums.hpp"
@@ -169,8 +169,7 @@ struct interface_container<
                   "impl_load requires edge_index type for IndexType");
 
     constexpr int nquad_intersection =
-        std::tuple_element_t<0,
-                             std::tuple<PointTypes...> >::n_quad_intersection;
+        std::tuple_element_t<0, std::tuple<PointTypes...>>::n_quad_intersection;
 
     const auto assign = [&](auto &point, auto &iquad, auto &index) {
       using PointType = std::decay_t<decltype(point)>;
@@ -204,7 +203,7 @@ struct interface_container<
                   "impl_load requires chunk_edge type for IndexType");
 
     constexpr int nquad_intersection =
-        std::tuple_element_t<0, std::tuple<EdgeTypes...> >::n_quad_intersection;
+        std::tuple_element_t<0, std::tuple<EdgeTypes...>>::n_quad_intersection;
 
     const auto factor_subview = [&]() {
       if constexpr (on_device) {

core/specfem/assembly/nonconforming_interfaces/dim3/data_access/load.hpp

@@ -0,0 +1,108 @@
+#pragma once
+
+#include "specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp"
+#include "specfem/data_access.hpp"
+
+namespace specfem::assembly {
+
+/**
+ * @defgroup CoupledInterfaceDataAccess
+ * @brief Data access functions for coupled interface computation data
+ *
+ */
+
+/**
+ * @brief Load interface data from container to point on host
+ *
+ * Loads coupled interface data using compile-time dispatch based on the point's
+ * template parameters (connection, interface, and boundary types).
+ *
+ * @ingroup CoupledInterfaceDataAccess
+ *
+ * @tparam IndexType Face index type
+ * @tparam ContainerType Coupled interfaces container type
+ * @tparam AccessorType Interface point (or face) type
+ *
+ * @param index Face index specifying the interface location
+ * @param container Coupled interfaces container holding interface data
+ * @param accessor Point or chunk_face object where loaded data will be stored
+ *
+ * @pre index refers to valid mesh face
+ * @pre container is properly initialized
+ * @pre point type matches supported interface combinations
+ *
+ * @note For host-side computations only. Use load_on_device for device code.
+ */
+template <typename IndexType, typename ContainerType, typename AccessorType,
+          typename std::enable_if_t<
+              ((specfem::data_access::is_face_index<IndexType>::value) &&
+               (specfem::data_access::is_nonconforming_interface<
+                   ContainerType>::value)),
+              int> = 0>
+inline void load_on_host(const IndexType &index, const ContainerType &container,
+                         AccessorType &accessor) {
+
+  static_assert(
+      specfem::data_access::CheckCompatibility<IndexType, ContainerType,
+                                               AccessorType>::value,
+      "Incompatible types in load_on_host");
+
+  using accessor_dispatch =
+      std::integral_constant<specfem::datatype::AccessorType,
+                             IndexType::accessor_type>;
+
+  container
+      .template get_interface_container<
+          AccessorType::interface_tag, AccessorType::boundary_tag,
+          AccessorType::connection_tag, AccessorType::flux_scheme_tag>()
+      .template impl_load<false>(accessor_dispatch(), index, accessor);
+}
+
+/**
+ * @brief Load interface data from container to face on device
+ *
+ * Loads coupled interface data using compile-time dispatch based on the face's
+ * template parameters (connection, interface, and boundary types).
+ *
+ * @ingroup CoupledInterfaceDataAccess
+ *
+ * @tparam IndexType Face index type
+ * @tparam ContainerType Coupled interfaces container type
+ * @tparam AccessorType Interface point (or face) type
+ *
+ * @param index Face index specifying the interface location
+ * @param container Coupled interfaces container holding interface data
+ * @param accessor Point or chunk_face object where loaded data will be stored
+ *
+ * @pre index refers to valid mesh face
+ * @pre container is properly initialized
+ * @pre point type matches supported interface combinations
+ *
+ * @note For device-side computations only. Use load_on_host for host code.
+ */
+template <typename IndexType, typename ContainerType, typename AccessorType,
+          typename std::enable_if_t<
+              ((specfem::data_access::is_face_index<IndexType>::value) &&
+               (specfem::data_access::is_nonconforming_interface<
+                   ContainerType>::value)),
+              int> = 0>
+KOKKOS_FORCEINLINE_FUNCTION void load_on_device(const IndexType &index,
+                                                const ContainerType &container,
+                                                AccessorType &accessor) {
+  static_assert(
+      specfem::data_access::CheckCompatibility<IndexType, ContainerType,
+                                               AccessorType>::value,
+      "Incompatible types in load_on_device");
+
+  using accessor_dispatch =
+      std::integral_constant<specfem::datatype::AccessorType,
+                             IndexType::accessor_type>;
+  container
+      .template get_interface_container<
+          AccessorType::interface_tag, AccessorType::boundary_tag,
+          AccessorType::connection_tag, AccessorType::flux_scheme_tag>()
+      .template impl_load<true>(accessor_dispatch(), index, accessor);
+
+  return;
+}
+} // namespace specfem::assembly