Commit 3edfb8f8 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#820 SolidDeltaResidual computes what is expected (with the 2. factor hack though).

parent 88238ac7
......@@ -94,6 +94,7 @@ namespace HappyHeart
const Unknown& vectorial_unknown,
unsigned int geom_mesh_region_dimension,
const HyperelasticLawT& hyperelastic_law,
const ParameterAtDof<ParameterNS::Type::scalar>::type& delta_fluid_mass,
QuadratureRulePerTopology::const_unique_ptr&& quadrature_rule_per_topology = nullptr);
//! Destructor.
......
......@@ -30,6 +30,7 @@ namespace HappyHeart
const Unknown& vectorial_unknown,
unsigned int geom_mesh_region_dimension,
const HyperelasticLawT& hyperelastic_law,
const ParameterAtDof<ParameterNS::Type::scalar>::type& delta_fluid_mass,
QuadratureRulePerTopology::const_unique_ptr&& quadrature_rule_per_topology)
: parent(felt_space,
vectorial_unknown,
......@@ -37,7 +38,8 @@ namespace HappyHeart
std::move(quadrature_rule_per_topology),
AllocateGradientFEltPhi::yes,
DoComputeProcessorWiseLocal2Global::yes,
hyperelastic_law)
hyperelastic_law,
delta_fluid_mass)
{
assert(vectorial_unknown.GetNature() == UnknownNS::Nature::vectorial);
}
......
......@@ -490,6 +490,9 @@ namespace HappyHeart
//! Constant accessor to the parameter that encapsulates delta velocity part of the current solution (in dH.edp) // \todo #820 Improve comment.
const ParameterAtDof<ParameterNS::Type::vector>::type& GetDeltaVelocitySolutionParam() const noexcept;
//! Constant accessor to the parameter that encapsulates delta fluidmass part of the current solution (in dH.edp) // \todo #820 Improve comment.
const ParameterAtDof<ParameterNS::Type::scalar>::type& GetDeltaFluidMassParam() const noexcept;
//! Parameter that encapsulates velSHalfVhf.
const ParameterAtDof<ParameterNS::Type::vector>::type& GetVelSHalfVhfParam() const noexcept;
......@@ -1002,12 +1005,14 @@ namespace HappyHeart
//! Parameter that encapsulates delta velocity part of the current solution #820 Improve comment
ParameterAtDof<ParameterNS::Type::vector>::type::const_unique_ptr delta_velocity_solution_param_ = nullptr;
//! Parameter that encapsulates delta fluid mass part of the current solution #820 Improve comment
ParameterAtDof<ParameterNS::Type::scalar>::type::const_unique_ptr delta_fluid_mass_param_ = nullptr;
//! Parameter that encapsulates pressure part of the current solution.
ParameterAtDof<ParameterNS::Type::scalar>::type::const_unique_ptr pressure_solution_param_ = nullptr;
//! Parameter that encapsulates velSHalfVhf.
ParameterAtDof<ParameterNS::Type::vector>::type::const_unique_ptr velSHalfVhf_param_ = nullptr;
//! Value of the fluid mass at the end of previous time step.
const GlobalVector& previous_time_step_fluid_mass_on_solid_mesh_;
......
......@@ -359,6 +359,18 @@ namespace HappyHeart
PETSC_VIEWER_ASCII_MATLAB);
}
{
auto& delta_fluid_mass = GetNonCstDeltaFluidMass();
delta_fluid_mass.ZeroEntries(__FILE__, __LINE__);
Wrappers::Petsc::MatMult(GetMonolithic2Mass().GetInterpolationMatrix(),
delta_mixt_variables,
delta_fluid_mass,
__FILE__, __LINE__);
}
{
auto& delta_residual = GetNonCstDeltaResidual();
......
......@@ -764,7 +764,8 @@ namespace HappyHeart
template<class HyperelasticLawT>
inline const GlobalVariationalOperatorNS::SolidDeltaResidual<HyperelasticLawT>& VariationalFormulation<HyperelasticLawT>
inline const GlobalVariationalOperatorNS::SolidDeltaResidual<HyperelasticLawT>&
VariationalFormulation<HyperelasticLawT>
::GetSolidDeltaResidualOperator() const noexcept
{
assert(!(!solid_delta_residual_operator_));
......@@ -787,7 +788,14 @@ namespace HappyHeart
}
template<class HyperelasticLawT>
inline const ParameterAtDof<ParameterNS::Type::scalar>::type& VariationalFormulation<HyperelasticLawT>
::GetDeltaFluidMassParam() const noexcept
{
assert(!(!delta_fluid_mass_param_));
return *delta_fluid_mass_param_;
}
} // namespace InnerLoopNS
......
......@@ -260,6 +260,14 @@ namespace HappyHeart
fluid_velocity,
GetCorrectedValues());
delta_fluid_mass_param_ =
std::make_unique<scalar_at_dof_type>("Delta fluid mass",
mesh,
fluid_felt_space,
fluid_mass,
GetDeltaFluidMass());
delta_velocity_solution_param_ = std::make_unique<vector_at_dof_type>("Fluid delta velocity solution",
mesh,
monolithic_felt_space,
......@@ -352,7 +360,8 @@ namespace HappyHeart
std::make_unique<type>(felt_space,
solid_disp,
mesh_dimension,
GetHyperelasticLaw());
GetHyperelasticLaw(),
GetDeltaFluidMassParam());
}
......
......@@ -127,7 +127,8 @@ namespace HappyHeart
*/
explicit SolidDeltaResidual(const ExtendedUnknown::vector_const_shared_ptr& unknown_list,
elementary_data_type&& elementary_data,
const HyperelasticLawT& hyperelastic_law);
const HyperelasticLawT& hyperelastic_law,
const ParameterAtDof<ParameterNS::Type::scalar>::type& delta_fluid_mass);
//! Destructor.
virtual ~SolidDeltaResidual() = default;
......@@ -260,7 +261,10 @@ namespace HappyHeart
//! Underlying hyperelasticity law.
const HyperelasticLawT& GetHyperelasticLaw() const noexcept;
//! Reference to delta_fluid_mass.
const ParameterAtDof<ParameterNS::Type::scalar>::type& GetDeltaFluidMass() const noexcept;
private:
......@@ -274,6 +278,9 @@ namespace HappyHeart
//! Underlying hyperelasticity law.
const HyperelasticLawT& hyperelastic_law_;
//! Reference to delta_fluid_mass.
const ParameterAtDof<ParameterNS::Type::scalar>::type& delta_fluid_mass_;
};
......
......@@ -36,11 +36,13 @@ namespace HappyHeart
SolidDeltaResidual<HyperelasticLawT>
::SolidDeltaResidual(const ExtendedUnknown::vector_const_shared_ptr& unknown_list,
elementary_data_type&& a_elementary_data,
const HyperelasticLawT& hyperelastic_law)
const HyperelasticLawT& hyperelastic_law,
const ParameterAtDof<ParameterNS::Type::scalar>::type& delta_fluid_mass)
: parent(unknown_list, std::move(a_elementary_data)),
matrix_parent(),
vector_parent(),
hyperelastic_law_(hyperelastic_law)
hyperelastic_law_(hyperelastic_law),
delta_fluid_mass_(delta_fluid_mass)
{
const auto& elementary_data = parent::GetElementaryData();
const unsigned int Ncomponent = elementary_data.GetGeomEltDimension();
......@@ -161,6 +163,10 @@ namespace HappyHeart
const auto& quad_pt = infos_at_quad_pt.GetQuadraturePoint();
const auto weight_meas = quad_pt.GetWeight() * infos_at_quad_pt.GetAbsoluteValueJacobianDeterminant();
decltype(auto) delta_fluid_mass = GetDeltaFluidMass();
const auto& current_geom_elt = elementary_data.GetCurrentGeomElt();
const auto factor = weight_meas * delta_fluid_mass.GetValue(quad_pt, current_geom_elt);
const auto& ref_felt = elementary_data.GetRefFElt(this->GetNthUnknown());
const int Nnode = static_cast<int>(elementary_data.Nnode());
......@@ -169,6 +175,8 @@ namespace HappyHeart
auto& vector_result = elementary_data.GetNonCstVectorResult();
for (unsigned int row_component = 0u; row_component < Ncomponent; ++row_component)
{
const auto dof_first_index = static_cast<int>(ref_felt.GetIndexFirstDofInElementaryData(row_component));
......@@ -183,7 +191,7 @@ namespace HappyHeart
for (int col = 0; col < int_Ncomponent; ++col)
value += dPhi(row_node, col) * rhs_part(col + component_first_index);
vector_result(dof_first_index + row_node) += value * weight_meas;
vector_result(dof_first_index + row_node) += value * factor;
}
}
}
......@@ -387,6 +395,14 @@ namespace HappyHeart
}
template <class HyperelasticLawT>
inline const ParameterAtDof<ParameterNS::Type::scalar>::type& SolidDeltaResidual<HyperelasticLawT>
::GetDeltaFluidMass() const noexcept
{
return delta_fluid_mass_;
}
} // namespace LocalVariationalOperatorNS
......
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