Commit 6a80f4d0 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1022 Reintroduce operator for block T21, which I removed too hastily. Doesn't...

#1022 Reintroduce operator for block T21, which I removed too hastily. Doesn't work yet: this raises issues as the block (scalar, scalar) is now a (vector, scalar) one...
parent d5079bde
......@@ -131,7 +131,6 @@ namespace HappyHeart
void Assemble(LinearAlgebraTupleT&& global_matrix_with_coeff_tuple,
const Domain& domain = Domain()) const;
};
......
......@@ -42,7 +42,7 @@ namespace HappyHeart
const QuadratureRulePerTopology* const quadrature_rule_per_topology)
: parent(felt_space,
fluid_mass,
std::move(quadrature_rule_per_topology),
quadrature_rule_per_topology,
AllocateGradientFEltPhi::yes,
DoComputeProcessorWiseLocal2Global::yes, // \todo #820 See if no ok
fluid_density,
......
......@@ -22,6 +22,8 @@
# include "ModelInstances/UnderDevelopment/Poromechanics/Crtp/Porosity.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/Crtp/CauchyGreenAccess.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/Crtp/Porosity.hpp"
namespace HappyHeart
{
......@@ -139,10 +141,7 @@ namespace HappyHeart
*/
void ComputeEltArray();
private:
//! Density of the fluid.
const ScalarParameter<>& GetFluidDensity() const noexcept;
......@@ -164,7 +163,6 @@ namespace HappyHeart
//! Hyperelastic law.
const HyperelasticLawT& hyperelastic_law_;
};
......
......@@ -100,8 +100,6 @@ namespace HappyHeart
const double density_value = density.GetValue(quad_pt, geom_elt);
const auto diff_value = difference_fluidmass.GetValue(quad_pt, geom_elt);
double hyperelastic_contrib;
......@@ -179,7 +177,6 @@ namespace HappyHeart
}
} // namespace LocalVariationalOperatorNS
......
......@@ -29,6 +29,7 @@
# include "ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/TransientSource.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/ScalarDivVectorial.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStep/ImplicitStepFluid/NewtonFixedPoint/GlobalVariationalOperatorInstances/T22.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStep/ImplicitStepFluid/NewtonFixedPoint/GlobalVariationalOperatorInstances/T21.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStep/ImplicitStepFluid/GlobalVariationalOperatorInstances/Darcy.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStep/ImplicitStepFluid/NewtonFixedPoint/GlobalVariationalOperatorInstances/HybridVector.hpp"
......@@ -151,6 +152,9 @@ namespace HappyHeart
//! Friendship to SnesInterface.
friend struct SolverNS::SnesInterface<self>;
//! Alias to t21 type.
using t21_type = GlobalVariationalOperatorNS::T21<HyperelasticLawT>;
//! Alias to fluid mass data crtp.
using fluid_mass_data_parent =
......@@ -271,14 +275,14 @@ namespace HappyHeart
//! Constant accessor to the darcy vector after Newton FP solve.
const GlobalVector& GetDarcyVector() const noexcept;
//! Constant accessor to the matrix obtained after convergence of Newton FP
//! Constant accessor to the matrix obtained after convergence of Newton FP with block T21 updated.
const GlobalMatrix& GetMonolithicMatrixAfterNewtonFP() const noexcept;
/*!
* \brief Constant accessor to the matrix obtained after convergence of Newton FP
* updated and without Dirichlet condition in T22 block
* \brief Constant accessor to the matrix obtained after convergence of Newton FP with block T21
* updated and without Dirichlet condition in T33 block
*/
const GlobalMatrix& GetMonolithicMatrixAfterNewtonFPWithoutDirichlet() const noexcept;
......@@ -311,14 +315,15 @@ namespace HappyHeart
private:
/*!
* \brief Non constant accessor to the matrix obtained after convergence of Newton FP
* \brief Non constant accessor to the matrix obtained after convergence of Newton FP with block T21
* updated.
*/
GlobalMatrix& GetNonCstMonolithicMatrixAfterNewtonFP() noexcept;
/*!
* \brief Non constant accessor to the matrix obtained after convergence of Newton FP
* without Dirichlet condition in T22 block
* \brief Non constant accessor to the matrix obtained after convergence of Newton FP with block T21
* updated and without Dirichlet condition in T33 block
*/
GlobalMatrix& GetNonCstMonolithicMatrixAfterNewtonFPWithoutDirichlet() noexcept;
......@@ -402,6 +407,10 @@ namespace HappyHeart
//! Constant accessor to the T22 operator.
const GlobalVariationalOperatorNS::T22& GetT22Operator() const noexcept;
//! Constant accessor to the t21 operator.
const GlobalVariationalOperatorNS::T21<HyperelasticLawT>& GetT21Operator() const noexcept;
//! Constant accessor to the mass operator for fluid pressure.
const ::HappyHeart::GlobalVariationalOperatorNS::Mass& GetMassOperatorFluidPressure() const noexcept;
......@@ -475,6 +484,8 @@ namespace HappyHeart
* numbering subsets).
*/
GlobalMatrix& GetNonCstT11InMonolithic() noexcept;
GlobalMatrix& GetNonCstT21InMonolithic() noexcept;
GlobalMatrix& GetNonCstT22InMonolithic() noexcept;
......@@ -508,6 +519,8 @@ namespace HappyHeart
*/
GlobalMatrix& GetNonCstFluidMassPressureMatrixOnSolidMesh() noexcept;
//! Non constant accessor to the t21 on the non monolithic FEltSpace.
GlobalMatrix& GetNonCstT21OnFluidMesh() noexcept;
///@}
......@@ -560,6 +573,12 @@ namespace HappyHeart
//! Compute T22.
void ComputeT22(DoApplyBoundaryCondition do_apply_boundary_condition,
DoRecomputeBlock do_recompute_block);
//! Compute T21.
void ComputeT21(DoRecomputeBlock do_recompute_block,
Wrappers::Petsc::DoReuseMatrix do_reuse_matrix = Wrappers::Petsc::DoReuseMatrix::yes);
private:
......@@ -610,6 +629,9 @@ namespace HappyHeart
//! Pressure div Velocity operator.
scalar_div_vectorial_type::const_unique_ptr pressure_div_velocity_ = nullptr;
//! T21 operator.
typename t21_type::const_unique_ptr t21_operator_ = nullptr;
///@}
......@@ -644,6 +666,16 @@ namespace HappyHeart
*/
GlobalMatrix::unique_ptr t11_in_monolithic_ = nullptr;
/*!
* \brief Work matrix that gets the exact same structure as system matrix.
*
* This matrix isn't allocated directly: the MatMatMult functions are actually writing its
* pattern (an assert is nonetheless there to check it does respect the pattern delimited by the
* numbering subsets).
*/
GlobalMatrix::unique_ptr t21_in_monolithic_ = nullptr;
/*!
* \brief Work matrix that gets the exact same structure as system matrix.
*
......@@ -657,9 +689,17 @@ namespace HappyHeart
GlobalMatrix::unique_ptr t12_block_ = nullptr;
//! T21 on the non monolithic \a FEltSpace.
GlobalMatrix::unique_ptr t21_on_fluid_mesh_ = nullptr;
GlobalMatrix::unique_ptr t22_in_monolithic_ = nullptr;
//! Intermediate work matrix to transfer T21 from solid to fluid.
GlobalMatrix::unique_ptr intermediate_t21_ = nullptr;
//! Intermediate work matrix to transfer T21 from block to monolithic.
GlobalMatrix::unique_ptr intermediate_t21_monolithic_ = nullptr;
// \todo #820 TMP for debug.
GlobalMatrix::unique_ptr t22_block_ = nullptr;
......@@ -667,13 +707,13 @@ namespace HappyHeart
///@}
/*!
* \brief Matrix obtained after convergence of Newton FP
* \brief Matrix obtained after convergence of Newton FP with block T21 updated.
*/
GlobalMatrix::unique_ptr monolithic_matrix_after_newton_FP_ = nullptr;
/*!
* \brief Matrix obtained after convergence of Newton FP without
* Dirichlet condition in T22 block.
* \brief Matrix obtained after convergence of Newton FP with block T21 updated and without
* Dirichlet condition in T33 block.
*/
GlobalMatrix::unique_ptr monolithic_matrix_after_newton_FP_without_dirichlet_ = nullptr;
......
......@@ -262,6 +262,7 @@ namespace HappyHeart
// IMPORTANT: T22 must be assembled first, as there is a Dirichlet boundary condition that deals only
// with its dofs.
ComputeT22(DoApplyBoundaryCondition::yes, DoRecomputeBlock::yes);
ComputeT21(DoRecomputeBlock::yes);
ComputeT12(DoRecomputeBlock::yes);
ComputeT11(DoRecomputeBlock::yes);
......@@ -273,6 +274,7 @@ namespace HappyHeart
// with its dofs.
ComputeT22(DoApplyBoundaryCondition::yes, DoRecomputeBlock::yes);
ComputeT12(DoRecomputeBlock::yes);
ComputeT21(DoRecomputeBlock::yes);
ComputeT11(DoRecomputeBlock::no);
break;
......@@ -284,6 +286,7 @@ namespace HappyHeart
ComputeT22(DoApplyBoundaryCondition::yes, DoRecomputeBlock::no);
ComputeT12(DoRecomputeBlock::no);
ComputeT11(DoRecomputeBlock::no);
ComputeT21(DoRecomputeBlock::yes);
break;
}
......@@ -292,6 +295,7 @@ namespace HappyHeart
ComputeT22(DoApplyBoundaryCondition::no, DoRecomputeBlock::yes);
ComputeT12(DoRecomputeBlock::no);
ComputeT11(DoRecomputeBlock::no);
ComputeT21(DoRecomputeBlock::no);
break;
}
......
......@@ -137,6 +137,16 @@ namespace HappyHeart
return *t22_in_monolithic_;
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>
::GetNonCstT21InMonolithic() noexcept
{
assert(!(!t21_in_monolithic_));
return *t21_in_monolithic_;
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>
......@@ -209,14 +219,33 @@ namespace HappyHeart
return *fluid_mass_pressure_matrix_on_solid_mesh_;
}
template<class HyperelasticLawT>
inline const GlobalVariationalOperatorNS::T21<HyperelasticLawT>& VariationalFormulation<HyperelasticLawT>
::GetT21Operator() const noexcept
{
assert(!(!t21_operator_));
return *t21_operator_;
}
template<class HyperelasticLawT>
inline const HyperelasticLawT& VariationalFormulation<HyperelasticLawT>
::GetHyperelasticLaw() const noexcept
{
return hyperelastic_law_;
}
template<class HyperelasticLawT>
inline GlobalMatrix& VariationalFormulation<HyperelasticLawT>::GetNonCstT21OnFluidMesh() noexcept
{
assert(!(!t21_on_fluid_mesh_));
return *t21_on_fluid_mesh_;
}
template<class HyperelasticLawT>
inline const GlobalVariationalOperatorNS::HybridVector<HyperelasticLawT>&
......
......@@ -250,6 +250,100 @@ namespace HappyHeart
system_matrix,
__FILE__, __LINE__);
}
template<class HyperelasticLawT>
void VariationalFormulation<HyperelasticLawT>
::ComputeT21(DoRecomputeBlock do_recompute_block,
Wrappers::Petsc::DoReuseMatrix do_reuse_matrix)
{
// \todo #820 Currently two transferts are done for dev purposes; however it is possible to reduce
// that to only one!
decltype(auto) god_of_dof = parent::GetGodOfDof();
decltype(auto) numbering_subset =
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_monolithic));
auto& t21_in_monolithic = *t21_in_monolithic_;
auto& system_matrix = parent::GetNonCstSystemMatrix(numbering_subset, numbering_subset);
decltype(auto) interpolator_holder = this->GetInterpolatorHolder();
switch(do_recompute_block)
{
case DoRecomputeBlock::no:
break;
case DoRecomputeBlock::yes:
{
auto& t21_on_solid_mesh = GetNonCstFluidMassPressureMatrixOnSolidMesh();
t21_on_solid_mesh.ZeroEntries(__FILE__, __LINE__);
{
GlobalMatrixWithCoefficient matrix_with_coeff(t21_on_solid_mesh, 1.);
GetT21Operator().Assemble(std::make_tuple(std::ref(matrix_with_coeff)));
}
auto& t21_on_fluid_mesh = GetNonCstT21OnFluidMesh();
Wrappers::Petsc::MatMatMult(t21_on_solid_mesh,
interpolator_holder.GetMatrixFluidmassFluidToSolid(),
*intermediate_t21_,
__FILE__, __LINE__,
do_reuse_matrix);
Wrappers::Petsc::MatMatMult(interpolator_holder.GetMatrixFluidmassSolidToFluid(),
*intermediate_t21_,
t21_on_fluid_mesh,
__FILE__, __LINE__,
do_reuse_matrix);
# ifndef NDEBUG
AssertMatrixRespectPattern(god_of_dof, t21_on_fluid_mesh, __FILE__, __LINE__);
# endif // NDEBUG
// Now do the transfert toward the monolithic matrix.
decltype(auto) monolithic_to_mass =
interpolator_holder.GetMatrixMonolithicToFluidmassOnFluid();
auto& work_var = *intermediate_t21_monolithic_;
Wrappers::Petsc::MatMatMult(t21_on_fluid_mesh,
monolithic_to_mass,
work_var,
__FILE__, __LINE__,
do_reuse_matrix);
Wrappers::Petsc::MatTransposeMatMult(monolithic_to_mass,
work_var,
t21_in_monolithic,
__FILE__, __LINE__,
do_reuse_matrix);
this->GetVariableHolder().template DevPrint<DevPhase::newton_fp>(t21_in_monolithic,
"t21");
# ifndef NDEBUG
AssertMatrixRespectPattern(god_of_dof,
t21_in_monolithic,
__FILE__, __LINE__);
# endif // NDEBUG
}
} // switch
// And now add t21 to system matrix.
Wrappers::Petsc::AXPY<NonZeroPattern::subset>(1.,
t21_in_monolithic,
system_matrix,
__FILE__, __LINE__);
}
......
......@@ -90,16 +90,21 @@ namespace HappyHeart
t11_in_monolithic_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
t12_in_monolithic_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
t21_in_monolithic_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset);
t22_in_monolithic_ = std::make_unique<GlobalMatrix>(parent::GetSystemMatrix(numbering_subset, numbering_subset));
intermediate_t12_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset); // Numbering subset inexact but we don't care here!
intermediate_t21_ = std::make_unique<GlobalMatrix>(numbering_subset, numbering_subset); // Numbering subset inexact but we don't care here!
t12_block_ =
std::make_unique<GlobalMatrix>(god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass)),
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)));
t21_on_fluid_mesh_ =
std::make_unique<GlobalMatrix>(god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)),
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_mass)));
t22_block_ =
std::make_unique<GlobalMatrix>(god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)),
god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity)));
......@@ -208,6 +213,25 @@ namespace HappyHeart
decltype(auto) solid_displacement_data = this->GetSolidDisplacementData();
decltype(auto) new_fluid_pressure_data = this->GetNonCstNewFluidPressureData();
{
decltype(auto) god_of_dof_manager = GodOfDofManager::GetInstance();
decltype(auto) solid_god_of_dof =
god_of_dof_manager.GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));
decltype(auto) solid_felt_space =
solid_god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::solid));
decltype(auto) hyperelastic_law = GetHyperelasticLaw();
t21_operator_ =
std::make_unique<t21_type>(solid_felt_space,
fluid_mass_unknown,
variable_holder.GetFluidDensity(),
hyperelastic_law,
this->GetCauchyGreenTensor());
}
darcy_operator_ =
std::make_unique<darcy_operator_type>(monolithic_felt_space,
fluid_velocity,
......@@ -318,6 +342,7 @@ namespace HappyHeart
// properly the pattern of intermediate matrixes.
ComputeT12(DoRecomputeBlock::yes, Wrappers::Petsc::DoReuseMatrix::no);
ComputeT11(DoRecomputeBlock::yes, Wrappers::Petsc::DoReuseMatrix::no);
ComputeT21(DoRecomputeBlock::yes, Wrappers::Petsc::DoReuseMatrix::no);
}
......
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