Commit 0cc9327c authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#820 Mimic the Dirichlet vector of Freefem script.

parent 431180ef
......@@ -43,21 +43,19 @@ NumberingSubset10 = {
} -- NumberingSubset10
-- NumberingSubset11 = {
--
-- -- 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_mass',
--
-- -- 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
--
-- } -- NumberingSubset11
NumberingSubset11 = {
-- 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_robin_interface',
-- 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
} -- NumberingSubset11
NumberingSubset13 = {
......@@ -484,12 +482,12 @@ Domain13 = {
-- dimensions.
-- Expected format: {VALUE1, VALUE2, ...}
-- Constraint: ops_in(v, {0, 1, 2, 3})
dimension_list = { },
dimension_list = { 1 },
-- 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 = { 1, 3 },
mesh_label_list = { 1, 3 },
-- 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
......@@ -910,7 +908,32 @@ FiniteElementSpace12 = {
numbering_subset_list = { 15, 15, 15, 13 }
}
-- FiniteElementSpace13 -- Fluid robin interface
FiniteElementSpace13 = {
-- Index of the god of dof into which the finite element space is defined.
-- Expected format: VALUE
god_of_dof_index = 1,
-- Index of the domain onto which the finite element space is defined. This domain must be unidimensional.
-- Expected format: VALUE
domain_index = 13,
-- 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 = { 'fluid_velocity' },
-- List of the shape function to use for each unknown;
-- Expected format: {"VALUE1", "VALUE2", ...}
-- shape_function_list = { 'P2', 'P1' },
shape_function_list = { 'P1' },
-- List of the numbering subset to use for each unknown;
-- Expected format: {VALUE1, VALUE2, ...}
numbering_subset_list = { 11 }
} -- FiniteElementSpace13
FiniteElementSpace20 = {
......
......@@ -374,6 +374,12 @@ namespace HappyHeart
//! Constant accessor to the interpolator to expand from fluid mass to monolithic (on fluid mesh).
const ConformInterpolatorNS::SubsetOrSuperset& GetMass2Monolithic() const noexcept;
/*!
* \brief Constant accessor to the interpolator that expands a velocity vector on the robin interface
* to the whole mesh.
*/
const ConformInterpolatorNS::SubsetOrSuperset& GetVelocityRobinInterfaceToFullMesh() const noexcept;
public:
//! Constant accessor to the interpolator to reduce from monolithic to only fluid mass (on fluid mesh).
......@@ -524,6 +530,22 @@ namespace HappyHeart
*/
GlobalVector& GetNonCstNonHomogeneousRobinVector() noexcept;
/*!
* \brief Constant accessor to the interpolator that limit a velocity vector to the dofs on the robin
* interface.
*/
const ConformInterpolatorNS::SubsetOrSuperset& GetVelocityToRobinInterface() const noexcept;
//! Constant accessor to the matrix to go from (whole mesh, monolithic) to (robin interface, velocity).
const GlobalMatrix& GetMonolithic2RobinVelocity() const noexcept;
/*!
* \brief Non constant accessor to the matrix to go from (whole mesh, monolithic) to (robin interface,
* velocity).
*/
GlobalMatrix& GetNonCstMonolithic2RobinVelocity() noexcept;
///@}
......@@ -564,6 +586,21 @@ namespace HappyHeart
GlobalVector& GetNonCstSolidDifferentialVelocity() noexcept;
/*!
* \brief Constant accessor to the differential velocity of the solid ('dispSdeltar' in dH function in
* Freefem script).
*/
const GlobalVector& GetSolidDifferentialVelocityOnRobinInterface() const noexcept;
/*!
* \brief Non constant accessor to the differential velocity of the solid ('dispSdeltar' in dH function
* in Freefem script).
*/
GlobalVector& GetNonCstSolidDifferentialVelocityOnRobinInterface() noexcept;
private:
......@@ -711,6 +748,12 @@ namespace HappyHeart
//! Test interpolator to expand from fluid velocity to monolithic.
ConformInterpolatorNS::SubsetOrSuperset::unique_ptr velocity_2_monolithic_ = nullptr;
//! Interpolator that limit a velocity vector to the dofs on the robin interface.
ConformInterpolatorNS::SubsetOrSuperset::unique_ptr velocity_to_robin_interface_ = nullptr;
//! Interpolator that expands a velocity vector on the robin interface to the whole mesh.
ConformInterpolatorNS::SubsetOrSuperset::unique_ptr velocity_robin_interface_to_full_mesh_ = nullptr;
///@}
......@@ -792,7 +835,8 @@ namespace HappyHeart
GlobalMatrix::unique_ptr t32_block_ = nullptr;
//! Matrix to go from (whole mesh, monolithic) to (robin interface, velocity).
GlobalMatrix::unique_ptr monolithic_2_robin_velocity_ = nullptr;
//! T21 on the non monolithic \a FEltSpace.
......@@ -887,6 +931,9 @@ namespace HappyHeart
//! Differential velocity of the solid ('dispSdeltar' in dH function in Freefem script).
GlobalVector::unique_ptr solid_differential_velocity_ = nullptr;
//! Differential velocity of the solid ('dispSdeltar' in dH function in Freefem script) on robin interface.
GlobalVector::unique_ptr solid_differential_velocity_on_robin_interface_ = nullptr;
};
......
......@@ -166,31 +166,48 @@ namespace HappyHeart
GetDifferentialDarcyOperator().Assemble(std::make_tuple(std::ref(contribution_with_coeff)));
}
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "rhs_" + GetIterationTag() + ".m";
std::cout << "FILE = " << filename << std::endl;
rhs.View(parent::MpiHappyHeart(),
filename,
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
{
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "rhs_" + GetIterationTag() + ".m";
rhs.View(parent::MpiHappyHeart(),
filename,
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
// auto& robin_vector = GetNonCstNonHomogeneousRobinVector();
// robin_vector.Copy(GetSolidDifferentialVelocity(), __FILE__, __LINE__);
//
decltype(auto) bc_manager = DirichletBoundaryConditionManager::GetInstance();
decltype(auto) non_homogeneous_bc =
bc_manager.GetNonCstDirichletBoundaryCondition(EnumUnderlyingType(BoundaryConditionIndex::fluid_robin_interface_dH));
non_homogeneous_bc.UpdateValues(GetSolidDifferentialVelocity());
auto& solid_differential_vel_on_robin_interface = GetNonCstSolidDifferentialVelocityOnRobinInterface();
Wrappers::Petsc::MatMult(GetMonolithic2RobinVelocity(),
GetSolidDifferentialVelocity(),
solid_differential_vel_on_robin_interface,
__FILE__, __LINE__);
non_homogeneous_bc.UpdateValues(solid_differential_vel_on_robin_interface);
auto& robin_vector = GetNonCstNonHomogeneousRobinVector();
robin_vector.ZeroEntries(__FILE__, __LINE__);
god_of_dof.template ApplyBoundaryCondition<BoundaryConditionMethod::penalization>(non_homogeneous_bc, robin_vector);
god_of_dof.template ApplyBoundaryCondition<BoundaryConditionMethod::pseudo_elimination>(non_homogeneous_bc, robin_vector);
robin_vector.View(this->MpiHappyHeart(), __FILE__, __LINE__);
{
const auto filename = parent::GetOutputDirectory(numbering_subset)
+ "/"
+ DifferentialPreffix(differential::yes)
+ "dirichletFunctionMixtVariables_" + GetIterationTag() + ".m";
robin_vector.View(parent::MpiHappyHeart(),
filename,
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
// \todo #820 At the moment, I accept my norm is not the same as Bruno's; see what Petsc gives when
......
......@@ -579,7 +579,8 @@ namespace HappyHeart
template<class HyperelasticLawT>
inline const GlobalVector& VariationalFormulation<HyperelasticLawT>::GetSolidDifferentialVelocity() const noexcept
inline const GlobalVector& VariationalFormulation<HyperelasticLawT>
::GetSolidDifferentialVelocity() const noexcept
{
assert(!(!solid_differential_velocity_));
return *solid_differential_velocity_;
......@@ -587,14 +588,34 @@ namespace HappyHeart
template<class HyperelasticLawT>
inline GlobalVector& VariationalFormulation<HyperelasticLawT>::GetNonCstSolidDifferentialVelocity() noexcept
inline GlobalVector& VariationalFormulation<HyperelasticLawT>
::GetNonCstSolidDifferentialVelocity() noexcept
{
return const_cast<GlobalVector&>(GetSolidDifferentialVelocity());
}
template<class HyperelasticLawT>
inline const GlobalVector& VariationalFormulation<HyperelasticLawT>::GetNonHomogeneousRobinVector() const noexcept
inline const GlobalVector& VariationalFormulation<HyperelasticLawT>
::GetSolidDifferentialVelocityOnRobinInterface() const noexcept
{
assert(!(!solid_differential_velocity_on_robin_interface_));
return *solid_differential_velocity_on_robin_interface_;
}
template<class HyperelasticLawT>
inline GlobalVector& VariationalFormulation<HyperelasticLawT>
::GetNonCstSolidDifferentialVelocityOnRobinInterface() noexcept
{
return const_cast<GlobalVector&>(GetSolidDifferentialVelocityOnRobinInterface());
}
template<class HyperelasticLawT>
inline const GlobalVector& VariationalFormulation<HyperelasticLawT>
::GetNonHomogeneousRobinVector() const noexcept
{
assert(!(!non_homogeneous_robin_vector_));
return *non_homogeneous_robin_vector_;
......@@ -602,12 +623,48 @@ namespace HappyHeart
template<class HyperelasticLawT>
inline GlobalVector& VariationalFormulation<HyperelasticLawT>::GetNonCstNonHomogeneousRobinVector() noexcept
inline GlobalVector& VariationalFormulation<HyperelasticLawT>
::GetNonCstNonHomogeneousRobinVector() noexcept
{
return const_cast<GlobalVector&>(GetNonHomogeneousRobinVector());
}
template<class HyperelasticLawT>
inline const ConformInterpolatorNS::SubsetOrSuperset& VariationalFormulation<HyperelasticLawT>
::GetVelocityToRobinInterface() const noexcept
{
assert(!(!velocity_to_robin_interface_));
return *velocity_to_robin_interface_;
}
template<class HyperelasticLawT>
inline const ConformInterpolatorNS::SubsetOrSuperset& VariationalFormulation<HyperelasticLawT>
::GetVelocityRobinInterfaceToFullMesh() const noexcept
{
assert(!(!velocity_robin_interface_to_full_mesh_));
return *velocity_robin_interface_to_full_mesh_;
}
template<class HyperelasticLawT>
inline const GlobalMatrix& VariationalFormulation<HyperelasticLawT>
::GetMonolithic2RobinVelocity() const noexcept
{
assert(!(!monolithic_2_robin_velocity_));
return *monolithic_2_robin_velocity_;
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>
::GetNonCstMonolithic2RobinVelocity() noexcept
{
return const_cast<GlobalMatrix&>(GetMonolithic2RobinVelocity());
}
} // namespace InnerLoopNS
......
......@@ -69,6 +69,10 @@ namespace HappyHeart
solid_differential_velocity_ = std::make_unique<GlobalVector>(parent::GetSystemSolution(numbering_subset));
non_homogeneous_robin_vector_ = std::make_unique<GlobalVector>(parent::GetSystemSolution(numbering_subset));
solid_differential_velocity_on_robin_interface_ =
std::make_unique<GlobalVector>(god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity_on_robin_interface)));
AllocateGlobalVector(god_of_dof, GetNonCstSolidDifferentialVelocityOnRobinInterface());
decltype(auto) system_matrix = parent::GetSystemMatrix(numbering_subset, numbering_subset);
mass_matrix_ = std::make_unique<GlobalMatrix>(system_matrix);
......@@ -483,6 +487,27 @@ namespace HappyHeart
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure)));
velocity_2_monolithic_->Init();
decltype(auto) robin_felt_space =
god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::fluid_robin));
velocity_to_robin_interface_ =
std::make_unique<type>(fluid_felt_space,
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)),
robin_felt_space,
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)));
velocity_to_robin_interface_->Init();
monolithic_2_robin_velocity_ =
std::make_unique<GlobalMatrix>(god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)),
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass_velocity_pressure)));
Wrappers::Petsc::MatMatMult(GetVelocityToRobinInterface().GetInterpolationMatrix(),
GetMonolithic2Velocity().GetInterpolationMatrix(),
*monolithic_2_robin_velocity_,
__FILE__, __LINE__);
}
}
......
......@@ -69,6 +69,7 @@ namespace HappyHeart
fluid = 10,
fluid_inlet_border = 11,
fluid_monolithic = 12,
fluid_robin = 13,
solid = 20,
solid_inlet_border = 21
......@@ -115,6 +116,7 @@ namespace HappyHeart
enum class NumberingSubsetIndex
{
fluid_velocity = 10,
fluid_velocity_on_robin_interface = 11,
porosity = 13,
solid_displacement_in_fluid_mesh = 14,
fluid_mass_velocity_pressure = 15,
......@@ -149,6 +151,7 @@ namespace HappyHeart
InputParameter::TimeManager,
InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)>,
InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity_on_robin_interface)>,
// InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::fluid_scalar)>,
InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::porosity)>,
InputParameter::NumberingSubset<EnumUnderlyingType(NumberingSubsetIndex::solid_velocity)>,
......@@ -192,6 +195,7 @@ namespace HappyHeart
InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::solid)>,
InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::fluid_inlet_border)>,
InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::fluid_monolithic)>,
InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::fluid_robin)>,
InputParameter::FEltSpace<EnumUnderlyingType(FEltSpaceIndex::solid_inlet_border)>,
InputParameter::Petsc<EnumUnderlyingType(SolverIndex::main_solver)>,
......
......@@ -35,6 +35,7 @@ namespace HappyHeart
//! Alias to unique pointer.
using unique_ptr = std::unique_ptr<self>;
private:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment