diff --git a/Data/Lua/demo_input_reaction_diffusion_1D_FHN.lua b/Data/Lua/demo_input_reaction_diffusion_1D_FHN.lua
index b4f142a2c7ab1f9562c526cfa3f50c9871105990..b743c824a16f8d5f63f3fce46aa56d2944cf712a 100644
--- a/Data/Lua/demo_input_reaction_diffusion_1D_FHN.lua
+++ b/Data/Lua/demo_input_reaction_diffusion_1D_FHN.lua
@@ -179,6 +179,39 @@ Domain3 = {
     geometric_element_type_list = { }
 }
 
+
+Domain4 = {
+    -- Index of the geometric mesh upon which the domain is defined (as defined in the present file). Might be
+    -- left empty if domain not limited to one mesh; at most one value is exxpected here.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    mesh_index = { 1 },
+    
+    -- List of dimensions encompassed by the domain. Might be left empty if no restriction at all upon
+    -- dimensions.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    -- Constraint: ops_in(v, {0, 1, 2, 3})
+    dimension_list = { },
+    
+    -- List of mesh labels encompassed by the domain. Might be left empty if no restriction at all upon mesh
+    -- labels. This parameter does not make sense if no mesh is defined for the domain.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    mesh_label_list = { },
+    
+    -- List of geometric element types considered in the domain. Might be left empty if no restriction upon
+    -- these. No constraint is applied at Ops level, as some geometric element types could be added after
+    -- generation of current input parameter file. Current list is below; if an incorrect value is put there it
+    -- will be detected a bit later when the domain object is built.
+    -- The known types when this file was generated are:
+    -- . Point1
+    -- . Segment2, Segment3
+    -- . Triangle3, Triangle6
+    -- . Quadrangle4, Quadrangle8, Quadrangle9
+    -- . Tetrahedron4, Tetrahedron10
+    -- . Hexahedron8, Hexahedron20, Hexahedron27.
+    -- Expected format: {"VALUE1", "VALUE2", ...}
+    geometric_element_type_list = { }
+}
+
 -- FiniteElementSpace1
 FiniteElementSpace1 = {
     -- Index of the god of dof into which the finite element space is defined.
diff --git a/Data/Lua/demo_input_reaction_diffusion_1D_MS.lua b/Data/Lua/demo_input_reaction_diffusion_1D_MS.lua
index bb42425e7b01ef2388d13b2ff3ab1d20d794253a..7dbd8a2938d7d81c651f1cecd8ee11b46847cbd5 100644
--- a/Data/Lua/demo_input_reaction_diffusion_1D_MS.lua
+++ b/Data/Lua/demo_input_reaction_diffusion_1D_MS.lua
@@ -179,6 +179,40 @@ Domain3 = {
     geometric_element_type_list = { }
 }
 
+
+Domain4 = {
+    -- Index of the geometric mesh upon which the domain is defined (as defined in the present file). Might be
+    -- left empty if domain not limited to one mesh; at most one value is exxpected here.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    mesh_index = { 1 },
+    
+    -- List of dimensions encompassed by the domain. Might be left empty if no restriction at all upon
+    -- dimensions.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    -- Constraint: ops_in(v, {0, 1, 2, 3})
+    dimension_list = { },
+    
+    -- List of mesh labels encompassed by the domain. Might be left empty if no restriction at all upon mesh
+    -- labels. This parameter does not make sense if no mesh is defined for the domain.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    mesh_label_list = { },
+    
+    -- List of geometric element types considered in the domain. Might be left empty if no restriction upon
+    -- these. No constraint is applied at Ops level, as some geometric element types could be added after
+    -- generation of current input parameter file. Current list is below; if an incorrect value is put there it
+    -- will be detected a bit later when the domain object is built.
+    -- The known types when this file was generated are:
+    -- . Point1
+    -- . Segment2, Segment3
+    -- . Triangle3, Triangle6
+    -- . Quadrangle4, Quadrangle8, Quadrangle9
+    -- . Tetrahedron4, Tetrahedron10
+    -- . Hexahedron8, Hexahedron20, Hexahedron27.
+    -- Expected format: {"VALUE1", "VALUE2", ...}
+    geometric_element_type_list = { }
+}
+
+
 -- FiniteElementSpace1
 FiniteElementSpace1 = {
     -- Index of the god of dof into which the finite element space is defined.
diff --git a/Data/Lua/demo_input_reaction_diffusion_2D_CRN.lua b/Data/Lua/demo_input_reaction_diffusion_2D_CRN.lua
index 4daa91d351a138821a97c7badc966fc42b6c2002..00d9511f3a4169b3611cdb9f913e6a854b9a713d 100644
--- a/Data/Lua/demo_input_reaction_diffusion_2D_CRN.lua
+++ b/Data/Lua/demo_input_reaction_diffusion_2D_CRN.lua
@@ -177,6 +177,40 @@ Domain3 = {
     geometric_element_type_list = { }
 }
 
+
+Domain4 = {
+    -- Index of the geometric mesh upon which the domain is defined (as defined in the present file). Might be
+    -- left empty if domain not limited to one mesh; at most one value is exxpected here.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    mesh_index = { 1 },
+    
+    -- List of dimensions encompassed by the domain. Might be left empty if no restriction at all upon
+    -- dimensions.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    -- Constraint: ops_in(v, {0, 1, 2, 3})
+    dimension_list = { },
+    
+    -- List of mesh labels encompassed by the domain. Might be left empty if no restriction at all upon mesh
+    -- labels. This parameter does not make sense if no mesh is defined for the domain.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    mesh_label_list = { },
+    
+    -- List of geometric element types considered in the domain. Might be left empty if no restriction upon
+    -- these. No constraint is applied at Ops level, as some geometric element types could be added after
+    -- generation of current input parameter file. Current list is below; if an incorrect value is put there it
+    -- will be detected a bit later when the domain object is built.
+    -- The known types when this file was generated are:
+    -- . Point1
+    -- . Segment2, Segment3
+    -- . Triangle3, Triangle6
+    -- . Quadrangle4, Quadrangle8, Quadrangle9
+    -- . Tetrahedron4, Tetrahedron10
+    -- . Hexahedron8, Hexahedron20, Hexahedron27.
+    -- Expected format: {"VALUE1", "VALUE2", ...}
+    geometric_element_type_list = { }
+}
+
+
 -- FiniteElementSpace1
 FiniteElementSpace1 = {
     -- Index of the god of dof into which the finite element space is defined.
diff --git a/Data/Lua/demo_input_reaction_diffusion_2D_FHN.lua b/Data/Lua/demo_input_reaction_diffusion_2D_FHN.lua
index 4e9708e229e4a2760895d07323de78b0b8b1f2bc..7d7a05615e0589024c7aec0d41855ec49e3e878a 100644
--- a/Data/Lua/demo_input_reaction_diffusion_2D_FHN.lua
+++ b/Data/Lua/demo_input_reaction_diffusion_2D_FHN.lua
@@ -179,6 +179,39 @@ Domain3 = {
     geometric_element_type_list = { }
 }
 
+
+Domain4 = {
+    -- Index of the geometric mesh upon which the domain is defined (as defined in the present file). Might be
+    -- left empty if domain not limited to one mesh; at most one value is exxpected here.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    mesh_index = { 1 },
+    
+    -- List of dimensions encompassed by the domain. Might be left empty if no restriction at all upon
+    -- dimensions.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    -- Constraint: ops_in(v, {0, 1, 2, 3})
+    dimension_list = { },
+    
+    -- List of mesh labels encompassed by the domain. Might be left empty if no restriction at all upon mesh
+    -- labels. This parameter does not make sense if no mesh is defined for the domain.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    mesh_label_list = { },
+    
+    -- List of geometric element types considered in the domain. Might be left empty if no restriction upon
+    -- these. No constraint is applied at Ops level, as some geometric element types could be added after
+    -- generation of current input parameter file. Current list is below; if an incorrect value is put there it
+    -- will be detected a bit later when the domain object is built.
+    -- The known types when this file was generated are:
+    -- . Point1
+    -- . Segment2, Segment3
+    -- . Triangle3, Triangle6
+    -- . Quadrangle4, Quadrangle8, Quadrangle9
+    -- . Tetrahedron4, Tetrahedron10
+    -- . Hexahedron8, Hexahedron20, Hexahedron27.
+    -- Expected format: {"VALUE1", "VALUE2", ...}
+    geometric_element_type_list = { }
+}
+
 -- FiniteElementSpace1
 FiniteElementSpace1 = {
     -- Index of the god of dof into which the finite element space is defined.
diff --git a/Data/Lua/demo_input_reaction_diffusion_2D_MS.lua b/Data/Lua/demo_input_reaction_diffusion_2D_MS.lua
index ae357ed1aff8c1c2428641ee98becdd67afe2d03..f0244a92ebccee185a10137e6cef22c6b1b4fd5d 100644
--- a/Data/Lua/demo_input_reaction_diffusion_2D_MS.lua
+++ b/Data/Lua/demo_input_reaction_diffusion_2D_MS.lua
@@ -179,6 +179,39 @@ Domain3 = {
     geometric_element_type_list = { }
 }
 
+
+Domain4 = {
+    -- Index of the geometric mesh upon which the domain is defined (as defined in the present file). Might be
+    -- left empty if domain not limited to one mesh; at most one value is exxpected here.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    mesh_index = { 1 },
+    
+    -- List of dimensions encompassed by the domain. Might be left empty if no restriction at all upon
+    -- dimensions.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    -- Constraint: ops_in(v, {0, 1, 2, 3})
+    dimension_list = { },
+    
+    -- List of mesh labels encompassed by the domain. Might be left empty if no restriction at all upon mesh
+    -- labels. This parameter does not make sense if no mesh is defined for the domain.
+    -- Expected format: {VALUE1, VALUE2, ...}
+    mesh_label_list = { },
+    
+    -- List of geometric element types considered in the domain. Might be left empty if no restriction upon
+    -- these. No constraint is applied at Ops level, as some geometric element types could be added after
+    -- generation of current input parameter file. Current list is below; if an incorrect value is put there it
+    -- will be detected a bit later when the domain object is built.
+    -- The known types when this file was generated are:
+    -- . Point1
+    -- . Segment2, Segment3
+    -- . Triangle3, Triangle6
+    -- . Quadrangle4, Quadrangle8, Quadrangle9
+    -- . Tetrahedron4, Tetrahedron10
+    -- . Hexahedron8, Hexahedron20, Hexahedron27.
+    -- Expected format: {"VALUE1", "VALUE2", ...}
+    geometric_element_type_list = { }
+}
+
 -- FiniteElementSpace1
 FiniteElementSpace1 = {
     -- Index of the god of dof into which the finite element space is defined.
diff --git a/Sources/ModelInstances/ReactionDiffusion/InputParameterList.hpp b/Sources/ModelInstances/ReactionDiffusion/InputParameterList.hpp
index 1d1a222f3a7474f5342968413aaf13ad2d8da0f8..3dfb3a51edd299bad69c81e8e31dadae1591e309 100644
--- a/Sources/ModelInstances/ReactionDiffusion/InputParameterList.hpp
+++ b/Sources/ModelInstances/ReactionDiffusion/InputParameterList.hpp
@@ -46,7 +46,8 @@ namespace HappyHeart
         {
             highest_dimension = 1,
             horizontal_source = 2,
-            spiral_source = 3
+            spiral_source = 3,
+            full_mesh = 4
         };
         
         
@@ -116,7 +117,8 @@ namespace HappyHeart
             InputParameter::Domain<EnumUnderlyingType(DomainIndex::highest_dimension)>,
             InputParameter::Domain<EnumUnderlyingType(DomainIndex::horizontal_source)>,
             InputParameter::Domain<EnumUnderlyingType(DomainIndex::spiral_source)>,
-                       
+            InputParameter::Domain<EnumUnderlyingType(DomainIndex::full_mesh)>,
+            
             InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::highest_dimension)>,
             InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::horizontal_source)>,
             InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::spiral_source)>,
diff --git a/Sources/ModelInstances/ReactionDiffusion/ReactionDiffusionVariationalFormulation.hpp b/Sources/ModelInstances/ReactionDiffusion/ReactionDiffusionVariationalFormulation.hpp
index fabca2f41ee337565c797401014cea902a0cbb85..1ff15aa9e949a23c1122f5a5032187ce04c82a77 100644
--- a/Sources/ModelInstances/ReactionDiffusion/ReactionDiffusionVariationalFormulation.hpp
+++ b/Sources/ModelInstances/ReactionDiffusion/ReactionDiffusionVariationalFormulation.hpp
@@ -18,6 +18,7 @@
 # include "HappyHeart/Utilities/Filesystem/Folder.hpp"
 
 # include "HappyHeart/Geometry/Domain/Domain.hpp"
+# include "HappyHeart/Geometry/Domain/DomainManager.hpp"
 
 # include "HappyHeart/FormulationSolver/VariationalFormulation.hpp"
 
diff --git a/Sources/ModelInstances/ReactionDiffusion/ReactionDiffusionVariationalFormulation.hxx b/Sources/ModelInstances/ReactionDiffusion/ReactionDiffusionVariationalFormulation.hxx
index 7d8508a65b2661570ffee90321018562e1156c48..7c1ea196ccb3f0a21ad44ac30e894c7dee5eec54 100644
--- a/Sources/ModelInstances/ReactionDiffusion/ReactionDiffusionVariationalFormulation.hxx
+++ b/Sources/ModelInstances/ReactionDiffusion/ReactionDiffusionVariationalFormulation.hxx
@@ -55,12 +55,14 @@ namespace HappyHeart
         template <Advanced::ReactionLawNS::ReactionLawName ReactionLawNameT>
         void ReactionDiffusionVariationalFormulation<ReactionLawNameT>::SupplInit(const InputParameterList& input_parameter_data)
         {
-            const auto& god_of_dof = this->GetGodOfDof();
-            const auto& mesh = god_of_dof.GetGeometricMeshRegion();
+            decltype(auto) domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::full_mesh));
             
             using Diffusion = InputParameter::Diffusion;
-            diffusion_tensor_ = InitScalarParameterFromInputData<Diffusion::Tensor<EnumUnderlyingType(TensorIndex::diffusion_tensor)>>("Diffusion tensor", mesh, input_parameter_data);
-            density_ = InitScalarParameterFromInputData<Diffusion::Density>("Density", mesh, input_parameter_data);
+            diffusion_tensor_ =
+                InitScalarParameterFromInputData<Diffusion::Tensor<EnumUnderlyingType(TensorIndex::diffusion_tensor)>>("Diffusion tensor",
+                                                                                                                       domain,
+                                                                                                                       input_parameter_data);
+            density_ = InitScalarParameterFromInputData<Diffusion::Density>("Density", domain, input_parameter_data);
             
             if (!GetDiffusionTensor().IsConstant())
                 throw Exception("Current Reaction-Diffusion model is restricted to a constant diffusion tensor.", __FILE__, __LINE__);
@@ -112,7 +114,9 @@ namespace HappyHeart
         void ReactionDiffusionVariationalFormulation<ReactionLawNameT>::DefineOperators(const InputParameterList& input_parameter_data)
         {
             const auto& god_of_dof = this->GetGodOfDof();
-            const auto& mesh = god_of_dof.GetGeometricMeshRegion();
+            
+            decltype(auto) domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::full_mesh));
+
             const auto& felt_space_highest_dimension = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::highest_dimension));
             const auto& felt_space_horizontal_source = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::horizontal_source));
             const auto& felt_space_spiral_source = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::spiral_source));
@@ -131,7 +135,7 @@ namespace HappyHeart
                 using parameter_type = InputParameter::TransientSource<EnumUnderlyingType(ForceIndexList::horizontal_source)>;
                 horizontal_source_parameter_ =
                     InitThreeDimensionalParameterFromInputData<parameter_type>("Horizontal source",
-                                                                  mesh,
+                                                                  domain,
                                                                   input_parameter_data);
                 
                 if (horizontal_source_parameter_ != nullptr)
@@ -146,7 +150,7 @@ namespace HappyHeart
                 using parameter_type = InputParameter::TransientSource<EnumUnderlyingType(ForceIndexList::spiral_source)>;
                 spiral_source_parameter_ =
                     InitThreeDimensionalParameterFromInputData<parameter_type>("Spiral source",
-                                                                  mesh,
+                                                                  domain,
                                                                   input_parameter_data);
                 
                 if (spiral_source_parameter_ != nullptr)
@@ -159,7 +163,7 @@ namespace HappyHeart
             }
             
             reaction_law_ = std::make_unique<reaction_law_type>(input_parameter_data,
-                                                                mesh,
+                                                                domain,
                                                                 this->GetTimeManager(),
                                                                 felt_space_highest_dimension.GetQuadratureRulePerTopology());
             
diff --git a/Sources/OperatorInstances/VariationalOperator/NonlinearForm/Local/NonLinearSource/ReactionLaw/Instantiations/CourtemancheRamirezNattel.hpp b/Sources/OperatorInstances/VariationalOperator/NonlinearForm/Local/NonLinearSource/ReactionLaw/Instantiations/CourtemancheRamirezNattel.hpp
index 155b6f0461b09486e605b04e7602582f79196f2c..c480e67db2bf0243677fff50bd275d0da15a5e4e 100644
--- a/Sources/OperatorInstances/VariationalOperator/NonlinearForm/Local/NonLinearSource/ReactionLaw/Instantiations/CourtemancheRamirezNattel.hpp
+++ b/Sources/OperatorInstances/VariationalOperator/NonlinearForm/Local/NonLinearSource/ReactionLaw/Instantiations/CourtemancheRamirezNattel.hpp
@@ -82,7 +82,7 @@ namespace HappyHeart
              */
             template <class InputParameterDataT>
             explicit ReactionLaw(const InputParameterDataT& input_parameter_data,
-                                 const GeometricMeshRegion& mesh,
+                                 const Domain& domain,
                                  const TimeManager& time_manager,
                                  const QuadratureRulePerTopology& default_quadrature_rule_set);
 
diff --git a/Sources/OperatorInstances/VariationalOperator/NonlinearForm/Local/NonLinearSource/ReactionLaw/Instantiations/CourtemancheRamirezNattel.hxx b/Sources/OperatorInstances/VariationalOperator/NonlinearForm/Local/NonLinearSource/ReactionLaw/Instantiations/CourtemancheRamirezNattel.hxx
index ef024b45c5f4cceb699adfd59bc0b230aa65ee78..a4047760657ddcca26a040226f5ca865bd4a7e1a 100644
--- a/Sources/OperatorInstances/VariationalOperator/NonlinearForm/Local/NonLinearSource/ReactionLaw/Instantiations/CourtemancheRamirezNattel.hxx
+++ b/Sources/OperatorInstances/VariationalOperator/NonlinearForm/Local/NonLinearSource/ReactionLaw/Instantiations/CourtemancheRamirezNattel.hxx
@@ -25,7 +25,7 @@ namespace HappyHeart
 
         template <class InputParameterDataT>
         ReactionLaw<ReactionLawName::CourtemancheRamirezNattel>::ReactionLaw(const InputParameterDataT& input_parameter_data,
-                                                                             const GeometricMeshRegion& mesh,
+                                                                             const Domain& domain,
                                                                              const TimeManager& time_manager,
                                                                              const QuadratureRulePerTopology& default_quadrature_rule_set)
         : time_manager_(time_manager)
@@ -64,213 +64,212 @@ namespace HappyHeart
             constexpr const double value_at_quad_jsr = 1.49;
 
 
-
             parameter_list_[EnumUnderlyingType(parameter_index::m)] =
             std::make_unique<ScalarParameterAtQuadPt>("m",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_m,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::h)] =
             std::make_unique<ScalarParameterAtQuadPt>("h",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_h,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::j)] =
             std::make_unique<ScalarParameterAtQuadPt>("j",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_j,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::ao)] =
             std::make_unique<ScalarParameterAtQuadPt>("ao",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_ao,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::io)] =
             std::make_unique<ScalarParameterAtQuadPt>("io",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_io,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::ua)] =
             std::make_unique<ScalarParameterAtQuadPt>("ua",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_ua,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::ui)] =
             std::make_unique<ScalarParameterAtQuadPt>("ui",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_ui,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::xr)] =
             std::make_unique<ScalarParameterAtQuadPt>("xr",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_xr,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::xs)] =
             std::make_unique<ScalarParameterAtQuadPt>("xs",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_xs,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::d)] =
             std::make_unique<ScalarParameterAtQuadPt>("d",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_d,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::f)] =
             std::make_unique<ScalarParameterAtQuadPt>("f",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_f,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::fca)] =
             std::make_unique<ScalarParameterAtQuadPt>("fca",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_fca,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::urel)] =
             std::make_unique<ScalarParameterAtQuadPt>("urel",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_urel,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::vrel)] =
             std::make_unique<ScalarParameterAtQuadPt>("vrel",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_vrel,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::wrel)] =
             std::make_unique<ScalarParameterAtQuadPt>("wrel",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_wrel,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::nai)] =
             std::make_unique<ScalarParameterAtQuadPt>("nai",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_nai,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::nao)] =
             std::make_unique<ScalarParameterAtQuadPt>("nao",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_nao,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::cao)] =
             std::make_unique<ScalarParameterAtQuadPt>("cao",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_cao,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::ki)] =
             std::make_unique<ScalarParameterAtQuadPt>("ki",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_ki,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::ko)] =
             std::make_unique<ScalarParameterAtQuadPt>("ko",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_ko,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::cai)] =
             std::make_unique<ScalarParameterAtQuadPt>("cai",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_cai,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::naiont)] =
             std::make_unique<ScalarParameterAtQuadPt>("naiont",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_naiont,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::kiont)] =
             std::make_unique<ScalarParameterAtQuadPt>("kiont",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_kiont,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::caiont)] =
             std::make_unique<ScalarParameterAtQuadPt>("caiont",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_caiont,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::ileak)] =
             std::make_unique<ScalarParameterAtQuadPt>("ileak",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_ileak,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::iup)] =
             std::make_unique<ScalarParameterAtQuadPt>("iup",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_iup,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::itr)] =
             std::make_unique<ScalarParameterAtQuadPt>("itr",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_itr,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::irel)] =
             std::make_unique<ScalarParameterAtQuadPt>("irel",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_irel,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::nsr)] =
             std::make_unique<ScalarParameterAtQuadPt>("nsr",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_nsr,
                                                       this->GetTimeManager());
 
             parameter_list_[EnumUnderlyingType(parameter_index::jsr)] =
             std::make_unique<ScalarParameterAtQuadPt>("jsr",
-                                                      mesh,
+                                                      domain,
                                                       default_quadrature_rule_set,
                                                       value_at_quad_jsr,
                                                       this->GetTimeManager());