#1022 Make fluid density parameter mesh dependent, to conform with new parameter/domain interface.

Data/Lua/demo_input_poromechanics.lua

@@ -248,6 +248,23 @@ NumberingSubset24 = {
 } -- NumberingSubset24
 
 
+NumberingSubset25 = {
+    
+    -- Name of the numbering subset (not really used; at the moment I just need one input parameter to ground
+    -- the possible values to choose elsewhere).
+    -- Expected format: "VALUE"
+    name = 'fluid_velocity_on_solid',
+    
+    -- Whether a vector defined on this numbering subset might be used to compute a movemesh. If true, a
+    -- FEltSpace featuring this numbering subset will compute additional quantities to enable fast computation.
+    -- This should be false for most numbering subsets, and when it's true the sole unknown involved should be a
+    -- displacement.
+    -- Expected format: 'true' or 'false' (without the quote)
+    do_move_mesh = false
+    
+} -- NumberingSubset25
+
+
 Unknown10 = {
 
 	-- Name of the unknown (used for displays in output).
@@ -1192,15 +1209,15 @@ FiniteElementSpace20 = {
     -- List of all unknowns defined in the finite element space. Unknowns here must be defined in this file as
     -- an 'Unknown' block; expected name/identifier is the name given there.
     -- Expected format: {"VALUE1", "VALUE2", ...}
-    unknown_list = { 'solid_displacement', 'solid_velocity', 'porosity', 'fluid_mass', 'fluid_pressure'  },
+    unknown_list = { 'solid_displacement', 'solid_velocity', 'porosity', 'fluid_mass', 'fluid_pressure', 'fluid_velocity'  },
 
     -- List of the shape function to use for each unknown;
     -- Expected format: {"VALUE1", "VALUE2", ...}
-    shape_function_list = { 'P1', 'P1', 'P1', 'P1', 'P1' },
+    shape_function_list = { 'P1', 'P1', 'P1', 'P1', 'P1', 'P1b' },
 
     -- List of the numbering subset to use for each unknown;
     -- Expected format: {VALUE1, VALUE2, ...}
-    numbering_subset_list = { 20, 21, 22, 23, 24 }
+    numbering_subset_list = { 20, 21, 22, 23, 24, 25 }
 
 } -- FiniteElementSpace20
 

Data/Lua/demo_input_poromechanics_cube.lua

@@ -247,6 +247,23 @@ NumberingSubset24 = {
 } -- NumberingSubset24
 
 
+NumberingSubset25 = {
+    
+    -- Name of the numbering subset (not really used; at the moment I just need one input parameter to ground
+    -- the possible values to choose elsewhere).
+    -- Expected format: "VALUE"
+    name = 'fluid_velocity_on_solid',
+    
+    -- Whether a vector defined on this numbering subset might be used to compute a movemesh. If true, a
+    -- FEltSpace featuring this numbering subset will compute additional quantities to enable fast computation.
+    -- This should be false for most numbering subsets, and when it's true the sole unknown involved should be a
+    -- displacement.
+    -- Expected format: 'true' or 'false' (without the quote)
+    do_move_mesh = false
+    
+} -- NumberingSubset25
+
+
 Unknown10 = {
 
 	-- Name of the unknown (used for displays in output).

Sources/FiniteElement/FiniteElementSpace/FEltSpace.cpp

@@ -246,7 +246,8 @@ namespace HappyHeart
                                }
                                );
         
-        assert(it != unknown_storage.cend());
+        assert(it != unknown_storage.cend() && "Your unknown is probably not present in the \a FEltSpace upon which "
+               "the operator is defined.");
         
         return *it;
     }

Sources/ModelInstances/UnderDevelopment/Poromechanics/ExplicitStep/VariationalFormulation.hxx

@@ -193,7 +193,7 @@ namespace HappyHeart
                 mass_matrix.ZeroEntries(__FILE__, __LINE__);
                 
                 assert(!NumericNS::IsZero(time_step));
-                const double coefficient = variable_holder.GetFluidDensity().GetConstantValue() / time_step;
+                const double coefficient = variable_holder.template GetFluidDensity<MeshIndex::fluid>().GetConstantValue() / time_step;
                 
                 GlobalMatrixWithCoefficient with_coeff_matrix(mass_matrix, coefficient);
                 GetMassOperatorWithPorosity().Assemble(std::make_tuple(std::ref(with_coeff_matrix)));
@@ -209,7 +209,7 @@ namespace HappyHeart
                 auto& mass_matrix = GetNonCstMassMatrixWithPressure();
                 mass_matrix.ZeroEntries(__FILE__, __LINE__);
                 
-                const double coefficient = variable_holder.GetFluidDensity().GetConstantValue() * GetFluidSourceTerm();
+                const double coefficient = variable_holder.template GetFluidDensity<MeshIndex::fluid>().GetConstantValue() * GetFluidSourceTerm();
                 
                 GlobalMatrixWithCoefficient with_coeff_matrix(mass_matrix, coefficient);
                 GetMassOperatorWithPressure().Assemble(std::make_tuple(std::ref(with_coeff_matrix)));
@@ -338,7 +338,7 @@ namespace HappyHeart
             decltype(auto) porosity_data = this->GetPorosityData();
             decltype(auto) porosity = porosity_data.GetNewOnFluidAsParam();
             decltype(auto) variable_holder = this->GetVariableHolder();
-            decltype(auto) fluid_density = variable_holder.GetFluidDensity();
+            decltype(auto) fluid_density = variable_holder.template GetFluidDensity<MeshIndex::fluid>();
             decltype(auto) new_fluid_pressure_data = this->GetNonCstNewFluidPressureData();
             
             mass_operator_with_porosity_ =

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStep/ImplicitStepFluid/NewtonFixedPoint/VariationalFormulationComputeMatrixBlock.hxx

@@ -155,7 +155,8 @@ namespace HappyHeart
                             auto& t12_block = GetNonCstT12Block();
                             t12_block.ZeroEntries(__FILE__, __LINE__);
 
-                            GlobalMatrixWithCoefficient matrix_with_coeff(t12_block, variable_holder.GetFluidDensity().GetConstantValue());
+                            GlobalMatrixWithCoefficient matrix_with_coeff(t12_block,
+                                                                          variable_holder.template GetFluidDensity<MeshIndex::fluid>().GetConstantValue());
                             GetMassDivVelocity().Assemble(std::make_tuple(std::ref(matrix_with_coeff)));
 
                             decltype(auto) monolithic_to_velocity = interpolator_holder.GetMatrixMonolithicToFluidVelocityOnFluid();

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStep/ImplicitStepFluid/NewtonFixedPoint/VariationalFormulationInit.hxx

@@ -28,6 +28,7 @@ namespace HappyHeart
             namespace NewtonFixedPointNS
             {
 
+                
                 template<class HyperelasticLawT>
                 const std::string& VariationalFormulation<HyperelasticLawT>::ClassName()
                 {
@@ -220,7 +221,7 @@ namespace HappyHeart
                             std::make_unique<t21_type>(solid_felt_space, // \todo #1040 Check adequate FEltSpace to use here.
                                                        fluid_mass_unknown,
                                                        fluid_velocity_unknown,
-                                                       variable_holder.GetFluidDensity(),
+                                                       variable_holder.template GetFluidDensity<MeshIndex::solid>(),
                                                        hyperelastic_law,
                                                        this->GetCauchyGreenTensor());
 
@@ -229,7 +230,7 @@ namespace HappyHeart
                     darcy_operator_ =
                         std::make_unique<darcy_operator_type>(monolithic_felt_space,
                                                               fluid_velocity_unknown,
-                                                              variable_holder.GetFluidDensity(),
+                                                              variable_holder.template GetFluidDensity<MeshIndex::fluid>(),
                                                               porosity,
                                                               fluid_velocity_data.GetNewAsParam(),
                                                               solid_velocity_data.GetHalfSumOnFluidAsParam(),
@@ -249,7 +250,7 @@ namespace HappyHeart
                         std::make_unique<GVO::HybridVector<HyperelasticLawT>>(monolithic_felt_space,
                                                                               fluid_mass_unknown,
                                                                               fluid_velocity_unknown,
-                                                                              variable_holder.GetFluidDensity(),
+                                                                              variable_holder.template GetFluidDensity<MeshIndex::fluid>(),
                                                                               porosity,
                                                                               new_fluid_pressure_data.GetNonCstUpdatePressureParamOnFluidCurrent(),
                                                                               fluid_source_term);
@@ -316,7 +317,7 @@ namespace HappyHeart
                                                                                monolithic_felt_space,
                                                                                fluid_velocity_unknown,
                                                                                porosity,
-                                                                               variable_holder.GetFluidDensity(),
+                                                                               variable_holder.template GetFluidDensity<MeshIndex::fluid>(),
                                                                                GetInternalFriction());
 
                         t22_operator_vel_felt_space_ =
@@ -324,7 +325,7 @@ namespace HappyHeart
                                                                                fluid_felt_space,
                                                                                fluid_velocity_unknown,
                                                                                porosity,
-                                                                               variable_holder.GetFluidDensity(),
+                                                                               variable_holder.template GetFluidDensity<MeshIndex::fluid>(),
                                                                                GetInternalFriction());
                     }
 

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStep/ImplicitStepFluid/dH/VariationalFormulationInit.hxx

@@ -170,14 +170,12 @@ namespace HappyHeart
                                                                                                   GetSolidDifferentialVelocity());
                     }
 
-
-
                     // \todo #1022 It was current on fluid here... See if 'usual' cauchy green tensor is ok or if
                     // a one computed on fluid should be used (the latter seems more logical).
                     delta_darcy_operator_ =
                         std::make_unique<darcy_operator_type>(monolithic_felt_space,
                                                           fluid_velocity,
-                                                          variable_holder.GetFluidDensity(),
+                                                          variable_holder.template GetFluidDensity<MeshIndex::fluid>(),
                                                           porosity,
                                                           fluid_velocity_data.GetDeltaAsParam(),
                                                           solid_velocity_data.GetHalfSumOnFluidAsParam(), // \todo #820 not elegant: it is not used as IsFullDarcy should be no... but if the user is wrong it might appear nonetheless! The IsFullDarcy option should be in constructor, even if that means creating more darcy objects.
@@ -199,7 +197,7 @@ namespace HappyHeart
                         std::make_unique<GVO::HybridVector>(monolithic_felt_space,
                                                             fluid_mass_unknown,
                                                             fluid_velocity,
-                                                            variable_holder.GetFluidDensity(),
+                                                            variable_holder.template GetFluidDensity<MeshIndex::fluid>(),
                                                             porosity);
 
                     {

Sources/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStep/VariationalFormulationInit.hxx

@@ -239,7 +239,7 @@ namespace HappyHeart
                                                                                    solid_displacement,
                                                                                    GetInletPressure(),
                                                                                    fluid_mass_data.GetCurrentOnInletBorderAsParam(),
-                                                                                   variable_holder.GetFluidDensity());
+                                                                                   variable_holder.template GetFluidDensity<MeshIndex::solid>());
                 }
                 
                 decltype(auto) solid_varf = GetSolidVarf();

Sources/ModelInstances/UnderDevelopment/Poromechanics/InputParameterList.hpp

@@ -143,7 +143,8 @@ namespace HappyHeart
             solid_velocity = 21,
             porosity_on_solid = 22,
             fluid_mass_on_solid = 23,
-            fluid_pressure_on_solid = 24
+            fluid_pressure_on_solid = 24,
+            fluid_velocity_on_solid = 25
         };
         
         
@@ -171,6 +172,7 @@ namespace HappyHeart
             InputParameter::TimeManager,
             
             InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)>,
+            InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity_on_solid)>,
             InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity_on_robin_interface)>,
             InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::solid_displacement_on_fluid)>,
             InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::porosity)>,

Sources/ModelInstances/UnderDevelopment/Poromechanics/Model.hpp

@@ -189,9 +189,6 @@ namespace HappyHeart
             ///{
             
             
-            //! Accessor to fluid density.
-            const ScalarParameter<>& GetFluidDensity() const noexcept;
-            
             //! Accessor to fluid viscosity.
             const ScalarParameter<>& GetFluidViscosity() const noexcept;
                      

Sources/ModelInstances/UnderDevelopment/Poromechanics/ModelAccessors.hxx

@@ -34,15 +34,7 @@ namespace HappyHeart
             return false; // ie no additional condition
         }
 
-             
-                
-        template<class SolidVariationalFormulationPolicyT>
-        inline const ScalarParameter<>& Model<SolidVariationalFormulationPolicyT>::GetFluidDensity() const noexcept
-        {
-            return GetVariableHolder().GetFluidDensity();
-        }
-        
-        
+               
         template<class SolidVariationalFormulationPolicyT>
         inline const ScalarParameter<>& Model<SolidVariationalFormulationPolicyT>::GetFluidViscosity() const noexcept
         {

Sources/ModelInstances/UnderDevelopment/Poromechanics/ModelInitialize.hxx

@@ -80,7 +80,7 @@ namespace HappyHeart
             {
                 hyperelastic_law_ =
                 std::make_unique<hyperelastic_law_type>(variable_holder.GetSolid(),
-                                                        variable_holder.GetFluidDensity(),
+                                                        variable_holder.template GetFluidDensity<MeshIndex::solid>(),
                                                         variable_holder.GetPenalizationPorosity(),
                                                         variable_holder.GetInitialPorosityValue(),
                                                         this->GetFluidmassData(),
@@ -109,7 +109,7 @@ namespace HappyHeart
                 // Set the initial monolithic data, with a non zero contribution from fluid mass.
                 const double phi_0 = variable_holder.GetInitialPorosityValue();
 
-                decltype(auto) fluid_density = GetFluidDensity().GetConstantValue();
+                decltype(auto) fluid_density = variable_holder.template GetFluidDensity<MeshIndex::solid>().GetConstantValue();
                 
                 decltype(auto) monolithic_data = GetNonCstMonolithicData();
                 monolithic_data.SetInitialValue(GetFluidmassData(),
@@ -125,7 +125,7 @@ namespace HappyHeart
             {
                 new_fluid_pressure_data_ =
                     std::make_unique<DataNS::NewFluidPressure<hyperelastic_law_type>>(GetHyperelasticLaw(),
-                                                                                      variable_holder.GetFluidDensity(),
+                                                                                      variable_holder.template GetFluidDensity<MeshIndex::solid>(),
                                                                                       this->GetCauchyGreenTensor(),
                                                                                       parent::GetTimeManager());
                
@@ -226,7 +226,7 @@ namespace HappyHeart
             
             const auto& time_manager = parent::GetTimeManager();
             const auto& mpi = parent::MpiHappyHeart();
-            decltype(auto) fluid_density = GetFluidDensity();
+            decltype(auto) fluid_density = this->GetVariableHolder().template GetFluidDensity<MeshIndex::solid>();
             decltype(auto) solid_displacement_data = GetSolidDisplacementData();
 
             decltype(auto) solid_god_of_dof = parent::GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));

Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.cpp

@@ -60,14 +60,32 @@ namespace HappyHeart
             }
             # endif // GENERATE_COMP_FILE_4_FREEFEM
             
-            decltype(auto) fluid_domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::fluid_full_mesh));
+            decltype(auto) domain_manager = DomainManager::GetInstance();
+            decltype(auto) fluid_domain = domain_manager.GetDomain(EnumUnderlyingType(DomainIndex::fluid_full_mesh));
+            decltype(auto) solid_domain = domain_manager.GetDomain(EnumUnderlyingType(DomainIndex::solid_full_mesh));
             
-            fluid_density_ = InitScalarParameterFromInputData<InputParameter::Fluid::Density>("Fluid density",
-                                                                                              fluid_domain,
-                                                                                              input_parameter_data);
+            {
+                constexpr auto fluid_mesh_index = EnumUnderlyingType(MeshIndex::fluid) - 1u;
+                assert(fluid_mesh_index < 2u);
+                
+                fluid_density_[fluid_mesh_index] =
+                    InitScalarParameterFromInputData<InputParameter::Fluid::Density>("Fluid density on fluid mesh",
+                                                                                     fluid_domain,
+                                                                                     input_parameter_data);
+
+                constexpr auto solid_mesh_index = EnumUnderlyingType(MeshIndex::solid) - 1u;
+                assert(solid_mesh_index < 2u);
+
+                
+                fluid_density_[solid_mesh_index] =
+                        InitScalarParameterFromInputData<InputParameter::Fluid::Density>("Fluid density on solid mesh",
+                                                                                         solid_domain,
+                                                                                         input_parameter_data);
+            }
             
             
-            if (!GetFluidDensity().IsConstant())
+            if (!GetFluidDensity<MeshIndex::fluid>().IsConstant()) // no need to test solid as well, as it is extracted
+                                                                   // from the very same data.
                 throw Exception("Current model is restricted to a constant fluid density! (in fact a density defined "
                                 "at dof would still be fine, but not at quadrature point...)",
                                 __FILE__, __LINE__);
@@ -80,8 +98,6 @@ namespace HappyHeart
                 throw Exception("Current model is restricted to a constant fluid viscosity!", __FILE__, __LINE__);
         
             
-            decltype(auto) solid_domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::solid_full_mesh));
-            
             {
                decltype(auto) felt_space = solid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::solid));
                 solid_ = std::make_unique<Solid>(input_parameter_data,

Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hpp

@@ -128,6 +128,7 @@ namespace HappyHeart
             
                       
             //! Accessor to fluid density.
+            template<MeshIndex MeshIndexT>
             const ScalarParameter<>& GetFluidDensity() const noexcept;
             
             //! Accessor to fluid viscosity.
@@ -149,8 +150,11 @@ namespace HappyHeart
             
         private:
             
+            static_assert(EnumUnderlyingType(MeshIndex::fluid) == 1, "Internal storage assumes that.");
+            static_assert(EnumUnderlyingType(MeshIndex::solid) == 2, "Internal storage assumes that.");
+            
             //! Density of the fluid.
-            ScalarParameter<>::unique_ptr fluid_density_ = nullptr;
+            ScalarParameter<>::array_unique_ptr<2> fluid_density_ = { { nullptr, nullptr } };
             
             //! Viscosity of the fluid.
             ScalarParameter<>::unique_ptr fluid_viscosity_ = nullptr;

Sources/ModelInstances/UnderDevelopment/Poromechanics/VariableHolder.hxx

@@ -20,12 +20,13 @@ namespace HappyHeart
     {
         
         
-             
-        
+        template<MeshIndex MeshIndexT>
         inline const ScalarParameter<>& VariableHolder::GetFluidDensity() const noexcept
         {
-            assert(!(!fluid_density_));
-            return *fluid_density_;
+            constexpr auto index = EnumUnderlyingType(MeshIndexT) - 1u;
+            static_assert(index == 0u || index == 1u, "Only two values are possible.");
+            assert(!(!fluid_density_[index]));
+            return *fluid_density_[index];
         }