diff --git a/src/language/utils/IntegrateCellArray.hpp b/src/language/utils/IntegrateCellArray.hpp
index 41450de3b9675fbe929fcaf0fde855ccb1be2fe4..380971323f28fd45faf05597558d9f5277ab6a8b 100644
--- a/src/language/utils/IntegrateCellArray.hpp
+++ b/src/language/utils/IntegrateCellArray.hpp
@@ -4,6 +4,7 @@
 #include <language/utils/IntegrateCellValue.hpp>
 #include <mesh/CellType.hpp>
 #include <mesh/ItemArray.hpp>
+#include <mesh/MeshTraits.hpp>
 
 template <typename T>
 class IntegrateCellArray;
@@ -13,7 +14,7 @@ class IntegrateCellArray<OutputType(InputType)>
   static constexpr size_t Dimension = OutputType::Dimension;
 
  public:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE static CellArray<OutputType>
   integrate(const std::vector<FunctionSymbolId>& function_symbol_id_list,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -33,7 +34,7 @@ class IntegrateCellArray<OutputType(InputType)>
     return cell_array;
   }
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE static Table<OutputType>
   integrate(const std::vector<FunctionSymbolId>& function_symbol_id_list,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -55,7 +56,7 @@ class IntegrateCellArray<OutputType(InputType)>
     return table;
   }
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE static Table<OutputType>
   integrate(const std::vector<FunctionSymbolId>& function_symbol_id_list,
             const IQuadratureDescriptor& quadrature_descriptor,
diff --git a/src/language/utils/IntegrateCellValue.hpp b/src/language/utils/IntegrateCellValue.hpp
index cde3bfacc6730cef86305090772a03e3c0074106..d76e9af6716ee2245691b91cefd27d355db1d8e4 100644
--- a/src/language/utils/IntegrateCellValue.hpp
+++ b/src/language/utils/IntegrateCellValue.hpp
@@ -15,7 +15,7 @@ template <typename OutputType, typename InputType>
 class IntegrateCellValue<OutputType(InputType)>
 {
  public:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE static CellValue<OutputType>
   integrate(const FunctionSymbolId& function_symbol_id,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -46,7 +46,7 @@ class IntegrateCellValue<OutputType(InputType)>
       mesh_v->variant());
   }
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE static Array<OutputType>
   integrate(const FunctionSymbolId& function_symbol_id,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -76,7 +76,7 @@ class IntegrateCellValue<OutputType(InputType)>
       mesh_v->variant());
   }
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE static Array<OutputType>
   integrate(const FunctionSymbolId& function_symbol_id,
             const IQuadratureDescriptor& quadrature_descriptor,
diff --git a/src/language/utils/IntegrateOnCells.hpp b/src/language/utils/IntegrateOnCells.hpp
index 79357544738109bc5c5861e01ca22a1ccaf95af8..2f359b265e670c16c621d9ba4fb9876362605a9d 100644
--- a/src/language/utils/IntegrateOnCells.hpp
+++ b/src/language/utils/IntegrateOnCells.hpp
@@ -14,6 +14,7 @@
 #include <mesh/CellType.hpp>
 #include <mesh/Connectivity.hpp>
 #include <mesh/ItemId.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <utils/Array.hpp>
 
 class FunctionSymbolId;
@@ -79,7 +80,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
     CellList(const Array<const CellId>& cell_list) : m_cell_list{cell_list} {}
   };
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT>
   static PUGS_INLINE void
   _tensorialIntegrateTo(const FunctionSymbolId& function_symbol_id,
                         const IQuadratureDescriptor& quadrature_descriptor,
@@ -340,7 +341,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
     // LCOV_EXCL_STOP
   }
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT>
   static PUGS_INLINE void
   _directIntegrateTo(const FunctionSymbolId& function_symbol_id,
                      const IQuadratureDescriptor& quadrature_descriptor,
@@ -586,7 +587,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
   }
 
  public:
-  template <typename MeshType, typename ArrayT>
+  template <MeshConcept MeshType, typename ArrayT>
   static PUGS_INLINE void
   integrateTo(const FunctionSymbolId& function_symbol_id,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -604,7 +605,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
     }
   }
 
-  template <typename MeshType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE ArrayT<OutputType>
   integrate(const FunctionSymbolId& function_symbol_id,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -623,7 +624,7 @@ class IntegrateOnCells<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
     return value;
   }
 
-  template <typename MeshType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE ArrayT<OutputType>
   integrate(const FunctionSymbolId& function_symbol_id,
             const IQuadratureDescriptor& quadrature_descriptor,
diff --git a/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp b/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
index bac8b5c544879688f081d8da0e09e08503d4f5ff..a386e81698b24b9d225cac6dd978e7e3aedbb32e 100644
--- a/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
+++ b/src/language/utils/ItemArrayVariantFunctionInterpoler.cpp
@@ -12,7 +12,7 @@
 
 #include <memory>
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 std::shared_ptr<ItemArrayVariant>
 ItemArrayVariantFunctionInterpoler::_interpolate() const
 {
@@ -52,7 +52,7 @@ ItemArrayVariantFunctionInterpoler::_interpolate() const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 std::shared_ptr<ItemArrayVariant>
 ItemArrayVariantFunctionInterpoler::_interpolate() const
 {
diff --git a/src/language/utils/ItemArrayVariantFunctionInterpoler.hpp b/src/language/utils/ItemArrayVariantFunctionInterpoler.hpp
index 8631212dc9b44e4bcac0275a6e389d83df9e1f7e..095614ffe1e3ff96f2cc47099d05a8777760f175 100644
--- a/src/language/utils/ItemArrayVariantFunctionInterpoler.hpp
+++ b/src/language/utils/ItemArrayVariantFunctionInterpoler.hpp
@@ -5,6 +5,7 @@
 #include <mesh/IZoneDescriptor.hpp>
 #include <mesh/ItemArrayVariant.hpp>
 #include <mesh/ItemType.hpp>
+#include <mesh/MeshTraits.hpp>
 
 class MeshVariant;
 
@@ -15,10 +16,10 @@ class ItemArrayVariantFunctionInterpoler
   const ItemType m_item_type;
   const std::vector<FunctionSymbolId> m_function_id_list;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   std::shared_ptr<ItemArrayVariant> _interpolate() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   std::shared_ptr<ItemArrayVariant> _interpolate() const;
 
  public:
diff --git a/src/language/utils/ItemValueVariantFunctionInterpoler.cpp b/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
index be073b907e1412c0de2445f2fa098af4cc83eb41..6759e5763dd66bf22b40cb771ae79b6f0541fe3d 100644
--- a/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
+++ b/src/language/utils/ItemValueVariantFunctionInterpoler.cpp
@@ -12,7 +12,7 @@
 
 #include <memory>
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 std::shared_ptr<ItemValueVariant>
 ItemValueVariantFunctionInterpoler::_interpolate() const
 {
@@ -52,7 +52,7 @@ ItemValueVariantFunctionInterpoler::_interpolate() const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 std::shared_ptr<ItemValueVariant>
 ItemValueVariantFunctionInterpoler::_interpolate() const
 {
diff --git a/src/language/utils/ItemValueVariantFunctionInterpoler.hpp b/src/language/utils/ItemValueVariantFunctionInterpoler.hpp
index c1c4daddd26e5819c2ee109b7c84846f559837ef..046d9bcf51aa8042016b09d2bff137b37fe44fa4 100644
--- a/src/language/utils/ItemValueVariantFunctionInterpoler.hpp
+++ b/src/language/utils/ItemValueVariantFunctionInterpoler.hpp
@@ -5,6 +5,7 @@
 #include <mesh/IZoneDescriptor.hpp>
 #include <mesh/ItemType.hpp>
 #include <mesh/ItemValueVariant.hpp>
+#include <mesh/MeshTraits.hpp>
 
 class MeshVariant;
 
@@ -15,10 +16,10 @@ class ItemValueVariantFunctionInterpoler
   const ItemType m_item_type;
   const FunctionSymbolId m_function_id;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   std::shared_ptr<ItemValueVariant> _interpolate() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   std::shared_ptr<ItemValueVariant> _interpolate() const;
 
  public:
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 43f5528a294f6f7b6912df989a9d0c064b128dcc..d542cee7f2d34c51bfd339d13c373a7008460a49 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -11,7 +11,7 @@
 #include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp>
 #include <utils/Stringify.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const MeshType& primal_mesh)
 {
diff --git a/src/mesh/DiamondDualMeshBuilder.hpp b/src/mesh/DiamondDualMeshBuilder.hpp
index d1c2b6dbd6f0eea21cfadaa8934edc50225e19a0..ae22c2667b9d8e43710ebc243ec53de39a37d34b 100644
--- a/src/mesh/DiamondDualMeshBuilder.hpp
+++ b/src/mesh/DiamondDualMeshBuilder.hpp
@@ -2,12 +2,13 @@
 #define DIAMOND_DUAL_MESH_BUILDER_HPP
 
 #include <mesh/MeshBuilderBase.hpp>
+#include <mesh/MeshTraits.hpp>
 
 class MeshVariant;
 class DiamondDualMeshBuilder : public MeshBuilderBase
 {
  private:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _buildDualDiamondMeshFrom(const MeshType&);
 
   friend class DualMeshManager;
diff --git a/src/mesh/Dual1DMeshBuilder.cpp b/src/mesh/Dual1DMeshBuilder.cpp
index b39cd000cd62c71817bdc1236dec441115504d65..64c48423ba0ecc70ef2d6c29018c9b5e3cd8c129 100644
--- a/src/mesh/Dual1DMeshBuilder.cpp
+++ b/src/mesh/Dual1DMeshBuilder.cpp
@@ -11,7 +11,7 @@
 #include <mesh/PrimalToDual1DConnectivityDataMapper.hpp>
 #include <utils/Stringify.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 Dual1DMeshBuilder::_buildDual1DMeshFrom(const MeshType& primal_mesh)
 {
diff --git a/src/mesh/Dual1DMeshBuilder.hpp b/src/mesh/Dual1DMeshBuilder.hpp
index 901a46fa5e1afbca55bd30a294b9f5ba2c4c0fe0..42810f5c089fb00a2aa7c834536b4cca5f6d4075 100644
--- a/src/mesh/Dual1DMeshBuilder.hpp
+++ b/src/mesh/Dual1DMeshBuilder.hpp
@@ -2,12 +2,13 @@
 #define DUAL_1D_MESH_BUILDER_HPP
 
 #include <mesh/MeshBuilderBase.hpp>
+#include <mesh/MeshTraits.hpp>
 
 class MeshVariant;
 class Dual1DMeshBuilder : public MeshBuilderBase
 {
  private:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _buildDual1DMeshFrom(const MeshType&);
 
   friend class DualMeshManager;
diff --git a/src/mesh/MedianDualMeshBuilder.cpp b/src/mesh/MedianDualMeshBuilder.cpp
index 411005ec94aff3b1c2f19f56054671ab03e070b2..016f729dc6ff39f8e089c7108744f148284cad17 100644
--- a/src/mesh/MedianDualMeshBuilder.cpp
+++ b/src/mesh/MedianDualMeshBuilder.cpp
@@ -11,7 +11,7 @@
 #include <mesh/PrimalToMedianDualConnectivityDataMapper.hpp>
 #include <utils/Stringify.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 MedianDualMeshBuilder::_buildMedianDualMeshFrom(const MeshType& primal_mesh)
 {
diff --git a/src/mesh/MedianDualMeshBuilder.hpp b/src/mesh/MedianDualMeshBuilder.hpp
index ad31b6aa1ad3583d28dee31c65443d8e9e8ae01b..eb33440c5379548c5282c4380fd37dd9b646b9a8 100644
--- a/src/mesh/MedianDualMeshBuilder.hpp
+++ b/src/mesh/MedianDualMeshBuilder.hpp
@@ -2,12 +2,13 @@
 #define MEDIAN_DUAL_MESH_BUILDER_HPP
 
 #include <mesh/MeshBuilderBase.hpp>
+#include <mesh/MeshTraits.hpp>
 
 class MeshVariant;
 class MedianDualMeshBuilder : public MeshBuilderBase
 {
  private:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _buildMedianDualMeshFrom(const MeshType&);
 
   friend class DualMeshManager;
diff --git a/src/mesh/MeshEdgeBoundary.cpp b/src/mesh/MeshEdgeBoundary.cpp
index 0ad15dbad0c229461ba1de4d3d7eaf3ca797610f..65d95e9f64ed53a60587fec05bffc085a9b541d6 100644
--- a/src/mesh/MeshEdgeBoundary.cpp
+++ b/src/mesh/MeshEdgeBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeBoundary::MeshEdgeBoundary(const MeshType&, const RefEdgeList& ref_edge_list) : m_ref_edge_list(ref_edge_list)
 {}
 
@@ -13,7 +13,7 @@ template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<1>&, const RefEdgeList&);
 template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<2>&, const RefEdgeList&);
 template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<3>&, const RefEdgeList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeBoundary::MeshEdgeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
 {
   constexpr size_t Dimension = MeshType::Dimension;
@@ -58,7 +58,7 @@ MeshEdgeBoundary::MeshEdgeBoundary(const MeshType& mesh, const RefFaceList& ref_
 template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<2>&, const RefFaceList&);
 template MeshEdgeBoundary::MeshEdgeBoundary(const Mesh<3>&, const RefFaceList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeBoundary
 getMeshEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshEdgeBoundary.hpp b/src/mesh/MeshEdgeBoundary.hpp
index c65106b2d3d0d7958f0c3010295f07f7a9e9b6d9..411de18857e3f64165538a8a0f5ce283dc46b3dc 100644
--- a/src/mesh/MeshEdgeBoundary.hpp
+++ b/src/mesh/MeshEdgeBoundary.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
@@ -12,41 +13,39 @@ class [[nodiscard]] MeshEdgeBoundary   // clazy:exclude=copyable-polymorphic
   RefEdgeList m_ref_edge_list;
 
  public:
-  template <typename MeshType>
+  template <MeshConcept 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:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshEdgeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list);
-  template <typename MeshType>
+  template <MeshConcept 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 <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeBoundary getMeshEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_EDGE_BOUNDARY_HPP
diff --git a/src/mesh/MeshEdgeInterface.cpp b/src/mesh/MeshEdgeInterface.cpp
index 55235d8c7b6ff2ff8e4a5a1d54e3beb46d21c4f4..b74964eeec64c09c8dcc7fdcbb201e7af2b2fac3 100644
--- a/src/mesh/MeshEdgeInterface.cpp
+++ b/src/mesh/MeshEdgeInterface.cpp
@@ -5,7 +5,7 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeInterface::MeshEdgeInterface(const MeshType&, const RefEdgeList& ref_edge_list) : m_ref_edge_list(ref_edge_list)
 {}
 
@@ -13,7 +13,7 @@ template MeshEdgeInterface::MeshEdgeInterface(const Mesh<1>&, const RefEdgeList&
 template MeshEdgeInterface::MeshEdgeInterface(const Mesh<2>&, const RefEdgeList&);
 template MeshEdgeInterface::MeshEdgeInterface(const Mesh<3>&, const RefEdgeList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeInterface::MeshEdgeInterface(const MeshType& mesh, const RefFaceList& ref_face_list)
 {
   constexpr size_t Dimension = MeshType::Dimension;
@@ -58,7 +58,7 @@ MeshEdgeInterface::MeshEdgeInterface(const MeshType& mesh, const RefFaceList& re
 template MeshEdgeInterface::MeshEdgeInterface(const Mesh<2>&, const RefFaceList&);
 template MeshEdgeInterface::MeshEdgeInterface(const Mesh<3>&, const RefFaceList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeInterface
 getMeshEdgeInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor)
 {
diff --git a/src/mesh/MeshEdgeInterface.hpp b/src/mesh/MeshEdgeInterface.hpp
index 45bcdbe84099303a48741697767887e5a9ca524f..b4a662d04424827fc9fc3bd0cc01a01570d2857b 100644
--- a/src/mesh/MeshEdgeInterface.hpp
+++ b/src/mesh/MeshEdgeInterface.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <mesh/IInterfaceDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
@@ -12,41 +13,39 @@ class [[nodiscard]] MeshEdgeInterface   // clazy:exclude=copyable-polymorphic
   RefEdgeList m_ref_edge_list;
 
  public:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   friend MeshEdgeInterface getMeshEdgeInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor);
 
   MeshEdgeInterface& operator=(const MeshEdgeInterface&) = default;
-  MeshEdgeInterface& operator=(MeshEdgeInterface&&)      = default;
+  MeshEdgeInterface& operator=(MeshEdgeInterface&&) = 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:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshEdgeInterface(const MeshType& mesh, const RefEdgeList& ref_edge_list);
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshEdgeInterface(const MeshType& mesh, const RefFaceList& ref_face_list);
 
  public:
   MeshEdgeInterface(const MeshEdgeInterface&) = default;   // LCOV_EXCL_LINE
-  MeshEdgeInterface(MeshEdgeInterface&&)      = default;   // LCOV_EXCL_LINE
+  MeshEdgeInterface(MeshEdgeInterface &&)     = default;   // LCOV_EXCL_LINE
 
   MeshEdgeInterface()          = default;
   virtual ~MeshEdgeInterface() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshEdgeInterface getMeshEdgeInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor);
 
 #endif   // MESH_EDGE_INTERFACE_HPP
diff --git a/src/mesh/MeshFaceBoundary.cpp b/src/mesh/MeshFaceBoundary.cpp
index d883ba63af186bf6cf63bcca36a8e89db9fc06cd..9d059861286d2a623a6ddff3d7f3040954f9a9c9 100644
--- a/src/mesh/MeshFaceBoundary.cpp
+++ b/src/mesh/MeshFaceBoundary.cpp
@@ -4,7 +4,7 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFaceBoundary::MeshFaceBoundary(const MeshType&, const RefFaceList& ref_face_list) : m_ref_face_list(ref_face_list)
 {}
 
@@ -12,7 +12,7 @@ template MeshFaceBoundary::MeshFaceBoundary(const Mesh<1>&, const RefFaceList&);
 template MeshFaceBoundary::MeshFaceBoundary(const Mesh<2>&, const RefFaceList&);
 template MeshFaceBoundary::MeshFaceBoundary(const Mesh<3>&, const RefFaceList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFaceBoundary
 getMeshFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshFaceBoundary.hpp b/src/mesh/MeshFaceBoundary.hpp
index c44aafa66c6f9f47dfc839a39025a42ae614b441..443f868c0520b6031221093a07d0af1c193538ed 100644
--- a/src/mesh/MeshFaceBoundary.hpp
+++ b/src/mesh/MeshFaceBoundary.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
@@ -12,39 +13,37 @@ class [[nodiscard]] MeshFaceBoundary   // clazy:exclude=copyable-polymorphic
   RefFaceList m_ref_face_list;
 
  public:
-  template <typename MeshTypeT>
+  template <MeshConcept 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:
-  template <typename MeshType>
+  template <MeshConcept 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 <typename MeshType>
+template <MeshConcept MeshType>
 MeshFaceBoundary getMeshFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_FACE_BOUNDARY_HPP
diff --git a/src/mesh/MeshFaceInterface.cpp b/src/mesh/MeshFaceInterface.cpp
index 663ac8d44027767f5318a424e83c9ca08f26cda8..8dfe8cb49011af58b2ab1e67a9170ba3a7f8e5e1 100644
--- a/src/mesh/MeshFaceInterface.cpp
+++ b/src/mesh/MeshFaceInterface.cpp
@@ -4,7 +4,7 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFaceInterface::MeshFaceInterface(const MeshType&, const RefFaceList& ref_face_list) : m_ref_face_list(ref_face_list)
 {}
 
@@ -12,7 +12,7 @@ template MeshFaceInterface::MeshFaceInterface(const Mesh<1>&, const RefFaceList&
 template MeshFaceInterface::MeshFaceInterface(const Mesh<2>&, const RefFaceList&);
 template MeshFaceInterface::MeshFaceInterface(const Mesh<3>&, const RefFaceList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFaceInterface
 getMeshFaceInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor)
 {
diff --git a/src/mesh/MeshFaceInterface.hpp b/src/mesh/MeshFaceInterface.hpp
index 0d8277e0bc9ee7a0504e01235a2ecfcfbcc749c3..64c75e14a758ad203f5773466c5010bb5938dda8 100644
--- a/src/mesh/MeshFaceInterface.hpp
+++ b/src/mesh/MeshFaceInterface.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <mesh/IInterfaceDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
@@ -12,39 +13,37 @@ class [[nodiscard]] MeshFaceInterface   // clazy:exclude=copyable-polymorphic
   RefFaceList m_ref_face_list;
 
  public:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   friend MeshFaceInterface getMeshFaceInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor);
 
   MeshFaceInterface& operator=(const MeshFaceInterface&) = default;
-  MeshFaceInterface& operator=(MeshFaceInterface&&)      = default;
+  MeshFaceInterface& operator=(MeshFaceInterface&&) = 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:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshFaceInterface(const MeshType& mesh, const RefFaceList& ref_face_list);
 
  public:
   MeshFaceInterface(const MeshFaceInterface&) = default;   // LCOV_EXCL_LINE
-  MeshFaceInterface(MeshFaceInterface&&)      = default;   // LCOV_EXCL_LINE
+  MeshFaceInterface(MeshFaceInterface &&)     = default;   // LCOV_EXCL_LINE
 
   MeshFaceInterface()          = default;
   virtual ~MeshFaceInterface() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFaceInterface getMeshFaceInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor);
 
 #endif   // MESH_FACE_INTERFACE_HPP
diff --git a/src/mesh/MeshFlatEdgeBoundary.cpp b/src/mesh/MeshFlatEdgeBoundary.cpp
index 7a2d4918d9c822303ba6debee91b5f4d707569c4..c20b9c1a5967fa511fdd4c9f99359696d7117cc2 100644
--- a/src/mesh/MeshFlatEdgeBoundary.cpp
+++ b/src/mesh/MeshFlatEdgeBoundary.cpp
@@ -4,7 +4,7 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshFlatNodeBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFlatEdgeBoundary<MeshType>
 getMeshFlatEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshFlatEdgeBoundary.hpp b/src/mesh/MeshFlatEdgeBoundary.hpp
index f395c4ebf3594c506795b537304ba970589f533c..10ab8edb835e8c7652665cae4b9db9cfe75c4418 100644
--- a/src/mesh/MeshFlatEdgeBoundary.hpp
+++ b/src/mesh/MeshFlatEdgeBoundary.hpp
@@ -3,7 +3,7 @@
 
 #include <mesh/MeshEdgeBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
@@ -24,7 +24,7 @@ class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=co
   MeshFlatEdgeBoundary& operator=(const MeshFlatEdgeBoundary&) = default;
   MeshFlatEdgeBoundary& operator=(MeshFlatEdgeBoundary&&)      = default;
 
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshFlatEdgeBoundary<MeshTypeT> getMeshFlatEdgeBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -40,7 +40,7 @@ class MeshFlatEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=co
   ~MeshFlatEdgeBoundary()                           = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFlatEdgeBoundary<MeshType> getMeshFlatEdgeBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshFlatFaceBoundary.cpp b/src/mesh/MeshFlatFaceBoundary.cpp
index 2ecdf6768eb865413386b10a1a419c6e30799def..632c6a033dee7fc5593363660cc179faa966fb80 100644
--- a/src/mesh/MeshFlatFaceBoundary.cpp
+++ b/src/mesh/MeshFlatFaceBoundary.cpp
@@ -4,7 +4,7 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshFlatNodeBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFlatFaceBoundary<MeshType>
 getMeshFlatFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshFlatFaceBoundary.hpp b/src/mesh/MeshFlatFaceBoundary.hpp
index 1aa2e854fbfd3816f006f1f2aaaf44abffd6376b..0c53a31e095251c6ffaa70b7f6d58bfce10a6864 100644
--- a/src/mesh/MeshFlatFaceBoundary.hpp
+++ b/src/mesh/MeshFlatFaceBoundary.hpp
@@ -3,7 +3,7 @@
 
 #include <mesh/MeshFaceBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
@@ -24,7 +24,7 @@ class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=co
   MeshFlatFaceBoundary& operator=(const MeshFlatFaceBoundary&) = default;
   MeshFlatFaceBoundary& operator=(MeshFlatFaceBoundary&&)      = default;
 
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshFlatFaceBoundary<MeshTypeT> getMeshFlatFaceBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -40,7 +40,7 @@ class MeshFlatFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=co
   ~MeshFlatFaceBoundary()                           = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFlatFaceBoundary<MeshType> getMeshFlatFaceBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
index 983f967385c0832ffaf5516a7b46410d22ffd74b..bdc7c20a2c854c2a9a6492d8f06f46a026b38b0b 100644
--- a/src/mesh/MeshFlatNodeBoundary.cpp
+++ b/src/mesh/MeshFlatNodeBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/MeshNodeBoundaryUtils.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 MeshFlatNodeBoundary<MeshType>::_checkBoundaryIsFlat(const TinyVector<MeshType::Dimension, double>& normal,
                                                      const TinyVector<MeshType::Dimension, double>& origin,
@@ -323,7 +323,7 @@ MeshFlatNodeBoundary<Mesh<3>>::_getOutgoingNormal(const Mesh<3>& mesh)
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFlatNodeBoundary<MeshType>
 getMeshFlatNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshFlatNodeBoundary.hpp b/src/mesh/MeshFlatNodeBoundary.hpp
index 453804724b306a0c114c0f8cd43ae99cf3430c93..2eb28d14899a1d70f47f8cd139b435c4a81390b7 100644
--- a/src/mesh/MeshFlatNodeBoundary.hpp
+++ b/src/mesh/MeshFlatNodeBoundary.hpp
@@ -2,8 +2,9 @@
 #define MESH_FLAT_NODE_BOUNDARY_HPP
 
 #include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshTraits.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
@@ -19,23 +20,21 @@ class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // cl
   Rd _getNormal(const MeshType& mesh);
 
   void _checkBoundaryIsFlat(const TinyVector<MeshType::Dimension, double>& normal,
-                            const TinyVector<MeshType::Dimension, double>& origin,
-                            const double length,
+                            const TinyVector<MeshType::Dimension, double>& origin, const double length,
                             const MeshType& mesh) const;
 
   Rd _getOutgoingNormal(const MeshType& mesh);
 
  public:
-  const Rd&
-  outgoingNormal() const
+  const Rd& outgoingNormal() const
   {
     return m_outgoing_normal;
   }
 
   MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
-  MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&)      = default;
+  MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
 
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshFlatNodeBoundary<MeshTypeT> getMeshFlatNodeBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -51,11 +50,11 @@ class [[nodiscard]] MeshFlatNodeBoundary final : public MeshNodeBoundary   // cl
  public:
   MeshFlatNodeBoundary()                            = default;
   MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
-  MeshFlatNodeBoundary(MeshFlatNodeBoundary&&)      = default;
+  MeshFlatNodeBoundary(MeshFlatNodeBoundary &&)     = default;
   ~MeshFlatNodeBoundary()                           = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshFlatNodeBoundary<MeshType> getMeshFlatNodeBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshLineEdgeBoundary.cpp b/src/mesh/MeshLineEdgeBoundary.cpp
index f1e0333a9afd8a642f3809560f0c16431e98e1e4..40de23d1e2dc86bbd82b8542b75e32d1cac40352 100644
--- a/src/mesh/MeshLineEdgeBoundary.cpp
+++ b/src/mesh/MeshLineEdgeBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/MeshLineNodeBoundary.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshLineEdgeBoundary<MeshType>
 getMeshLineEdgeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshLineEdgeBoundary.hpp b/src/mesh/MeshLineEdgeBoundary.hpp
index 39c97a194f756019d9f80fe7ec67b7baf28261fb..0fc5db0ba6b73960cd1a64c5a0626585c857a4b8 100644
--- a/src/mesh/MeshLineEdgeBoundary.hpp
+++ b/src/mesh/MeshLineEdgeBoundary.hpp
@@ -4,7 +4,7 @@
 #include <algebra/TinyMatrix.hpp>
 #include <mesh/MeshEdgeBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class [[nodiscard]] MeshLineEdgeBoundary final : public MeshEdgeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
@@ -18,7 +18,7 @@ class [[nodiscard]] MeshLineEdgeBoundary final : public MeshEdgeBoundary   // cl
   const Rd m_direction;
 
  public:
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshLineEdgeBoundary<MeshTypeT> getMeshLineEdgeBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -44,7 +44,7 @@ class [[nodiscard]] MeshLineEdgeBoundary final : public MeshEdgeBoundary   // cl
   ~MeshLineEdgeBoundary()                           = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshLineEdgeBoundary<MeshType> getMeshLineEdgeBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshLineFaceBoundary.cpp b/src/mesh/MeshLineFaceBoundary.cpp
index 81ea991ed6dcf999e3144a9b4d9d48f080a1499e..3608e7aef37dfeddf9eb2c6ad7d207d1acae3f07 100644
--- a/src/mesh/MeshLineFaceBoundary.cpp
+++ b/src/mesh/MeshLineFaceBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/MeshLineNodeBoundary.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshLineFaceBoundary<MeshType>
 getMeshLineFaceBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshLineFaceBoundary.hpp b/src/mesh/MeshLineFaceBoundary.hpp
index 203faac8ef170dc1a15e6dc92d00793bc878fee5..206961a08f88f9a4008a14a7a223b9d90ff323ed 100644
--- a/src/mesh/MeshLineFaceBoundary.hpp
+++ b/src/mesh/MeshLineFaceBoundary.hpp
@@ -4,7 +4,7 @@
 #include <algebra/TinyMatrix.hpp>
 #include <mesh/MeshFaceBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class [[nodiscard]] MeshLineFaceBoundary final : public MeshFaceBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
@@ -17,7 +17,7 @@ class [[nodiscard]] MeshLineFaceBoundary final : public MeshFaceBoundary   // cl
   const Rd m_direction;
 
  public:
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshLineFaceBoundary<MeshTypeT> getMeshLineFaceBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -43,7 +43,7 @@ class [[nodiscard]] MeshLineFaceBoundary final : public MeshFaceBoundary   // cl
   ~MeshLineFaceBoundary()                           = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshLineFaceBoundary<MeshType> getMeshLineFaceBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshLineNodeBoundary.cpp b/src/mesh/MeshLineNodeBoundary.cpp
index 1228fa5aaaa21f59f7d76c89f41682787426bc96..1cd3bb64d749d49de82d925bff9eaa3ecf88ead9 100644
--- a/src/mesh/MeshLineNodeBoundary.cpp
+++ b/src/mesh/MeshLineNodeBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/MeshNodeBoundaryUtils.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 MeshLineNodeBoundary<MeshType>::_checkBoundaryIsLine(const TinyVector<MeshType::Dimension, double>& direction,
                                                      const TinyVector<MeshType::Dimension, double>& origin,
@@ -105,7 +105,7 @@ MeshLineNodeBoundary<Mesh<3>>::_getDirection(const Mesh<3>& mesh)
   return direction;
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshLineNodeBoundary<MeshType>
 getMeshLineNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshLineNodeBoundary.hpp b/src/mesh/MeshLineNodeBoundary.hpp
index bb6fc4b0424fbdf1f96a21821bc003825a6a773b..fb75ef9a49d77095c47d1a45db42e636b905c572 100644
--- a/src/mesh/MeshLineNodeBoundary.hpp
+++ b/src/mesh/MeshLineNodeBoundary.hpp
@@ -4,7 +4,7 @@
 #include <algebra/TinyMatrix.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class [[nodiscard]] MeshLineNodeBoundary final : public MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 {
  public:
@@ -24,7 +24,7 @@ class [[nodiscard]] MeshLineNodeBoundary final : public MeshNodeBoundary   // cl
                             const MeshType& mesh) const;
 
  public:
-  template <typename MeshTypeT>
+  template <MeshConcept MeshTypeT>
   friend MeshLineNodeBoundary<MeshTypeT> getMeshLineNodeBoundary(const MeshTypeT& mesh,
                                                                  const IBoundaryDescriptor& boundary_descriptor);
 
@@ -58,7 +58,7 @@ class [[nodiscard]] MeshLineNodeBoundary final : public MeshNodeBoundary   // cl
   ~MeshLineNodeBoundary()                           = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshLineNodeBoundary<MeshType> getMeshLineNodeBoundary(const MeshType& mesh,
                                                        const IBoundaryDescriptor& boundary_descriptor);
 
diff --git a/src/mesh/MeshNodeBoundary.cpp b/src/mesh/MeshNodeBoundary.cpp
index e750580ce1298e4ea110847d89ddf385ca893c38..8cdcad2fb91700aaeed9bd8cd9301ac4173c7294 100644
--- a/src/mesh/MeshNodeBoundary.cpp
+++ b/src/mesh/MeshNodeBoundary.cpp
@@ -5,7 +5,7 @@
 #include <mesh/Mesh.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeBoundary::MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
 {
   constexpr size_t Dimension = MeshType::Dimension;
@@ -52,7 +52,7 @@ MeshNodeBoundary::MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_
   const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_node_list);
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeBoundary::MeshNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list)
 {
   constexpr size_t Dimension = MeshType::Dimension;
@@ -98,7 +98,7 @@ MeshNodeBoundary::MeshNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_
   const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_node_list);
 }
 
-template <typename MeshType>
+template <MeshConcept 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) {
@@ -121,7 +121,7 @@ template MeshNodeBoundary::MeshNodeBoundary(const Mesh<1>&, const RefNodeList&);
 template MeshNodeBoundary::MeshNodeBoundary(const Mesh<2>&, const RefNodeList&);
 template MeshNodeBoundary::MeshNodeBoundary(const Mesh<3>&, const RefNodeList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeBoundary
 getMeshNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
 {
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 0e869e2eef5b9783aba66840e050071cc3ddf57d..36d79a308e8edd730d0fab33d2ced9a6af0e4e95 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
@@ -12,43 +13,41 @@ class [[nodiscard]] MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
   RefNodeList m_ref_node_list;
 
  public:
-  template <typename MeshTypeT>
+  template <MeshConcept 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:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list);
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list);
-  template <typename MeshType>
+  template <MeshConcept 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 <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeBoundary getMeshNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshNodeBoundaryUtils.hpp b/src/mesh/MeshNodeBoundaryUtils.hpp
index a62f9d8f3affa0cf2a840e69776b91c88fa8ec0e..f187a210d2c3db4ec40fbafc3ad276454df3cbe5 100644
--- a/src/mesh/MeshNodeBoundaryUtils.hpp
+++ b/src/mesh/MeshNodeBoundaryUtils.hpp
@@ -2,11 +2,12 @@
 #define MESH_NODE_BOUNDARY_UTILS_HPP
 
 #include <algebra/TinyVector.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 
 #include <array>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 std::array<TinyVector<MeshType::Dimension>, MeshType::Dimension*(MeshType::Dimension - 1)> getBounds(
   const MeshType& mesh,
   const RefNodeList& ref_node_list);
diff --git a/src/mesh/MeshNodeInterface.cpp b/src/mesh/MeshNodeInterface.cpp
index 54da43007219d7b6de8c68f90c8d2b80aecb0697..57eed14cb4f90037cf9a382a62f5a6917886bd88 100644
--- a/src/mesh/MeshNodeInterface.cpp
+++ b/src/mesh/MeshNodeInterface.cpp
@@ -5,7 +5,7 @@
 #include <mesh/MeshTraits.hpp>
 #include <utils/Messenger.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeInterface::MeshNodeInterface(const MeshType& mesh, const RefFaceList& ref_face_list)
 {
   static_assert(is_polygonal_mesh_v<MeshType>);
@@ -53,7 +53,7 @@ MeshNodeInterface::MeshNodeInterface(const MeshType& mesh, const RefFaceList& re
   const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_node_list);
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeInterface::MeshNodeInterface(const MeshType& mesh, const RefEdgeList& ref_edge_list)
 {
   static_assert(is_polygonal_mesh_v<MeshType>);
@@ -100,7 +100,7 @@ MeshNodeInterface::MeshNodeInterface(const MeshType& mesh, const RefEdgeList& re
   const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_node_list);
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeInterface::MeshNodeInterface(const MeshType&, const RefNodeList& ref_node_list) : m_ref_node_list(ref_node_list)
 {
   if (ref_node_list.type() != RefItemListBase::Type::interface) {
@@ -123,7 +123,7 @@ template MeshNodeInterface::MeshNodeInterface(const Mesh<1>&, const RefNodeList&
 template MeshNodeInterface::MeshNodeInterface(const Mesh<2>&, const RefNodeList&);
 template MeshNodeInterface::MeshNodeInterface(const Mesh<3>&, const RefNodeList&);
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeInterface
 getMeshNodeInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor)
 {
diff --git a/src/mesh/MeshNodeInterface.hpp b/src/mesh/MeshNodeInterface.hpp
index 4eaafd06b56dfcb07d9e1bc5c0b29d661c0386e1..bb2105649bb5acee027821a9c227b7323789cc99 100644
--- a/src/mesh/MeshNodeInterface.hpp
+++ b/src/mesh/MeshNodeInterface.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <mesh/IInterfaceDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
 
@@ -12,43 +13,41 @@ class [[nodiscard]] MeshNodeInterface   // clazy:exclude=copyable-polymorphic
   RefNodeList m_ref_node_list;
 
  public:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   friend MeshNodeInterface getMeshNodeInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor);
 
   MeshNodeInterface& operator=(const MeshNodeInterface&) = default;
-  MeshNodeInterface& operator=(MeshNodeInterface&&)      = default;
+  MeshNodeInterface& operator=(MeshNodeInterface&&) = 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:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshNodeInterface(const MeshType& mesh, const RefFaceList& ref_face_list);
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshNodeInterface(const MeshType& mesh, const RefEdgeList& ref_edge_list);
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshNodeInterface(const MeshType&, const RefNodeList& ref_node_list);
 
  public:
   MeshNodeInterface(const MeshNodeInterface&) = default;
-  MeshNodeInterface(MeshNodeInterface&&)      = default;
+  MeshNodeInterface(MeshNodeInterface &&)     = default;
 
   MeshNodeInterface()          = default;
   virtual ~MeshNodeInterface() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 MeshNodeInterface getMeshNodeInterface(const MeshType& mesh, const IInterfaceDescriptor& interface_descriptor);
 
 #endif   // MESH_NODE_INTERFACE_HPP
diff --git a/src/mesh/MeshRandomizer.cpp b/src/mesh/MeshRandomizer.cpp
index d6b8c2f1b43616955db4f48ec67210738f582e4c..8c94807db37ad95e4f905e992f9a1c000d14c59e 100644
--- a/src/mesh/MeshRandomizer.cpp
+++ b/src/mesh/MeshRandomizer.cpp
@@ -16,7 +16,7 @@
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <utils/RandomEngine.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshRandomizerHandler::MeshRandomizer
 {
  private:
@@ -274,7 +274,7 @@ class MeshRandomizerHandler::MeshRandomizer
   ~MeshRandomizer() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshRandomizerHandler::MeshRandomizer<MeshType>::AxisBoundaryCondition
 {
  public:
@@ -313,7 +313,7 @@ class MeshRandomizerHandler::MeshRandomizer<Mesh<1>>::AxisBoundaryCondition
   ~AxisBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshRandomizerHandler::MeshRandomizer<MeshType>::FixedBoundaryCondition
 {
  private:
@@ -331,7 +331,7 @@ class MeshRandomizerHandler::MeshRandomizer<MeshType>::FixedBoundaryCondition
   ~FixedBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshRandomizerHandler::MeshRandomizer<MeshType>::SymmetryBoundaryCondition
 {
  public:
diff --git a/src/mesh/MeshRandomizer.hpp b/src/mesh/MeshRandomizer.hpp
index 048e04bfdad1780fc41362595ca0378400a6266e..6db719cb8fac9082624be806e75a6e4a59091f25 100644
--- a/src/mesh/MeshRandomizer.hpp
+++ b/src/mesh/MeshRandomizer.hpp
@@ -1,6 +1,8 @@
 #ifndef MESH_RANDOMIZER_HPP
 #define MESH_RANDOMIZER_HPP
 
+#include <mesh/MeshTraits.hpp>
+
 #include <memory>
 #include <vector>
 
@@ -11,7 +13,7 @@ class FunctionSymbolId;
 class MeshRandomizerHandler
 {
  private:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   class MeshRandomizer;
 
  public:
diff --git a/src/mesh/MeshRelaxer.cpp b/src/mesh/MeshRelaxer.cpp
index 7dcc86bc0df00f33b199935ca5268b2d060e81f6..65302649d995851ced9a485788f8e1fc2b9ef14b 100644
--- a/src/mesh/MeshRelaxer.cpp
+++ b/src/mesh/MeshRelaxer.cpp
@@ -5,12 +5,12 @@
 #include <mesh/MeshTraits.hpp>
 #include <mesh/MeshVariant.hpp>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 std::shared_ptr<const MeshVariant>
 MeshRelaxer::_relax(const MeshType& source_mesh, const MeshType& destination_mesh, const double& theta) const
 {
   static_assert(is_polygonal_mesh_v<MeshType>);
-  if (source_mesh.shared_connectivity() == destination_mesh.shared_connectivity()) {
+  if (source_mesh.connectivity().id() == destination_mesh.connectivity().id()) {
     using ConnectivityType               = typename MeshType::Connectivity;
     constexpr size_t Dimension           = MeshType::Dimension;
     const ConnectivityType& connectivity = source_mesh.connectivity();
diff --git a/src/mesh/MeshRelaxer.hpp b/src/mesh/MeshRelaxer.hpp
index 1e2b082064dadb1f85d1c1e2c93151b331993c1c..c771842cf947fb4cf9b00c05a0a6cb7fd301fe83 100644
--- a/src/mesh/MeshRelaxer.hpp
+++ b/src/mesh/MeshRelaxer.hpp
@@ -1,6 +1,8 @@
 #ifndef MESH_RELAXER_HPP
 #define MESH_RELAXER_HPP
 
+#include <mesh/MeshTraits.hpp>
+
 class MeshVariant;
 
 #include <memory>
@@ -8,7 +10,7 @@ class MeshVariant;
 class MeshRelaxer
 {
  private:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   std::shared_ptr<const MeshVariant> _relax(const MeshType& source_mesh,
                                             const MeshType& destination_mesh,
                                             const double& theta) const;
diff --git a/src/mesh/MeshSmoother.cpp b/src/mesh/MeshSmoother.cpp
index eafc9b94bbd9a5a61edaf2d7440fc025b38ecee3..7950c732abe7a07961e95fbd7da1728a8f4ba937 100644
--- a/src/mesh/MeshSmoother.cpp
+++ b/src/mesh/MeshSmoother.cpp
@@ -20,7 +20,7 @@
 
 #include <variant>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshSmootherHandler::MeshSmoother
 {
  private:
@@ -321,7 +321,7 @@ class MeshSmootherHandler::MeshSmoother
   ~MeshSmoother() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshSmootherHandler::MeshSmoother<MeshType>::AxisBoundaryCondition
 {
  public:
@@ -360,7 +360,7 @@ class MeshSmootherHandler::MeshSmoother<Mesh<1>>::AxisBoundaryCondition
   ~AxisBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshSmootherHandler::MeshSmoother<MeshType>::FixedBoundaryCondition
 {
  private:
@@ -378,7 +378,7 @@ class MeshSmootherHandler::MeshSmoother<MeshType>::FixedBoundaryCondition
   ~FixedBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class MeshSmootherHandler::MeshSmoother<MeshType>::SymmetryBoundaryCondition
 {
  public:
diff --git a/src/mesh/MeshSmoother.hpp b/src/mesh/MeshSmoother.hpp
index c5875561a5bce80257dbf6a118517c2b8443c021..d43e380a97987e495c785a509f85f61743e2f50f 100644
--- a/src/mesh/MeshSmoother.hpp
+++ b/src/mesh/MeshSmoother.hpp
@@ -1,6 +1,8 @@
 #ifndef MESH_SMOOTHER_HPP
 #define MESH_SMOOTHER_HPP
 
+#include <mesh/MeshTraits.hpp>
+
 #include <memory>
 #include <vector>
 
@@ -13,7 +15,7 @@ class DiscreteFunctionVariant;
 class MeshSmootherHandler
 {
  private:
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   class MeshSmoother;
 
  public:
diff --git a/src/mesh/MeshTransformer.cpp b/src/mesh/MeshTransformer.cpp
index d44d3251c92dcc521e993906cea021467b9391d7..81f723a6cee92a43ab4b06f52b4dffc7a3b92095 100644
--- a/src/mesh/MeshTransformer.cpp
+++ b/src/mesh/MeshTransformer.cpp
@@ -1,6 +1,7 @@
 #include <mesh/MeshTransformer.hpp>
 
 #include <mesh/Connectivity.hpp>
+#include <mesh/ItemArray.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshTraits.hpp>
 #include <mesh/MeshVariant.hpp>
@@ -13,15 +14,15 @@ class MeshTransformer::MeshTransformation<OutputType(InputType)>
   static constexpr size_t Dimension = OutputType::Dimension;
 
  public:
-  template <typename MeshType>
-  static std::shared_ptr<const MeshType>
-  transform(const FunctionSymbolId& function_symbol_id, const MeshType& mesh)
+  template <size_t Dimension>
+  static std::shared_ptr<const Mesh<Dimension>>
+  transform(const FunctionSymbolId& function_symbol_id, const Mesh<Dimension>& mesh)
   {
     NodeValue<OutputType> xr(mesh.connectivity());
     NodeValue<const InputType> given_xr = mesh.xr();
     EvaluateAtPoints<OutputType(InputType)>::evaluateTo(function_symbol_id, given_xr, xr);
 
-    return std::make_shared<const MeshType>(mesh.shared_connectivity(), xr);
+    return std::make_shared<const Mesh<Dimension>>(mesh.shared_connectivity(), xr);
   }
 };
 
diff --git a/src/mesh/MeshVariant.hpp b/src/mesh/MeshVariant.hpp
index 7278878a9d7be440a179e913cb6491c5175d7094..c37d2dd01343f0b4e58f785b05200a45a0c2cb44 100644
--- a/src/mesh/MeshVariant.hpp
+++ b/src/mesh/MeshVariant.hpp
@@ -40,7 +40,7 @@ class MeshVariant
   const IConnectivity& connectivity() const;
   std::shared_ptr<const IConnectivity> shared_connectivity() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   PUGS_INLINE std::shared_ptr<const MeshType>
   get() const
   {
@@ -66,7 +66,7 @@ class MeshVariant
 
   MeshVariant() = delete;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   MeshVariant(const std::shared_ptr<const MeshType>& p_mesh) : m_p_mesh_variant{p_mesh}
   {}
 
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index 5415d59f829c1662bc41d37dd455ffa862d4ee94..991ec74537fffaaa2bef855f5c4af88f7acfa03f 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -148,7 +148,7 @@ GnuplotWriter::_writeData(const ItemDataT& item_data,
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 GnuplotWriter::_writeDataAtNodes(const MeshType& mesh,
                                  const OutputNamedItemDataSet& output_named_item_data_set,
@@ -240,7 +240,7 @@ GnuplotWriter::_writeNodeData(const NodeArray<DataType>& node_array, NodeId node
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 GnuplotWriter::_write(const MeshType& mesh,
                       const OutputNamedItemDataSet& output_named_item_data_set,
diff --git a/src/output/GnuplotWriter.hpp b/src/output/GnuplotWriter.hpp
index 9cf22a7efe20878379a57f9025b2f537dd2830bc..3069b79da716d2640e8446822e23fcab37b98316 100644
--- a/src/output/GnuplotWriter.hpp
+++ b/src/output/GnuplotWriter.hpp
@@ -5,6 +5,7 @@
 
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <output/OutputNamedItemValueSet.hpp>
 
 class MeshVariant;
@@ -40,12 +41,12 @@ class GnuplotWriter final : public WriterBase
                   [[maybe_unused]] NodeId node_id,
                   std::ostream& fout) const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _writeDataAtNodes(const MeshType& mesh,
                          const OutputNamedItemDataSet& output_named_item_data_set,
                          std::ostream& fout) const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _write(const MeshType& mesh,
               const OutputNamedItemDataSet& output_named_item_data_set,
               std::optional<double> time) const;
diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
index a31d5432f825ac2611d1029814a685c0e740e632..84e8aa080c6802eb546a4d0183fb446201e31ac8 100644
--- a/src/output/GnuplotWriter1D.cpp
+++ b/src/output/GnuplotWriter1D.cpp
@@ -142,7 +142,7 @@ GnuplotWriter1D::_itemDataNbRow(const ItemArray<DataType, item_type>& item_array
   return item_array.sizeOfArrays();
 }
 
-template <typename MeshType, ItemType item_type>
+template <MeshConcept MeshType, ItemType item_type>
 void
 GnuplotWriter1D::_writeItemDatas(const MeshType& mesh,
                                  const OutputNamedItemDataSet& output_named_item_data_set,
@@ -275,7 +275,7 @@ GnuplotWriter1D::_writeItemDatas(const MeshType& mesh,
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 GnuplotWriter1D::_write(const MeshType& mesh,
                         const OutputNamedItemDataSet& output_named_item_data_set,
diff --git a/src/output/GnuplotWriter1D.hpp b/src/output/GnuplotWriter1D.hpp
index 61b697dd1e83f61252aa81459d0673d3e4036920..66661918252856d8b6a069773c931c705a9b9330 100644
--- a/src/output/GnuplotWriter1D.hpp
+++ b/src/output/GnuplotWriter1D.hpp
@@ -5,6 +5,7 @@
 
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <output/OutputNamedItemValueSet.hpp>
 
 class MeshVariant;
@@ -37,14 +38,14 @@ class GnuplotWriter1D final : public WriterBase
   template <typename DataType, ItemType item_type>
   size_t _itemDataNbRow(const ItemArray<DataType, item_type>&) const;
 
-  template <typename MeshType, ItemType item_type>
+  template <MeshConcept MeshType, ItemType item_type>
   void _writeItemDatas(const MeshType& mesh,
                        const OutputNamedItemDataSet& output_named_item_data_set,
                        std::ostream& fout) const;
 
   void _writePreamble(const OutputNamedItemDataSet& output_named_item_value_set, std::ostream& fout) const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _write(const MeshType& mesh,
               const OutputNamedItemDataSet& output_named_item_value_set,
               std::optional<double> time) const;
diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index 889b6b554932d4a93322eab5b0c1ac40548ee982..cd7834d713d2f9ca94d6dd65aa301ed444c6340f 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -355,7 +355,7 @@ VTKWriter::_write_node_data(std::ofstream& os,
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 VTKWriter::_write(const MeshType& mesh,
                   const OutputNamedItemDataSet& given_output_named_item_data_set,
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 1141daec8535f7d4c18af65a28b4816e8ea7b775..6c40f1005952d1649dcbd989f639907826a312fd 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -5,6 +5,7 @@
 
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <output/OutputNamedItemValueSet.hpp>
 
 #include <optional>
@@ -93,7 +94,7 @@ class VTKWriter final : public WriterBase
                         const ItemDataT& item_data,
                         SerializedDataList& serialized_data_list) const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   void _write(const MeshType& mesh,
               const OutputNamedItemDataSet& output_named_item_data_set,
               std::optional<double> time) const;
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 43476df571d4ac14e7fcfb8e6fd21593d6049896..90998d3aba2e9ecf0a9b5c4e8961c7b87889eaa5 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -48,7 +48,7 @@ acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
     c.meshVariant()->variant());
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler::IAcousticSolver
 {
  private:
@@ -514,7 +514,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   ~AcousticSolver()                = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 AcousticSolverHandler::AcousticSolver<MeshType>::_applyPressureBC(const BoundaryConditionList& bc_list,
                                                                   const MeshType& mesh,
@@ -579,7 +579,7 @@ AcousticSolverHandler::AcousticSolver<MeshType>::_applyPressureBC(const Boundary
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 AcousticSolverHandler::AcousticSolver<MeshType>::_applySymmetryBC(const BoundaryConditionList& bc_list,
                                                                   NodeValue<Rdxd>& Ar,
@@ -610,7 +610,7 @@ AcousticSolverHandler::AcousticSolver<MeshType>::_applySymmetryBC(const Boundary
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 AcousticSolverHandler::AcousticSolver<MeshType>::_applyVelocityBC(const BoundaryConditionList& bc_list,
                                                                   NodeValue<Rdxd>& Ar,
@@ -647,7 +647,7 @@ AcousticSolverHandler::AcousticSolver<MeshType>::_applyVelocityBC(const Boundary
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 AcousticSolverHandler::AcousticSolver<MeshType>::_applyExternalVelocityBC(const BoundaryConditionList& bc_list,
                                                                           const DiscreteScalarFunction& p,
@@ -676,7 +676,7 @@ AcousticSolverHandler::AcousticSolver<MeshType>::_applyExternalVelocityBC(const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class AcousticSolverHandler::AcousticSolver<MeshType>::FixedBoundaryCondition
 {
  private:
@@ -694,7 +694,7 @@ class AcousticSolverHandler::AcousticSolver<MeshType>::FixedBoundaryCondition
   ~FixedBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class AcousticSolverHandler::AcousticSolver<MeshType>::PressureBoundaryCondition
 {
  private:
@@ -748,7 +748,7 @@ class AcousticSolverHandler::AcousticSolver<Mesh<1>>::PressureBoundaryCondition
   ~PressureBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class AcousticSolverHandler::AcousticSolver<MeshType>::ExternalFSIVelocityBoundaryCondition
 {
  private:
@@ -797,7 +797,7 @@ class AcousticSolverHandler::AcousticSolver<MeshType>::ExternalFSIVelocityBounda
   ~ExternalFSIVelocityBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class AcousticSolverHandler::AcousticSolver<MeshType>::VelocityBoundaryCondition
 {
  private:
@@ -826,7 +826,7 @@ class AcousticSolverHandler::AcousticSolver<MeshType>::VelocityBoundaryCondition
   ~VelocityBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class AcousticSolverHandler::AcousticSolver<MeshType>::SymmetryBoundaryCondition
 {
  public:
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 75e08981dfb50993233eb80ddfdf67266a24cf06..73f9522ba12ebb5f2570486a1bcb3a108b6b4270 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -1,6 +1,8 @@
 #ifndef ACOUSTIC_SOLVER_HPP
 #define ACOUSTIC_SOLVER_HPP
 
+#include <mesh/MeshTraits.hpp>
+
 #include <memory>
 #include <tuple>
 #include <vector>
@@ -59,14 +61,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 <typename MeshType>
+  template <MeshConcept MeshType>
   class AcousticSolver;
 
   std::unique_ptr<IAcousticSolver> m_acoustic_solver;
diff --git a/src/scheme/CellIntegrator.hpp b/src/scheme/CellIntegrator.hpp
index 26544b7c4f742dee16295922d9a9206320847844..33241f9b8683700d24507000a50ed0f0ae507517 100644
--- a/src/scheme/CellIntegrator.hpp
+++ b/src/scheme/CellIntegrator.hpp
@@ -75,7 +75,7 @@ class CellIntegrator
     CellList(const ArrayT<CellId>& cell_list) : m_cell_list{cell_list} {}
   };
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   _tensorialIntegrateTo(std::function<OutputType(InputType)> function,
                         const IQuadratureDescriptor& quadrature_descriptor,
@@ -303,7 +303,7 @@ class CellIntegrator
     // LCOV_EXCL_STOP
   }
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   _directIntegrateTo(std::function<OutputType(InputType)> function,
                      const IQuadratureDescriptor& quadrature_descriptor,
@@ -514,7 +514,7 @@ class CellIntegrator
   }
 
  public:
-  template <typename MeshType, typename OutputArrayT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   integrateTo(const std::function<OutputType(InputType)>& function,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -532,7 +532,7 @@ class CellIntegrator
     }
   }
 
-  template <typename MeshType, typename OutputArrayT, typename FunctionType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename FunctionType>
   static PUGS_INLINE void
   integrateTo(const FunctionType& function,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -542,7 +542,7 @@ class CellIntegrator
     integrateTo(std::function(function), quadrature_descriptor, mesh, value);
   }
 
-  template <typename MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE ArrayT<OutputType>
   integrate(const std::function<OutputType(InputType)>& function,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -561,7 +561,7 @@ class CellIntegrator
     return value;
   }
 
-  template <typename MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE auto
   integrate(const FunctionType& function,
             const IQuadratureDescriptor& quadrature_descriptor,
diff --git a/src/scheme/DiscreteFunctionIntegrator.cpp b/src/scheme/DiscreteFunctionIntegrator.cpp
index 92dd526d20acb01fb8605e3595dbac0491ac7570..5caf8e5012a835417a99017c3625e9934092e748 100644
--- a/src/scheme/DiscreteFunctionIntegrator.cpp
+++ b/src/scheme/DiscreteFunctionIntegrator.cpp
@@ -8,7 +8,7 @@
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <utils/Exceptions.hpp>
 
-template <typename MeshType, typename DataType, typename ValueType>
+template <MeshConcept MeshType, typename DataType, typename ValueType>
 DiscreteFunctionVariant
 DiscreteFunctionIntegrator::_integrateOnZoneList() const
 {
@@ -67,7 +67,7 @@ DiscreteFunctionIntegrator::_integrateOnZoneList() const
   return DiscreteFunctionP0<ValueType>(m_mesh, cell_value);
 }
 
-template <typename MeshType, typename DataType, typename ValueType>
+template <MeshConcept MeshType, typename DataType, typename ValueType>
 DiscreteFunctionVariant
 DiscreteFunctionIntegrator::_integrateGlobally() const
 {
@@ -82,7 +82,7 @@ DiscreteFunctionIntegrator::_integrateGlobally() const
                                          template integrate<MeshType>(m_function_id, *m_quadrature_descriptor, *mesh));
 }
 
-template <typename MeshType, typename DataType, typename ValueType>
+template <MeshConcept MeshType, typename DataType, typename ValueType>
 DiscreteFunctionVariant
 DiscreteFunctionIntegrator::_integrate() const
 {
@@ -93,7 +93,7 @@ DiscreteFunctionIntegrator::_integrate() const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 DiscreteFunctionVariant
 DiscreteFunctionIntegrator::_integrate() const
 {
diff --git a/src/scheme/DiscreteFunctionIntegrator.hpp b/src/scheme/DiscreteFunctionIntegrator.hpp
index b832ee6f4dd2ece42628f7768d04e84c8fd3e16b..f294d39e1c75cd10090cde9414bed8288c9fee48 100644
--- a/src/scheme/DiscreteFunctionIntegrator.hpp
+++ b/src/scheme/DiscreteFunctionIntegrator.hpp
@@ -4,6 +4,7 @@
 #include <analysis/IQuadratureDescriptor.hpp>
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IZoneDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 
 #include <memory>
 
@@ -18,16 +19,16 @@ class DiscreteFunctionIntegrator
   std::shared_ptr<const IQuadratureDescriptor> m_quadrature_descriptor;
   const FunctionSymbolId m_function_id;
 
-  template <typename MeshType, typename DataType, typename ValueType = DataType>
+  template <MeshConcept MeshType, typename DataType, typename ValueType = DataType>
   DiscreteFunctionVariant _integrateOnZoneList() const;
 
-  template <typename MeshType, typename DataType, typename ValueType = DataType>
+  template <MeshConcept MeshType, typename DataType, typename ValueType = DataType>
   DiscreteFunctionVariant _integrateGlobally() const;
 
-  template <typename MeshType, typename DataType, typename ValueType = DataType>
+  template <MeshConcept MeshType, typename DataType, typename ValueType = DataType>
   DiscreteFunctionVariant _integrate() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionVariant _integrate() const;
 
  public:
diff --git a/src/scheme/DiscreteFunctionInterpoler.cpp b/src/scheme/DiscreteFunctionInterpoler.cpp
index f8e5a58ca5dd4c736fff870e9d58468535ef4325..8f2e0eef738402c67275fdc0a66ff96fa1de744e 100644
--- a/src/scheme/DiscreteFunctionInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionInterpoler.cpp
@@ -8,7 +8,7 @@
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <utils/Exceptions.hpp>
 
-template <typename MeshType, typename DataType, typename ValueType>
+template <MeshConcept MeshType, typename DataType, typename ValueType>
 DiscreteFunctionVariant
 DiscreteFunctionInterpoler::_interpolateOnZoneList() const
 {
@@ -68,7 +68,7 @@ DiscreteFunctionInterpoler::_interpolateOnZoneList() const
   return DiscreteFunctionP0<ValueType>{m_mesh, cell_value};
 }
 
-template <typename MeshType, typename DataType, typename ValueType>
+template <MeshConcept MeshType, typename DataType, typename ValueType>
 DiscreteFunctionVariant
 DiscreteFunctionInterpoler::_interpolateGlobally() const
 {
@@ -98,7 +98,7 @@ DiscreteFunctionInterpoler::_interpolateGlobally() const
   }
 }
 
-template <typename MeshType, typename DataType, typename ValueType>
+template <MeshConcept MeshType, typename DataType, typename ValueType>
 DiscreteFunctionVariant
 DiscreteFunctionInterpoler::_interpolate() const
 {
@@ -109,7 +109,7 @@ DiscreteFunctionInterpoler::_interpolate() const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 DiscreteFunctionVariant
 DiscreteFunctionInterpoler::_interpolate() const
 {
diff --git a/src/scheme/DiscreteFunctionInterpoler.hpp b/src/scheme/DiscreteFunctionInterpoler.hpp
index 73c3e55f2590b3b241b7425a316afd28b48890f6..95350e34c42f851fedb2c3bbf7420140522d36ef 100644
--- a/src/scheme/DiscreteFunctionInterpoler.hpp
+++ b/src/scheme/DiscreteFunctionInterpoler.hpp
@@ -3,6 +3,7 @@
 
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IZoneDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
 class DiscreteFunctionVariant;
@@ -18,16 +19,16 @@ class DiscreteFunctionInterpoler
   std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
   const FunctionSymbolId m_function_id;
 
-  template <typename MeshType, typename DataType, typename ValueType = DataType>
+  template <MeshConcept MeshType, typename DataType, typename ValueType = DataType>
   DiscreteFunctionVariant _interpolateOnZoneList() const;
 
-  template <typename MeshType, typename DataType, typename ValueType = DataType>
+  template <MeshConcept MeshType, typename DataType, typename ValueType = DataType>
   DiscreteFunctionVariant _interpolateGlobally() const;
 
-  template <typename MeshType, typename DataType, typename ValueType = DataType>
+  template <MeshConcept MeshType, typename DataType, typename ValueType = DataType>
   DiscreteFunctionVariant _interpolate() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionVariant _interpolate() const;
 
  public:
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index 56642754bf97f93dba8e16c0be51ab9bdd0f92cd..df43f52905837036d0ef843fbc5635f0d1274c6d 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -735,12 +735,12 @@ class DiscreteFunctionP0
     Assert(mesh_v->connectivity().id() == cell_value.connectivity_ptr()->id());
   }
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionP0(const std::shared_ptr<const MeshType>& p_mesh)
     : m_mesh_v{p_mesh->meshVariant()}, m_cell_values{m_mesh_v->connectivity()}
   {}
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionP0(const std::shared_ptr<const MeshType>& p_mesh, const CellValue<DataType>& cell_value)
     : m_mesh_v{p_mesh->meshVariant()}, m_cell_values{cell_value}
   {
diff --git a/src/scheme/DiscreteFunctionP0Vector.hpp b/src/scheme/DiscreteFunctionP0Vector.hpp
index aaa549173a7933d6670c2ae621d7caba94397d59..da3904b2ad9853a05a521920c3a3a9434aac14a9 100644
--- a/src/scheme/DiscreteFunctionP0Vector.hpp
+++ b/src/scheme/DiscreteFunctionP0Vector.hpp
@@ -240,12 +240,12 @@ class DiscreteFunctionP0Vector
     Assert(m_mesh->connectivity().id() == cell_arrays.connectivity_ptr()->id());
   }
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionP0Vector(const std::shared_ptr<const MeshType>& p_mesh, size_t size)
     : m_mesh{p_mesh->meshVariant()}, m_cell_arrays{m_mesh->connectivity(), size}
   {}
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionP0Vector(const std::shared_ptr<const MeshType>& p_mesh, const CellArray<DataType>& cell_arrays)
     : m_mesh{p_mesh->meshVariant()}, m_cell_arrays{cell_arrays}
   {
diff --git a/src/scheme/DiscreteFunctionVectorIntegrator.cpp b/src/scheme/DiscreteFunctionVectorIntegrator.cpp
index 74bb43d663cb776d50abce97fe1a34e96e2ceb42..fccc6f79db41dd61cd6feea80baae7c8021487d3 100644
--- a/src/scheme/DiscreteFunctionVectorIntegrator.cpp
+++ b/src/scheme/DiscreteFunctionVectorIntegrator.cpp
@@ -7,7 +7,7 @@
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <utils/Exceptions.hpp>
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorIntegrator::_integrateOnZoneList() const
 {
@@ -68,7 +68,7 @@ DiscreteFunctionVectorIntegrator::_integrateOnZoneList() const
   return DiscreteFunctionP0Vector<DataType>(m_mesh, cell_array);
 }
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorIntegrator::_integrateGlobally() const
 {
@@ -83,7 +83,7 @@ DiscreteFunctionVectorIntegrator::_integrateGlobally() const
                                                                                       *m_quadrature_descriptor, *mesh));
 }
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorIntegrator::_integrate() const
 {
@@ -94,7 +94,7 @@ DiscreteFunctionVectorIntegrator::_integrate() const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorIntegrator::_integrate() const
 {
diff --git a/src/scheme/DiscreteFunctionVectorIntegrator.hpp b/src/scheme/DiscreteFunctionVectorIntegrator.hpp
index e89943cfb949b3ea154bb243feb781ef88cf9674..b9b78525fda0d3ea12c71270aed85e8371c9a3a8 100644
--- a/src/scheme/DiscreteFunctionVectorIntegrator.hpp
+++ b/src/scheme/DiscreteFunctionVectorIntegrator.hpp
@@ -4,6 +4,7 @@
 #include <analysis/IQuadratureDescriptor.hpp>
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IZoneDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
 #include <memory>
@@ -21,16 +22,16 @@ class DiscreteFunctionVectorIntegrator
   std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
   const std::vector<FunctionSymbolId> m_function_id_list;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   DiscreteFunctionVariant _integrateOnZoneList() const;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   DiscreteFunctionVariant _integrateGlobally() const;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   DiscreteFunctionVariant _integrate() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionVariant _integrate() const;
 
  public:
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.cpp b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
index fffead074b14f08acff58d83a29d8a356ffb8801..661eb483c2da2bb7a30ea1e7adb02d770482fa32 100644
--- a/src/scheme/DiscreteFunctionVectorInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.cpp
@@ -7,7 +7,7 @@
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <utils/Exceptions.hpp>
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorInterpoler::_interpolateOnZoneList() const
 {
@@ -65,7 +65,7 @@ DiscreteFunctionVectorInterpoler::_interpolateOnZoneList() const
   return DiscreteFunctionP0Vector<DataType>(p_mesh, cell_array);
 }
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorInterpoler::_interpolateGlobally() const
 {
@@ -81,7 +81,7 @@ DiscreteFunctionVectorInterpoler::_interpolateGlobally() const
                                               ItemType::cell>(m_function_id_list, mesh_data.xj()));
 }
 
-template <typename MeshType, typename DataType>
+template <MeshConcept MeshType, typename DataType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorInterpoler::_interpolate() const
 {
@@ -92,7 +92,7 @@ DiscreteFunctionVectorInterpoler::_interpolate() const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 DiscreteFunctionVariant
 DiscreteFunctionVectorInterpoler::_interpolate() const
 {
diff --git a/src/scheme/DiscreteFunctionVectorInterpoler.hpp b/src/scheme/DiscreteFunctionVectorInterpoler.hpp
index 2999f1220ce4294effea168c9c76cae5c88fb507..ee8536d2d3655548b92d95c45e5d247390eaf9d9 100644
--- a/src/scheme/DiscreteFunctionVectorInterpoler.hpp
+++ b/src/scheme/DiscreteFunctionVectorInterpoler.hpp
@@ -3,6 +3,7 @@
 
 #include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/IZoneDescriptor.hpp>
+#include <mesh/MeshTraits.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 
 class DiscreteFunctionVariant;
@@ -19,16 +20,16 @@ class DiscreteFunctionVectorInterpoler
   std::shared_ptr<const IDiscreteFunctionDescriptor> m_discrete_function_descriptor;
   const std::vector<FunctionSymbolId> m_function_id_list;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   DiscreteFunctionVariant _interpolateOnZoneList() const;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   DiscreteFunctionVariant _interpolateGlobally() const;
 
-  template <typename MeshType, typename DataType>
+  template <MeshConcept MeshType, typename DataType>
   DiscreteFunctionVariant _interpolate() const;
 
-  template <typename MeshType>
+  template <MeshConcept MeshType>
   DiscreteFunctionVariant _interpolate() const;
 
  public:
diff --git a/src/scheme/EdgeIntegrator.hpp b/src/scheme/EdgeIntegrator.hpp
index 2bf287c3d6f61cb83ef3619124334a2124898f8c..7557a7ed851fa999e129e99f68bdd7812dbafd3d 100644
--- a/src/scheme/EdgeIntegrator.hpp
+++ b/src/scheme/EdgeIntegrator.hpp
@@ -67,7 +67,7 @@ class EdgeIntegrator
     EdgeList(const ArrayT<EdgeId>& edge_list) : m_edge_list{edge_list} {}
   };
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   _integrateTo(std::function<OutputType(InputType)> function,
                const IQuadratureDescriptor& quadrature_descriptor,
@@ -116,7 +116,7 @@ class EdgeIntegrator
   }
 
  public:
-  template <typename MeshType, typename OutputArrayT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   integrateTo(const std::function<OutputType(InputType)>& function,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -129,7 +129,7 @@ class EdgeIntegrator
                                          AllEdges<Dimension>{mesh.connectivity()}, value);
   }
 
-  template <typename MeshType, typename OutputArrayT, typename FunctionType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename FunctionType>
   static PUGS_INLINE void
   integrateTo(const FunctionType& function,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -139,7 +139,7 @@ class EdgeIntegrator
     integrateTo(std::function(function), quadrature_descriptor, mesh, value);
   }
 
-  template <typename MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE ArrayT<OutputType>
   integrate(const std::function<OutputType(InputType)>& function,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -152,7 +152,7 @@ class EdgeIntegrator
     return value;
   }
 
-  template <typename MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE auto
   integrate(const FunctionType& function,
             const IQuadratureDescriptor& quadrature_descriptor,
diff --git a/src/scheme/FaceIntegrator.hpp b/src/scheme/FaceIntegrator.hpp
index 657ca7a88a3a959f2856bc34ae8938ab91561785..3fc6b25f7b1bf99419e99e22c3873083832b8a7c 100644
--- a/src/scheme/FaceIntegrator.hpp
+++ b/src/scheme/FaceIntegrator.hpp
@@ -69,7 +69,7 @@ class FaceIntegrator
     FaceList(const ArrayT<FaceId>& face_list) : m_face_list{face_list} {}
   };
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   _tensorialIntegrateTo(std::function<OutputType(InputType)> function,
                         const IQuadratureDescriptor& quadrature_descriptor,
@@ -167,7 +167,7 @@ class FaceIntegrator
     }
   }
 
-  template <typename MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename ListTypeT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   _directIntegrateTo(std::function<OutputType(InputType)> function,
                      const IQuadratureDescriptor& quadrature_descriptor,
@@ -265,7 +265,7 @@ class FaceIntegrator
   }
 
  public:
-  template <typename MeshType, typename OutputArrayT, typename OutputType, typename InputType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename OutputType, typename InputType>
   static PUGS_INLINE void
   integrateTo(const std::function<OutputType(InputType)>& function,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -283,7 +283,7 @@ class FaceIntegrator
     }
   }
 
-  template <typename MeshType, typename OutputArrayT, typename FunctionType>
+  template <MeshConcept MeshType, typename OutputArrayT, typename FunctionType>
   static PUGS_INLINE void
   integrateTo(const FunctionType& function,
               const IQuadratureDescriptor& quadrature_descriptor,
@@ -293,7 +293,7 @@ class FaceIntegrator
     integrateTo(std::function(function), quadrature_descriptor, mesh, value);
   }
 
-  template <typename MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, typename OutputType, typename InputType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE ArrayT<OutputType>
   integrate(const std::function<OutputType(InputType)>& function,
             const IQuadratureDescriptor& quadrature_descriptor,
@@ -312,7 +312,7 @@ class FaceIntegrator
     return value;
   }
 
-  template <typename MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
+  template <MeshConcept MeshType, typename FunctionType, template <typename DataType> typename ArrayT>
   static PUGS_INLINE auto
   integrate(const FunctionType& function,
             const IQuadratureDescriptor& quadrature_descriptor,
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index 509156b1ecbb8eadc3de8c21d69b1d1faa5a9402..a2ebb2b256a3957169cf52457fd87a3d3236e8f3 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -27,7 +27,7 @@
 #include <variant>
 #include <vector>
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class FluxingAdvectionSolver
 {
  private:
@@ -207,7 +207,7 @@ class FluxingAdvectionSolver
   ~FluxingAdvectionSolver() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 FluxingAdvectionSolver<MeshType>::_computeDonorCells(FaceValue<const double> algebraic_fluxing_volumes)
 {
@@ -385,7 +385,7 @@ FluxingAdvectionSolver<Mesh<3>>::_computeAlgebraicFluxingVolume()
   return algebraic_fluxing_volume;
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 FaceValue<double>
 FluxingAdvectionSolver<MeshType>::_computeFluxingVolume(FaceValue<double> algebraic_fluxing_volumes)
 {
@@ -400,7 +400,7 @@ FluxingAdvectionSolver<MeshType>::_computeFluxingVolume(FaceValue<double> algebr
   return algebraic_fluxing_volumes;
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 FluxingAdvectionSolver<MeshType>::_computeCycleNumber(FaceValue<double> fluxing_volumes)
 {
@@ -443,7 +443,7 @@ FluxingAdvectionSolver<MeshType>::_computeCycleNumber(FaceValue<double> fluxing_
   m_cycle_fluxing_volume = fluxing_volumes;
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 FluxingAdvectionSolver<MeshType>::_computeGeometricalData()
 {
@@ -453,7 +453,7 @@ FluxingAdvectionSolver<MeshType>::_computeGeometricalData()
   this->_computeCycleNumber(fluxing_volumes);
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 template <typename CellDataType>
 void
 FluxingAdvectionSolver<MeshType>::_remapOne(const CellValue<const double>& step_Vj,
@@ -633,7 +633,7 @@ FluxingAdvectionSolver<MeshType>::_remapOne(const CellValue<const double>& step_
   old_q = new_q;
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 FluxingAdvectionSolver<MeshType>::_remapAllQuantities()
 {
@@ -696,7 +696,7 @@ FluxingAdvectionSolver<MeshType>::_remapAllQuantities()
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 std::vector<std::shared_ptr<const DiscreteFunctionVariant>>
 FluxingAdvectionSolver<MeshType>::remap(
   const std::vector<std::shared_ptr<const VariableBCDescriptor>>& quantity_list_with_bc)
@@ -748,7 +748,7 @@ FluxingAdvectionSolver<MeshType>::remap(
   return new_variables;
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 template <typename DataType>
 class FluxingAdvectionSolver<MeshType>::InflowValueBoundaryCondition
 {
@@ -779,7 +779,7 @@ class FluxingAdvectionSolver<MeshType>::InflowValueBoundaryCondition
   ~InflowValueBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class FluxingAdvectionSolver<MeshType>::InflowArrayBoundaryCondition
 {
  private:
@@ -809,7 +809,7 @@ class FluxingAdvectionSolver<MeshType>::InflowArrayBoundaryCondition
   ~InflowArrayBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class FluxingAdvectionSolver<MeshType>::OutflowBoundaryCondition
 {
  private:
@@ -830,7 +830,7 @@ class FluxingAdvectionSolver<MeshType>::OutflowBoundaryCondition
   ~OutflowBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class FluxingAdvectionSolver<MeshType>::SymmetryBoundaryCondition
 {
  private:
diff --git a/src/scheme/HyperelasticSolver.cpp b/src/scheme/HyperelasticSolver.cpp
index 5309385ecb210c20dec2995c3b6a603a21ca40a7..d88067f982e244a53cb6f87d062e404392188a69 100644
--- a/src/scheme/HyperelasticSolver.cpp
+++ b/src/scheme/HyperelasticSolver.cpp
@@ -48,7 +48,7 @@ hyperelastic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
     c.meshVariant()->variant());
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticSolverHandler::IHyperelasticSolver
 {
  private:
@@ -549,7 +549,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
   ~HyperelasticSolver()                    = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyPressureBC(const BoundaryConditionList& bc_list,
                                                                           const MeshType& mesh,
@@ -614,7 +614,7 @@ HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyPressureBC(const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyNormalStressBC(const BoundaryConditionList& bc_list,
                                                                               const MeshType& mesh,
@@ -679,7 +679,7 @@ HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyNormalStressBC(co
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applySymmetryBC(const BoundaryConditionList& bc_list,
                                                                           NodeValue<Rdxd>& Ar,
@@ -710,7 +710,7 @@ HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applySymmetryBC(const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 void
 HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyVelocityBC(const BoundaryConditionList& bc_list,
                                                                           NodeValue<Rdxd>& Ar,
@@ -747,7 +747,7 @@ HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyVelocityBC(const
   }
 }
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::FixedBoundaryCondition
 {
  private:
@@ -765,7 +765,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::FixedBoundaryCond
   ~FixedBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::PressureBoundaryCondition
 {
  private:
@@ -819,7 +819,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<Mesh<1>>::PressureBoundaryCo
   ~PressureBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::NormalStressBoundaryCondition
 {
  private:
@@ -873,7 +873,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<Mesh<1>>::NormalStressBounda
   ~NormalStressBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::VelocityBoundaryCondition
 {
  private:
@@ -902,7 +902,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::VelocityBoundaryC
   ~VelocityBoundaryCondition() = default;
 };
 
-template <typename MeshType>
+template <MeshConcept MeshType>
 class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::SymmetryBoundaryCondition
 {
  public:
diff --git a/src/scheme/HyperelasticSolver.hpp b/src/scheme/HyperelasticSolver.hpp
index 75963a51971e623717bb363ec7e6c24fa4cced3e..0fe58e8572aac1a4dabc66da20f2634211f500de 100644
--- a/src/scheme/HyperelasticSolver.hpp
+++ b/src/scheme/HyperelasticSolver.hpp
@@ -1,6 +1,8 @@
 #ifndef HYPERELASTIC_SOLVER_HPP
 #define HYPERELASTIC_SOLVER_HPP
 
+#include <mesh/MeshTraits.hpp>
+
 #include <memory>
 #include <tuple>
 #include <vector>
@@ -65,14 +67,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 <typename MeshType>
+  template <MeshConcept MeshType>
   class HyperelasticSolver;
 
   std::unique_ptr<IHyperelasticSolver> m_hyperelastic_solver;