Commit f6e9a0d8 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1049 Fix poromechanics model.

parent c8825b39
......@@ -563,6 +563,41 @@ Domain14 = {
} -- Domain14
-- Dirichlet z -- Solely for 3D; use inexistant mesh label index...
Domain15 = {
-- 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 expected 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 = {}
} -- Domain15
......@@ -749,6 +784,43 @@ Domain24 = {
} -- Domain24
Domain25 = {
-- 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 expected here.
-- Expected format: {VALUE1, VALUE2, ...}
mesh_index = { 2 },
-- 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 = {}
} -- Domain25
EssentialBoundaryCondition10 = {
-- Name of the boundary condition (must be unique).
......
......@@ -560,7 +560,40 @@ Domain14 = {
} -- Domain14
-- Dirichlet z -- Solely for 3D; use inexistant mesh label index...
Domain15 = {
-- 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 expected 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 = {}
} -- Domain15
Domain20 = {
......@@ -746,6 +779,42 @@ Domain24 = {
} -- Domain24
Domain25 = {
-- 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 expected here.
-- Expected format: {VALUE1, VALUE2, ...}
mesh_index = { 2 },
-- 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 = {}
} -- Domain25
EssentialBoundaryCondition10 = {
-- Name of the boundary condition (must be unique).
......
......@@ -64,8 +64,8 @@ namespace HappyHeart
{
decltype(auto) unknown_manager = UnknownManager::GetInstance();
using vector_at_dof_type = ParameterAtDof<ParameterNS::Type::vector>;
decltype(auto) mesh = fluid_god_of_dof.GetGeometricMeshRegion();
decltype(auto) domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::fluid_full_mesh));
decltype(auto) monolithic_felt_space =
fluid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::fluid_monolithic));
......@@ -74,13 +74,13 @@ namespace HappyHeart
unknown_manager.GetUnknown(EnumUnderlyingType(UnknownIndex::fluid_velocity));
new_as_param_ = std::make_unique<vector_at_dof_type>("Fluid velocity solution",
mesh,
domain,
monolithic_felt_space,
fluid_velocity_unknown,
monolithic_data.GetNew());
delta_as_param_ = std::make_unique<vector_at_dof_type>("Fluid delta velocity solution",
mesh,
domain,
monolithic_felt_space,
fluid_velocity_unknown,
monolithic_data.GetDelta());
......@@ -90,7 +90,7 @@ namespace HappyHeart
velFtr_as_param_ = std::make_unique<vector_at_dof_type>("velFtr",
mesh,
domain,
fluid_felt_space,
fluid_velocity_unknown,
GetVelFtr());
......
......@@ -48,8 +48,6 @@ namespace HappyHeart
decltype(auto) solid_numbering_subset =
solid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_on_solid));
decltype(auto) solid_mesh = solid_god_of_dof.GetGeometricMeshRegion();
constexpr auto solid_index = EnumUnderlyingType(InternalMeshIndex::solid);
constexpr auto fluid_index = EnumUnderlyingType(InternalMeshIndex::fluid);
......@@ -95,34 +93,36 @@ namespace HappyHeart
const auto& fluid_mass_unknown =
unknown_manager.GetUnknown(EnumUnderlyingType(UnknownIndex::fluid_mass));
new_as_param_[solid_index] = std::make_unique<scalar_at_dof_type>("Fluid mass",
solid_mesh,
felt_space,
fluid_mass_unknown,
decltype(auto) solid_domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::solid_full_mesh));
new_as_param_[solid_index] = std::make_unique<scalar_at_dof_type>("Fluid mass",
solid_domain,
felt_space,
fluid_mass_unknown,
GetNew());
current_value_as_param_[solid_index] = std::make_unique<scalar_at_dof_type>("Current fluid mass",
solid_mesh,
felt_space,
fluid_mass_unknown,
GetCurrent());
solid_domain,
felt_space,
fluid_mass_unknown,
GetCurrent());
former_value_as_param_[solid_index] = std::make_unique<scalar_at_dof_type>("Former fluid mass",
solid_mesh,
felt_space,
fluid_mass_unknown,
GetFormer());
solid_domain,
felt_space,
fluid_mass_unknown,
GetFormer());
new_minus_current_as_param_ =
std::make_unique<scalar_at_dof_type>("New - current fluid mass",
solid_mesh,
solid_domain,
felt_space,
fluid_mass_unknown,
GetNewMinusCurrentOnSolid());
current_minus_former_as_param_ =
std::make_unique<scalar_at_dof_type>("Current - former fluid mass",
solid_mesh,
solid_domain,
felt_space,
fluid_mass_unknown,
GetCurrentMinusFormerOnSolid());
......@@ -132,14 +132,14 @@ namespace HappyHeart
current_value_on_inlet_border_as_param_ =
std::make_unique<scalar_at_dof_type>("Current fluid mass on inlet border",
solid_mesh,
solid_domain,
inlet_border_felt_space,
fluid_mass_unknown,
GetCurrent());
delta_on_solid_as_param_ =
std::make_unique<scalar_at_dof_type>("Delta fluid mass",
solid_mesh,
solid_domain,
felt_space,
fluid_mass_unknown,
GetDeltaOnSolid());
......@@ -218,7 +218,6 @@ namespace HappyHeart
return *current_minus_former_;
}
} // namespace DataNS
......
......@@ -58,7 +58,14 @@ namespace HappyHeart
decltype(auto) hyperelastic_law = GetHyperelasticLaw();
decltype(auto) domain_manager = DomainManager::GetInstance();
decltype(auto) solid_domain =
domain_manager.GetDomain(EnumUnderlyingType(DomainIndex::solid_full_mesh));
decltype(auto) fluid_domain =
domain_manager.GetDomain(EnumUnderlyingType(DomainIndex::fluid_full_mesh));
decltype(auto) fluid_felt_space =
fluid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::fluid));
......@@ -89,9 +96,10 @@ namespace HappyHeart
constexpr auto time_label_type = TimeLabel2::new_value;
constexpr auto index = EnumUnderlyingType(time_label_type);
pressure_parameter_on_solid_[index] =
std::make_unique<pressure_parameter_type>("Fluid pressure on solid",
solid_god_of_dof.GetGeometricMeshRegion(),
solid_domain,
solid_felt_space.GetQuadratureRulePerTopology(),
0.,
time_manager);
......@@ -105,9 +113,10 @@ namespace HappyHeart
hyperelastic_law,
this->GetCauchyGreenTensor());
pressure_parameter_on_fluid_[index] =
std::make_unique<pressure_parameter_type>("Fluid pressure on fluid",
fluid_god_of_dof.GetGeometricMeshRegion(),
fluid_domain,
fluid_felt_space.GetQuadratureRulePerTopology(),
0.,
time_manager);
......@@ -128,7 +137,7 @@ namespace HappyHeart
pressure_parameter_on_solid_[index] =
std::make_unique<pressure_parameter_type>("Fluid pressure on solid",
solid_god_of_dof.GetGeometricMeshRegion(),
solid_domain,
solid_felt_space.GetQuadratureRulePerTopology(),
0.,
time_manager);
......@@ -144,7 +153,7 @@ namespace HappyHeart
pressure_parameter_on_fluid_[index] =
std::make_unique<pressure_parameter_type>("Fluid pressure on fluid",
fluid_god_of_dof.GetGeometricMeshRegion(),
fluid_domain,
fluid_felt_space.GetQuadratureRulePerTopology(),
0.,
time_manager);
......
......@@ -58,20 +58,22 @@ namespace HappyHeart
{
using scalar_at_dof_type = ParameterAtDof<ParameterNS::Type::scalar>;
decltype(auto) mesh = solid_god_of_dof.GetGeometricMeshRegion();
decltype(auto) felt_space =
solid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::solid));
decltype(auto) solid_domain =
DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::solid_full_mesh));
new_on_solid_as_param_ = std::make_unique<scalar_at_dof_type>("Porosity solution",
mesh,
felt_space,
porosity_unknown,
new_value);
solid_domain,
felt_space,
porosity_unknown,
new_value);
current_value_on_solid_as_param_ =
std::make_unique<scalar_at_dof_type>("Porosity solution",
mesh,
solid_domain,
felt_space,
porosity_unknown,
new_value);
......@@ -91,24 +93,24 @@ namespace HappyHeart
current_value_on_fluid_ = std::make_unique<GlobalVector>(new_value);
{
decltype(auto) mesh = fluid_god_of_dof.GetGeometricMeshRegion();
decltype(auto) felt_space =
fluid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::fluid));
decltype(auto) felt_space_dim_minus_1 =
fluid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::fluid_inlet_border));
decltype(auto) fluid_domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::fluid_full_mesh));
new_on_fluid_as_param_ = std::make_unique<fluid_param_at_dof_type>("Porosity solution on fluid",
mesh,
felt_space,
felt_space_dim_minus_1,
porosity_unknown,
new_value);
fluid_domain,
felt_space,
felt_space_dim_minus_1,
porosity_unknown,
new_value);
current_value_on_fluid_as_param_ =
std::make_unique<fluid_param_at_dof_type>("Porosity solution",
mesh,
fluid_domain,
felt_space,
felt_space_dim_minus_1,
porosity_unknown,
......
......@@ -70,8 +70,6 @@ namespace HappyHeart
{
decltype(auto) unknown_manager = UnknownManager::GetInstance();
using vector_at_dof_type = ParameterAtDof<ParameterNS::Type::vector>;
decltype(auto) mesh = fluid_god_of_dof.GetGeometricMeshRegion();
decltype(auto) fluid_felt_space =
fluid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::fluid));
......@@ -79,8 +77,10 @@ namespace HappyHeart
const auto& fluid_velocity_unknown =
unknown_manager.GetUnknown(EnumUnderlyingType(UnknownIndex::fluid_velocity));
decltype(auto) fluid_domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::fluid_full_mesh));
half_sum_on_fluid_as_param_ = std::make_unique<vector_at_dof_type>("velSHalfVhf",
mesh,
fluid_domain,
fluid_felt_space,
fluid_velocity_unknown,
GetHalfSumOnFluid());
......
......@@ -147,17 +147,22 @@ namespace HappyHeart
{
using parameter_type = ::HappyHeart::Internal::ParameterNS::ParameterInstance
<
ParameterNS::Type::scalar,
ParameterNS::Policy::Constant,
ParameterNS::TimeDependencyNS::None
ParameterNS::Type::scalar,
ParameterNS::Policy::Constant,
ParameterNS::TimeDependencyNS::None
>;
decltype(auto) domain_manager = DomainManager::GetInstance();
decltype(auto) fluid_domain =
domain_manager.GetDomain(EnumUnderlyingType(DomainIndex::fluid_full_mesh));
pseudo_young_modulus_ = std::make_unique<parameter_type>("Pseudo young modulus",
mesh,
fluid_domain,
2. * variable_holder.GetFluidViscosity().GetConstantValue());
pseudo_poisson_ratio_ = std::make_unique<parameter_type>("Zero poisson ratio",
mesh,
fluid_domain,
0.);
}
......
......@@ -136,6 +136,7 @@ namespace HappyHeart
unknown_manager.GetUnknown(EnumUnderlyingType(UnknownIndex::fluid_mass));
decltype(auto) mesh = god_of_dof.GetGeometricMeshRegion();
assert(mesh.GetUniqueId() == EnumUnderlyingType(MeshIndex::fluid));
namespace GVO = GlobalVariationalOperatorNS;
......@@ -156,12 +157,14 @@ namespace HappyHeart
decltype(auto) solid_velocity_data = this->GetSolidVelocityData();
decltype(auto) solid_displacement_data = this->GetSolidDisplacementData();
decltype(auto) fluid_domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::fluid_full_mesh));
{
using vector_at_dof_type = ParameterAtDof<ParameterNS::Type::vector>;
solid_differential_velocity_param_ = std::make_unique<vector_at_dof_type>("solid_differential_velocity",
mesh,
fluid_domain,
monolithic_felt_space,
fluid_velocity,
GetSolidDifferentialVelocity());
......
......@@ -61,15 +61,17 @@ namespace HappyHeart
{
fluid = 10,
fluid_dirichlet_y = 11,
fluid_dirichlet_z = 14,
fluid_inlet_border = 12,
fluid_robin_interface = 13,
fluid_dirichlet_z = 14,
fluid_full_mesh = 15,
solid = 20,
solid_inlet_border = 21,
solid_dirichlet_x = 22, // \todo #963 CHECK IF IT IS REALLY USEFUL (ALMOST same as 21)
solid_dirichlet_y = 23,
solid_dirichlet_z = 24
solid_dirichlet_z = 24,
solid_full_mesh = 25
};
......@@ -197,11 +199,13 @@ namespace HappyHeart
InputParameter::Domain<EnumUnderlyingType(DomainIndex::fluid_dirichlet_z)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::fluid_inlet_border)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::fluid_robin_interface)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::fluid_full_mesh)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::solid)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::solid_inlet_border)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::solid_dirichlet_x)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::solid_dirichlet_y)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::solid_dirichlet_z)>,
InputParameter::Domain<EnumUnderlyingType(DomainIndex::solid_full_mesh)>,
InputParameter::DirichletBoundaryCondition<EnumUnderlyingType(BoundaryConditionIndex::fluid_dirichlet_y)>,
InputParameter::DirichletBoundaryCondition<EnumUnderlyingType(BoundaryConditionIndex::fluid_dirichlet_z)>,
......
......@@ -147,9 +147,6 @@ namespace HappyHeart
template<class SolidVariationalFormulationPolicyT>
void Model<SolidVariationalFormulationPolicyT>::InitParameters(const InputParameterList& input_parameter_data)
{
decltype(auto) fluid_god_of_dof = parent::GetGodOfDof(EnumUnderlyingType(MeshIndex::fluid));
decltype(auto) fluid_mesh = fluid_god_of_dof.GetGeometricMeshRegion();
decltype(auto) solid_god_of_dof = parent::GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));
decltype(auto) solid_mesh = solid_god_of_dof.GetGeometricMeshRegion();
......@@ -157,6 +154,9 @@ namespace HappyHeart
decltype(auto) time_manager = parent::GetTimeManager();
decltype(auto) solid_domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::solid_full_mesh));
decltype(auto) fluid_domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::fluid_full_mesh));
{
LocalVector initial_value;
......@@ -179,7 +179,7 @@ namespace HappyHeart
cauchy_green_tensor_ =
std::make_unique<cauchy_green_tensor_type>("Cauchy-Green tensor",
solid_mesh,
solid_domain,
quadrature_rule_per_topology,
initial_value,
time_manager);
......@@ -202,7 +202,7 @@ namespace HappyHeart
parameter_type,
ParameterNS::TimeDependencyFunctor
>("inlet_pressure",
fluid_mesh,
fluid_domain,
input_parameter_data);
// \todo 820 Put that in input parameter file?
......
......@@ -76,8 +76,6 @@ namespace HappyHeart
const auto& ref_felt = elementary_data.GetRefFElt(GetNthUnknown());
const unsigned int Ncomponent = elementary_data.GetMeshDimension();
decltype(auto) cauchy_green_tensor = this->GetCauchyGreenTensor();
decltype(auto) invariant_holder = GetNonCstInvariantHolder();
......
......@@ -38,10 +38,7 @@ namespace HappyHeart
{
decltype(auto) god_of_dof_manager = GodOfDofManager::GetInstance();
decltype(auto) fluid_god_of_dof = god_of_dof_manager.GetGodOfDof(EnumUnderlyingType(MeshIndex::fluid));
decltype(auto) solid_god_of_dof = god_of_dof_manager.GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));
decltype(auto) fluid_mesh = fluid_god_of_dof.GetGeometricMeshRegion();
decltype(auto) solid_mesh = solid_god_of_dof.GetGeometricMeshRegion();
# ifdef GENERATE_COMP_FILE_4_FREEFEM
{
......@@ -63,9 +60,11 @@ namespace HappyHeart
}
# endif // GENERATE_COMP_FILE_4_FREEFEM
decltype(auto) fluid_domain = DomainManager::GetInstance().GetDomain(EnumUnderlyingType(DomainIndex::fluid_full_mesh));