Forgot adding files that I added for convenience

include/IO/mesh/impl/fortran/dim3/interface.hpp

@@ -0,0 +1,92 @@
+#pragma once
+
+#include "IO/fortranio/interface.hpp"
+#include "mesh/coordinates/coordinates.hpp"
+#include "mesh/mapping/mapping.hpp"
+#include "mesh/mesh.hpp"
+#include "mesh/parameters/parameters.hpp"
+#include "mesh/partial_derivatives/partial_derivatives.hpp"
+#include "specfem_mpi/interface.hpp"
+
+namespace specfem {
+namespace IO {
+namespace mesh {
+namespace impl {
+namespace fortran {
+namespace dim3 {
+
+/*
+ * @brief Read paramters from 3D mesh database
+ *
+ * @param stream Input stream
+ * @param mpi MPI object
+ * @return specfem::mesh::parameters<specfem::dimension::type::dim2> Mesh
+ * parameters
+ */
+specfem::mesh::parameters<specfem::dimension::type::dim3>
+read_mesh_parameters(std::ifstream &stream, const specfem::MPI::MPI *mpi);
+
+/*
+ * @brief Read mapping from 3D mesh database
+ *
+ * @param stream Input stream
+ * @param mapping Mapping object
+ * @param mpi MPI object
+ */
+void read_ibool(std::ifstream &stream,
+                specfem::mesh::mapping<specfem::dimension::type::dim3> &mapping,
+                const specfem::MPI::MPI *mpi);
+
+/*
+ * @brief Read coordinates from 3D mesh database
+ *
+ * @param stream Input stream
+ * @param coordinates Coordinates object
+ * @param mpi MPI object
+ */
+void read_xyz(
+    std::ifstream &stream,
+    specfem::mesh::coordinates<specfem::dimension::type::dim3> &coordinates,
+    const specfem::MPI::MPI *mpi);
+
+/*
+ * @brief Read Jacobian from 3D mesh database
+ *
+ * @param stream Input stream
+ * @param coordinates Coordinates object
+ * @param mpi MPI object
+ */
+void read_partial_derivatives(
+    std::ifstream &stream,
+    specfem::mesh::partial_derivatives<specfem::dimension::type::dim3>
+        &partial_derivatives,
+    const specfem::MPI::MPI *mpi);
+
+template <typename T> using View1D = Kokkos::View<T *, Kokkos::HostSpace>;
+
+template <typename T> using View4D = Kokkos::View<T ****, Kokkos::HostSpace>;
+
+template <typename T> using View5D = Kokkos::View<T *****, Kokkos::HostSpace>;
+
+template <typename T> void read_array(std::ifstream &stream, View1D<T> &array);
+
+template <typename T> void read_array(std::ifstream &stream, View4D<T> &array);
+
+template <typename T> void read_array(std::ifstream &stream, View5D<T> &array);
+
+// Read index array will subtract 1 from each value when reading to account for
+// Fortran 1-based indexing
+template <typename T>
+void read_index_array(std::ifstream &stream, View1D<T> &array);
+
+template <typename T>
+void read_index_array(std::ifstream &stream, View4D<T> &array);
+
+} // namespace dim3
+} // namespace fortran
+} // namespace impl
+} // namespace mesh
+} // namespace IO
+} // namespace specfem
+
+#include "IO/mesh/impl/fortran/dim3/utilities.hpp"

include/IO/mesh/impl/fortran/dim3/utilities.hpp

@@ -0,0 +1,173 @@
+#include "IO/fortranio/interface.hpp"
+#include "IO/mesh/impl/fortran/dim3/interface.hpp"
+#include "enumerations/dimension.hpp"
+#include "mesh/parameters/parameters.hpp"
+#include "specfem_mpi/interface.hpp"
+#include <Kokkos_Core.hpp>
+
+template <typename T> using View1D = Kokkos::View<T *, Kokkos::HostSpace>;
+
+template <typename T> using View4D = Kokkos::View<T ****, Kokkos::HostSpace>;
+
+template <typename T> using View5D = Kokkos::View<T *****, Kokkos::HostSpace>;
+
+template <typename T>
+void specfem::IO::mesh::impl::fortran::dim3::read_array(std::ifstream &stream,
+                                                        View1D<T> &array) {
+  const int n = array.extent(0);
+  std::vector<T> dummy_T(n, -999999);
+
+  try {
+    // Read into dummy vector
+    specfem::IO::fortran_read_line(stream, &dummy_T);
+    // Assign to KokkosView
+    for (int i = 0; i < n; i++) {
+      array(i) = dummy_T[i];
+    }
+  } catch (std::runtime_error &e) {
+    std::ostringstream error_message;
+    error_message << "Error reading array from database file:\n"
+                  << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")";
+    throw std::runtime_error(error_message.str());
+  }
+}
+
+template <typename T>
+void specfem::IO::mesh::impl::fortran::dim3::read_array(std::ifstream &stream,
+                                                        View4D<T> &array) {
+  const int nspec = array.extent(0);
+  const int ngllx = array.extent(1);
+  const int nglly = array.extent(2);
+  const int ngllz = array.extent(3);
+
+  std::vector<T> dummy_T(ngllx * nglly * ngllz, -999999);
+
+  try {
+    for (int ispec = 0; ispec < nspec; ispec++) {
+      // Read into dummy vector for each ispec
+      specfem::IO::fortran_read_line(stream, &dummy_T);
+
+      // Assign to KokkosView
+      int counter = 0;
+      for (int igllx = 0; igllx < ngllx; igllx++) {
+        for (int iglly = 0; iglly < nglly; iglly++) {
+          for (int igllz = 0; igllz < ngllz; igllz++) {
+            array(ispec, igllx, iglly, igllz) = dummy_T[counter];
+            counter++;
+          }
+        }
+      }
+      if ((ispec == 0) || (ispec == nspec - 1)) {
+        std::cout << "ispec=" << ispec
+                  << " array(ispec,0,0,0)=" << array(ispec, 0, 0, 0)
+                  << std::endl;
+      }
+    }
+
+  } catch (std::runtime_error &e) {
+    std::ostringstream error_message;
+    error_message << "Error reading array from database file:\n"
+                  << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")";
+    throw std::runtime_error(error_message.str());
+  }
+}
+
+template <typename T>
+void specfem::IO::mesh::impl::fortran::dim3::read_array(std::ifstream &stream,
+                                                        View5D<T> &array) {
+  const int nspec = array.extent(0);
+  const int ngllx = array.extent(1);
+  const int nglly = array.extent(2);
+  const int ngllz = array.extent(3);
+  const int ncomp = array.extent(4);
+
+  std::vector<T> dummy_T(ngllx * nglly * ngllz * ncomp, -9999999);
+
+  try {
+    for (int ispec = 0; ispec < nspec; ispec++) {
+      // Read into dummy vector for each ispec
+      specfem::IO::fortran_read_line(stream, &dummy_T);
+
+      // Assign to KokkosView
+      int counter = 0;
+      for (int igllx = 0; igllx < ngllx; igllx++) {
+        for (int iglly = 0; iglly < nglly; iglly++) {
+          for (int igllz = 0; igllz < ngllz; igllz++) {
+            for (int icomp = 0; icomp < ncomp; icomp++) {
+              array(ispec, igllx, iglly, igllz, icomp) = dummy_T[counter];
+              counter++;
+            }
+          }
+        }
+      }
+    }
+  } catch (std::runtime_error &e) {
+    std::ostringstream error_message;
+    error_message << "Error reading array from database file:\n"
+                  << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")";
+    throw std::runtime_error(error_message.str());
+  }
+}
+
+template <typename T>
+void specfem::IO::mesh::impl::fortran::dim3::read_index_array(
+    std::ifstream &stream, View1D<T> &array) {
+  const int n = array.extent(0);
+  std::vector<T> dummy_T(n, -999999);
+
+  try {
+    // Read into dummy vector
+    specfem::IO::fortran_read_line(stream, &dummy_T);
+
+    // Assign to KokkosView
+    for (int i = 0; i < n; i++) {
+      array(i) = dummy_T[i] - 1;
+    }
+
+  } catch (std::runtime_error &e) {
+    std::ostringstream error_message;
+    error_message << "Error reading 1D index_array from database file:\n"
+                  << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")";
+    throw std::runtime_error(error_message.str());
+  }
+}
+
+template <typename T>
+void specfem::IO::mesh::impl::fortran::dim3::read_index_array(
+    std::ifstream &stream, View4D<T> &array) {
+  const int nspec = array.extent(0);
+  const int ngllx = array.extent(1);
+  const int nglly = array.extent(2);
+  const int ngllz = array.extent(3);
+
+  std::vector<T> dummy_T(ngllx * nglly * ngllz, -9999.0);
+
+  try {
+    for (int ispec = 0; ispec < nspec; ispec++) {
+      // Read into dummy vector for each ispec
+      specfem::IO::fortran_read_line(stream, &dummy_T);
+
+      // Assign to KokkosView
+      int counter = 0;
+      for (int igllx = 0; igllx < ngllx; igllx++) {
+        for (int iglly = 0; iglly < nglly; iglly++) {
+          for (int igllz = 0; igllz < ngllz; igllz++) {
+            array(ispec, igllx, iglly, igllz) = dummy_T[counter] - 1;
+            counter++;
+          }
+        }
+      }
+      if ((ispec == 0) || (ispec == nspec - 1)) {
+        std::cout << "ispec=" << ispec
+                  << " array(ispec,0,0,0)=" << array(ispec, 0, 0, 0)
+                  << std::endl;
+      }
+    }
+
+  } catch (std::runtime_error &e) {
+    std::ostringstream error_message;
+    error_message << "Error reading array from database file:\n"
+                  << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")";
+    throw std::runtime_error(error_message.str());
+  }
+}

include/mesh/partial_derivatives/partial_derivatives.hpp

@@ -0,0 +1,56 @@
+#pragma once
+#include "enumerations/dimension.hpp"
+#include "specfem_setup.hpp"
+#include <Kokkos_Core.hpp>
+
+namespace specfem {
+namespace mesh {
+
+template <specfem::dimension::type DimensionType> struct partial_derivatives;
+
+template <> struct partial_derivatives<specfem::dimension::type::dim3> {
+  constexpr static auto dimension = specfem::dimension::type::dim3;
+
+  using LocalView = Kokkos::View<type_real ****, Kokkos::HostSpace>;
+
+  // Parameters
+  int nspec;
+  int ngllx;
+  int nglly;
+  int ngllz;
+
+  LocalView xix;
+  LocalView xiy;
+  LocalView xiz;
+  LocalView etax;
+  LocalView etay;
+  LocalView etaz;
+  LocalView gammax;
+  LocalView gammay;
+  LocalView gammaz;
+  LocalView jacobian;
+
+  // Constructors
+  partial_derivatives(){}; // Default constructor
+
+  // Constructor to initialize the coordinates
+  partial_derivatives(int nspec, int ngllx, int nglly, int ngllz)
+      : nspec(nspec), ngllx(ngllx), nglly(nglly), ngllz(ngllz),
+        xix("xix", nspec, ngllx, nglly, ngllz),
+        xiy("xiy", nspec, ngllx, nglly, ngllz),
+        xiz("xiz", nspec, ngllx, nglly, ngllz),
+        etax("etax", nspec, ngllx, nglly, ngllz),
+        etay("etay", nspec, ngllx, nglly, ngllz),
+        etaz("etaz", nspec, ngllx, nglly, ngllz),
+        gammax("gammax", nspec, ngllx, nglly, ngllz),
+        gammay("gammay", nspec, ngllx, nglly, ngllz),
+        gammaz("gammaz", nspec, ngllx, nglly, ngllz),
+        jacobian("jacobian", nspec, ngllx, nglly, ngllz){};
+
+  void print() const;
+
+  void print(int ispec, int igllx, int iglly, int igllz) const;
+};
+
+} // namespace mesh
+} // namespace specfem