diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6aa1b867..b7b88919 100644
--- a/CMakeLists.txt
+++ b/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(
diff --git a/docs/parameter_documentation/index.rst b/docs/parameter_documentation/index.rst
index ecac2253..515b68bf 100644
--- a/docs/parameter_documentation/index.rst
+++ b/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
diff --git a/docs/parameter_documentation/mesh_modifiers.rst b/docs/parameter_documentation/mesh_modifiers.rst
new file mode 100644
index 00000000..5a2c5475
--- /dev/null
+++ b/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
diff --git a/include/mesh/mesh.hpp b/include/mesh/mesh.hpp
index cd64b19d..9f3f910e 100644
--- a/include/mesh/mesh.hpp
+++ b/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;
diff --git a/include/mesh/modifiers/modifiers.hpp b/include/mesh/modifiers/modifiers.hpp
new file mode 100644
index 00000000..bd7ba8a3
--- /dev/null
+++ b/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
diff --git a/include/parameter_parser/mesh_modifiers.hpp b/include/parameter_parser/mesh_modifiers.hpp
new file mode 100644
index 00000000..a636ba52
--- /dev/null
+++ b/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"
diff --git a/include/parameter_parser/mesh_modifiers.tpp b/include/parameter_parser/mesh_modifiers.tpp
new file mode 100644
index 00000000..9c5d73d6
--- /dev/null
+++ b/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));
+    }
+  }
+}
diff --git a/include/parameter_parser/setup.hpp b/include/parameter_parser/setup.hpp
index a82f42e0..8528cdbf 100644
--- a/include/parameter_parser/setup.hpp
+++ b/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
diff --git a/src/execute.cpp b/src/execute.cpp
index 138d5190..6c1a4b3a 100644
--- a/src/execute.cpp
+++ b/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);
   // --------------------------------------------------------------
 
   // --------------------------------------------------------------
diff --git a/src/mesh/modifiers/apply2d.cpp b/src/mesh/modifiers/apply2d.cpp
new file mode 100644
index 00000000..532efea7
--- /dev/null
+++ b/src/mesh/modifiers/apply2d.cpp
@@ -0,0 +1,301 @@
+#include "jacobian/shape_functions.hpp"
+#include "kokkos_abstractions.h"
+#include "mesh/boundaries/acoustic_free_surface.hpp"
+#include "mesh/control_nodes/control_nodes.hpp"
+#include "mesh/elements/axial_elements.hpp"
+#include "mesh/materials/materials.hpp"
+#include "mesh/modifiers/modifiers.hpp"
+#include <stdexcept>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+/**
+ * @brief Subdivides the boundary by adding to (index_mapping,type) the
+ * sub-elements along the given edge in the given ispec_old.
+ *
+ * @param ispec_matspec - reference to the original ispec->material mapping to
+ * reference
+ * @param modifiers - reference to the modifiers object (to access subx,subz)
+ * @param index_mapping - the index_mapping to append to
+ * @param type - the edge type array to append to
+ * @param ispec_old - ispec of the original element
+ * @param edge_type - edge of the original element
+ * @param ispec_old_to_new_offset - the first subelement ispec, used to compute
+ * isub -> ispec_new.
+ */
+static void subdivide_inherit_edgecond(
+    const specfem::kokkos::HostView1d<
+        specfem::mesh::materials::material_specification> &ispec_matspec,
+    const specfem::mesh::modifiers<specfem::dimension::type::dim2> &modifiers,
+    std::vector<int> &index_mapping,
+    std::vector<specfem::enums::boundaries::type> &type, const int ispec_old,
+    const specfem::enums::boundaries::type edge_type,
+    const int ispec_old_to_new_offset) {
+  int subz, subx;
+  const auto mat = ispec_matspec(ispec_old);
+  std::tie(subx, subz) = modifiers.get_subdivision(mat.database_index);
+
+  // number of subdivided elements to set
+  const int nsub = (edge_type == specfem::enums::boundaries::type::TOP ||
+                    edge_type == specfem::enums::boundaries::type::BOTTOM)
+                       ? subx
+                       : subz;
+  for (int isub = 0; isub < nsub; isub++) {
+    int isubz, isubx;
+    switch (edge_type) {
+    case specfem::enums::boundaries::type::RIGHT:
+      isubz = isub;
+      isubx = subx - 1;
+      break;
+    case specfem::enums::boundaries::type::TOP:
+      isubz = subz - 1;
+      isubx = isub;
+      break;
+    case specfem::enums::boundaries::type::LEFT:
+      isubz = subz - 1;
+      isubx = isub;
+      break;
+    case specfem::enums::boundaries::type::BOTTOM:
+      isubz = subz - 1;
+      isubx = isub;
+      break;
+    default:
+      throw std::runtime_error("Edge type NONE found in specfem::mesh::mesh.");
+      break;
+    }
+    const int ispec_new = ispec_old_to_new_offset + (subx * isubz + isubx);
+    index_mapping.push_back(ispec_new);
+    type.push_back(edge_type);
+  }
+}
+
+static void subdivide(
+    specfem::mesh::mesh<specfem::dimension::type::dim2> &mesh,
+    const specfem::mesh::modifiers<specfem::dimension::type::dim2> &modifiers) {
+  constexpr auto DimensionType = specfem::dimension::type::dim2;
+  const int nspec_old = mesh.nspec;
+  const int ngnod = mesh.control_nodes.ngnod;
+  const int ngnod_new = 9;
+  const auto material_index_mapping_old = mesh.materials.material_index_mapping;
+
+  /* we will expand coords, but we don't know to what size yet; use dynamic
+   * resize of vector.
+   */
+  std::vector<std::pair<type_real, type_real> > coord_new(mesh.npgeo);
+  // transfer coords to vector
+  for (int i = 0; i < mesh.npgeo; i++) {
+    coord_new[i] = std::make_pair(mesh.control_nodes.coord(0, i),
+                                  mesh.control_nodes.coord(1, i));
+  }
+
+  // this is how we map from old ispec to new ispec
+  std::vector<int> ispec_old_to_new_offsets(nspec_old);
+  // const auto ispec_old_to_new = [&](const int ispec, const int isubz,
+  //                                   const int isubx) {
+  //   const auto mat = material_index_mapping_old(ispec);
+  //   return ispec_old_to_new_offsets[ispec] +
+  //          std::get<0>(modifiers.get_subdivision(mat.database_index)) * isubz
+  //          + isubx;
+  // };
+
+  int nspec_new = 0; // update as we iterate
+  int subz, subx;
+  std::vector<int> knods_new;
+  std::vector<specfem::mesh::materials::material_specification>
+      material_index_mapping_new;
+  for (int ispec = 0; ispec < nspec_old; ispec++) {
+    const auto mat = material_index_mapping_old(ispec);
+    std::tie(subx, subz) = modifiers.get_subdivision(mat.database_index);
+
+    // form a matrix of node indices that will be used, making new ones where
+    // necessary.
+    const int ispec_ngnod_z = 2 * subz + 1;
+    const int ispec_ngnod_x = 2 * subx + 1;
+    const int total_knods = ispec_ngnod_x * ispec_ngnod_z;
+    std::vector<int> ispec_knod(total_knods);
+    for (int iz = 0; iz < ispec_ngnod_z; iz++) {
+      for (int ix = 0; ix < ispec_ngnod_x; ix++) {
+        const int ispec_inod = ispec_ngnod_x * iz + ix;
+        if (iz == 0) {
+          if (ix == 0) {
+            ispec_knod[ispec_inod] = mesh.control_nodes.knods(0, ispec);
+            continue;
+          } else if (ix == ispec_ngnod_x - 1) {
+            ispec_knod[ispec_inod] = mesh.control_nodes.knods(1, ispec);
+            continue;
+          } else if (ix == subx && ngnod == 9) {
+            ispec_knod[ispec_inod] = mesh.control_nodes.knods(4, ispec);
+            continue;
+          }
+        } else if (iz == ispec_ngnod_z - 1) {
+          if (ix == 0) {
+            ispec_knod[ispec_inod] = mesh.control_nodes.knods(3, ispec);
+            continue;
+          } else if (ix == ispec_ngnod_x - 1) {
+            ispec_knod[ispec_inod] = mesh.control_nodes.knods(2, ispec);
+            continue;
+          } else if (ix == subx && ngnod == 9) {
+            ispec_knod[ispec_inod] = mesh.control_nodes.knods(6, ispec);
+            continue;
+          }
+        } else if (iz == subz && ngnod == 9) {
+          if (ix == 0) {
+            ispec_knod[ispec_inod] = mesh.control_nodes.knods(7, ispec);
+            continue;
+          } else if (ix == ispec_ngnod_x - 1) {
+            ispec_knod[ispec_inod] = mesh.control_nodes.knods(5, ispec);
+            continue;
+          } else if (ix == subx && ngnod == 9) {
+            ispec_knod[ispec_inod] = mesh.control_nodes.knods(8, ispec);
+            continue;
+          }
+        }
+        // nod does not already exist. Make it exist
+        auto shape_functions = specfem::jacobian::define_shape_functions(
+            ix / ((type_real)subx) - 1, iz / ((type_real)subz) - 1, ngnod);
+
+        double xcor = 0.0;
+        double zcor = 0.0;
+
+        for (int in = 0; in < ngnod; in++) {
+          const int knod = mesh.control_nodes.knods(in, ispec);
+          xcor += mesh.control_nodes.coord(0, knod) * shape_functions[in];
+          zcor += mesh.control_nodes.coord(1, knod) * shape_functions[in];
+        }
+        ispec_knod[ispec_inod] = coord_new.size();
+        coord_new.push_back(std::make_pair(xcor, zcor));
+      }
+    }
+    ispec_old_to_new_offsets[ispec] = nspec_new;
+    for (int isubz = 0; isubz < subz; isubz++) {
+      for (int isubx = 0; isubx < subx; isubx++) {
+        // new element (based on new ngnod=9)
+        const int inod = ispec_ngnod_x * (isubz * 2) + isubx * 2;
+        knods_new.push_back(ispec_knod[inod]);
+        knods_new.push_back(ispec_knod[inod + 2]);
+        knods_new.push_back(ispec_knod[inod + 2 + ispec_ngnod_x * 2]);
+        knods_new.push_back(ispec_knod[inod + ispec_ngnod_x * 2]);
+        knods_new.push_back(ispec_knod[inod + 1]);
+        knods_new.push_back(ispec_knod[inod + 2 + ispec_ngnod_x]);
+        knods_new.push_back(ispec_knod[inod + 1 + ispec_ngnod_x * 2]);
+        knods_new.push_back(ispec_knod[inod + ispec_ngnod_x]);
+        knods_new.push_back(ispec_knod[inod + 1 + ispec_ngnod_x]);
+        material_index_mapping_new.push_back(mat);
+        nspec_new++;
+      }
+    }
+  }
+
+  mesh.nspec = nspec_new;
+
+  // skip coupled interface subdivision, since we don't need it.
+  // instead force it at assembly time.
+  mesh.requires_coupled_interface_recalculation = true;
+
+  // finalize control nodes struct and materials struct
+  const int npgeo_new = coord_new.size();
+  mesh.npgeo = npgeo_new;
+  mesh.materials.material_index_mapping = specfem::kokkos::HostView1d<
+      specfem::mesh::materials::material_specification>(
+      "specfem::mesh::material_index_mapping", nspec_new);
+  mesh.control_nodes =
+      specfem::mesh::control_nodes<specfem::dimension::type::dim2>(
+          2, nspec_new, ngnod_new, coord_new.size());
+
+  for (int igeo = 0; igeo < npgeo_new; igeo++) {
+    std::tie(mesh.control_nodes.coord(0, igeo),
+             mesh.control_nodes.coord(1, igeo)) = coord_new[igeo];
+  }
+  for (int ispec = 0; ispec < nspec_new; ispec++) {
+    for (int inod = 0; inod < ngnod_new; inod++) {
+      mesh.control_nodes.knods(inod, ispec) =
+          knods_new[ispec * ngnod_new + inod];
+    }
+    mesh.materials.material_index_mapping(ispec) =
+        material_index_mapping_new[ispec];
+  }
+
+  // boundaries
+  { // absorbing
+    std::vector<int> index_mapping;
+    std::vector<specfem::enums::boundaries::type> type;
+    auto &bdry = mesh.boundaries.absorbing_boundary;
+    for (int ielem = 0; ielem < bdry.nelements; ielem++) {
+      const int ispec = bdry.index_mapping(ielem);
+      const specfem::enums::boundaries::type edge = bdry.type(ielem);
+      subdivide_inherit_edgecond(material_index_mapping_old, modifiers,
+                                 index_mapping, type, ispec, edge,
+                                 ispec_old_to_new_offsets[ispec]);
+    }
+    bdry =
+        specfem::mesh::absorbing_boundary<DimensionType>(index_mapping.size());
+    for (int ielem = 0; ielem < bdry.nelements; ielem++) {
+      bdry.index_mapping(ielem) = index_mapping[ielem];
+      bdry.type(ielem) = type[ielem];
+    }
+  }
+  { // acoustic free surface
+    std::vector<int> index_mapping;
+    std::vector<specfem::enums::boundaries::type> type;
+    auto &bdry = mesh.boundaries.acoustic_free_surface;
+    for (int ielem = 0; ielem < bdry.nelem_acoustic_surface; ielem++) {
+      const int ispec = bdry.index_mapping(ielem);
+      const specfem::enums::boundaries::type edge = bdry.type(ielem);
+      subdivide_inherit_edgecond(material_index_mapping_old, modifiers,
+                                 index_mapping, type, ispec, edge,
+                                 ispec_old_to_new_offsets[ispec]);
+    }
+    bdry = specfem::mesh::acoustic_free_surface<DimensionType>(
+        index_mapping.size());
+    for (int ielem = 0; ielem < bdry.nelem_acoustic_surface; ielem++) {
+      bdry.index_mapping(ielem) = index_mapping[ielem];
+      bdry.type(ielem) = type[ielem];
+    }
+  }
+
+  { // axial nodes
+    auto &axi = mesh.axial_nodes;
+    auto axiflag = axi.is_on_the_axis;
+    axi = specfem::mesh::elements::axial_elements<DimensionType>(nspec_new);
+    for (int ispec = 0; ispec < nspec_old; ispec++) {
+      if (!axiflag(ispec)) {
+        continue;
+      }
+      auto mat = material_index_mapping_old(ispec);
+      std::tie(subx, subz) = modifiers.get_subdivision(mat.database_index);
+      if (subz != 1 || subx != 1) {
+        throw std::runtime_error(
+            "Subdividing elements on the axis not yet supported");
+      }
+      axi.is_on_the_axis(ispec_old_to_new_offsets[ispec]) = true;
+    }
+  }
+
+  // parameters
+  mesh.parameters.ngnod = ngnod_new;
+  mesh.parameters.nspec = nspec_new;
+  mesh.parameters.nelemabs = mesh.boundaries.absorbing_boundary.nelements;
+  mesh.parameters.nelem_acoustic_surface =
+      mesh.boundaries.acoustic_free_surface.nelem_acoustic_surface;
+  // mesh.parameters.nelem_acforcing = 0;
+  mesh.parameters.num_fluid_solid_edges =
+      mesh.coupled_interfaces.elastic_acoustic.num_interfaces;
+  mesh.parameters.num_fluid_poro_edges =
+      mesh.coupled_interfaces.acoustic_poroelastic.num_interfaces;
+  mesh.parameters.num_solid_poro_edges =
+      mesh.coupled_interfaces.elastic_poroelastic.num_interfaces;
+  // mesh.parameters.nelem_on_the_axis = 0;
+  // mesh.parameters.nnodes_tangential_curve = 0;
+
+  mesh.tags = specfem::mesh::tags<specfem::dimension::type::dim2>(
+      mesh.materials, mesh.boundaries);
+}
+
+template <>
+void specfem::mesh::modifiers<specfem::dimension::type::dim2>::apply(
+    specfem::mesh::mesh<specfem::dimension::type::dim2> &mesh) const {
+  subdivide(mesh, *this);
+
+  // forcing boundary skipped
+}
diff --git a/src/mesh/modifiers/modifiers.cpp b/src/mesh/modifiers/modifiers.cpp
new file mode 100644
index 00000000..24601887
--- /dev/null
+++ b/src/mesh/modifiers/modifiers.cpp
@@ -0,0 +1,52 @@
+#include "mesh/modifiers/modifiers.hpp"
+#include <cstdio>
+#include <string>
+
+// apply() handled in apply*.cpp in this directory
+
+//===== display / debug / info =====
+template <specfem::dimension::type DimensionType>
+std::string
+specfem::mesh::modifiers<DimensionType>::subdivisions_to_string() const {
+  std::string repr =
+      "subdivisions (set: " + std::to_string(subdivisions.size()) + "):";
+#define BUFSIZE 50
+  char buf[BUFSIZE];
+  for (const auto &[matID, subs] : subdivisions) {
+    repr += "\n  - material %d: " + dimtuple<int, dim>::subdiv_str(subs) +
+            " cell subdivision";
+  }
+#undef BUFSIZE
+  return repr;
+}
+
+template <specfem::dimension::type DimensionType>
+std::string specfem::mesh::modifiers<DimensionType>::to_string() const {
+  std::string repr = "mesh modifiers: \n";
+  repr += subdivisions_to_string();
+
+  return repr;
+}
+
+//===== setting modifiers =====
+template <specfem::dimension::type DimensionType>
+void specfem::mesh::modifiers<DimensionType>::set_subdivision(
+    const int material, subdiv_tuple subdivs) {
+  subdivisions.insert(std::make_pair(0, subdivs));
+}
+//===== getting modifiers =====
+template <specfem::dimension::type DimensionType>
+typename specfem::mesh::modifiers<DimensionType>::subdiv_tuple
+specfem::mesh::modifiers<DimensionType>::get_subdivision(
+    const int material) const {
+  auto got = subdivisions.find(material);
+  if (got == subdivisions.end()) {
+    // default: no subdividing (1 subdiv in z, 1 in x)
+    return {};
+  } else {
+    return got->second;
+  }
+}
+
+template class specfem::mesh::modifiers<specfem::dimension::type::dim2>;
+template class specfem::mesh::modifiers<specfem::dimension::type::dim3>;
diff --git a/src/parameter_parser/setup.cpp b/src/parameter_parser/setup.cpp
index 4e7e8840..12317df8 100644
--- a/src/parameter_parser/setup.cpp
+++ b/src/parameter_parser/setup.cpp
@@ -77,6 +77,13 @@ specfem::runtime_configuration::setup::setup(const YAML::Node &parameter_dict,
     throw std::runtime_error("Error reading specfem runtime configuration.");
   }
 
+  // mesh modifiers: not required
+  if (const YAML::Node &n_mesh_modifiers = runtime_config["mesh-modifiers"]) {
+    this->mesh_modifiers =
+        std::make_unique<specfem::runtime_configuration::mesh_modifiers>(
+            n_mesh_modifiers);
+  }
+
   try {
     this->databases = std::make_unique<
         specfem::runtime_configuration::database_configuration>(n_databases);
diff --git a/tests/unit-tests/CMakeLists.txt b/tests/unit-tests/CMakeLists.txt
index c2ead6bd..21ebb772 100644
--- a/tests/unit-tests/CMakeLists.txt
+++ b/tests/unit-tests/CMakeLists.txt
@@ -356,6 +356,31 @@ target_link_libraries(
   -lpthread -lm
 )
 
+
+add_executable(
+  mesh_modifiers_tests
+  mesh_modifiers/runner.cpp
+  mesh_modifiers/subdivision/test_subdivide.cpp
+)
+
+target_link_libraries(
+  mesh_modifiers_tests
+  quadrature
+  mesh
+  IO
+  # material_class
+  yaml-cpp
+  kokkos_environment
+  mpi_environment
+  compute
+  parameter_reader
+  timescheme
+  coupled_interface
+  periodic_tasks
+  -lpthread -lm
+)
+
+
 # add_executable(
 #   seismogram_elastic_tests
 #   seismogram/elastic/seismogram_tests.cpp
diff --git a/tests/unit-tests/mesh_modifiers/runner.cpp b/tests/unit-tests/mesh_modifiers/runner.cpp
new file mode 100644
index 00000000..7ff9188a
--- /dev/null
+++ b/tests/unit-tests/mesh_modifiers/runner.cpp
@@ -0,0 +1,10 @@
+#include "../Kokkos_Environment.hpp"
+#include "../MPI_environment.hpp"
+#include "gtest/gtest.h"
+
+int main(int argc, char *argv[]) {
+  ::testing::InitGoogleTest(&argc, argv);
+  ::testing::AddGlobalTestEnvironment(new MPIEnvironment);
+  ::testing::AddGlobalTestEnvironment(new KokkosEnvironment);
+  return RUN_ALL_TESTS();
+}
diff --git a/tests/unit-tests/mesh_modifiers/subdivision/test_subdivide.cpp b/tests/unit-tests/mesh_modifiers/subdivision/test_subdivide.cpp
new file mode 100644
index 00000000..2ff3359d
--- /dev/null
+++ b/tests/unit-tests/mesh_modifiers/subdivision/test_subdivide.cpp
@@ -0,0 +1,516 @@
+#include "../../Kokkos_Environment.hpp"
+#include "../../MPI_environment.hpp"
+#include "IO/interface.hpp"
+#include "kokkos_abstractions.h"
+#include "mesh/mesh.hpp"
+#include "mesh/modifiers/modifiers.hpp"
+#include <limits>
+#include <stdexcept>
+
+namespace specfem {
+namespace test {
+namespace mesh_modifiers {
+namespace subdivisions {
+
+// we may want to look into implementing a kd tree. Instead, just bin things
+// into a quadtree, using heap-stored nodes
+template <typename T, int num_children> struct tree_node {
+  tree_node<T, num_children> *parent;
+  T data;
+  std::vector<tree_node<T, num_children> > children;
+  int depth;
+  bool leaf;
+
+  tree_node() : depth(-1), leaf(true) {}
+
+  void gen_children() {
+    if (leaf) {
+      children = std::vector<tree_node<T, num_children> >(num_children);
+      for (int i = 0; i < num_children; i++) {
+        children[i] = tree_node<T, num_children>();
+        children[i].parent = this;
+        children[i].depth = depth + 1;
+      }
+      leaf = false;
+    }
+  }
+};
+
+// utility to store control nodes. contains interpolation functions
+struct elem_coords {
+  type_real pts[3][3][2];
+  type_real cx, cz;
+  unsigned char mode;
+
+  elem_coords() = default;
+  elem_coords(const Kokkos::View<int **, Kokkos::HostSpace> &knods,
+              const Kokkos::View<type_real **, Kokkos::HostSpace> &coord,
+              const int ispec) {
+    const int ngnod = knods.extent(0);
+    switch (ngnod) {
+    case 4: {
+      for (int idim = 0; idim < 2; idim++) {
+        pts[0][0][idim] = coord(idim, knods(0, ispec));
+        pts[0][1][idim] = coord(idim, knods(1, ispec));
+        pts[1][1][idim] = coord(idim, knods(2, ispec));
+        pts[1][0][idim] = coord(idim, knods(3, ispec));
+      }
+      mode = 0;
+      std::tie(cx, cz) = interpolate(0, 0);
+      return;
+    }
+    case 9: {
+      for (int idim = 0; idim < 2; idim++) {
+        pts[0][0][idim] = coord(idim, knods(0, ispec));
+        pts[0][2][idim] = coord(idim, knods(1, ispec));
+        pts[2][2][idim] = coord(idim, knods(2, ispec));
+        pts[2][0][idim] = coord(idim, knods(3, ispec));
+        pts[0][1][idim] = coord(idim, knods(4, ispec));
+        pts[1][2][idim] = coord(idim, knods(5, ispec));
+        pts[2][1][idim] = coord(idim, knods(6, ispec));
+        pts[1][0][idim] = coord(idim, knods(7, ispec));
+        pts[1][1][idim] = coord(idim, knods(8, ispec));
+      }
+      mode = 1;
+      cx = pts[1][1][0];
+      cz = pts[1][1][1];
+      return;
+    }
+    default:
+      throw std::runtime_error(
+          "number of mesh control nodes invalid: ngnod = " +
+          std::to_string(ngnod));
+      return;
+    }
+  }
+  elem_coords(const specfem::mesh::mesh<specfem::dimension::type::dim2> &mesh,
+              const int ispec)
+      : elem_coords(mesh.control_nodes.knods, mesh.control_nodes.coord, ispec) {
+  }
+
+  std::pair<type_real, type_real> interpolate(type_real xi, type_real gamma) {
+    switch (mode) {
+    case 0: {
+      type_real bot, top;
+      type_real out[2];
+      for (int idim = 0; idim < 2; idim++) {
+        bot = (pts[0][0][idim] * (1 - xi) + pts[0][1][idim] * (1 + xi)) / 2;
+        top = (pts[1][0][idim] * (1 - xi) + pts[1][1][idim] * (1 + xi)) / 2;
+        out[idim] = (bot * (1 - gamma) + top * (1 + gamma)) / 2;
+      }
+      return std::make_pair(out[0], out[1]);
+    }
+    case 1: {
+      type_real xavgs[3];
+      type_real out[2];
+      type_real xi2 = xi * xi;
+      type_real ga2 = gamma * gamma;
+      type_real coefxi1 = (xi2 - xi) / 2;
+      type_real coefxi2 = (1 - xi2);
+      type_real coefxi3 = (xi + xi2) / 2;
+
+      type_real coefga1 = (ga2 - gamma) / 2;
+      type_real coefga2 = (1 - ga2);
+      type_real coefga3 = (gamma + ga2) / 2;
+      for (int idim = 0; idim < 2; idim++) {
+        for (int iz = 0; iz < 3; iz++) {
+          xavgs[iz] = pts[iz][0][idim] * coefxi1 + pts[iz][1][idim] * coefxi2 +
+                      pts[iz][2][idim] * coefxi3;
+        }
+        out[idim] =
+            xavgs[0] * coefga1 + xavgs[1] * coefga2 + xavgs[2] * coefga3;
+      }
+      return std::make_pair(out[0], out[1]);
+    }
+    default:
+      return std::make_pair(-1, -1);
+    }
+  }
+  type_real small_len2() {
+    type_real x1, z1, x2, z2;
+    std::tie(x1, z1) = interpolate(-1, 0);
+    std::tie(x2, z2) = interpolate(1, 0);
+    type_real len = (x2 - x1) * (x2 - x1) + (z2 - z1) * (z2 - z1);
+    std::tie(x1, z1) = interpolate(0, -1);
+    std::tie(x2, z2) = interpolate(0, 1);
+    return std::min(len, (x2 - x1) * (x2 - x1) + (z2 - z1) * (z2 - z1));
+  }
+};
+
+// bins to store elements into.
+struct elem_coord_bin2D {
+  std::vector<int> indices; // this will only be nonempty for leaf nodes
+  type_real xmin, xmax, ymin, ymax, cx, cy; // coordinates for bin (aabb)
+  tree_node<elem_coord_bin2D, 4> *parent;   // parent->data == *this
+
+  elem_coord_bin2D() = default;
+  elem_coord_bin2D(tree_node<elem_coord_bin2D, 4> &treenode, type_real xmin,
+                   type_real xmax, type_real ymin, type_real ymax)
+      : xmin(xmin), xmax(xmax), ymin(ymin), ymax(ymax), parent(&treenode),
+        cx((xmin + xmax) / 2), cy((ymin + ymax) / 2) {}
+  elem_coord_bin2D(tree_node<elem_coord_bin2D, 4> &treenode, int index)
+      : parent(&treenode) {
+    tree_node<elem_coord_bin2D, 4> &parent = *treenode.parent;
+    type_real cx = (parent.data.xmax + parent.data.xmin) / 2;
+    type_real cy = (parent.data.ymax + parent.data.ymin) / 2;
+    // inherit aabb from quad-subdivided parent's aabb
+    switch (index) {
+    case 0:
+      // lower left
+      xmin = parent.data.xmin;
+      ymin = parent.data.ymin;
+      xmax = cx;
+      ymax = cy;
+      break;
+    case 1:
+      // lower right
+      xmin = cx;
+      ymin = parent.data.ymin;
+      xmax = parent.data.xmax;
+      ymax = cy;
+      break;
+    case 2:
+      // upper left
+      xmin = parent.data.xmin;
+      ymin = cy;
+      xmax = cx;
+      ymax = parent.data.ymax;
+      break;
+    case 3:
+      // upper right
+      xmin = cx;
+      ymin = cy;
+      xmax = parent.data.ymax;
+      ymax = parent.data.ymax;
+      break;
+    }
+    this->cx = (xmin + xmax) / 2;
+    this->cy = (ymin + ymax) / 2;
+  }
+
+  /**
+   * @brief Subdivides this node. every index in this node is then passed to one
+   * of the children.
+   *
+   * @param coordref - ispec -> elemdata mapping
+   */
+  void subdivide(const std::vector<elem_coords> &coordref) {
+    if (!parent->leaf) {
+      // subdivide underlying tree structure
+      parent->gen_children();
+      for (int i = 0; i < 4; i++) {
+        parent->children[i].data = elem_coord_bin2D(parent->children[i], i);
+      }
+    }
+
+    // insert() will properly insert each index into the correct child
+    for (const int &i : indices) {
+      insert(coordref, i);
+    }
+    // we're moving, not copying
+    indices.clear();
+  }
+
+  /**
+   * @brief Inserts the ispec into the tree at this node. If we are a leaf,
+   * places into this index. Otherwise, it is inserted into the correct child.
+   *
+   * @param coordref - ispec -> elemdata mapping
+   * @param ispec - element to store
+   * @param max_storage - threshold for when a leaf is large enough to warrant
+   * subdivision.
+   * @param maxdepth - the furthest subdivision allowed before max_storage is
+   * ignored.
+   */
+  void insert(const std::vector<elem_coords> &coordref, int ispec,
+              int max_storage = 20, int maxdepth = 10) {
+    if (parent->leaf || parent->depth >= maxdepth) {
+      indices.push_back(ispec);
+      return;
+    } else {
+      const int sub = ((coordref[ispec].cx > cx) ? 1 : 0) |
+                      ((coordref[ispec].cz > cy) ? 2 : 0);
+      auto &child = parent->children[sub].data;
+      if (child.indices.size() >= max_storage) {
+        child.subdivide(coordref);
+      }
+      child.insert(coordref, ispec, max_storage, maxdepth);
+    }
+  }
+
+  /**
+   * @brief Finds the k nearest elements to the given coordinates.
+   *
+   * @param coordref - element coord vector
+   * @param x - x coordinate to find elements near
+   * @param y - y (z) coordinate to find elements near
+   * @param k - number of elements to search for
+   * @param nearest_inds a vector of <ispec, dist^2> for the solutions. This
+   * should start empty, but will be set, in order from closest to farthest.
+   */
+  void knearest(const std::vector<elem_coords> &coordref, const type_real x,
+                const type_real y, const int k,
+                std::vector<std::pair<int, type_real> > &nearest_inds) {
+    // trace down to leaf containing (x,y), unless we collected everything, and
+    // we're too far
+
+    int ncollected = nearest_inds.size();
+    if (ncollected >= k) {
+      // nearest_inds is full. return if the this box cannot beat the last of
+      // the closest
+      type_real boxclosestx = std::min(xmax, std::max(xmin, x));
+      type_real boxclosesty = std::min(ymax, std::max(ymin, y));
+      type_real dist2 = (x - boxclosestx) * (x - boxclosestx) +
+                        (y - boxclosesty) * (y - boxclosesty);
+      if (dist2 > nearest_inds[k - 1].second) {
+        return;
+      }
+    }
+    if (parent->leaf) {
+      // linsearch
+      int nelems = indices.size();
+
+      type_real maxdist = (ncollected >= k)
+                              ? nearest_inds[k - 1].second
+                              : std::numeric_limits<type_real>::max();
+      for (int ielem = 0; ielem < nelems; ielem++) {
+        int ispec = indices[ielem];
+        type_real dist = (x - coordref[ispec].cx) * (x - coordref[ispec].cx) +
+                         (y - coordref[ispec].cz) * (y - coordref[ispec].cz);
+        if (dist < maxdist) {
+          if (ncollected < k) {
+            ncollected++;
+            // the last entry will always be overwritten
+            nearest_inds.push_back(std::make_pair(0, 0));
+          }
+
+          // make it so ibubble is the open space where the collecteds are in
+          // order.
+          int ibubble = ncollected - 1;
+          while (ibubble > 0 && nearest_inds[ibubble - 1].second > dist) {
+            nearest_inds[ibubble] = nearest_inds[ibubble - 1];
+            ibubble--;
+          }
+          nearest_inds[ibubble] = std::make_pair(ispec, dist);
+          // update maxdist in case we need to
+          maxdist = (ncollected >= k) ? nearest_inds[k - 1].second
+                                      : std::numeric_limits<type_real>::max();
+        }
+      }
+    } else {
+      // try the "closest" sub-boxes in order. no complicated math; just simple
+      // heuristics
+      int xprio = (x > cx) ? 1 : 0;
+      int yprio = (y > cy) ? 2 : 0;
+      parent->children[yprio + xprio].data.knearest(coordref, x, y, k,
+                                                    nearest_inds);
+      parent->children[(2 - yprio) + xprio].data.knearest(coordref, x, y, k,
+                                                          nearest_inds);
+      parent->children[yprio + (1 - xprio)].data.knearest(coordref, x, y, k,
+                                                          nearest_inds);
+      parent->children[(2 - yprio) + (1 - xprio)].data.knearest(
+          coordref, x, y, k, nearest_inds);
+    }
+  }
+};
+
+void test_mesh_unifsubdivisions2D(
+    specfem::mesh::mesh<specfem::dimension::type::dim2> &mesh,
+    int global_subdiv_factor_x, int global_subdiv_factor_z) {
+  const int ngnod = mesh.control_nodes.ngnod;
+
+  type_real xmin, xmax, zmin, zmax;
+  xmin = xmax = mesh.control_nodes.coord(0, 0);
+  zmin = zmax = mesh.control_nodes.coord(1, 0);
+
+  for (int ispec = 0; ispec < mesh.nspec; ispec++) {
+    xmin = std::min(xmin, mesh.control_nodes.coord(0, ispec));
+    zmin = std::min(zmin, mesh.control_nodes.coord(1, ispec));
+    xmax = std::max(xmax, mesh.control_nodes.coord(0, ispec));
+    zmax = std::max(zmax, mesh.control_nodes.coord(1, ispec));
+  }
+
+  specfem::mesh::mesh<specfem::dimension::type::dim2> meshcopy = mesh;
+  specfem::mesh::modifiers<specfem::dimension::type::dim2> modifiers;
+  for (int imat = 0; imat < mesh.materials.n_materials; imat++) {
+    modifiers.set_subdivision(
+        imat, { global_subdiv_factor_x, global_subdiv_factor_z });
+  }
+  modifiers.apply(meshcopy);
+
+  // check meshes
+  specfem::kokkos::HostView3d<bool> hits("test_mesh_unifsubdivisions2D hits",
+                                         mesh.nspec, global_subdiv_factor_z,
+                                         global_subdiv_factor_x);
+  // storing control nodes
+  std::vector<elem_coords> base_elems;
+  // quadtree bins
+  tree_node<elem_coord_bin2D, 4> base_bins;
+  base_bins.data = elem_coord_bin2D(base_bins, xmin, xmax, zmin, zmax);
+  base_bins.depth = 0;
+  for (int ispec = 0; ispec < mesh.nspec; ispec++) {
+    base_elems.push_back(elem_coords(mesh, ispec));
+    base_bins.data.insert(base_elems, ispec);
+  }
+
+  for (int ispec = 0; ispec < meshcopy.nspec; ispec++) {
+    elem_coords elem(meshcopy, ispec);
+    // find the tuple (ispec, isubx, isubz) that this ispec_new should be in
+    std::vector<std::pair<int, type_real> > nearest_inds;
+    constexpr int num_nearest_search = 9;
+    base_bins.data.knearest(base_elems, elem.cx, elem.cz, num_nearest_search,
+                            nearest_inds);
+    int ispec_near, isubx_near, isubz_near;
+    type_real min_err = std::numeric_limits<type_real>::max();
+    type_real equality_threshold = elem.small_len2() * 1e-3;
+    type_real inv_subfact_x = 1.0 / global_subdiv_factor_x;
+    type_real inv_subfact_z = 1.0 / global_subdiv_factor_z;
+    for (int inear = 0; inear < num_nearest_search; inear++) {
+      int ispec_orig = nearest_inds[inear].first;
+      elem_coords &elem_orig = base_elems[ispec_orig];
+      // the center of elem *should* have parameter 2*[ (isub + 0.5)/sub ] - 1
+      for (int isubx = 0; isubx < global_subdiv_factor_x; isubx++) {
+        for (int isubz = 0; isubz < global_subdiv_factor_z; isubz++) {
+          type_real x, z;
+          std::tie(x, z) =
+              elem_orig.interpolate((2 * isubx + 1) * inv_subfact_x - 1,
+                                    (2 * isubz + 1) * inv_subfact_z - 1);
+          type_real err =
+              (x - elem.cx) * (x - elem.cx) + (z - elem.cz) * (z - elem.cz);
+          if (err < min_err) {
+            min_err = err;
+            ispec_near = ispec_orig;
+            isubx_near = isubx;
+            isubz_near = isubz;
+          }
+        }
+      }
+      if (min_err < equality_threshold) {
+        // we can be sure we found it.
+        break;
+      }
+    }
+    // now check: are the control nodes in good places?
+    type_real pt_errs[3][3];
+    elem_coords &elem_near = base_elems[ispec_near];
+    type_real xi[3] = { (2 * isubx_near) * inv_subfact_x - 1,
+                        (2 * isubx_near + 1) * inv_subfact_x - 1,
+                        (2 * isubx_near + 2) * inv_subfact_x - 1 };
+    type_real ga[3] = { (2 * isubz_near) * inv_subfact_z - 1,
+                        (2 * isubz_near + 1) * inv_subfact_z - 1,
+                        (2 * isubz_near + 2) * inv_subfact_z - 1 };
+    bool fail = false;
+    for (int ix = 0; ix < 3; ix++) {
+      for (int iz = 0; iz < 3; iz++) {
+        type_real x1, z1, x2, z2;
+        std::tie(x1, z1) = elem.interpolate(ix - 1, iz - 1);
+        std::tie(x2, z2) = elem_near.interpolate(xi[ix], ga[iz]);
+        pt_errs[ix][iz] = (x2 - x1) * (x2 - x1) + (z2 - z1) * (z2 - z1);
+        if (pt_errs[ix][iz] > equality_threshold) {
+          fail = true;
+        }
+      }
+    }
+
+    if (fail) {
+      char pt_char[3][3];
+      for (int ix = 0; ix < 3; ix++) {
+        for (int iz = 0; iz < 3; iz++) {
+          pt_char[ix][iz] = (pt_errs[ix][iz] > equality_threshold) ? 'x' : ' ';
+        }
+      }
+
+      FAIL() << "--------------------------------------------------\n"
+             << "\033[0;31m[FAILED]\033[0m Test failed\n"
+             << " - global subdivision: " << global_subdiv_factor_x << "x"
+             << global_subdiv_factor_z << "\n"
+             << " - subdivided element (" << isubx_near << "," << isubz_near
+             << ") of element " << ispec_near << "\n"
+             << " - 9 node chart:\n"
+             << "     " << pt_char[0][2] << pt_char[1][2] << pt_char[2][2]
+             << "\n"
+             << "     " << pt_char[0][1] << pt_char[1][1] << pt_char[2][1]
+             << "\n"
+             << "     " << pt_char[0][0] << pt_char[1][0] << pt_char[2][0]
+             << "\n\n"
+             << "     " << pt_errs[0][2] << "\t" << pt_errs[1][2] << "\t"
+             << pt_errs[2][2] << "\n"
+             << "     " << pt_errs[0][1] << "\t" << pt_errs[1][1] << "\t"
+             << pt_errs[2][1] << "\n"
+             << "     " << pt_errs[0][0] << "\t" << pt_errs[1][0] << "\t"
+             << pt_errs[2][0] << "\n"
+             << "--------------------------------------------------\n\n"
+             << std::endl;
+    }
+    // this is precisely the subdivision. Ensure the material is correct
+    if (mesh.materials.material_index_mapping(ispec_near).database_index !=
+        meshcopy.materials.material_index_mapping(ispec).database_index) {
+      FAIL() << "--------------------------------------------------\n"
+             << "\033[0;31m[FAILED]\033[0m Test failed\n"
+             << " - global subdivision: " << global_subdiv_factor_x << "x"
+             << global_subdiv_factor_z << "\n"
+             << " - subdivided element (" << isubx_near << "," << isubz_near
+             << ") of element " << ispec_near << "\n"
+             << " - Incorrect material index mapping:\n"
+             << "     should be "
+             << mesh.materials.material_index_mapping(ispec_near).database_index
+             << "\n"
+             << "     found "
+             << meshcopy.materials.material_index_mapping(ispec).database_index
+             << "\n"
+             << "--------------------------------------------------\n\n"
+             << std::endl;
+    }
+    if (hits(ispec_near, isubz_near, isubx_near)) {
+      // we already hit this element; something went wrong.
+      FAIL() << "--------------------------------------------------\n"
+             << "\033[0;31m[FAILED]\033[0m Test failed\n"
+             << " - global subdivision: " << global_subdiv_factor_x << "x"
+             << global_subdiv_factor_z << "\n"
+             << " - subdivided element (" << isubx_near << "," << isubz_near
+             << ") of element " << ispec_near << "\n"
+             << " - This subelement was hit multiple times!\n"
+             << "     second ispec_new that hit: " << ispec << "\n"
+             << "--------------------------------------------------\n\n"
+             << std::endl;
+    }
+    // mark this element as hit
+    hits(ispec_near, isubz_near, isubx_near) = true;
+  }
+
+  // ensure every element was hit
+  for (int ispec = 0; ispec < mesh.nspec; ispec++) {
+    for (int isubz = 0; isubz < global_subdiv_factor_z; isubz++) {
+      for (int isubx = 0; isubx < global_subdiv_factor_x; isubx++) {
+        if (!hits(ispec, isubz, isubx)) {
+          FAIL() << "--------------------------------------------------\n"
+                 << "\033[0;31m[FAILED]\033[0m Test failed\n"
+                 << " - global subdivision: " << global_subdiv_factor_x << "x"
+                 << global_subdiv_factor_z << "\n"
+                 << " - subdivided element (" << isubx << "," << isubz
+                 << ") of element " << ispec << "\n"
+                 << " - This subelement was never hit!\n"
+                 << "--------------------------------------------------\n\n"
+                 << std::endl;
+        }
+      }
+    }
+  }
+}
+} // namespace subdivisions
+} // namespace mesh_modifiers
+} // namespace test
+} // namespace specfem
+
+TEST(MESH_MODIFIERS, subdivision) {
+  auto mpi = MPIEnvironment::get_mpi();
+  auto mesh =
+      specfem::IO::read_mesh("../../../tests/unit-tests/displacement_tests/"
+                             "Newmark/serial/test1/database.bin",
+                             mpi);
+
+  specfem::test::mesh_modifiers::subdivisions::test_mesh_unifsubdivisions2D(
+      mesh, 2, 3);
+  specfem::test::mesh_modifiers::subdivisions::test_mesh_unifsubdivisions2D(
+      mesh, 4, 2);
+}