Iterators and Access functions for edges

CMakeLists.txt

@@ -355,6 +355,7 @@ add_library(
 target_link_libraries(
         enumerations
         ${BOOST_LIBS}
+        Kokkos::kokkos
 )
 
 

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

@@ -217,4 +217,47 @@ load_after_simd_dispatch(const std::true_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_edge<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_index();
+        const auto iedge = point_index.get_policy_index();
+        const int iglob = field.template get_iglob<on_device>(
+            point_index.ispec, point_index.iz, point_index.ix, MediumTag);
+        for (int icomp = 0; icomp < ncomponents; ++icomp) {
+          (specfem::assembly::fields_impl::base_load_accessor<
+               on_device, AccessorTypes::data_class>(
+               iglob, icomp, current_field,
+               accessors(iedge, point_index.ipoint, icomp)),
+           ...);
+        }
+      });
+
+  return;
+}
+
 } // namespace specfem::assembly::simulation_field_impl

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

@@ -36,30 +36,9 @@ KOKKOS_FORCEINLINE_FUNCTION void load_on_device(const IndexType &index,
 template <typename IndexType, typename ContainerType, typename... AccessorTypes,
           typename std::enable_if_t<
               ((specfem::data_access::is_index_type<IndexType>::value) &&
-               (specfem::data_access::is_point<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;
-
-  // Check that all accessors have the same medium tag
-  fields_impl::check_accessor_compatibility<AccessorTypes...>();
-
-  using simd_accessor_type =
-      std::integral_constant<bool, IndexType::using_simd>;
-  simulation_field_impl::load_after_simd_dispatch<true>(
-      simd_accessor_type(), index, field, accessors...);
-  return;
-}
-
-template <typename IndexType, typename ContainerType, typename... AccessorTypes,
-          typename std::enable_if_t<
-              ((specfem::data_access::is_chunk_element<IndexType>::value) &&
-               (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_field<AccessorTypes>::value && ...)),
               int> = 0>
 KOKKOS_FORCEINLINE_FUNCTION void load_on_device(const IndexType &index,

core/specfem/chunk_edge/index.hpp

@@ -0,0 +1,52 @@
+#pragma once
+
+#include "enumerations/interface.hpp"
+#include "execution/chunked_edge_iterator.hpp"
+#include "specfem/data_access.hpp"
+
+namespace specfem::chunk_edge {
+
+template <specfem::dimension::type DimensionTag, typename ViewType,
+          typename TeamMemberType>
+class Index
+    : public specfem::execution::ChunkEdgeIndex<DimensionTag, TeamMemberType,
+                                                ViewType>,
+      public specfem::data_access::Accessor<
+          specfem::data_access::AccessorType::chunk_edge,
+          specfem::data_access::DataClassType::index, DimensionTag, false> {
+private:
+  using base_type =
+      specfem::execution::ChunkEdgeIndex<DimensionTag, TeamMemberType,
+                                         ViewType>;
+
+public:
+  /// @brief Iterator type for traversing elements in the chunk
+  using iterator_type = typename base_type::iterator_type;
+
+  /**
+   * @brief Construct index from existing chunk element index base.
+   *
+   * Creates a chunk element index by wrapping an existing ChunkElementIndex
+   * object. This constructor is useful when you already have a base chunk
+   * element index and want to add the data access layer functionality.
+   *
+   * @param base The base chunk element index to wrap
+   *
+   * @code{.cpp}
+   * // Create base chunk element index
+   * auto base_index = specfem::execution::ChunkElementIndex<...>(...);
+   *
+   * // Wrap it with data access layer
+   * IndexType chunk_index(base_index);
+   * @endcode
+   */
+  KOKKOS_INLINE_FUNCTION
+  Index(const base_type &base) : base_type(base) {}
+
+  KOKKOS_INLINE_FUNCTION
+  Index(const ViewType indices, const int &ngllz, const int &ngllx,
+        const TeamMemberType &kokkos_index)
+      : base_type(indices, ngllz, ngllx, kokkos_index) {}
+};
+
+} // namespace specfem::chunk_edge

core/specfem/point.hpp

@@ -7,6 +7,7 @@
 #include "point/displacement.hpp"
 #include "point/field_derivatives.hpp"
 #include "point/index.hpp"
+#include "point/interface_index.hpp"
 #include "point/jacobian_matrix.hpp"
 #include "point/kernels.hpp"
 #include "point/mapped_index.hpp"

core/specfem/point/interface_index.hpp

@@ -0,0 +1,15 @@
+
+#pragma once
+
+#include "enumerations/interface.hpp"
+
+namespace specfem::point {
+
+template <specfem::dimension::type DimensionTag> class interface_index {
+public:
+  specfem::point::index<DimensionTag> self_index;
+  specfem::point::index<DimensionTag> coupled_index;
+  int ipoint;
+};
+
+} // namespace specfem::point

core/specfem/point/mapped_index.hpp

@@ -29,18 +29,10 @@ template <specfem::dimension::type DimensionTag, bool UseSIMD>
 struct mapped_index : public index<DimensionTag, UseSIMD> {
 private:
   using base_type = index<DimensionTag, UseSIMD>;
-  using accessor_type = specfem::data_access::Accessor<
-      specfem::data_access::AccessorType::point,
-      specfem::data_access::DataClassType::mapped_index, DimensionTag,
-      UseSIMD>; ///< Accessor type for
-                ///< mapped index
 
 public:
   int imap; ///< Index of the mapped element
 
-  constexpr static auto data_class =
-      accessor_type::data_class; ///< Data class of the mapped index
-
   /**
    * @brief Constructor for the mapped index
    *

include/enumerations/mesh_entities.hpp

@@ -1,5 +1,6 @@
 #pragma once
 
+#include <Kokkos_Core.hpp>
 #include <algorithm>
 #include <list>
 
@@ -124,4 +125,19 @@ std::list<type> edges_of_corner(const type &corner);
  */
 std::list<type> corners_of_edge(const type &edge);
 
+struct edge {
+  specfem::mesh_entity::type edge_type;
+  int ispec;
+  bool reverse_orientation;
+
+  KOKKOS_INLINE_FUNCTION
+  edge(const int ispec, const specfem::mesh_entity::type edge_type,
+       const bool reverse_orientation = false)
+      : edge_type(edge_type), ispec(ispec),
+        reverse_orientation(reverse_orientation) {}
+
+  KOKKOS_INLINE_FUNCTION
+  edge() = default;
+};
+
 } // namespace specfem::mesh_entity

include/execution/chunked_domain_iterator.hpp

@@ -358,6 +358,15 @@ class ChunkedDomainIterator : public TeamPolicy<ParallelConfig> {
                           ///< specfem::execution::for_each_level
   using execution_space =
       typename base_type::execution_space; ///< Execution space type.
+  using base_index_type = specfem::point::index<
+      ParallelConfig::dimension,
+      ParallelConfig::simd::using_simd>; ///< Index type
+                                         ///< to be used
+                                         ///< when calling
+                                         ///< @ref
+                                         ///< specfem::execution::for_all
+                                         ///< with this
+                                         ///< iterator.
 
   /**
    * @brief Construct a new Chunked Domain Iterator object