diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 99d9680a55263e0de7434676fb5478df7d6438cb..78b72bcafc640ecc888c0997eeb3ff8f3056db58 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -35,6 +35,7 @@ add_library(
   MeshLineFaceBoundary.cpp
   MeshLineNodeBoundary.cpp
   MeshNodeBoundary.cpp
+  MeshNodeBoundaryUtils.cpp
   MeshNodeInterface.cpp
   MeshRandomizer.cpp
   MeshRelaxer.cpp
diff --git a/src/mesh/MeshEdgeBoundary.cpp b/src/mesh/MeshEdgeBoundary.cpp
index 4c60b3424e93470e0135ee34fb81e870058e2fa9..0a4aed04bb9e54e0f35dda08e3b9ce5e8f584559 100644
--- a/src/mesh/MeshEdgeBoundary.cpp
+++ b/src/mesh/MeshEdgeBoundary.cpp
@@ -5,19 +5,19 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <size_t Dimension>
-MeshEdgeBoundary<Dimension>::MeshEdgeBoundary(const Mesh<Connectivity<Dimension>>&, const RefEdgeList& ref_edge_list)
-  : m_ref_edge_list(ref_edge_list)
+template <typename MeshType>
+MeshEdgeBoundary::MeshEdgeBoundary(const MeshType&, const RefEdgeList& ref_edge_list) : m_ref_edge_list(ref_edge_list)
 {}
 
-template MeshEdgeBoundary<1>::MeshEdgeBoundary(const Mesh<Connectivity<1>>&, const RefEdgeList&);
-template MeshEdgeBoundary<2>::MeshEdgeBoundary(const Mesh<Connectivity<2>>&, const RefEdgeList&);
-template MeshEdgeBoundary<3>::MeshEdgeBoundary(const Mesh<Connectivity<3>>&, const RefEdgeList&);
+template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<Connectivity<1>>&, const RefEdgeList&);
+template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<Connectivity<2>>&, const RefEdgeList&);
+template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<Connectivity<3>>&, const RefEdgeList&);
 
-template <size_t Dimension>
-MeshEdgeBoundary<Dimension>::MeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
-                                              const RefFaceList& ref_face_list)
+template <typename MeshType>
+MeshEdgeBoundary::MeshEdgeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
 {
+  constexpr size_t Dimension = MeshType::Dimension;
+
   const Array<const FaceId>& face_list = ref_face_list.list();
   static_assert(Dimension > 1, "conversion from to edge from face is valid in dimension > 1");
 
@@ -55,13 +55,15 @@ MeshEdgeBoundary<Dimension>::MeshEdgeBoundary(const Mesh<Connectivity<Dimension>
   const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_edge_list);
 }
 
-template MeshEdgeBoundary<2>::MeshEdgeBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&);
-template MeshEdgeBoundary<3>::MeshEdgeBoundary(const Mesh<Connectivity<3>>&, const RefFaceList&);
+template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&);
+template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<Connectivity<3>>&, const RefFaceList&);
 
-template <size_t Dimension>
-MeshEdgeBoundary<Dimension>
-getMeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+template <typename MeshType>
+MeshEdgeBoundary
+getMeshEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
+  constexpr size_t Dimension = MeshType::Dimension;
+
   for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<ItemType::edge>();
        ++i_ref_edge_list) {
     const auto& ref_edge_list = mesh.connectivity().template refItemList<ItemType::edge>(i_ref_edge_list);
@@ -76,7 +78,7 @@ getMeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDe
         throw NormalError(ost.str());
       }
 
-      return MeshEdgeBoundary<Dimension>{mesh, ref_edge_list};
+      return MeshEdgeBoundary{mesh, ref_edge_list};
     }
   }
   if constexpr (Dimension > 1) {
@@ -94,7 +96,7 @@ getMeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDe
           throw NormalError(ost.str());
         }
 
-        return MeshEdgeBoundary<Dimension>{mesh, ref_face_list};
+        return MeshEdgeBoundary{mesh, ref_face_list};
       }
     }
   }
@@ -105,6 +107,6 @@ getMeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDe
   throw NormalError(ost.str());
 }
 
-template MeshEdgeBoundary<1> getMeshEdgeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
-template MeshEdgeBoundary<2> getMeshEdgeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
-template MeshEdgeBoundary<3> getMeshEdgeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
+template MeshEdgeBoundary getMeshEdgeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
+template MeshEdgeBoundary getMeshEdgeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
+template MeshEdgeBoundary getMeshEdgeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshEdgeBoundary.hpp b/src/mesh/MeshEdgeBoundary.hpp
index e37c4fdffefdc494c97bd8ddbea61f656ef938b7..c65106b2d3d0d7958f0c3010295f07f7a9e9b6d9 100644
--- a/src/mesh/MeshEdgeBoundary.hpp
+++ b/src/mesh/MeshEdgeBoundary.hpp
@@ -6,52 +6,47 @@
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
-template <size_t Dimension>
-class Connectivity;
-
-template <typename ConnectivityType>
-class Mesh;
-
-template <size_t Dimension>
 class [[nodiscard]] MeshEdgeBoundary   // clazy:exclude=copyable-polymorphic
 {
  protected:
   RefEdgeList m_ref_edge_list;
 
  public:
-  template <size_t MeshDimension>
-  friend MeshEdgeBoundary<MeshDimension> getMeshEdgeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
-                                                             const IBoundaryDescriptor& boundary_descriptor);
+  template <typename MeshType>
+  friend MeshEdgeBoundary getMeshEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
   MeshEdgeBoundary& operator=(const MeshEdgeBoundary&) = default;
-  MeshEdgeBoundary& operator=(MeshEdgeBoundary&&) = default;
+  MeshEdgeBoundary& operator=(MeshEdgeBoundary&&)      = default;
 
   PUGS_INLINE
-  const RefEdgeList& refEdgeList() const
+  const RefEdgeList&
+  refEdgeList() const
   {
     return m_ref_edge_list;
   }
 
   PUGS_INLINE
-  const Array<const EdgeId>& edgeList() const
+  const Array<const EdgeId>&
+  edgeList() const
   {
     return m_ref_edge_list.list();
   }
 
  protected:
-  MeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefEdgeList& ref_edge_list);
-  MeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list);
+  template <typename MeshType>
+  MeshEdgeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list);
+  template <typename MeshType>
+  MeshEdgeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list);
 
  public:
   MeshEdgeBoundary(const MeshEdgeBoundary&) = default;   // LCOV_EXCL_LINE
-  MeshEdgeBoundary(MeshEdgeBoundary &&)     = default;   // LCOV_EXCL_LINE
+  MeshEdgeBoundary(MeshEdgeBoundary&&)      = default;   // LCOV_EXCL_LINE
 
   MeshEdgeBoundary()          = default;
   virtual ~MeshEdgeBoundary() = default;
 };
 
-template <size_t Dimension>
-MeshEdgeBoundary<Dimension> getMeshEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
-                                                const IBoundaryDescriptor& boundary_descriptor);
+template <typename MeshType>
+MeshEdgeBoundary getMeshEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_EDGE_BOUNDARY_HPP
diff --git a/src/mesh/MeshFaceBoundary.cpp b/src/mesh/MeshFaceBoundary.cpp
index ad77830f5cbcb10070cfb7cb4989a7065a650f16..3a578bb7893040ca3838581b0b7c5ad5d297d9a3 100644
--- a/src/mesh/MeshFaceBoundary.cpp
+++ b/src/mesh/MeshFaceBoundary.cpp
@@ -4,18 +4,17 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <size_t Dimension>
-MeshFaceBoundary<Dimension>::MeshFaceBoundary(const Mesh<Connectivity<Dimension>>&, const RefFaceList& ref_face_list)
-  : m_ref_face_list(ref_face_list)
+template <typename MeshType>
+MeshFaceBoundary::MeshFaceBoundary(const MeshType&, const RefFaceList& ref_face_list) : m_ref_face_list(ref_face_list)
 {}
 
-template MeshFaceBoundary<1>::MeshFaceBoundary(const Mesh<Connectivity<1>>&, const RefFaceList&);
-template MeshFaceBoundary<2>::MeshFaceBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&);
-template MeshFaceBoundary<3>::MeshFaceBoundary(const Mesh<Connectivity<3>>&, const RefFaceList&);
+template MeshFaceBoundary::MeshFaceBoundary(const Mesh<Connectivity<1>>&, const RefFaceList&);
+template MeshFaceBoundary::MeshFaceBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&);
+template MeshFaceBoundary::MeshFaceBoundary(const Mesh<Connectivity<3>>&, const RefFaceList&);
 
-template <size_t Dimension>
-MeshFaceBoundary<Dimension>
-getMeshFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+template <typename MeshType>
+MeshFaceBoundary
+getMeshFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
   for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
        ++i_ref_face_list) {
@@ -31,7 +30,7 @@ getMeshFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDe
         throw NormalError(ost.str());
       }
 
-      return MeshFaceBoundary<Dimension>{mesh, ref_face_list};
+      return MeshFaceBoundary{mesh, ref_face_list};
     }
   }
 
@@ -41,6 +40,6 @@ getMeshFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDe
   throw NormalError(ost.str());
 }
 
-template MeshFaceBoundary<1> getMeshFaceBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
-template MeshFaceBoundary<2> getMeshFaceBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
-template MeshFaceBoundary<3> getMeshFaceBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
+template MeshFaceBoundary getMeshFaceBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
+template MeshFaceBoundary getMeshFaceBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
+template MeshFaceBoundary getMeshFaceBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshFaceBoundary.hpp b/src/mesh/MeshFaceBoundary.hpp
index 98a63774138e60bbb01437ee4fb6d88bd1400e40..c44aafa66c6f9f47dfc839a39025a42ae614b441 100644
--- a/src/mesh/MeshFaceBoundary.hpp
+++ b/src/mesh/MeshFaceBoundary.hpp
@@ -6,51 +6,45 @@
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
-template <size_t Dimension>
-class Connectivity;
-
-template <typename ConnectivityType>
-class Mesh;
-
-template <size_t Dimension>
 class [[nodiscard]] MeshFaceBoundary   // clazy:exclude=copyable-polymorphic
 {
  protected:
   RefFaceList m_ref_face_list;
 
  public:
-  template <size_t MeshDimension>
-  friend MeshFaceBoundary<MeshDimension> getMeshFaceBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
-                                                             const IBoundaryDescriptor& boundary_descriptor);
+  template <typename MeshTypeT>
+  friend MeshFaceBoundary getMeshFaceBoundary(const MeshTypeT& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
   MeshFaceBoundary& operator=(const MeshFaceBoundary&) = default;
-  MeshFaceBoundary& operator=(MeshFaceBoundary&&) = default;
+  MeshFaceBoundary& operator=(MeshFaceBoundary&&)      = default;
 
   PUGS_INLINE
-  const RefFaceList& refFaceList() const
+  const RefFaceList&
+  refFaceList() const
   {
     return m_ref_face_list;
   }
 
   PUGS_INLINE
-  const Array<const FaceId>& faceList() const
+  const Array<const FaceId>&
+  faceList() const
   {
     return m_ref_face_list.list();
   }
 
  protected:
-  MeshFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list);
+  template <typename MeshType>
+  MeshFaceBoundary(const MeshType& mesh, const RefFaceList& ref_face_list);
 
  public:
   MeshFaceBoundary(const MeshFaceBoundary&) = default;   // LCOV_EXCL_LINE
-  MeshFaceBoundary(MeshFaceBoundary &&)     = default;   // LCOV_EXCL_LINE
+  MeshFaceBoundary(MeshFaceBoundary&&)      = default;   // LCOV_EXCL_LINE
 
   MeshFaceBoundary()          = default;
   virtual ~MeshFaceBoundary() = default;
 };
 
-template <size_t Dimension>
-MeshFaceBoundary<Dimension> getMeshFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh,
-                                                const IBoundaryDescriptor& boundary_descriptor);
+template <typename MeshType>
+MeshFaceBoundary getMeshFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_FACE_BOUNDARY_HPP
diff --git a/src/mesh/MeshFlatEdgeBoundary.cpp b/src/mesh/MeshFlatEdgeBoundary.cpp
index d4c1ec155a0e8c4d863facbf6280d71265944933..3783049ff2422e899f90210c1ffba722f71a4767 100644
--- a/src/mesh/MeshFlatEdgeBoundary.cpp
+++ b/src/mesh/MeshFlatEdgeBoundary.cpp
@@ -4,17 +4,20 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshFlatNodeBoundary.hpp>
 
-template <size_t Dimension>
-MeshFlatEdgeBoundary<Dimension>
-getMeshFlatEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+template <typename MeshType>
+MeshFlatEdgeBoundary<const MeshType>
+getMeshFlatEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
-  MeshEdgeBoundary<Dimension> mesh_edge_boundary          = getMeshEdgeBoundary(mesh, boundary_descriptor);
-  MeshFlatNodeBoundary<Dimension> mesh_flat_node_boundary = getMeshFlatNodeBoundary(mesh, boundary_descriptor);
+  MeshEdgeBoundary mesh_edge_boundary          = getMeshEdgeBoundary(mesh, boundary_descriptor);
+  MeshFlatNodeBoundary mesh_flat_node_boundary = getMeshFlatNodeBoundary(mesh, boundary_descriptor);
 
-  return MeshFlatEdgeBoundary<Dimension>{mesh, mesh_edge_boundary.refEdgeList(),
-                                         mesh_flat_node_boundary.outgoingNormal()};
+  return MeshFlatEdgeBoundary<const MeshType>{mesh, mesh_edge_boundary.refEdgeList(),
+                                              mesh_flat_node_boundary.outgoingNormal()};
 }
 
-template MeshFlatEdgeBoundary<1> getMeshFlatEdgeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
-template MeshFlatEdgeBoundary<2> getMeshFlatEdgeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
-template MeshFlatEdgeBoundary<3> getMeshFlatEdgeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
+template MeshFlatEdgeBoundary<const Mesh<Connectivity<1>>> getMeshFlatEdgeBoundary(const Mesh<Connectivity<1>>&,
+                                                                                   const IBoundaryDescriptor&);
+template MeshFlatEdgeBoundary<const Mesh<Connectivity<2>>> getMeshFlatEdgeBoundary(const Mesh<Connectivity<2>>&,
+                                                                                   const IBoundaryDescriptor&);
+template MeshFlatEdgeBoundary<const Mesh<Connectivity<3>>> getMeshFlatEdgeBoundary(const Mesh<Connectivity<3>>&,
+                                                                                   const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshFlatEdgeBoundary.hpp b/src/mesh/MeshFlatEdgeBoundary.hpp
index a104d1413a7319f08b74bf1ef103a5157b803f18..60b54257fb8b038a110ba4860c4b5c54586d0307 100644
--- a/src/mesh/MeshFlatEdgeBoundary.hpp
+++ b/src/mesh/MeshFlatEdgeBoundary.hpp
@@ -3,11 +3,11 @@
 
 #include <mesh/MeshEdgeBoundary.hpp>
 
-template <size_t Dimension>
-class MeshFlatEdgeBoundary final : public MeshEdgeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
+template <typename MeshType>
+class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
-  using Rd = TinyVector<Dimension, double>;
+  using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
   const Rd m_outgoing_normal;
@@ -20,16 +20,15 @@ class MeshFlatEdgeBoundary final : public MeshEdgeBoundary<Dimension>   // clazy
   }
 
   MeshFlatEdgeBoundary& operator=(const MeshFlatEdgeBoundary&) = default;
-  MeshFlatEdgeBoundary& operator=(MeshFlatEdgeBoundary&&) = default;
+  MeshFlatEdgeBoundary& operator=(MeshFlatEdgeBoundary&&)      = default;
 
-  template <size_t MeshDimension>
-  friend MeshFlatEdgeBoundary<MeshDimension> getMeshFlatEdgeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
-                                                                     const IBoundaryDescriptor& boundary_descriptor);
+  template <typename MeshTypeT>
+  friend MeshFlatEdgeBoundary<const MeshTypeT> getMeshFlatEdgeBoundary(const MeshTypeT& mesh,
+                                                                       const IBoundaryDescriptor& boundary_descriptor);
 
  private:
-  template <typename MeshType>
   MeshFlatEdgeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list, const Rd& outgoing_normal)
-    : MeshEdgeBoundary<Dimension>(mesh, ref_edge_list), m_outgoing_normal(outgoing_normal)
+    : MeshEdgeBoundary(mesh, ref_edge_list), m_outgoing_normal(outgoing_normal)
   {}
 
  public:
@@ -39,8 +38,8 @@ class MeshFlatEdgeBoundary final : public MeshEdgeBoundary<Dimension>   // clazy
   ~MeshFlatEdgeBoundary()                           = default;
 };
 
-template <size_t Dimension>
-MeshFlatEdgeBoundary<Dimension> getMeshFlatEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
-                                                        const IBoundaryDescriptor& boundary_descriptor);
+template <typename MeshType>
+MeshFlatEdgeBoundary<const MeshType> getMeshFlatEdgeBoundary(const MeshType& mesh,
+                                                             const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_FLAT_EDGE_BOUNDARY_HPP
diff --git a/src/mesh/MeshFlatFaceBoundary.cpp b/src/mesh/MeshFlatFaceBoundary.cpp
index 87d353e25404c3b5996a63945a9e220e49e97e55..e05373897dce3421af35fa94a06bb89dd760ec35 100644
--- a/src/mesh/MeshFlatFaceBoundary.cpp
+++ b/src/mesh/MeshFlatFaceBoundary.cpp
@@ -4,17 +4,20 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshFlatNodeBoundary.hpp>
 
-template <size_t Dimension>
-MeshFlatFaceBoundary<Dimension>
-getMeshFlatFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+template <typename MeshType>
+MeshFlatFaceBoundary<const MeshType>
+getMeshFlatFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
-  MeshFaceBoundary<Dimension> mesh_face_boundary          = getMeshFaceBoundary(mesh, boundary_descriptor);
-  MeshFlatNodeBoundary<Dimension> mesh_flat_node_boundary = getMeshFlatNodeBoundary(mesh, boundary_descriptor);
+  MeshFaceBoundary mesh_face_boundary          = getMeshFaceBoundary(mesh, boundary_descriptor);
+  MeshFlatNodeBoundary mesh_flat_node_boundary = getMeshFlatNodeBoundary(mesh, boundary_descriptor);
 
-  return MeshFlatFaceBoundary<Dimension>{mesh, mesh_face_boundary.refFaceList(),
-                                         mesh_flat_node_boundary.outgoingNormal()};
+  return MeshFlatFaceBoundary<const MeshType>{mesh, mesh_face_boundary.refFaceList(),
+                                              mesh_flat_node_boundary.outgoingNormal()};
 }
 
-template MeshFlatFaceBoundary<1> getMeshFlatFaceBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
-template MeshFlatFaceBoundary<2> getMeshFlatFaceBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
-template MeshFlatFaceBoundary<3> getMeshFlatFaceBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
+template MeshFlatFaceBoundary<const Mesh<Connectivity<1>>> getMeshFlatFaceBoundary(const Mesh<Connectivity<1>>&,
+                                                                                   const IBoundaryDescriptor&);
+template MeshFlatFaceBoundary<const Mesh<Connectivity<2>>> getMeshFlatFaceBoundary(const Mesh<Connectivity<2>>&,
+                                                                                   const IBoundaryDescriptor&);
+template MeshFlatFaceBoundary<const Mesh<Connectivity<3>>> getMeshFlatFaceBoundary(const Mesh<Connectivity<3>>&,
+                                                                                   const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshFlatFaceBoundary.hpp b/src/mesh/MeshFlatFaceBoundary.hpp
index f3d34e16c5d52af38033e82dcd75e6897f392a3d..3620d0f3db39b0754fa85c72fa2b441327d7012b 100644
--- a/src/mesh/MeshFlatFaceBoundary.hpp
+++ b/src/mesh/MeshFlatFaceBoundary.hpp
@@ -3,11 +3,11 @@
 
 #include <mesh/MeshFaceBoundary.hpp>
 
-template <size_t Dimension>
-class MeshFlatFaceBoundary final : public MeshFaceBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
+template <typename MeshType>
+class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
-  using Rd = TinyVector<Dimension, double>;
+  using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
   const Rd m_outgoing_normal;
@@ -20,16 +20,15 @@ class MeshFlatFaceBoundary final : public MeshFaceBoundary<Dimension>   // clazy
   }
 
   MeshFlatFaceBoundary& operator=(const MeshFlatFaceBoundary&) = default;
-  MeshFlatFaceBoundary& operator=(MeshFlatFaceBoundary&&) = default;
+  MeshFlatFaceBoundary& operator=(MeshFlatFaceBoundary&&)      = default;
 
-  template <size_t MeshDimension>
-  friend MeshFlatFaceBoundary<MeshDimension> getMeshFlatFaceBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
-                                                                     const IBoundaryDescriptor& boundary_descriptor);
+  template <typename MeshTypeT>
+  friend MeshFlatFaceBoundary<const MeshTypeT> getMeshFlatFaceBoundary(const MeshTypeT& mesh,
+                                                                       const IBoundaryDescriptor& boundary_descriptor);
 
  private:
-  template <typename MeshType>
   MeshFlatFaceBoundary(const MeshType& mesh, const RefFaceList& ref_face_list, const Rd& outgoing_normal)
-    : MeshFaceBoundary<Dimension>(mesh, ref_face_list), m_outgoing_normal(outgoing_normal)
+    : MeshFaceBoundary(mesh, ref_face_list), m_outgoing_normal(outgoing_normal)
   {}
 
  public:
@@ -39,8 +38,8 @@ class MeshFlatFaceBoundary final : public MeshFaceBoundary<Dimension>   // clazy
   ~MeshFlatFaceBoundary()                           = default;
 };
 
-template <size_t Dimension>
-MeshFlatFaceBoundary<Dimension> getMeshFlatFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh,
-                                                        const IBoundaryDescriptor& boundary_descriptor);
+template <typename MeshType>
+MeshFlatFaceBoundary<const MeshType> getMeshFlatFaceBoundary(const MeshType& mesh,
+                                                             const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_FLAT_FACE_BOUNDARY_HPP
diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
index f5e60af046b49903304762aec86a237d7b2a868a..3ee564332b57c1bdf5a1f219af6417d9b9c0844c 100644
--- a/src/mesh/MeshFlatNodeBoundary.cpp
+++ b/src/mesh/MeshFlatNodeBoundary.cpp
@@ -2,14 +2,15 @@
 
 #include <mesh/Connectivity.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/MeshNodeBoundaryUtils.hpp>
 #include <utils/Messenger.hpp>
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-MeshFlatNodeBoundary<Dimension>::_checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
-                                                      const TinyVector<Dimension, double>& origin,
-                                                      const double length,
-                                                      const Mesh<Connectivity<Dimension>>& mesh) const
+MeshFlatNodeBoundary<MeshType>::_checkBoundaryIsFlat(const TinyVector<MeshType::Dimension, double>& normal,
+                                                     const TinyVector<MeshType::Dimension, double>& origin,
+                                                     const double length,
+                                                     const MeshType& mesh) const
 {
   const NodeValue<const Rd>& xr = mesh.xr();
 
@@ -33,7 +34,7 @@ MeshFlatNodeBoundary<Dimension>::_checkBoundaryIsFlat(const TinyVector<Dimension
 
 template <>
 TinyVector<1, double>
-MeshFlatNodeBoundary<1>::_getNormal(const Mesh<Connectivity<1>>& mesh)
+MeshFlatNodeBoundary<const Mesh<Connectivity<1>>>::_getNormal(const Mesh<Connectivity<1>>& mesh)
 {
   using R = TinyVector<1, double>;
 
@@ -59,11 +60,11 @@ MeshFlatNodeBoundary<1>::_getNormal(const Mesh<Connectivity<1>>& mesh)
 
 template <>
 TinyVector<2, double>
-MeshFlatNodeBoundary<2>::_getNormal(const Mesh<Connectivity<2>>& mesh)
+MeshFlatNodeBoundary<const Mesh<Connectivity<2>>>::_getNormal(const Mesh<Connectivity<2>>& mesh)
 {
   using R2 = TinyVector<2, double>;
 
-  std::array<R2, 2> bounds = this->_getBounds(mesh);
+  std::array<R2, 2> bounds = getBounds(mesh, m_ref_node_list);
 
   const R2& xmin = bounds[0];
   const R2& xmax = bounds[1];
@@ -87,7 +88,9 @@ MeshFlatNodeBoundary<2>::_getNormal(const Mesh<Connectivity<2>>& mesh)
 
 template <>
 TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getFarestNode(const Mesh<Connectivity<3>>& mesh, const Rd& x0, const Rd& x1)
+MeshFlatNodeBoundary<const Mesh<Connectivity<3>>>::_getFarestNode(const Mesh<Connectivity<3>>& mesh,
+                                                                  const Rd& x0,
+                                                                  const Rd& x1)
 {
   const NodeValue<const Rd>& xr = mesh.xr();
   const auto node_number        = mesh.connectivity().nodeNumber();
@@ -138,7 +141,7 @@ MeshFlatNodeBoundary<3>::_getFarestNode(const Mesh<Connectivity<3>>& mesh, const
 
 template <>
 TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getNormal(const Mesh<Connectivity<3>>& mesh)
+MeshFlatNodeBoundary<const Mesh<Connectivity<3>>>::_getNormal(const Mesh<Connectivity<3>>& mesh)
 {
   using R3 = TinyVector<3, double>;
 
@@ -159,7 +162,7 @@ MeshFlatNodeBoundary<3>::_getNormal(const Mesh<Connectivity<3>>& mesh)
     }
 
     return std::array<R3, 2>{bounds[max_i], bounds[max_j]};
-  }(this->_getBounds(mesh));
+  }(getBounds(mesh, m_ref_node_list));
 
   const R3& x0 = diagonal[0];
   const R3& x1 = diagonal[1];
@@ -195,7 +198,7 @@ MeshFlatNodeBoundary<3>::_getNormal(const Mesh<Connectivity<3>>& mesh)
 
 template <>
 TinyVector<1, double>
-MeshFlatNodeBoundary<1>::_getOutgoingNormal(const Mesh<Connectivity<1>>& mesh)
+MeshFlatNodeBoundary<const Mesh<Connectivity<1>>>::_getOutgoingNormal(const Mesh<Connectivity<1>>& mesh)
 {
   using R = TinyVector<1, double>;
 
@@ -238,7 +241,7 @@ MeshFlatNodeBoundary<1>::_getOutgoingNormal(const Mesh<Connectivity<1>>& mesh)
 
 template <>
 TinyVector<2, double>
-MeshFlatNodeBoundary<2>::_getOutgoingNormal(const Mesh<Connectivity<2>>& mesh)
+MeshFlatNodeBoundary<const Mesh<Connectivity<2>>>::_getOutgoingNormal(const Mesh<Connectivity<2>>& mesh)
 {
   using R2 = TinyVector<2, double>;
 
@@ -281,7 +284,7 @@ MeshFlatNodeBoundary<2>::_getOutgoingNormal(const Mesh<Connectivity<2>>& mesh)
 
 template <>
 TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getOutgoingNormal(const Mesh<Connectivity<3>>& mesh)
+MeshFlatNodeBoundary<const Mesh<Connectivity<3>>>::_getOutgoingNormal(const Mesh<Connectivity<3>>& mesh)
 {
   using R3 = TinyVector<3, double>;
 
@@ -322,16 +325,16 @@ MeshFlatNodeBoundary<3>::_getOutgoingNormal(const Mesh<Connectivity<3>>& mesh)
   }
 }
 
-template <size_t Dimension>
-MeshFlatNodeBoundary<Dimension>
-getMeshFlatNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+template <typename MeshType>
+MeshFlatNodeBoundary<const MeshType>
+getMeshFlatNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
   for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
        ++i_ref_node_list) {
     const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
     const RefId& ref          = ref_node_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshFlatNodeBoundary<Dimension>{mesh, ref_node_list};
+      return MeshFlatNodeBoundary<const MeshType>{mesh, ref_node_list};
     }
   }
   for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
@@ -339,7 +342,7 @@ getMeshFlatNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBounda
     const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
     const RefId& ref          = ref_face_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshFlatNodeBoundary<Dimension>{mesh, ref_face_list};
+      return MeshFlatNodeBoundary<const MeshType>{mesh, ref_face_list};
     }
   }
 
@@ -349,6 +352,9 @@ getMeshFlatNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBounda
   throw NormalError(ost.str());
 }
 
-template MeshFlatNodeBoundary<1> getMeshFlatNodeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
-template MeshFlatNodeBoundary<2> getMeshFlatNodeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
-template MeshFlatNodeBoundary<3> getMeshFlatNodeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
+template MeshFlatNodeBoundary<const Mesh<Connectivity<1>>> getMeshFlatNodeBoundary(const Mesh<Connectivity<1>>&,
+                                                                                   const IBoundaryDescriptor&);
+template MeshFlatNodeBoundary<const Mesh<Connectivity<2>>> getMeshFlatNodeBoundary(const Mesh<Connectivity<2>>&,
+                                                                                   const IBoundaryDescriptor&);
+template MeshFlatNodeBoundary<const Mesh<Connectivity<3>>> getMeshFlatNodeBoundary(const Mesh<Connectivity<3>>&,
+                                                                                   const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshFlatNodeBoundary.hpp b/src/mesh/MeshFlatNodeBoundary.hpp
index 54a10ae1c1412005440cd9e6cc49973c017265e0..df3c9d28851c635d6079f87082356717edfa0b52 100644
--- a/src/mesh/MeshFlatNodeBoundary.hpp
+++ b/src/mesh/MeshFlatNodeBoundary.hpp
@@ -3,26 +3,25 @@
 
 #include <mesh/MeshNodeBoundary.hpp>
 
-template <size_t Dimension>
-class [[nodiscard]] MeshFlatNodeBoundary final
-  : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
+template <typename MeshType>
+class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
-  using Rd = TinyVector<Dimension, double>;
+  using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
   const Rd m_outgoing_normal;
 
-  Rd _getFarestNode(const Mesh<Connectivity<Dimension>>& mesh, const Rd& x0, const Rd& x1);
+  Rd _getFarestNode(const MeshType& mesh, const Rd& x0, const Rd& x1);
 
-  Rd _getNormal(const Mesh<Connectivity<Dimension>>& mesh);
+  Rd _getNormal(const MeshType& mesh);
 
-  void _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
-                            const TinyVector<Dimension, double>& origin,
+  void _checkBoundaryIsFlat(const TinyVector<MeshType::Dimension, double>& normal,
+                            const TinyVector<MeshType::Dimension, double>& origin,
                             const double length,
-                            const Mesh<Connectivity<Dimension>>& mesh) const;
+                            const MeshType& mesh) const;
 
-  Rd _getOutgoingNormal(const Mesh<Connectivity<Dimension>>& mesh);
+  Rd _getOutgoingNormal(const MeshType& mesh);
 
  public:
   const Rd&
@@ -34,19 +33,17 @@ class [[nodiscard]] MeshFlatNodeBoundary final
   MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
   MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&)      = default;
 
-  template <size_t MeshDimension>
-  friend MeshFlatNodeBoundary<MeshDimension> getMeshFlatNodeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
-                                                                     const IBoundaryDescriptor& boundary_descriptor);
+  template <typename MeshTypeT>
+  friend MeshFlatNodeBoundary<const MeshTypeT> getMeshFlatNodeBoundary(const MeshTypeT& mesh,
+                                                                       const IBoundaryDescriptor& boundary_descriptor);
 
  private:
-  template <typename MeshType>
   MeshFlatNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
-    : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh))
+    : MeshNodeBoundary(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh))
   {}
 
-  template <typename MeshType>
   MeshFlatNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
-    : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_outgoing_normal(_getOutgoingNormal(mesh))
+    : MeshNodeBoundary(mesh, ref_node_list), m_outgoing_normal(_getOutgoingNormal(mesh))
   {}
 
  public:
@@ -56,8 +53,8 @@ class [[nodiscard]] MeshFlatNodeBoundary final
   ~MeshFlatNodeBoundary()                           = default;
 };
 
-template <size_t Dimension>
-MeshFlatNodeBoundary<Dimension> getMeshFlatNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
-                                                        const IBoundaryDescriptor& boundary_descriptor);
+template <typename MeshType>
+MeshFlatNodeBoundary<const MeshType> getMeshFlatNodeBoundary(const MeshType& mesh,
+                                                             const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_FLAT_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshLineEdgeBoundary.cpp b/src/mesh/MeshLineEdgeBoundary.cpp
index 16e2e3e9764d1465cd767c130fbf27b0e7589afd..31c4bb52465870eb4c6def597d1faa89e8df9e73 100644
--- a/src/mesh/MeshLineEdgeBoundary.cpp
+++ b/src/mesh/MeshLineEdgeBoundary.cpp
@@ -5,15 +5,18 @@
 #include <mesh/MeshLineNodeBoundary.hpp>
 #include <utils/Messenger.hpp>
 
-template <size_t Dimension>
-MeshLineEdgeBoundary<Dimension>
-getMeshLineEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+template <typename MeshType>
+MeshLineEdgeBoundary<const MeshType>
+getMeshLineEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
-  MeshEdgeBoundary<Dimension> mesh_edge_boundary          = getMeshEdgeBoundary(mesh, boundary_descriptor);
-  MeshLineNodeBoundary<Dimension> mesh_line_node_boundary = getMeshLineNodeBoundary(mesh, boundary_descriptor);
+  MeshEdgeBoundary mesh_edge_boundary          = getMeshEdgeBoundary(mesh, boundary_descriptor);
+  MeshLineNodeBoundary mesh_line_node_boundary = getMeshLineNodeBoundary(mesh, boundary_descriptor);
 
-  return MeshLineEdgeBoundary<Dimension>{mesh, mesh_edge_boundary.refEdgeList(), mesh_line_node_boundary.direction()};
+  return MeshLineEdgeBoundary<const MeshType>{mesh, mesh_edge_boundary.refEdgeList(),
+                                              mesh_line_node_boundary.direction()};
 }
 
-template MeshLineEdgeBoundary<2> getMeshLineEdgeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
-template MeshLineEdgeBoundary<3> getMeshLineEdgeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
+template MeshLineEdgeBoundary<const Mesh<Connectivity<2>>> getMeshLineEdgeBoundary(const Mesh<Connectivity<2>>&,
+                                                                                   const IBoundaryDescriptor&);
+template MeshLineEdgeBoundary<const Mesh<Connectivity<3>>> getMeshLineEdgeBoundary(const Mesh<Connectivity<3>>&,
+                                                                                   const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshLineEdgeBoundary.hpp b/src/mesh/MeshLineEdgeBoundary.hpp
index 160444648524f266941a8d2903124bf91b77ee8b..b120a6882896e354f899a8c59b178c48f8934610 100644
--- a/src/mesh/MeshLineEdgeBoundary.hpp
+++ b/src/mesh/MeshLineEdgeBoundary.hpp
@@ -4,46 +4,46 @@
 #include <algebra/TinyMatrix.hpp>
 #include <mesh/MeshEdgeBoundary.hpp>
 
-template <size_t Dimension>
-class [[nodiscard]] MeshLineEdgeBoundary final
-  : public MeshEdgeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
+template <typename MeshType>
+class [[nodiscard]] MeshLineEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
-  static_assert(Dimension > 1, "MeshLineEdgeBoundary makes only sense in dimension 1");
+  static_assert(MeshType::Dimension > 1, "MeshLineEdgeBoundary makes only sense in dimension 1");
 
-  using Rd = TinyVector<Dimension, double>;
+  using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
   const Rd m_direction;
 
  public:
-  template <size_t MeshDimension>
-  friend MeshLineEdgeBoundary<MeshDimension> getMeshLineEdgeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
-                                                                     const IBoundaryDescriptor& boundary_descriptor);
+  template <typename MeshTypeT>
+  friend MeshLineEdgeBoundary<const MeshTypeT> getMeshLineEdgeBoundary(const MeshTypeT& mesh,
+                                                                       const IBoundaryDescriptor& boundary_descriptor);
 
   PUGS_INLINE
-  const Rd& direction() const
+  const Rd&
+  direction() const
   {
     return m_direction;
   }
 
   MeshLineEdgeBoundary& operator=(const MeshLineEdgeBoundary&) = default;
-  MeshLineEdgeBoundary& operator=(MeshLineEdgeBoundary&&) = default;
+  MeshLineEdgeBoundary& operator=(MeshLineEdgeBoundary&&)      = default;
 
  private:
-  MeshLineEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefEdgeList& ref_edge_list, const Rd& direction)
-    : MeshEdgeBoundary<Dimension>(mesh, ref_edge_list), m_direction(direction)
+  MeshLineEdgeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list, const Rd& direction)
+    : MeshEdgeBoundary(mesh, ref_edge_list), m_direction(direction)
   {}
 
  public:
   MeshLineEdgeBoundary()                            = default;
   MeshLineEdgeBoundary(const MeshLineEdgeBoundary&) = default;
-  MeshLineEdgeBoundary(MeshLineEdgeBoundary &&)     = default;
+  MeshLineEdgeBoundary(MeshLineEdgeBoundary&&)      = default;
   ~MeshLineEdgeBoundary()                           = default;
 };
 
-template <size_t Dimension>
-MeshLineEdgeBoundary<Dimension> getMeshLineEdgeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
-                                                        const IBoundaryDescriptor& boundary_descriptor);
+template <typename MeshType>
+MeshLineEdgeBoundary<const MeshType> getMeshLineEdgeBoundary(const MeshType& mesh,
+                                                             const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_LINE_EDGE_BOUNDARY_HPP
diff --git a/src/mesh/MeshLineFaceBoundary.cpp b/src/mesh/MeshLineFaceBoundary.cpp
index cd67176349d0edb475001c58acc99d9cac9792ec..e81aef554062fe97509a2ebd770652828743d265 100644
--- a/src/mesh/MeshLineFaceBoundary.cpp
+++ b/src/mesh/MeshLineFaceBoundary.cpp
@@ -5,14 +5,16 @@
 #include <mesh/MeshLineNodeBoundary.hpp>
 #include <utils/Messenger.hpp>
 
-template <size_t Dimension>
-MeshLineFaceBoundary<Dimension>
-getMeshLineFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+template <typename MeshType>
+MeshLineFaceBoundary<const MeshType>
+getMeshLineFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
-  MeshFaceBoundary<Dimension> mesh_face_boundary          = getMeshFaceBoundary(mesh, boundary_descriptor);
-  MeshLineNodeBoundary<Dimension> mesh_line_node_boundary = getMeshLineNodeBoundary(mesh, boundary_descriptor);
+  MeshFaceBoundary mesh_face_boundary          = getMeshFaceBoundary(mesh, boundary_descriptor);
+  MeshLineNodeBoundary mesh_line_node_boundary = getMeshLineNodeBoundary(mesh, boundary_descriptor);
 
-  return MeshLineFaceBoundary<Dimension>{mesh, mesh_face_boundary.refFaceList(), mesh_line_node_boundary.direction()};
+  return MeshLineFaceBoundary<const MeshType>{mesh, mesh_face_boundary.refFaceList(),
+                                              mesh_line_node_boundary.direction()};
 }
 
-template MeshLineFaceBoundary<2> getMeshLineFaceBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
+template MeshLineFaceBoundary<const Mesh<Connectivity<2>>> getMeshLineFaceBoundary(const Mesh<Connectivity<2>>&,
+                                                                                   const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshLineFaceBoundary.hpp b/src/mesh/MeshLineFaceBoundary.hpp
index 5631361f5f995883d80c5317a178eab122ddcda9..f95611bb8fd7f2f02ad11e505bb11df45395e32b 100644
--- a/src/mesh/MeshLineFaceBoundary.hpp
+++ b/src/mesh/MeshLineFaceBoundary.hpp
@@ -4,46 +4,46 @@
 #include <algebra/TinyMatrix.hpp>
 #include <mesh/MeshFaceBoundary.hpp>
 
-template <size_t Dimension>
-class [[nodiscard]] MeshLineFaceBoundary final
-  : public MeshFaceBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
+template <typename MeshType>
+class [[nodiscard]] MeshLineFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
-  static_assert(Dimension == 2, "MeshLineFaceBoundary makes only sense in dimension 2");
+  static_assert(MeshType::Dimension == 2, "MeshLineFaceBoundary makes only sense in dimension 2");
 
-  using Rd = TinyVector<Dimension, double>;
+  using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
   const Rd m_direction;
 
  public:
-  template <size_t MeshDimension>
-  friend MeshLineFaceBoundary<MeshDimension> getMeshLineFaceBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
-                                                                     const IBoundaryDescriptor& boundary_descriptor);
+  template <typename MeshTypeT>
+  friend MeshLineFaceBoundary<const MeshTypeT> getMeshLineFaceBoundary(const MeshTypeT& mesh,
+                                                                       const IBoundaryDescriptor& boundary_descriptor);
 
   PUGS_INLINE
-  const Rd& direction() const
+  const Rd&
+  direction() const
   {
     return m_direction;
   }
 
   MeshLineFaceBoundary& operator=(const MeshLineFaceBoundary&) = default;
-  MeshLineFaceBoundary& operator=(MeshLineFaceBoundary&&) = default;
+  MeshLineFaceBoundary& operator=(MeshLineFaceBoundary&&)      = default;
 
  private:
-  MeshLineFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list, const Rd& direction)
-    : MeshFaceBoundary<Dimension>(mesh, ref_face_list), m_direction(direction)
+  MeshLineFaceBoundary(const MeshType& mesh, const RefFaceList& ref_face_list, const Rd& direction)
+    : MeshFaceBoundary(mesh, ref_face_list), m_direction(direction)
   {}
 
  public:
   MeshLineFaceBoundary()                            = default;
   MeshLineFaceBoundary(const MeshLineFaceBoundary&) = default;
-  MeshLineFaceBoundary(MeshLineFaceBoundary &&)     = default;
+  MeshLineFaceBoundary(MeshLineFaceBoundary&&)      = default;
   ~MeshLineFaceBoundary()                           = default;
 };
 
-template <size_t Dimension>
-MeshLineFaceBoundary<Dimension> getMeshLineFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh,
-                                                        const IBoundaryDescriptor& boundary_descriptor);
+template <typename MeshType>
+MeshLineFaceBoundary<const MeshType> getMeshLineFaceBoundary(const MeshType& mesh,
+                                                             const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_LINE_FACE_BOUNDARY_HPP
diff --git a/src/mesh/MeshLineNodeBoundary.cpp b/src/mesh/MeshLineNodeBoundary.cpp
index 1b126f6207719d648415edbc41be6bb83c36abe7..a03c6102ba1ceeb37c35b75f3cd02040313b569d 100644
--- a/src/mesh/MeshLineNodeBoundary.cpp
+++ b/src/mesh/MeshLineNodeBoundary.cpp
@@ -2,16 +2,17 @@
 
 #include <mesh/Connectivity.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/MeshNodeBoundaryUtils.hpp>
 #include <utils/Messenger.hpp>
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-MeshLineNodeBoundary<Dimension>::_checkBoundaryIsLine(const TinyVector<Dimension, double>& direction,
-                                                      const TinyVector<Dimension, double>& origin,
-                                                      const double length,
-                                                      const Mesh<Connectivity<Dimension>>& mesh) const
+MeshLineNodeBoundary<MeshType>::_checkBoundaryIsLine(const TinyVector<MeshType::Dimension, double>& direction,
+                                                     const TinyVector<MeshType::Dimension, double>& origin,
+                                                     const double length,
+                                                     const MeshType& mesh) const
 {
-  using Rdxd = TinyMatrix<Dimension>;
+  using Rdxd = TinyMatrix<MeshType::Dimension>;
 
   const NodeValue<const Rd>& xr = mesh.xr();
 
@@ -35,14 +36,13 @@ MeshLineNodeBoundary<Dimension>::_checkBoundaryIsLine(const TinyVector<Dimension
   }
 }
 
-template <>
 template <>
 TinyVector<2>
-MeshLineNodeBoundary<2>::_getDirection(const Mesh<Connectivity<2>>& mesh)
+MeshLineNodeBoundary<const Mesh<Connectivity<2>>>::_getDirection(const Mesh<Connectivity<2>>& mesh)
 {
   using R2 = TinyVector<2, double>;
 
-  std::array<R2, 2> bounds = this->_getBounds(mesh);
+  std::array<R2, 2> bounds = getBounds(mesh, m_ref_node_list);
 
   const R2& xmin = bounds[0];
   const R2& xmax = bounds[1];
@@ -63,14 +63,13 @@ MeshLineNodeBoundary<2>::_getDirection(const Mesh<Connectivity<2>>& mesh)
   return direction;
 }
 
-template <>
 template <>
 TinyVector<3>
-MeshLineNodeBoundary<3>::_getDirection(const Mesh<Connectivity<3>>& mesh)
+MeshLineNodeBoundary<const Mesh<Connectivity<3>>>::_getDirection(const Mesh<Connectivity<3>>& mesh)
 {
   using R3 = TinyVector<3, double>;
 
-  std::array<R3, 6> bounds = this->_getBounds(mesh);
+  std::array<R3, 6> bounds = getBounds(mesh, m_ref_node_list);
 
   const R3& xmin = bounds[0];
   const R3& ymin = bounds[1];
@@ -106,16 +105,16 @@ MeshLineNodeBoundary<3>::_getDirection(const Mesh<Connectivity<3>>& mesh)
   return direction;
 }
 
-template <size_t Dimension>
-MeshLineNodeBoundary<Dimension>
-getMeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+template <typename MeshType>
+MeshLineNodeBoundary<const MeshType>
+getMeshLineNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
   for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
        ++i_ref_node_list) {
     const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
     const RefId& ref          = ref_node_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshLineNodeBoundary<Dimension>{mesh, ref_node_list};
+      return MeshLineNodeBoundary<const MeshType>{mesh, ref_node_list};
     }
   }
   for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<ItemType::edge>();
@@ -123,7 +122,7 @@ getMeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBounda
     const auto& ref_edge_list = mesh.connectivity().template refItemList<ItemType::edge>(i_ref_edge_list);
     const RefId& ref          = ref_edge_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshLineNodeBoundary<Dimension>{mesh, ref_edge_list};
+      return MeshLineNodeBoundary<const MeshType>{mesh, ref_edge_list};
     }
   }
   for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
@@ -131,7 +130,7 @@ getMeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBounda
     const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
     const RefId& ref          = ref_face_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshLineNodeBoundary<Dimension>{mesh, ref_face_list};
+      return MeshLineNodeBoundary<const MeshType>{mesh, ref_face_list};
     }
   }
 
@@ -141,5 +140,7 @@ getMeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBounda
   throw NormalError(ost.str());
 }
 
-template MeshLineNodeBoundary<2> getMeshLineNodeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
-template MeshLineNodeBoundary<3> getMeshLineNodeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
+template MeshLineNodeBoundary<const Mesh<Connectivity<2>>> getMeshLineNodeBoundary(const Mesh<Connectivity<2>>&,
+                                                                                   const IBoundaryDescriptor&);
+template MeshLineNodeBoundary<const Mesh<Connectivity<3>>> getMeshLineNodeBoundary(const Mesh<Connectivity<3>>&,
+                                                                                   const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshLineNodeBoundary.hpp b/src/mesh/MeshLineNodeBoundary.hpp
index 730b9f536bbf5563d656b59c0d246ff71193ea57..50dbfa00b7e1bad64882562876ff40bb37adff77 100644
--- a/src/mesh/MeshLineNodeBoundary.hpp
+++ b/src/mesh/MeshLineNodeBoundary.hpp
@@ -4,60 +4,61 @@
 #include <algebra/TinyMatrix.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
 
-template <size_t Dimension>
-class [[nodiscard]] MeshLineNodeBoundary final
-  : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
+template <typename MeshType>
+class [[nodiscard]] MeshLineNodeBoundary final : public MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
-  static_assert(Dimension > 1, "MeshLineNodeBoundary makes only sense in dimension greater than 1");
+  static_assert(MeshType::Dimension > 1, "MeshLineNodeBoundary makes only sense in dimension greater than 1");
 
-  using Rd = TinyVector<Dimension, double>;
+  using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
   const Rd m_direction;
 
-  template <size_t MeshDimension>
-  TinyVector<MeshDimension> _getDirection(const Mesh<Connectivity<MeshDimension>>&);
+  TinyVector<MeshType::Dimension> _getDirection(const MeshType&);
 
-  void _checkBoundaryIsLine(const TinyVector<Dimension, double>& direction, const TinyVector<Dimension, double>& origin,
-                            const double length, const Mesh<Connectivity<Dimension>>& mesh) const;
+  void _checkBoundaryIsLine(const TinyVector<MeshType::Dimension, double>& direction,
+                            const TinyVector<MeshType::Dimension, double>& origin,
+                            const double length,
+                            const MeshType& mesh) const;
 
  public:
-  template <size_t MeshDimension>
-  friend MeshLineNodeBoundary<MeshDimension> getMeshLineNodeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
-                                                                     const IBoundaryDescriptor& boundary_descriptor);
+  template <typename MeshTypeT>
+  friend MeshLineNodeBoundary<const MeshTypeT> getMeshLineNodeBoundary(const MeshTypeT& mesh,
+                                                                       const IBoundaryDescriptor& boundary_descriptor);
 
   PUGS_INLINE
-  const Rd& direction() const
+  const Rd&
+  direction() const
   {
     return m_direction;
   }
 
   MeshLineNodeBoundary& operator=(const MeshLineNodeBoundary&) = default;
-  MeshLineNodeBoundary& operator=(MeshLineNodeBoundary&&) = default;
+  MeshLineNodeBoundary& operator=(MeshLineNodeBoundary&&)      = default;
 
  private:
-  MeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefEdgeList& ref_edge_list)
-    : MeshNodeBoundary<Dimension>(mesh, ref_edge_list), m_direction(_getDirection(mesh))
+  MeshLineNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list)
+    : MeshNodeBoundary(mesh, ref_edge_list), m_direction(_getDirection(mesh))
   {}
 
-  MeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list)
-    : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_direction(_getDirection(mesh))
+  MeshLineNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
+    : MeshNodeBoundary(mesh, ref_face_list), m_direction(_getDirection(mesh))
   {}
 
-  MeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefNodeList& ref_node_list)
-    : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_direction(_getDirection(mesh))
+  MeshLineNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
+    : MeshNodeBoundary(mesh, ref_node_list), m_direction(_getDirection(mesh))
   {}
 
  public:
   MeshLineNodeBoundary()                            = default;
   MeshLineNodeBoundary(const MeshLineNodeBoundary&) = default;
-  MeshLineNodeBoundary(MeshLineNodeBoundary &&)     = default;
+  MeshLineNodeBoundary(MeshLineNodeBoundary&&)      = default;
   ~MeshLineNodeBoundary()                           = default;
 };
 
-template <size_t Dimension>
-MeshLineNodeBoundary<Dimension> getMeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
-                                                        const IBoundaryDescriptor& boundary_descriptor);
+template <typename MeshType>
+MeshLineNodeBoundary<const MeshType> getMeshLineNodeBoundary(const MeshType& mesh,
+                                                             const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_LINE_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshNodeBoundary.cpp b/src/mesh/MeshNodeBoundary.cpp
index e883aa4eeb80b40962218208fc044f964a7d4a8d..97f32e9d433553e4d1826f21be08a46e48a91a78 100644
--- a/src/mesh/MeshNodeBoundary.cpp
+++ b/src/mesh/MeshNodeBoundary.cpp
@@ -5,175 +5,11 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <>
-std::array<TinyVector<2>, 2>
-MeshNodeBoundary<2>::_getBounds(const Mesh<Connectivity<2>>& mesh) const
+template <typename MeshType>
+MeshNodeBoundary::MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
 {
-  using R2 = TinyVector<2, double>;
+  constexpr size_t Dimension = MeshType::Dimension;
 
-  const NodeValue<const R2>& xr = mesh.xr();
-
-  std::array<R2, 2> bounds;
-  R2& xmin = bounds[0];
-  R2& xmax = bounds[1];
-
-  xmin = R2{std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-  xmax = R2{-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()};
-
-  auto update_xmin = [](const R2& x, R2& x_min) {
-    if ((x[0] < x_min[0]) or ((x[0] == x_min[0]) and (x[1] < x_min[1]))) {
-      x_min = x;
-    }
-  };
-
-  auto update_xmax = [](const R2& x, R2& x_max) {
-    if ((x[0] > x_max[0]) or ((x[0] == x_max[0]) and (x[1] > x_max[1]))) {
-      x_max = x;
-    }
-  };
-
-  auto node_list = m_ref_node_list.list();
-  for (size_t r = 0; r < node_list.size(); ++r) {
-    const R2& x = xr[node_list[r]];
-    update_xmin(x, xmin);
-    update_xmax(x, xmax);
-  }
-
-  if (parallel::size() > 1) {
-    Array<R2> xmin_array = parallel::allGather(xmin);
-    Array<R2> xmax_array = parallel::allGather(xmax);
-    for (size_t i = 0; i < xmin_array.size(); ++i) {
-      update_xmin(xmin_array[i], xmin);
-    }
-    for (size_t i = 0; i < xmax_array.size(); ++i) {
-      update_xmax(xmax_array[i], xmax);
-    }
-  }
-
-  return bounds;
-}
-
-template <>
-std::array<TinyVector<3>, 6>
-MeshNodeBoundary<3>::_getBounds(const Mesh<Connectivity<3>>& mesh) const
-{
-  using R3 = TinyVector<3, double>;
-
-  auto update_xmin = [](const R3& x, R3& xmin) {
-    // XMIN: X.xmin X.ymax X.zmax
-    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] > xmin[1])) or
-        ((x[0] == xmin[0]) and (x[1] == xmin[1]) and (x[2] > xmin[2]))) {
-      xmin = x;
-    }
-  };
-
-  auto update_xmax = [](const R3& x, R3& xmax) {
-    // XMAX: X.xmax X.ymin X.zmin
-    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] < xmax[1])) or
-        ((x[0] == xmax[0]) and (x[1] == xmax[1]) and (x[2] < xmax[2]))) {
-      xmax = x;
-    }
-  };
-
-  auto update_ymin = [](const R3& x, R3& ymin) {
-    // YMIN: X.ymin X.zmax X.xmin
-    if ((x[1] < ymin[1]) or ((x[1] == ymin[1]) and (x[2] > ymin[2])) or
-        ((x[1] == ymin[1]) and (x[2] == ymin[2]) and (x[0] < ymin[0]))) {
-      ymin = x;
-    }
-  };
-
-  auto update_ymax = [](const R3& x, R3& ymax) {
-    // YMAX: X.ymax X.zmin X.xmax
-    if ((x[1] > ymax[1]) or ((x[1] == ymax[1]) and (x[2] < ymax[2])) or
-        ((x[1] == ymax[1]) and (x[2] == ymax[2]) and (x[0] > ymax[0]))) {
-      ymax = x;
-    }
-  };
-
-  auto update_zmin = [](const R3& x, R3& zmin) {
-    // ZMIN: X.zmin X.xmin X.ymin
-    if ((x[2] < zmin[2]) or ((x[2] == zmin[2]) and (x[0] < zmin[0])) or
-        ((x[2] == zmin[2]) and (x[0] == zmin[0]) and (x[1] < zmin[1]))) {
-      zmin = x;
-    }
-  };
-
-  auto update_zmax = [](const R3& x, R3& zmax) {
-    // ZMAX: X.zmax X.xmax X.ymax
-    if ((x[2] > zmax[2]) or ((x[2] == zmax[2]) and (x[0] > zmax[0])) or
-        ((x[2] == zmax[2]) and (x[0] == zmax[0]) and (x[1] > zmax[1]))) {
-      zmax = x;
-    }
-  };
-
-  std::array<R3, 6> bounds;
-  R3& xmin = bounds[0];
-  R3& ymin = bounds[1];
-  R3& zmin = bounds[2];
-  R3& xmax = bounds[3];
-  R3& ymax = bounds[4];
-  R3& zmax = bounds[5];
-
-  xmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-  ymin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-  zmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-
-  xmax =
-    -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-  ymax =
-    -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-  zmax =
-    -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-
-  const NodeValue<const R3>& xr = mesh.xr();
-
-  auto node_list = m_ref_node_list.list();
-  for (size_t r = 0; r < node_list.size(); ++r) {
-    const R3& x = xr[node_list[r]];
-    update_xmin(x, xmin);
-    update_ymin(x, ymin);
-    update_zmin(x, zmin);
-    update_xmax(x, xmax);
-    update_ymax(x, ymax);
-    update_zmax(x, zmax);
-  }
-
-  if (parallel::size() > 1) {
-    Array<const R3> xmin_array = parallel::allGather(xmin);
-    Array<const R3> ymin_array = parallel::allGather(ymin);
-    Array<const R3> zmin_array = parallel::allGather(zmin);
-    Array<const R3> xmax_array = parallel::allGather(xmax);
-    Array<const R3> ymax_array = parallel::allGather(ymax);
-    Array<const R3> zmax_array = parallel::allGather(zmax);
-
-    for (size_t i = 0; i < xmin_array.size(); ++i) {
-      update_xmin(xmin_array[i], xmin);
-    }
-    for (size_t i = 0; i < ymin_array.size(); ++i) {
-      update_ymin(ymin_array[i], ymin);
-    }
-    for (size_t i = 0; i < zmin_array.size(); ++i) {
-      update_zmin(zmin_array[i], zmin);
-    }
-    for (size_t i = 0; i < xmax_array.size(); ++i) {
-      update_xmax(xmax_array[i], xmax);
-    }
-    for (size_t i = 0; i < ymax_array.size(); ++i) {
-      update_ymax(ymax_array[i], ymax);
-    }
-    for (size_t i = 0; i < zmax_array.size(); ++i) {
-      update_zmax(zmax_array[i], zmax);
-    }
-  }
-
-  return bounds;
-}
-
-template <size_t Dimension>
-MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
-                                              const RefFaceList& ref_face_list)
-{
   const Array<const FaceId>& face_list = ref_face_list.list();
   if (ref_face_list.type() != RefItemListBase::Type::boundary) {
     std::ostringstream ost;
@@ -216,10 +52,11 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>
   const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_node_list);
 }
 
-template <size_t Dimension>
-MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
-                                              const RefEdgeList& ref_edge_list)
+template <typename MeshType>
+MeshNodeBoundary::MeshNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list)
 {
+  constexpr size_t Dimension = MeshType::Dimension;
+
   const Array<const EdgeId>& edge_list = ref_edge_list.list();
   if (ref_edge_list.type() != RefItemListBase::Type::boundary) {
     std::ostringstream ost;
@@ -261,9 +98,8 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>
   const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_node_list);
 }
 
-template <size_t Dimension>
-MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>>&, const RefNodeList& ref_node_list)
-  : m_ref_node_list(ref_node_list)
+template <typename MeshType>
+MeshNodeBoundary::MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list) : m_ref_node_list(ref_node_list)
 {
   if (ref_node_list.type() != RefItemListBase::Type::boundary) {
     std::ostringstream ost;
@@ -273,28 +109,28 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>
   }
 }
 
-template MeshNodeBoundary<1>::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefFaceList&);
-template MeshNodeBoundary<2>::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&);
-template MeshNodeBoundary<3>::MeshNodeBoundary(const Mesh<Connectivity<3>>&, const RefFaceList&);
+template MeshNodeBoundary::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefFaceList&);
+template MeshNodeBoundary::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&);
+template MeshNodeBoundary::MeshNodeBoundary(const Mesh<Connectivity<3>>&, const RefFaceList&);
 
-template MeshNodeBoundary<1>::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefEdgeList&);
-template MeshNodeBoundary<2>::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefEdgeList&);
-template MeshNodeBoundary<3>::MeshNodeBoundary(const Mesh<Connectivity<3>>&, const RefEdgeList&);
+template MeshNodeBoundary::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefEdgeList&);
+template MeshNodeBoundary::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefEdgeList&);
+template MeshNodeBoundary::MeshNodeBoundary(const Mesh<Connectivity<3>>&, const RefEdgeList&);
 
-template MeshNodeBoundary<1>::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefNodeList&);
-template MeshNodeBoundary<2>::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefNodeList&);
-template MeshNodeBoundary<3>::MeshNodeBoundary(const Mesh<Connectivity<3>>&, const RefNodeList&);
+template MeshNodeBoundary::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefNodeList&);
+template MeshNodeBoundary::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefNodeList&);
+template MeshNodeBoundary::MeshNodeBoundary(const Mesh<Connectivity<3>>&, const RefNodeList&);
 
-template <size_t Dimension>
-MeshNodeBoundary<Dimension>
-getMeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+template <typename MeshType>
+MeshNodeBoundary
+getMeshNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
   for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
        ++i_ref_node_list) {
     const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
     const RefId& ref          = ref_node_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshNodeBoundary<Dimension>{mesh, ref_node_list};
+      return MeshNodeBoundary{mesh, ref_node_list};
     }
   }
   for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<ItemType::edge>();
@@ -302,7 +138,7 @@ getMeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDe
     const auto& ref_edge_list = mesh.connectivity().template refItemList<ItemType::edge>(i_ref_edge_list);
     const RefId& ref          = ref_edge_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshNodeBoundary<Dimension>{mesh, ref_edge_list};
+      return MeshNodeBoundary{mesh, ref_edge_list};
     }
   }
   for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
@@ -310,7 +146,7 @@ getMeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDe
     const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
     const RefId& ref          = ref_face_list.refId();
     if (ref == boundary_descriptor) {
-      return MeshNodeBoundary<Dimension>{mesh, ref_face_list};
+      return MeshNodeBoundary{mesh, ref_face_list};
     }
   }
 
@@ -320,6 +156,6 @@ getMeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDe
   throw NormalError(ost.str());
 }
 
-template MeshNodeBoundary<1> getMeshNodeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
-template MeshNodeBoundary<2> getMeshNodeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
-template MeshNodeBoundary<3> getMeshNodeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
+template MeshNodeBoundary getMeshNodeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
+template MeshNodeBoundary getMeshNodeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
+template MeshNodeBoundary getMeshNodeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index ce56579e853171e4b5be0d1f195936a75b838f78..0e869e2eef5b9783aba66840e050071cc3ddf57d 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -6,56 +6,49 @@
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
-template <size_t Dimension>
-class Connectivity;
-
-template <typename ConnectivityType>
-class Mesh;
-
-template <size_t Dimension>
 class [[nodiscard]] MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 {
  protected:
   RefNodeList m_ref_node_list;
 
-  std::array<TinyVector<Dimension>, Dimension*(Dimension - 1)> _getBounds(const Mesh<Connectivity<Dimension>>& mesh)
-    const;
-
  public:
-  template <size_t MeshDimension>
-  friend MeshNodeBoundary<MeshDimension> getMeshNodeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
-                                                             const IBoundaryDescriptor& boundary_descriptor);
+  template <typename MeshTypeT>
+  friend MeshNodeBoundary getMeshNodeBoundary(const MeshTypeT& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
   MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
-  MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default;
+  MeshNodeBoundary& operator=(MeshNodeBoundary&&)      = default;
 
   PUGS_INLINE
-  const RefNodeList& refNodeList() const
+  const RefNodeList&
+  refNodeList() const
   {
     return m_ref_node_list;
   }
 
   PUGS_INLINE
-  const Array<const NodeId>& nodeList() const
+  const Array<const NodeId>&
+  nodeList() const
   {
     return m_ref_node_list.list();
   }
 
  protected:
-  MeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list);
-  MeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefEdgeList& ref_edge_list);
-  MeshNodeBoundary(const Mesh<Connectivity<Dimension>>&, const RefNodeList& ref_node_list);
+  template <typename MeshType>
+  MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list);
+  template <typename MeshType>
+  MeshNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list);
+  template <typename MeshType>
+  MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list);
 
  public:
   MeshNodeBoundary(const MeshNodeBoundary&) = default;
-  MeshNodeBoundary(MeshNodeBoundary &&)     = default;
+  MeshNodeBoundary(MeshNodeBoundary&&)      = default;
 
   MeshNodeBoundary()          = default;
   virtual ~MeshNodeBoundary() = default;
 };
 
-template <size_t Dimension>
-MeshNodeBoundary<Dimension> getMeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
-                                                const IBoundaryDescriptor& boundary_descriptor);
+template <typename MeshType>
+MeshNodeBoundary getMeshNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshNodeBoundaryUtils.cpp b/src/mesh/MeshNodeBoundaryUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..63df3e36475d2d60029fe5cd27dba9f37eac21cc
--- /dev/null
+++ b/src/mesh/MeshNodeBoundaryUtils.cpp
@@ -0,0 +1,170 @@
+#include <mesh/MeshNodeBoundaryUtils.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/Messenger.hpp>
+
+template <>
+std::array<TinyVector<2>, 2>
+getBounds(const Mesh<Connectivity<2>>& mesh, const RefNodeList& ref_node_list)
+{
+  using R2 = TinyVector<2, double>;
+
+  const NodeValue<const R2>& xr = mesh.xr();
+
+  std::array<R2, 2> bounds;
+  R2& xmin = bounds[0];
+  R2& xmax = bounds[1];
+
+  xmin = R2{std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  xmax = R2{-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()};
+
+  auto update_xmin = [](const R2& x, R2& x_min) {
+    if ((x[0] < x_min[0]) or ((x[0] == x_min[0]) and (x[1] < x_min[1]))) {
+      x_min = x;
+    }
+  };
+
+  auto update_xmax = [](const R2& x, R2& x_max) {
+    if ((x[0] > x_max[0]) or ((x[0] == x_max[0]) and (x[1] > x_max[1]))) {
+      x_max = x;
+    }
+  };
+
+  auto node_list = ref_node_list.list();
+  for (size_t r = 0; r < node_list.size(); ++r) {
+    const R2& x = xr[node_list[r]];
+    update_xmin(x, xmin);
+    update_xmax(x, xmax);
+  }
+
+  if (parallel::size() > 1) {
+    Array<R2> xmin_array = parallel::allGather(xmin);
+    Array<R2> xmax_array = parallel::allGather(xmax);
+    for (size_t i = 0; i < xmin_array.size(); ++i) {
+      update_xmin(xmin_array[i], xmin);
+    }
+    for (size_t i = 0; i < xmax_array.size(); ++i) {
+      update_xmax(xmax_array[i], xmax);
+    }
+  }
+
+  return bounds;
+}
+
+template <>
+std::array<TinyVector<3>, 6>
+getBounds(const Mesh<Connectivity<3>>& mesh, const RefNodeList& ref_node_list)
+{
+  using R3 = TinyVector<3, double>;
+
+  auto update_xmin = [](const R3& x, R3& xmin) {
+    // XMIN: X.xmin X.ymax X.zmax
+    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] > xmin[1])) or
+        ((x[0] == xmin[0]) and (x[1] == xmin[1]) and (x[2] > xmin[2]))) {
+      xmin = x;
+    }
+  };
+
+  auto update_xmax = [](const R3& x, R3& xmax) {
+    // XMAX: X.xmax X.ymin X.zmin
+    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] < xmax[1])) or
+        ((x[0] == xmax[0]) and (x[1] == xmax[1]) and (x[2] < xmax[2]))) {
+      xmax = x;
+    }
+  };
+
+  auto update_ymin = [](const R3& x, R3& ymin) {
+    // YMIN: X.ymin X.zmax X.xmin
+    if ((x[1] < ymin[1]) or ((x[1] == ymin[1]) and (x[2] > ymin[2])) or
+        ((x[1] == ymin[1]) and (x[2] == ymin[2]) and (x[0] < ymin[0]))) {
+      ymin = x;
+    }
+  };
+
+  auto update_ymax = [](const R3& x, R3& ymax) {
+    // YMAX: X.ymax X.zmin X.xmax
+    if ((x[1] > ymax[1]) or ((x[1] == ymax[1]) and (x[2] < ymax[2])) or
+        ((x[1] == ymax[1]) and (x[2] == ymax[2]) and (x[0] > ymax[0]))) {
+      ymax = x;
+    }
+  };
+
+  auto update_zmin = [](const R3& x, R3& zmin) {
+    // ZMIN: X.zmin X.xmin X.ymin
+    if ((x[2] < zmin[2]) or ((x[2] == zmin[2]) and (x[0] < zmin[0])) or
+        ((x[2] == zmin[2]) and (x[0] == zmin[0]) and (x[1] < zmin[1]))) {
+      zmin = x;
+    }
+  };
+
+  auto update_zmax = [](const R3& x, R3& zmax) {
+    // ZMAX: X.zmax X.xmax X.ymax
+    if ((x[2] > zmax[2]) or ((x[2] == zmax[2]) and (x[0] > zmax[0])) or
+        ((x[2] == zmax[2]) and (x[0] == zmax[0]) and (x[1] > zmax[1]))) {
+      zmax = x;
+    }
+  };
+
+  std::array<R3, 6> bounds;
+  R3& xmin = bounds[0];
+  R3& ymin = bounds[1];
+  R3& zmin = bounds[2];
+  R3& xmax = bounds[3];
+  R3& ymax = bounds[4];
+  R3& zmax = bounds[5];
+
+  xmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  ymin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  zmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+
+  xmax =
+    -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  ymax =
+    -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  zmax =
+    -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+
+  const NodeValue<const R3>& xr = mesh.xr();
+
+  auto node_list = ref_node_list.list();
+  for (size_t r = 0; r < node_list.size(); ++r) {
+    const R3& x = xr[node_list[r]];
+    update_xmin(x, xmin);
+    update_ymin(x, ymin);
+    update_zmin(x, zmin);
+    update_xmax(x, xmax);
+    update_ymax(x, ymax);
+    update_zmax(x, zmax);
+  }
+
+  if (parallel::size() > 1) {
+    Array<const R3> xmin_array = parallel::allGather(xmin);
+    Array<const R3> ymin_array = parallel::allGather(ymin);
+    Array<const R3> zmin_array = parallel::allGather(zmin);
+    Array<const R3> xmax_array = parallel::allGather(xmax);
+    Array<const R3> ymax_array = parallel::allGather(ymax);
+    Array<const R3> zmax_array = parallel::allGather(zmax);
+
+    for (size_t i = 0; i < xmin_array.size(); ++i) {
+      update_xmin(xmin_array[i], xmin);
+    }
+    for (size_t i = 0; i < ymin_array.size(); ++i) {
+      update_ymin(ymin_array[i], ymin);
+    }
+    for (size_t i = 0; i < zmin_array.size(); ++i) {
+      update_zmin(zmin_array[i], zmin);
+    }
+    for (size_t i = 0; i < xmax_array.size(); ++i) {
+      update_xmax(xmax_array[i], xmax);
+    }
+    for (size_t i = 0; i < ymax_array.size(); ++i) {
+      update_ymax(ymax_array[i], ymax);
+    }
+    for (size_t i = 0; i < zmax_array.size(); ++i) {
+      update_zmax(zmax_array[i], zmax);
+    }
+  }
+
+  return bounds;
+}
diff --git a/src/mesh/MeshNodeBoundaryUtils.hpp b/src/mesh/MeshNodeBoundaryUtils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a62f9d8f3affa0cf2a840e69776b91c88fa8ec0e
--- /dev/null
+++ b/src/mesh/MeshNodeBoundaryUtils.hpp
@@ -0,0 +1,14 @@
+#ifndef MESH_NODE_BOUNDARY_UTILS_HPP
+#define MESH_NODE_BOUNDARY_UTILS_HPP
+
+#include <algebra/TinyVector.hpp>
+#include <mesh/RefItemList.hpp>
+
+#include <array>
+
+template <typename MeshType>
+std::array<TinyVector<MeshType::Dimension>, MeshType::Dimension*(MeshType::Dimension - 1)> getBounds(
+  const MeshType& mesh,
+  const RefNodeList& ref_node_list);
+
+#endif   // MESH_NODE_BOUNDARY_UTILS_HPP
diff --git a/src/mesh/MeshRandomizer.cpp b/src/mesh/MeshRandomizer.cpp
index 57b826dff5ba409e69816b352c572b67547217e1..9a2a56d5a754d96a195cb88a6e70ac5f0d5c6dab 100644
--- a/src/mesh/MeshRandomizer.cpp
+++ b/src/mesh/MeshRandomizer.cpp
@@ -16,14 +16,15 @@
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <utils/RandomEngine.hpp>
 
-template <size_t Dimension>
+template <typename MeshType>
 class MeshRandomizerHandler::MeshRandomizer
 {
  private:
+  static constexpr size_t Dimension = MeshType::Dimension;
+
   using Rd               = TinyVector<Dimension>;
   using Rdxd             = TinyMatrix<Dimension>;
-  using ConnectivityType = Connectivity<Dimension>;
-  using MeshType         = Mesh<ConnectivityType>;
+  using ConnectivityType = typename MeshType::Connectivity;
 
   const MeshType& m_given_mesh;
 
@@ -273,14 +274,14 @@ class MeshRandomizerHandler::MeshRandomizer
   ~MeshRandomizer() = default;
 };
 
-template <size_t Dimension>
-class MeshRandomizerHandler::MeshRandomizer<Dimension>::AxisBoundaryCondition
+template <typename MeshType>
+class MeshRandomizerHandler::MeshRandomizer<MeshType>::AxisBoundaryCondition
 {
  public:
   using Rd = TinyVector<Dimension, double>;
 
  private:
-  const MeshLineNodeBoundary<Dimension> m_mesh_line_node_boundary;
+  const MeshLineNodeBoundary<const MeshType> m_mesh_line_node_boundary;
 
  public:
   const Rd&
@@ -295,7 +296,7 @@ class MeshRandomizerHandler::MeshRandomizer<Dimension>::AxisBoundaryCondition
     return m_mesh_line_node_boundary.nodeList();
   }
 
-  AxisBoundaryCondition(MeshLineNodeBoundary<Dimension>&& mesh_line_node_boundary)
+  AxisBoundaryCondition(MeshLineNodeBoundary<const MeshType>&& mesh_line_node_boundary)
     : m_mesh_line_node_boundary(mesh_line_node_boundary)
   {
     ;
@@ -305,18 +306,18 @@ class MeshRandomizerHandler::MeshRandomizer<Dimension>::AxisBoundaryCondition
 };
 
 template <>
-class MeshRandomizerHandler::MeshRandomizer<1>::AxisBoundaryCondition
+class MeshRandomizerHandler::MeshRandomizer<Mesh<Connectivity<1>>>::AxisBoundaryCondition
 {
  public:
   AxisBoundaryCondition()  = default;
   ~AxisBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class MeshRandomizerHandler::MeshRandomizer<Dimension>::FixedBoundaryCondition
+template <typename MeshType>
+class MeshRandomizerHandler::MeshRandomizer<MeshType>::FixedBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
 
  public:
   const Array<const NodeId>&
@@ -325,19 +326,19 @@ class MeshRandomizerHandler::MeshRandomizer<Dimension>::FixedBoundaryCondition
     return m_mesh_node_boundary.nodeList();
   }
 
-  FixedBoundaryCondition(MeshNodeBoundary<Dimension>&& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
+  FixedBoundaryCondition(MeshNodeBoundary&& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
 
   ~FixedBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class MeshRandomizerHandler::MeshRandomizer<Dimension>::SymmetryBoundaryCondition
+template <typename MeshType>
+class MeshRandomizerHandler::MeshRandomizer<MeshType>::SymmetryBoundaryCondition
 {
  public:
-  using Rd = TinyVector<Dimension, double>;
+  using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
-  const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+  const MeshFlatNodeBoundary<const MeshType> m_mesh_flat_node_boundary;
 
  public:
   const Rd&
@@ -352,7 +353,7 @@ class MeshRandomizerHandler::MeshRandomizer<Dimension>::SymmetryBoundaryConditio
     return m_mesh_flat_node_boundary.nodeList();
   }
 
-  SymmetryBoundaryCondition(MeshFlatNodeBoundary<Dimension>&& mesh_flat_node_boundary)
+  SymmetryBoundaryCondition(MeshFlatNodeBoundary<const MeshType>&& mesh_flat_node_boundary)
     : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
   {
     ;
diff --git a/src/mesh/MeshRandomizer.hpp b/src/mesh/MeshRandomizer.hpp
index b71d99a676c8c18201871efa11902eeaab63f742..048e04bfdad1780fc41362595ca0378400a6266e 100644
--- a/src/mesh/MeshRandomizer.hpp
+++ b/src/mesh/MeshRandomizer.hpp
@@ -11,7 +11,7 @@ class FunctionSymbolId;
 class MeshRandomizerHandler
 {
  private:
-  template <size_t Dimension>
+  template <typename MeshType>
   class MeshRandomizer;
 
  public:
diff --git a/src/mesh/MeshSmoother.cpp b/src/mesh/MeshSmoother.cpp
index 281548bf1401c3c33f209a582d2d8d33a9b99bb2..d90b95a4d9356e441155b844bc33eeb754aaee62 100644
--- a/src/mesh/MeshSmoother.cpp
+++ b/src/mesh/MeshSmoother.cpp
@@ -20,14 +20,15 @@
 
 #include <variant>
 
-template <size_t Dimension>
+template <typename MeshType>
 class MeshSmootherHandler::MeshSmoother
 {
  private:
+  static constexpr size_t Dimension = MeshType::Dimension;
+
   using Rd               = TinyVector<Dimension>;
   using Rdxd             = TinyMatrix<Dimension>;
-  using ConnectivityType = Connectivity<Dimension>;
-  using MeshType         = Mesh<ConnectivityType>;
+  using ConnectivityType = typename MeshType::Connectivity;
 
   const MeshType& m_given_mesh;
 
@@ -320,14 +321,14 @@ class MeshSmootherHandler::MeshSmoother
   ~MeshSmoother() = default;
 };
 
-template <size_t Dimension>
-class MeshSmootherHandler::MeshSmoother<Dimension>::AxisBoundaryCondition
+template <typename MeshType>
+class MeshSmootherHandler::MeshSmoother<MeshType>::AxisBoundaryCondition
 {
  public:
   using Rd = TinyVector<Dimension, double>;
 
  private:
-  const MeshLineNodeBoundary<Dimension> m_mesh_line_node_boundary;
+  const MeshLineNodeBoundary<const MeshType> m_mesh_line_node_boundary;
 
  public:
   const Rd&
@@ -342,7 +343,7 @@ class MeshSmootherHandler::MeshSmoother<Dimension>::AxisBoundaryCondition
     return m_mesh_line_node_boundary.nodeList();
   }
 
-  AxisBoundaryCondition(MeshLineNodeBoundary<Dimension>&& mesh_line_node_boundary)
+  AxisBoundaryCondition(MeshLineNodeBoundary<const MeshType>&& mesh_line_node_boundary)
     : m_mesh_line_node_boundary(mesh_line_node_boundary)
   {
     ;
@@ -352,18 +353,18 @@ class MeshSmootherHandler::MeshSmoother<Dimension>::AxisBoundaryCondition
 };
 
 template <>
-class MeshSmootherHandler::MeshSmoother<1>::AxisBoundaryCondition
+class MeshSmootherHandler::MeshSmoother<Mesh<Connectivity<1>>>::AxisBoundaryCondition
 {
  public:
   AxisBoundaryCondition()  = default;
   ~AxisBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class MeshSmootherHandler::MeshSmoother<Dimension>::FixedBoundaryCondition
+template <typename MeshType>
+class MeshSmootherHandler::MeshSmoother<MeshType>::FixedBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
 
  public:
   const Array<const NodeId>&
@@ -372,19 +373,19 @@ class MeshSmootherHandler::MeshSmoother<Dimension>::FixedBoundaryCondition
     return m_mesh_node_boundary.nodeList();
   }
 
-  FixedBoundaryCondition(MeshNodeBoundary<Dimension>&& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
+  FixedBoundaryCondition(MeshNodeBoundary&& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
 
   ~FixedBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class MeshSmootherHandler::MeshSmoother<Dimension>::SymmetryBoundaryCondition
+template <typename MeshType>
+class MeshSmootherHandler::MeshSmoother<MeshType>::SymmetryBoundaryCondition
 {
  public:
-  using Rd = TinyVector<Dimension, double>;
+  using Rd = TinyVector<MeshType::Dimension, double>;
 
  private:
-  const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+  const MeshFlatNodeBoundary<const MeshType> m_mesh_flat_node_boundary;
 
  public:
   const Rd&
@@ -399,7 +400,7 @@ class MeshSmootherHandler::MeshSmoother<Dimension>::SymmetryBoundaryCondition
     return m_mesh_flat_node_boundary.nodeList();
   }
 
-  SymmetryBoundaryCondition(MeshFlatNodeBoundary<Dimension>&& mesh_flat_node_boundary)
+  SymmetryBoundaryCondition(MeshFlatNodeBoundary<const MeshType>&& mesh_flat_node_boundary)
     : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
   {
     ;
diff --git a/src/mesh/MeshSmoother.hpp b/src/mesh/MeshSmoother.hpp
index 209be6c1339495203df264b967e7995d9485bdb3..c5875561a5bce80257dbf6a118517c2b8443c021 100644
--- a/src/mesh/MeshSmoother.hpp
+++ b/src/mesh/MeshSmoother.hpp
@@ -13,7 +13,7 @@ class DiscreteFunctionVariant;
 class MeshSmootherHandler
 {
  private:
-  template <size_t Dimension>
+  template <typename MeshType>
   class MeshSmoother;
 
  public:
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 575a281c19fddd3a74448c725080ca2025c0c4b5..e60372f0c46b1d3eed0852c5fc60486b6f1c611d 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -61,14 +61,15 @@ acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c)
   }
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler::IAcousticSolver
 {
  private:
+  static constexpr size_t Dimension = MeshType::Dimension;
+
   using Rdxd = TinyMatrix<Dimension>;
   using Rd   = TinyVector<Dimension>;
 
-  using MeshType     = Mesh<Connectivity<Dimension>>;
   using MeshDataType = MeshData<Dimension>;
 
   using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, const double>;
@@ -271,8 +272,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
         const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
           dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
         if (dirichlet_bc_descriptor.name() == "velocity") {
-          MeshNodeBoundary<Dimension> mesh_node_boundary =
-            getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
+          MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
 
           Array<const Rd> value_list =
             InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
@@ -284,17 +284,15 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
           const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
 
           if constexpr (Dimension == 1) {
-            MeshNodeBoundary<Dimension> mesh_node_boundary =
-              getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+            MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
 
             Array<const double> node_values =
               InterpolateItemValue<double(Rd)>::template interpolate<ItemType::node>(pressure_id, mesh.xr(),
                                                                                      mesh_node_boundary.nodeList());
-
-            bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
+            PressureBoundaryCondition p(mesh_node_boundary, node_values);
+            bc_list.emplace_back(p);
           } else {
-            MeshFaceBoundary<Dimension> mesh_face_boundary =
-              getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+            MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
 
             MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
             Array<const double> face_values =
@@ -529,11 +527,11 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   ~AcousticSolver()                = default;
 };
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-AcousticSolverHandler::AcousticSolver<Dimension>::_applyPressureBC(const BoundaryConditionList& bc_list,
-                                                                   const MeshType& mesh,
-                                                                   NodeValue<Rd>& br) const
+AcousticSolverHandler::AcousticSolver<MeshType>::_applyPressureBC(const BoundaryConditionList& bc_list,
+                                                                  const MeshType& mesh,
+                                                                  NodeValue<Rd>& br) const
 {
   for (const auto& boundary_condition : bc_list) {
     std::visit(
@@ -594,11 +592,11 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applyPressureBC(const Boundar
   }
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-AcousticSolverHandler::AcousticSolver<Dimension>::_applySymmetryBC(const BoundaryConditionList& bc_list,
-                                                                   NodeValue<Rdxd>& Ar,
-                                                                   NodeValue<Rd>& br) const
+AcousticSolverHandler::AcousticSolver<MeshType>::_applySymmetryBC(const BoundaryConditionList& bc_list,
+                                                                  NodeValue<Rdxd>& Ar,
+                                                                  NodeValue<Rd>& br) const
 {
   for (const auto& boundary_condition : bc_list) {
     std::visit(
@@ -625,11 +623,11 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applySymmetryBC(const Boundar
   }
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-AcousticSolverHandler::AcousticSolver<Dimension>::_applyVelocityBC(const BoundaryConditionList& bc_list,
-                                                                   NodeValue<Rdxd>& Ar,
-                                                                   NodeValue<Rd>& br) const
+AcousticSolverHandler::AcousticSolver<MeshType>::_applyVelocityBC(const BoundaryConditionList& bc_list,
+                                                                  NodeValue<Rdxd>& Ar,
+                                                                  NodeValue<Rd>& br) const
 {
   for (const auto& boundary_condition : bc_list) {
     std::visit(
@@ -662,12 +660,12 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applyVelocityBC(const Boundar
   }
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-AcousticSolverHandler::AcousticSolver<Dimension>::_applyExternalVelocityBC(const BoundaryConditionList& bc_list,
-                                                                           const DiscreteScalarFunction& p,
-                                                                           NodeValue<Rdxd>& Ar,
-                                                                           NodeValue<Rd>& br) const
+AcousticSolverHandler::AcousticSolver<MeshType>::_applyExternalVelocityBC(const BoundaryConditionList& bc_list,
+                                                                          const DiscreteScalarFunction& p,
+                                                                          NodeValue<Rdxd>& Ar,
+                                                                          NodeValue<Rd>& br) const
 {
   for (const auto& boundary_condition : bc_list) {
     std::visit(
@@ -691,11 +689,11 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applyExternalVelocityBC(const
   }
 }
 
-template <size_t Dimension>
-class AcousticSolverHandler::AcousticSolver<Dimension>::FixedBoundaryCondition
+template <typename MeshType>
+class AcousticSolverHandler::AcousticSolver<MeshType>::FixedBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
 
  public:
   const Array<const NodeId>&
@@ -704,18 +702,16 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::FixedBoundaryCondition
     return m_mesh_node_boundary.nodeList();
   }
 
-  FixedBoundaryCondition(const MeshNodeBoundary<Dimension> mesh_node_boundary)
-    : m_mesh_node_boundary{mesh_node_boundary}
-  {}
+  FixedBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
 
   ~FixedBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class AcousticSolverHandler::AcousticSolver<Dimension>::PressureBoundaryCondition
+template <typename MeshType>
+class AcousticSolverHandler::AcousticSolver<MeshType>::PressureBoundaryCondition
 {
  private:
-  const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+  const MeshFaceBoundary m_mesh_face_boundary;
   const Array<const double> m_value_list;
 
  public:
@@ -731,8 +727,7 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::PressureBoundaryConditio
     return m_value_list;
   }
 
-  PressureBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary,
-                            const Array<const double>& value_list)
+  PressureBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, const Array<const double>& value_list)
     : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
   {}
 
@@ -740,10 +735,10 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::PressureBoundaryConditio
 };
 
 template <>
-class AcousticSolverHandler::AcousticSolver<1>::PressureBoundaryCondition
+class AcousticSolverHandler::AcousticSolver<const Mesh<Connectivity<1>>>::PressureBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<1> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
   const Array<const double> m_value_list;
 
  public:
@@ -759,19 +754,19 @@ class AcousticSolverHandler::AcousticSolver<1>::PressureBoundaryCondition
     return m_value_list;
   }
 
-  PressureBoundaryCondition(const MeshNodeBoundary<1>& mesh_node_boundary, const Array<const double>& value_list)
+  PressureBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const Array<const double>& value_list)
     : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
   {}
 
   ~PressureBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class AcousticSolverHandler::AcousticSolver<Dimension>::ExternalFSIVelocityBoundaryCondition
+template <typename MeshType>
+class AcousticSolverHandler::AcousticSolver<MeshType>::ExternalFSIVelocityBoundaryCondition
 {
  private:
   const ItemToItemMatrix<ItemType::node, ItemType::cell> m_node_to_cell_matrix;
-  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
   const Array<TinyVector<Dimension>> m_value_list;
   const std::shared_ptr<const Socket> m_socket;
 
@@ -800,7 +795,7 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::ExternalFSIVelocityBound
   }
 
   ExternalFSIVelocityBoundaryCondition(const Mesh<Connectivity<Dimension>>& mesh,
-                                       const MeshNodeBoundary<Dimension>& mesh_node_boundary,
+                                       const MeshNodeBoundary& mesh_node_boundary,
                                        const std::shared_ptr<const Socket>& socket)
     : m_node_to_cell_matrix{mesh.connectivity().nodeToCellMatrix()},
       m_mesh_node_boundary{mesh_node_boundary},
@@ -815,11 +810,11 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::ExternalFSIVelocityBound
   ~ExternalFSIVelocityBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class AcousticSolverHandler::AcousticSolver<Dimension>::VelocityBoundaryCondition
+template <typename MeshType>
+class AcousticSolverHandler::AcousticSolver<MeshType>::VelocityBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
 
   const Array<const TinyVector<Dimension>> m_value_list;
 
@@ -836,7 +831,7 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::VelocityBoundaryConditio
     return m_value_list;
   }
 
-  VelocityBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary,
+  VelocityBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
                             const Array<const TinyVector<Dimension>>& value_list)
     : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
   {}
@@ -844,14 +839,14 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::VelocityBoundaryConditio
   ~VelocityBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class AcousticSolverHandler::AcousticSolver<Dimension>::SymmetryBoundaryCondition
+template <typename MeshType>
+class AcousticSolverHandler::AcousticSolver<MeshType>::SymmetryBoundaryCondition
 {
  public:
   using Rd = TinyVector<Dimension, double>;
 
  private:
-  const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+  const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
 
  public:
   const Rd&
@@ -872,7 +867,7 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::SymmetryBoundaryConditio
     return m_mesh_flat_node_boundary.nodeList();
   }
 
-  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
+  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary)
     : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
   {
     ;
@@ -891,7 +886,7 @@ AcousticSolverHandler::AcousticSolverHandler(const std::shared_ptr<const MeshVar
     [&](auto&& mesh) {
       using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
       if constexpr (is_polygonal_mesh<MeshType>) {
-        m_acoustic_solver = std::make_unique<AcousticSolver<MeshType::Dimension>>();
+        m_acoustic_solver = std::make_unique<AcousticSolver<MeshType>>();
       } else {
         throw NormalError("unexpected mesh type");
       }
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 85d829c1b126cf4ae1aebdc7db598de82fade217..75e08981dfb50993233eb80ddfdf67266a24cf06 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -59,14 +59,14 @@ class AcousticSolverHandler
           const std::shared_ptr<const DiscreteFunctionVariant>& p,
           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
 
-    IAcousticSolver()                  = default;
-    IAcousticSolver(IAcousticSolver&&) = default;
+    IAcousticSolver()                             = default;
+    IAcousticSolver(IAcousticSolver&&)            = default;
     IAcousticSolver& operator=(IAcousticSolver&&) = default;
 
     virtual ~IAcousticSolver() = default;
   };
 
-  template <size_t Dimension>
+  template <typename MeshType>
   class AcousticSolver;
 
   std::unique_ptr<IAcousticSolver> m_acoustic_solver;
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index 01ff407bd35a597a4285ce62ffec530bc77e81fe..25d9cf3de756c2b21f6a85447b4a69359cea4767 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -26,13 +26,13 @@
 #include <variant>
 #include <vector>
 
-template <size_t Dimension>
+template <typename MeshType>
 class FluxingAdvectionSolver
 {
  private:
-  using Rd = TinyVector<Dimension>;
+  static constexpr size_t Dimension = MeshType::Dimension;
 
-  using MeshType     = Mesh<Connectivity<Dimension>>;
+  using Rd           = TinyVector<Dimension>;
   using MeshDataType = MeshData<Dimension>;
 
   const std::shared_ptr<const MeshType> m_old_mesh;
@@ -114,8 +114,7 @@ class FluxingAdvectionSolver
         const InflowBoundaryConditionDescriptor& inflow_bc_descriptor =
           dynamic_cast<const InflowBoundaryConditionDescriptor&>(*i_bc);
 
-        MeshFaceBoundary<Dimension> mesh_face_boundary =
-          getMeshFaceBoundary(mesh, inflow_bc_descriptor.boundaryDescriptor());
+        MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, inflow_bc_descriptor.boundaryDescriptor());
 
         FaceValue<const Rd> xl = MeshDataManager::instance().getMeshData(*m_old_mesh).xl();
         Array<const DataType> value_list =
@@ -158,8 +157,7 @@ class FluxingAdvectionSolver
         const InflowBoundaryConditionDescriptor& inflow_bc_descriptor =
           dynamic_cast<const InflowBoundaryConditionDescriptor&>(*i_bc);
 
-        MeshFaceBoundary<Dimension> mesh_face_boundary =
-          getMeshFaceBoundary(mesh, inflow_bc_descriptor.boundaryDescriptor());
+        MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, inflow_bc_descriptor.boundaryDescriptor());
 
         FaceValue<const Rd> xl = MeshDataManager::instance().getMeshData(*m_old_mesh).xl();
         Table<const double> array_list =
@@ -209,9 +207,9 @@ class FluxingAdvectionSolver
   ~FluxingAdvectionSolver() = default;
 };
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-FluxingAdvectionSolver<Dimension>::_computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes)
+FluxingAdvectionSolver<MeshType>::_computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes)
 {
   m_donnor_cell = [&] {
     const FaceValuePerCell<const bool> cell_face_is_reversed = m_new_mesh->connectivity().cellFaceIsReversed();
@@ -242,7 +240,8 @@ FluxingAdvectionSolver<Dimension>::_computeDonorCells(FaceValue<const double> al
 
 template <>
 void
-FluxingAdvectionSolver<1>::_computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes)
+FluxingAdvectionSolver<const Mesh<Connectivity<1>>>::_computeDonorCells(
+  FaceValue<const double> algebraic_fluxing_volumes)
 {
   m_donnor_cell = [&] {
     const auto face_to_cell_matrix = m_new_mesh->connectivity().faceToCellMatrix();
@@ -276,7 +275,7 @@ FluxingAdvectionSolver<1>::_computeDonorCells(FaceValue<const double> algebraic_
 
 template <>
 FaceValue<double>
-FluxingAdvectionSolver<1>::_computeAlgebraicFluxingVolume()
+FluxingAdvectionSolver<const Mesh<Connectivity<1>>>::_computeAlgebraicFluxingVolume()
 {
   Array<double> fluxing_volumes{m_new_mesh->numberOfNodes()};
   NodeValue<double> nodal_fluxing_volume(m_new_mesh->connectivity(), fluxing_volumes);
@@ -295,7 +294,7 @@ FluxingAdvectionSolver<1>::_computeAlgebraicFluxingVolume()
 
 template <>
 FaceValue<double>
-FluxingAdvectionSolver<2>::_computeAlgebraicFluxingVolume()
+FluxingAdvectionSolver<const Mesh<Connectivity<2>>>::_computeAlgebraicFluxingVolume()
 {
   const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix();
   FaceValue<double> algebraic_fluxing_volume(m_new_mesh->connectivity());
@@ -322,7 +321,7 @@ FluxingAdvectionSolver<2>::_computeAlgebraicFluxingVolume()
 
 template <>
 FaceValue<double>
-FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume()
+FluxingAdvectionSolver<const Mesh<Connectivity<3>>>::_computeAlgebraicFluxingVolume()
 {
   // due to the z component of the jacobian determinant, degree 3
   // polynomials must be exactly integrated
@@ -387,9 +386,9 @@ FluxingAdvectionSolver<3>::_computeAlgebraicFluxingVolume()
   return algebraic_fluxing_volume;
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 FaceValue<double>
-FluxingAdvectionSolver<Dimension>::_computeFluxingVolume(FaceValue<double> algebraic_fluxing_volumes)
+FluxingAdvectionSolver<MeshType>::_computeFluxingVolume(FaceValue<double> algebraic_fluxing_volumes)
 {
   Assert(m_donnor_cell.isBuilt());
   // Now that donnor cells are clearly defined, we consider the
@@ -402,9 +401,9 @@ FluxingAdvectionSolver<Dimension>::_computeFluxingVolume(FaceValue<double> algeb
   return algebraic_fluxing_volumes;
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-FluxingAdvectionSolver<Dimension>::_computeCycleNumber(FaceValue<double> fluxing_volumes)
+FluxingAdvectionSolver<MeshType>::_computeCycleNumber(FaceValue<double> fluxing_volumes)
 {
   const auto cell_to_face_matrix = m_old_mesh->connectivity().cellToFaceMatrix();
 
@@ -445,9 +444,9 @@ FluxingAdvectionSolver<Dimension>::_computeCycleNumber(FaceValue<double> fluxing
   m_cycle_fluxing_volume = fluxing_volumes;
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-FluxingAdvectionSolver<Dimension>::_computeGeometricalData()
+FluxingAdvectionSolver<MeshType>::_computeGeometricalData()
 {
   auto fluxing_volumes = this->_computeAlgebraicFluxingVolume();
   this->_computeDonorCells(fluxing_volumes);
@@ -455,12 +454,12 @@ FluxingAdvectionSolver<Dimension>::_computeGeometricalData()
   this->_computeCycleNumber(fluxing_volumes);
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 template <typename CellDataType>
 void
-FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step_Vj,
-                                             const std::vector<BoundaryConditionVariant>& q_bc_list,
-                                             CellDataType& old_q)
+FluxingAdvectionSolver<MeshType>::_remapOne(const CellValue<const double>& step_Vj,
+                                            const std::vector<BoundaryConditionVariant>& q_bc_list,
+                                            CellDataType& old_q)
 {
   using DataType = std::decay_t<typename CellDataType::data_type>;
 
@@ -635,9 +634,9 @@ FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step
   old_q = new_q;
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-FluxingAdvectionSolver<Dimension>::_remapAllQuantities()
+FluxingAdvectionSolver<MeshType>::_remapAllQuantities()
 {
   const auto cell_to_face_matrix              = m_new_mesh->connectivity().cellToFaceMatrix();
   const auto face_local_number_in_their_cells = m_new_mesh->connectivity().faceLocalNumbersInTheirCells();
@@ -698,9 +697,9 @@ FluxingAdvectionSolver<Dimension>::_remapAllQuantities()
   }
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
-FluxingAdvectionSolver<Dimension>::remap(
+FluxingAdvectionSolver<MeshType>::remap(
   const std::vector<std::shared_ptr<const VariableBCDescriptor>>& quantity_list_with_bc)
 {
   m_boundary_condition_list.resize(quantity_list_with_bc.size());
@@ -758,12 +757,12 @@ FluxingAdvectionSolver<Dimension>::remap(
   return new_variables;
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 template <typename DataType>
-class FluxingAdvectionSolver<Dimension>::InflowValueBoundaryCondition
+class FluxingAdvectionSolver<MeshType>::InflowValueBoundaryCondition
 {
  private:
-  const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+  const MeshFaceBoundary m_mesh_face_boundary;
 
   Array<const DataType> m_value_list;
 
@@ -780,7 +779,7 @@ class FluxingAdvectionSolver<Dimension>::InflowValueBoundaryCondition
     return m_value_list;
   }
 
-  InflowValueBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary, Array<const DataType> value_list)
+  InflowValueBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, Array<const DataType> value_list)
     : m_mesh_face_boundary(mesh_face_boundary), m_value_list(value_list)
   {
     ;
@@ -789,11 +788,11 @@ class FluxingAdvectionSolver<Dimension>::InflowValueBoundaryCondition
   ~InflowValueBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class FluxingAdvectionSolver<Dimension>::InflowArrayBoundaryCondition
+template <typename MeshType>
+class FluxingAdvectionSolver<MeshType>::InflowArrayBoundaryCondition
 {
  private:
-  const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+  const MeshFaceBoundary m_mesh_face_boundary;
 
   Table<const double> m_array_list;
 
@@ -810,7 +809,7 @@ class FluxingAdvectionSolver<Dimension>::InflowArrayBoundaryCondition
     return m_array_list;
   }
 
-  InflowArrayBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary, Table<const double> array_list)
+  InflowArrayBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, Table<const double> array_list)
     : m_mesh_face_boundary(mesh_face_boundary), m_array_list(array_list)
   {
     ;
@@ -819,11 +818,11 @@ class FluxingAdvectionSolver<Dimension>::InflowArrayBoundaryCondition
   ~InflowArrayBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class FluxingAdvectionSolver<Dimension>::OutflowBoundaryCondition
+template <typename MeshType>
+class FluxingAdvectionSolver<MeshType>::OutflowBoundaryCondition
 {
  private:
-  const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+  const MeshFaceBoundary m_mesh_face_boundary;
 
  public:
   const Array<const FaceId>&
@@ -832,8 +831,7 @@ class FluxingAdvectionSolver<Dimension>::OutflowBoundaryCondition
     return m_mesh_face_boundary.faceList();
   }
 
-  OutflowBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary)
-    : m_mesh_face_boundary(mesh_face_boundary)
+  OutflowBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary) : m_mesh_face_boundary(mesh_face_boundary)
   {
     ;
   }
@@ -841,11 +839,11 @@ class FluxingAdvectionSolver<Dimension>::OutflowBoundaryCondition
   ~OutflowBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class FluxingAdvectionSolver<Dimension>::SymmetryBoundaryCondition
+template <typename MeshType>
+class FluxingAdvectionSolver<MeshType>::SymmetryBoundaryCondition
 {
  private:
-  const MeshFlatFaceBoundary<Dimension> m_mesh_flat_face_boundary;
+  const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;
 
  public:
   const Array<const FaceId>&
@@ -854,7 +852,7 @@ class FluxingAdvectionSolver<Dimension>::SymmetryBoundaryCondition
     return m_mesh_flat_face_boundary.faceList();
   }
 
-  SymmetryBoundaryCondition(const MeshFlatFaceBoundary<Dimension>& mesh_flat_face_boundary)
+  SymmetryBoundaryCondition(const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
     : m_mesh_flat_face_boundary(mesh_flat_face_boundary)
   {
     ;
@@ -884,8 +882,7 @@ advectByFluxing(const std::shared_ptr<const MeshVariant> new_mesh_v,
       using NewMeshType = typename std::decay_t<decltype(new_mesh)>::element_type;
 
       if constexpr (std::is_same_v<OldMeshType, NewMeshType>) {
-        constexpr size_t Dimension = OldMeshType::Dimension;
-        return FluxingAdvectionSolver<Dimension>{old_mesh, new_mesh}.remap(remapped_variables_with_bc);
+        return FluxingAdvectionSolver<OldMeshType>{old_mesh, new_mesh}.remap(remapped_variables_with_bc);
       } else {
         throw NormalError("incompatible mesh types");
       }
diff --git a/src/scheme/HyperelasticSolver.cpp b/src/scheme/HyperelasticSolver.cpp
index a85ae5b2bebd4cece7f81d8e486d0f360863f73b..1e0654aebedcdc00108963f7ca05416752422008 100644
--- a/src/scheme/HyperelasticSolver.cpp
+++ b/src/scheme/HyperelasticSolver.cpp
@@ -55,14 +55,15 @@ hyperelastic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c)
     mesh_v->meshPointer());
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticSolverHandler::IHyperelasticSolver
 {
  private:
+  static constexpr size_t Dimension = MeshType::Dimension;
+
   using Rdxd = TinyMatrix<Dimension>;
   using Rd   = TinyVector<Dimension>;
 
-  using MeshType     = Mesh<Connectivity<Dimension>>;
   using MeshDataType = MeshData<Dimension>;
 
   using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, const double>;
@@ -263,8 +264,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
         const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
           dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
         if (dirichlet_bc_descriptor.name() == "velocity") {
-          MeshNodeBoundary<Dimension> mesh_node_boundary =
-            getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
+          MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
 
           Array<const Rd> value_list =
             InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
@@ -276,8 +276,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
           const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
 
           if constexpr (Dimension == 1) {
-            MeshNodeBoundary<Dimension> mesh_node_boundary =
-              getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+            MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
 
             Array<const double> node_values =
               InterpolateItemValue<double(Rd)>::template interpolate<ItemType::node>(pressure_id, mesh.xr(),
@@ -285,8 +284,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
 
             bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
           } else {
-            MeshFaceBoundary<Dimension> mesh_face_boundary =
-              getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+            MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
 
             MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
             Array<const double> face_values =
@@ -299,8 +297,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
           const FunctionSymbolId normal_stress_id = dirichlet_bc_descriptor.rhsSymbolId();
 
           if constexpr (Dimension == 1) {
-            MeshNodeBoundary<Dimension> mesh_node_boundary =
-              getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+            MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
 
             Array<const Rdxd> node_values =
               InterpolateItemValue<Rdxd(Rd)>::template interpolate<ItemType::node>(normal_stress_id, mesh.xr(),
@@ -308,8 +305,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
 
             bc_list.emplace_back(NormalStressBoundaryCondition{mesh_node_boundary, node_values});
           } else {
-            MeshFaceBoundary<Dimension> mesh_face_boundary =
-              getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+            MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
 
             MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
             Array<const Rdxd> face_values =
@@ -560,11 +556,11 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
   ~HyperelasticSolver()                    = default;
 };
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyPressureBC(const BoundaryConditionList& bc_list,
-                                                                           const MeshType& mesh,
-                                                                           NodeValue<Rd>& br) const
+HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyPressureBC(const BoundaryConditionList& bc_list,
+                                                                          const MeshType& mesh,
+                                                                          NodeValue<Rd>& br) const
 {
   for (const auto& boundary_condition : bc_list) {
     std::visit(
@@ -625,11 +621,11 @@ HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyPressureBC(const
   }
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyNormalStressBC(const BoundaryConditionList& bc_list,
-                                                                               const MeshType& mesh,
-                                                                               NodeValue<Rd>& br) const
+HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyNormalStressBC(const BoundaryConditionList& bc_list,
+                                                                              const MeshType& mesh,
+                                                                              NodeValue<Rd>& br) const
 {
   for (const auto& boundary_condition : bc_list) {
     std::visit(
@@ -690,11 +686,11 @@ HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyNormalStressBC(c
   }
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applySymmetryBC(const BoundaryConditionList& bc_list,
-                                                                           NodeValue<Rdxd>& Ar,
-                                                                           NodeValue<Rd>& br) const
+HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applySymmetryBC(const BoundaryConditionList& bc_list,
+                                                                          NodeValue<Rdxd>& Ar,
+                                                                          NodeValue<Rd>& br) const
 {
   for (const auto& boundary_condition : bc_list) {
     std::visit(
@@ -721,11 +717,11 @@ HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applySymmetryBC(const
   }
 }
 
-template <size_t Dimension>
+template <typename MeshType>
 void
-HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyVelocityBC(const BoundaryConditionList& bc_list,
-                                                                           NodeValue<Rdxd>& Ar,
-                                                                           NodeValue<Rd>& br) const
+HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyVelocityBC(const BoundaryConditionList& bc_list,
+                                                                          NodeValue<Rdxd>& Ar,
+                                                                          NodeValue<Rd>& br) const
 {
   for (const auto& boundary_condition : bc_list) {
     std::visit(
@@ -758,11 +754,11 @@ HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyVelocityBC(const
   }
 }
 
-template <size_t Dimension>
-class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::FixedBoundaryCondition
+template <typename MeshType>
+class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::FixedBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
 
  public:
   const Array<const NodeId>&
@@ -771,18 +767,16 @@ class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::FixedBoundaryCon
     return m_mesh_node_boundary.nodeList();
   }
 
-  FixedBoundaryCondition(const MeshNodeBoundary<Dimension> mesh_node_boundary)
-    : m_mesh_node_boundary{mesh_node_boundary}
-  {}
+  FixedBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
 
   ~FixedBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::PressureBoundaryCondition
+template <typename MeshType>
+class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::PressureBoundaryCondition
 {
  private:
-  const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+  const MeshFaceBoundary m_mesh_face_boundary;
   const Array<const double> m_value_list;
 
  public:
@@ -798,8 +792,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::PressureBoundary
     return m_value_list;
   }
 
-  PressureBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary,
-                            const Array<const double>& value_list)
+  PressureBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, const Array<const double>& value_list)
     : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
   {}
 
@@ -807,10 +800,10 @@ class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::PressureBoundary
 };
 
 template <>
-class HyperelasticSolverHandler::HyperelasticSolver<1>::PressureBoundaryCondition
+class HyperelasticSolverHandler::HyperelasticSolver<const Mesh<Connectivity<1>>>::PressureBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<1> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
   const Array<const double> m_value_list;
 
  public:
@@ -826,18 +819,18 @@ class HyperelasticSolverHandler::HyperelasticSolver<1>::PressureBoundaryConditio
     return m_value_list;
   }
 
-  PressureBoundaryCondition(const MeshNodeBoundary<1>& mesh_node_boundary, const Array<const double>& value_list)
+  PressureBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const Array<const double>& value_list)
     : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
   {}
 
   ~PressureBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::NormalStressBoundaryCondition
+template <typename MeshType>
+class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::NormalStressBoundaryCondition
 {
  private:
-  const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+  const MeshFaceBoundary m_mesh_face_boundary;
   const Array<const Rdxd> m_value_list;
 
  public:
@@ -853,8 +846,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::NormalStressBoun
     return m_value_list;
   }
 
-  NormalStressBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary,
-                                const Array<const Rdxd>& value_list)
+  NormalStressBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, const Array<const Rdxd>& value_list)
     : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
   {}
 
@@ -862,10 +854,10 @@ class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::NormalStressBoun
 };
 
 template <>
-class HyperelasticSolverHandler::HyperelasticSolver<1>::NormalStressBoundaryCondition
+class HyperelasticSolverHandler::HyperelasticSolver<const Mesh<Connectivity<1>>>::NormalStressBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<1> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
   const Array<const Rdxd> m_value_list;
 
  public:
@@ -881,18 +873,18 @@ class HyperelasticSolverHandler::HyperelasticSolver<1>::NormalStressBoundaryCond
     return m_value_list;
   }
 
-  NormalStressBoundaryCondition(const MeshNodeBoundary<1>& mesh_node_boundary, const Array<const Rdxd>& value_list)
+  NormalStressBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const Array<const Rdxd>& value_list)
     : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
   {}
 
   ~NormalStressBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::VelocityBoundaryCondition
+template <typename MeshType>
+class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::VelocityBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
 
   const Array<const TinyVector<Dimension>> m_value_list;
 
@@ -909,7 +901,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::VelocityBoundary
     return m_value_list;
   }
 
-  VelocityBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary,
+  VelocityBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
                             const Array<const TinyVector<Dimension>>& value_list)
     : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
   {}
@@ -917,14 +909,14 @@ class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::VelocityBoundary
   ~VelocityBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::SymmetryBoundaryCondition
+template <typename MeshType>
+class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::SymmetryBoundaryCondition
 {
  public:
   using Rd = TinyVector<Dimension, double>;
 
  private:
-  const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+  const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
 
  public:
   const Rd&
@@ -945,7 +937,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::SymmetryBoundary
     return m_mesh_flat_node_boundary.nodeList();
   }
 
-  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
+  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary)
     : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
   {
     ;
@@ -964,7 +956,7 @@ HyperelasticSolverHandler::HyperelasticSolverHandler(const std::shared_ptr<const
     [&](auto&& mesh) {
       using MeshType = typename std::decay_t<decltype(mesh)>::element_type;
       if constexpr (is_polygonal_mesh<MeshType>) {
-        m_hyperelastic_solver = std::make_unique<HyperelasticSolver<MeshType::Dimension>>();
+        m_hyperelastic_solver = std::make_unique<HyperelasticSolver<MeshType>>();
       } else {
         throw NormalError("unexpected mesh type");
       }
diff --git a/src/scheme/HyperelasticSolver.hpp b/src/scheme/HyperelasticSolver.hpp
index 37a20acbf372706656962718ab80230de399cd7a..75963a51971e623717bb363ec7e6c24fa4cced3e 100644
--- a/src/scheme/HyperelasticSolver.hpp
+++ b/src/scheme/HyperelasticSolver.hpp
@@ -65,14 +65,14 @@ class HyperelasticSolverHandler
           const std::shared_ptr<const DiscreteFunctionVariant>& p,
           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
 
-    IHyperelasticSolver()                      = default;
-    IHyperelasticSolver(IHyperelasticSolver&&) = default;
+    IHyperelasticSolver()                                 = default;
+    IHyperelasticSolver(IHyperelasticSolver&&)            = default;
     IHyperelasticSolver& operator=(IHyperelasticSolver&&) = default;
 
     virtual ~IHyperelasticSolver() = default;
   };
 
-  template <size_t Dimension>
+  template <typename MeshType>
   class HyperelasticSolver;
 
   std::unique_ptr<IHyperelasticSolver> m_hyperelastic_solver;