dG Dependecies: Mesh Modifiers and Subdivision

CMakeLists.txt

@@ -305,6 +305,8 @@ add_library(
         src/mesh/coupled_interfaces/coupled_interfaces.cpp
         src/mesh/tags/tags.cpp
         src/mesh/mesh.cpp
+        src/mesh/modifiers/modifiers.cpp
+        src/mesh/modifiers/apply2d.cpp
 )
 
 target_link_libraries(

docs/parameter_documentation/index.rst

@@ -14,3 +14,4 @@ The run-time behaviour of the code is controlled by YAML parameter files. The pa
     receivers
     run_setup
     databases
+    mesh_modifiers

docs/parameter_documentation/mesh_modifiers.rst

@@ -0,0 +1,72 @@
+Mesh Modifiers
+##############
+
+Mesh modifiers are rules to apply to a mesh loaded in from a database file.
+
+Parameter definitions
+=====================
+
+**Parameter name** : ``mesh-modifiers``
+---------------------------------------
+
+**default value** : None
+
+**possible values** : [YAML Node]
+
+**documentation** : Define databases section
+
+
+**Parameter name** : ``mesh-modifiers.subdivisions``
+****************************************************
+
+**default value**: None
+
+**possible values**: [YAML list]
+
+**documentation**: Collection of subdivision rules.
+
+
+**Parameter name** : ``mesh-modifiers.subdivisions.material``
+*************************************************************
+
+**default value**: None
+
+**possible values**: [int]
+
+**documentation**: The database index of the material (specified in the meshfem parfile) to be subdivided.
+
+
+**Parameter name** : ``mesh-modifiers.subdivisions.x``
+******************************************************
+
+**default value**: 1
+
+**possible values**: [int] (must be positive)
+
+**documentation**: The number of subdivisions along the x-axis. 1 keeps the axis unmodified.
+
+**Parameter name** : ``mesh-modifiers.subdivisions.z``
+******************************************************
+
+**default value**: 1
+
+**possible values**: [int] (must be positive)
+
+**documentation**: The number of subdivisions along the z-axis. 1 keeps the axis unmodified.
+
+
+
+
+
+.. admonition:: Example of mesh-modifiers section
+
+    .. code-block:: yaml
+
+        mesh-modifiers:
+            subdivisions:
+              - material: 1
+                x: 2
+                z: 2
+              - material: 2
+                x: 3
+                z: 3

include/mesh/mesh.hpp

@@ -27,6 +27,8 @@ template <> struct mesh<specfem::dimension::type::dim2> {
   constexpr static auto dimension =
       specfem::dimension::type::dim2; ///< Dimension
 
+  bool requires_coupled_interface_recalculation; ///< False by default; set to
+                                                 ///< true by mesh modifiers
   int npgeo; ///< Total number of spectral element control nodes
   int nspec; ///< Total number of spectral elements
   int nproc; ///< Total number of processors
@@ -93,7 +95,8 @@ template <> struct mesh<specfem::dimension::type::dim2> {
       : npgeo(npgeo), nspec(nspec), nproc(nproc), control_nodes(control_nodes),
         parameters(parameters), coupled_interfaces(coupled_interfaces),
         boundaries(boundaries), tags(tags), tangential_nodes(tangential_nodes),
-        axial_nodes(axial_nodes), materials(materials){};
+        axial_nodes(axial_nodes), materials(materials),
+        requires_coupled_interface_recalculation(false){};
   ///@}
 
   std::string print() const;

include/mesh/modifiers/modifiers.hpp

@@ -0,0 +1,55 @@
+#pragma once
+
+#include "enumerations/dimension.hpp"
+#include "mesh/mesh.hpp"
+
+#include <string>
+#include <unordered_map>
+
+namespace specfem {
+namespace mesh {
+
+// dimtuple<T, dim> gives tuple<T,...> (size = dim)
+// dimtuple<T, toset, set> gives tuple<T,...> (size = toset + sizeof(set))
+template <typename T, int toset, typename... set> struct dimtuple {
+  // using type = typename std::conditional<toset <= 0, std::tuple<set...>,
+  // typename dimtuple<T,toset-1,T,set...>::type>::type;
+  using unravelstruct = dimtuple<T, toset - 1, T, set...>;
+  using type = typename unravelstruct::type;
+  static inline std::string subdiv_str(type tup) {
+    return unravelstruct::subdiv_str(tup) + ((toset > 1) ? "x" : "") +
+           std::to_string(std::get<toset - 1>(tup));
+  }
+};
+template <typename T, typename... set> struct dimtuple<T, 0, set...> {
+  using type = std::tuple<set...>;
+  static inline std::string subdiv_str(type tup) { return ""; }
+};
+
+template <specfem::dimension::type DimensionType> class modifiers {
+public:
+  static constexpr int dim = specfem::dimension::dimension<DimensionType>::dim;
+  using subdiv_tuple = typename dimtuple<int, dim>::type;
+  modifiers() {}
+
+  //===== application =====
+  void apply(specfem::mesh::mesh<DimensionType> &mesh) const;
+
+  //===== display / debug / info =====
+  std::string to_string() const;
+  std::string subdivisions_to_string() const;
+
+  //===== setting modifiers =====
+  void set_subdivision(const int material, subdiv_tuple subdivisions);
+  //===== getting modifiers =====
+  specfem::mesh::modifiers<DimensionType>::subdiv_tuple
+  get_subdivision(const int material) const;
+
+private:
+  std::unordered_map<int, subdiv_tuple> subdivisions; ///< map
+                                                      ///< materialID ->
+                                                      ///< (subdivide_z,
+                                                      ///< subdivide_x)
+};
+} // namespace mesh
+} // namespace specfem

include/parameter_parser/mesh_modifiers.hpp

@@ -0,0 +1,32 @@
+#pragma once
+
+#include "enumerations/dimension.hpp"
+#include "mesh/modifiers/modifiers.hpp"
+
+#include "yaml-cpp/yaml.h"
+#include <memory>
+
+namespace specfem {
+namespace runtime_configuration {
+/**
+ * @brief class to read mesh modifier information
+ *
+ */
+class mesh_modifiers {
+public:
+  mesh_modifiers(const YAML::Node &Node) : mesh_modifiers_node(Node) {}
+
+  template <specfem::dimension::type DimensionType>
+  std::shared_ptr<specfem::mesh::modifiers<DimensionType> >
+  instantiate_mesh_modifiers();
+  template <specfem::dimension::type DimensionType>
+  void load_subdivisions(specfem::mesh::modifiers<DimensionType> &modifiers);
+
+private:
+  YAML::Node mesh_modifiers_node; /// Node that contains receiver information
+};
+
+} // namespace runtime_configuration
+} // namespace specfem
+
+#include "parameter_parser/mesh_modifiers.tpp"

include/parameter_parser/mesh_modifiers.tpp

@@ -0,0 +1,91 @@
+#pragma once
+#include "parameter_parser/mesh_modifiers.hpp"
+
+// #include <fstream>
+// #include <iostream>
+#include <ostream>
+
+template<specfem::dimension::type DimensionType>
+std::shared_ptr<specfem::mesh::modifiers<DimensionType>>
+specfem::runtime_configuration::mesh_modifiers::instantiate_mesh_modifiers() {
+  std::shared_ptr<specfem::mesh::modifiers<DimensionType>> modifiers =
+      std::make_shared<specfem::mesh::modifiers<DimensionType>>();
+
+  load_subdivisions(*modifiers);
+
+  return modifiers;
+}
+
+template<specfem::dimension::type DimensionType>
+void specfem::runtime_configuration::mesh_modifiers::load_subdivisions(
+    specfem::mesh::modifiers<DimensionType> &modifiers) {
+  if (const YAML::Node &subdivisions_node =
+          mesh_modifiers_node["subdivisions"]) {
+    if (!subdivisions_node.IsSequence()) {
+      std::ostringstream message;
+
+      message << "Error reading specfem mesh modifiers.\n";
+      message << "\"subdivisions\" must be a sequence.\n";
+
+      std::runtime_error(message.str());
+    }
+    for (YAML::Node subdiv_entry : subdivisions_node) {
+      int subx = 1;
+      int subz = 1;
+      int materialID;
+      if (const YAML::Node &materialentry = subdiv_entry["material"]) {
+        try {
+          materialID = materialentry.as<int>();
+        } catch (YAML::ParserException &e) {
+          std::ostringstream message;
+          message << "Error reading specfem mesh modifiers: subdivisions.\n";
+          message << "\"material\" must be an integer.\n";
+          message << e.what();
+          std::runtime_error(message.str());
+        }
+      } else {
+        std::ostringstream message;
+        message << "Error reading specfem mesh modifiers: subdivisions.\n";
+        message
+            << "\"material\" must be specified for each subdivision entry.\n";
+        std::runtime_error(message.str());
+      }
+      if (const YAML::Node &subdiv_z = subdiv_entry["z"]) {
+        try {
+          subz = subdiv_z.as<int>();
+          if (subz < 1) {
+            std::ostringstream message;
+            message << "Error reading specfem mesh modifiers: subdivisions.\n";
+            message << "\"z\" must be a positive integer.\n";
+            std::runtime_error(message.str());
+          }
+        } catch (YAML::ParserException &e) {
+          std::ostringstream message;
+          message << "Error reading specfem mesh modifiers: subdivisions.\n";
+          message << "\"z\" must be a positive integer.\n";
+          message << e.what();
+          std::runtime_error(message.str());
+        }
+      }
+      if (const YAML::Node &subdiv_x = subdiv_entry["x"]) {
+        try {
+          subx = subdiv_x.as<int>();
+          if (subx < 1) {
+            std::ostringstream message;
+            message << "Error reading specfem mesh modifiers: subdivisions.\n";
+            message << "\"x\" must be a positive integer.\n";
+            std::runtime_error(message.str());
+          }
+        } catch (YAML::ParserException &e) {
+          std::ostringstream message;
+          message << "Error reading specfem mesh modifiers: subdivisions.\n";
+          message << "\"x\" must be a positive integer.\n";
+          message << e.what();
+          std::runtime_error(message.str());
+        }
+      }
+      // minus 1 for zero-indexing.
+      modifiers.set_subdivision(materialID - 1, std::make_tuple(subx, subz));
+    }
+  }
+}

include/parameter_parser/setup.hpp

@@ -4,6 +4,7 @@
 #include "IO/reader.hpp"
 #include "database_configuration.hpp"
 #include "header.hpp"
+#include "mesh_modifiers.hpp"
 #include "parameter_parser/solver/interface.hpp"
 #include "quadrature.hpp"
 #include "receivers.hpp"
@@ -218,6 +219,17 @@ class setup {
     return this->solver->instantiate<NGLL>(dt, assembly, time_scheme, tasks);
   }
 
+  template <specfem::dimension::type DimensionType>
+  std::shared_ptr<specfem::mesh::modifiers<DimensionType> >
+  instantiate_mesh_modifiers() const {
+    if (this->mesh_modifiers) {
+      return this->mesh_modifiers->instantiate_mesh_modifiers<DimensionType>();
+    } else {
+      // default modifiers object; not nullptr
+      return std::make_shared<specfem::mesh::modifiers<DimensionType> >();
+    }
+  }
+
   int get_nsteps() const { return this->time_scheme->get_nsteps(); }
 
 private:
@@ -255,6 +267,8 @@ class setup {
       databases; ///< Get database filenames
   std::unique_ptr<specfem::runtime_configuration::solver::solver>
       solver; ///< Pointer to solver object
+  std::unique_ptr<specfem::runtime_configuration::mesh_modifiers>
+      mesh_modifiers; ///< Pointer to mesh modifiers container
 };
 } // namespace runtime_configuration
 } // namespace specfem

src/execute.cpp

@@ -42,7 +42,10 @@ void execute(
   //                   Read mesh and materials
   // --------------------------------------------------------------
   const auto quadrature = setup.instantiate_quadrature();
-  const auto mesh = specfem::IO::read_mesh(database_filename, mpi);
+  const auto mesh_modifiers =
+      setup.instantiate_mesh_modifiers<specfem::dimension::type::dim2>();
+  auto mesh = specfem::IO::read_mesh(database_filename, mpi);
+  mesh_modifiers->apply(mesh);
   // --------------------------------------------------------------
 
   // --------------------------------------------------------------