Commit 8ae41dab authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#820 Two different quantities were represented by the same attribute (the...

#820 Two different quantities were represented by the same attribute (the unique_ptr was overwritten...); corrected now by renaming one of them.
parent 24b01992
......@@ -31,6 +31,7 @@ namespace HappyHeart
const ScalarParameter<>& density,
const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>::type& porosity,
const ParameterAtDof<ParameterNS::Type::vector>::type& velocity_solution,
const ParameterAtDof<ParameterNS::Type::vector>::type& velSHalfVhfr,
const ParameterAtDof<ParameterNS::Type::scalar>::type& pressure_solution,
const double internal_friction,
const TimeManager& time_manager,
......@@ -46,6 +47,7 @@ namespace HappyHeart
density,
porosity,
velocity_solution,
velSHalfVhfr,
pressure_solution,
internal_friction),
velocity_solution_(velocity_solution)
......
......@@ -89,6 +89,7 @@ namespace HappyHeart
const ScalarParameter<>& density,
const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>::type& porosity,
const ParameterAtDof<ParameterNS::Type::vector>::type& velocity_solution,
const ParameterAtDof<ParameterNS::Type::vector>::type& velSHalfVhfr,
const ParameterAtDof<ParameterNS::Type::scalar>::type& pressure_solution,
double internal_friction,
const TimeManager& time_manager,
......
......@@ -448,6 +448,7 @@ namespace HappyHeart
//! Constant accessor to the parameter that encapsulates velocity part of the current solution.
const ParameterAtDof<ParameterNS::Type::vector>::type& GetVelocitySolutionParam() const noexcept;
//! 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;
......@@ -457,6 +458,8 @@ namespace HappyHeart
//! Parameter that encapsulates velSHalfVhf.
const ParameterAtDof<ParameterNS::Type::vector>::type& GetVelSHalfVhfParam() const noexcept;
const ParameterAtDof<ParameterNS::Type::vector>::type& GetSolidDifferentialVelocityParam() const noexcept;
//! Constant accessor to the parameter that encapsulates pressure part of the current solution.
const ParameterAtDof<ParameterNS::Type::scalar>::type& GetPressureSolutionParam() const noexcept;
......@@ -910,6 +913,10 @@ namespace HappyHeart
//! Parameter that encapsulates velSHalfVhf.
ParameterAtDof<ParameterNS::Type::vector>::type::const_unique_ptr velSHalfVhf_param_ = nullptr;
//! Parameter that encapsulates solid_differential_velocity_param_ (also called velSHalfVhf in Freefem,
// but not the same as the one above...
ParameterAtDof<ParameterNS::Type::vector>::type::const_unique_ptr solid_differential_velocity_param_ = nullptr;
//! Value of the fluid mass at the end of previous time step.
const GlobalVector& previous_time_step_fluid_mass_on_solid_mesh_;
......
......@@ -506,6 +506,9 @@ namespace HappyHeart
// \todo #820 I assume here Darcy is the first thing assembled into rhs. Anyway this line will
// disappear once model is validated.
variable_holder.template DevPrint<DevPhase::newton_fp>(GetVelSHalfVhfParam().GetGlobalVector(),
"velSHalfVhf_before_darcy");
variable_holder.template DevPrint<DevPhase::newton_fp>(rhs,
"darcy");
......
......@@ -387,20 +387,28 @@ namespace HappyHeart
template<class HyperelasticLawT>
inline const ParameterAtDof<ParameterNS::Type::vector>::type& VariationalFormulation<HyperelasticLawT>
::GetDeltaVelocitySolutionParam() const noexcept
::GetVelSHalfVhfParam() const noexcept
{
assert(!(!delta_velocity_solution_param_));
return *delta_velocity_solution_param_;
assert(!(!velSHalfVhf_param_));
return *velSHalfVhf_param_;
}
template<class HyperelasticLawT>
inline const ParameterAtDof<ParameterNS::Type::vector>::type& VariationalFormulation<HyperelasticLawT>
::GetVelSHalfVhfParam() const noexcept
::GetSolidDifferentialVelocityParam() const noexcept
{
assert(!(!velSHalfVhf_param_));
return *velSHalfVhf_param_;
assert(!(!solid_differential_velocity_param_));
return *solid_differential_velocity_param_;
}
template<class HyperelasticLawT>
inline const ParameterAtDof<ParameterNS::Type::vector>::type& VariationalFormulation<HyperelasticLawT>
::GetDeltaVelocitySolutionParam() const noexcept
{
assert(!(!delta_velocity_solution_param_));
return *delta_velocity_solution_param_;
}
......
......@@ -253,6 +253,8 @@ namespace HappyHeart
delta_solid_varf_solver_ = InitSolver<EnumUnderlyingType(SolverIndex::main_solver)>(parent::MpiHappyHeart(),
input_parameter_data);
decltype(auto) variable_holder = this->GetVariableHolder();
{
......@@ -265,6 +267,13 @@ namespace HappyHeart
fluid_velocity,
GetCorrectedValues());
velSHalfVhf_param_ = std::make_unique<vector_at_dof_type>("velSHalfVhf",
mesh,
fluid_felt_space,
fluid_velocity,
variable_holder.GetVelSHalfVhf());
delta_fluid_mass_param_ =
std::make_unique<scalar_at_dof_type>("Delta fluid mass",
mesh,
......@@ -279,11 +288,11 @@ namespace HappyHeart
fluid_velocity,
GetDeltaMixtVariables());
velSHalfVhf_param_ = std::make_unique<vector_at_dof_type>("velSHalfVhf",
mesh,
monolithic_felt_space,
fluid_velocity,
GetSolidDifferentialVelocity());
solid_differential_velocity_param_ = std::make_unique<vector_at_dof_type>("solid_differential_velocity",
mesh,
monolithic_felt_space,
fluid_velocity,
GetSolidDifferentialVelocity());
pressure_solution_param_ = std::make_unique<scalar_at_dof_type>("Fluid pressure solution",
mesh,
......@@ -300,6 +309,7 @@ namespace HappyHeart
GetFluidDensity(),
porosity,
GetVelocitySolutionParam(),
GetVelSHalfVhfParam(),
GetPressureSolutionParam(),
GetInternalFriction(),
parent::GetTimeManager(),
......@@ -312,6 +322,7 @@ namespace HappyHeart
GetFluidDensity(),
porosity,
GetDeltaVelocitySolutionParam(),
GetVelSHalfVhfParam(), // \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.
GetPressureSolutionParam(),
GetInternalFriction(),
parent::GetTimeManager(),
......@@ -324,7 +335,7 @@ namespace HappyHeart
fluid_velocity,
mesh_dimension,
porosity,
GetVelSHalfVhfParam(),
GetSolidDifferentialVelocityParam(),
GetInternalFriction(),
DoComputeProcessorWiseLocal2Global::yes); // \todo #820 Check if no is ok.
......
......@@ -29,12 +29,14 @@ namespace HappyHeart
const ScalarParameter<>& density,
const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>::type& porosity,
const ParameterAtDof<ParameterNS::Type::vector>::type& velocity_solution,
const ParameterAtDof<ParameterNS::Type::vector>::type& velSHalfVhfr,
const ParameterAtDof<ParameterNS::Type::scalar>::type& pressure_solution,
const double internal_friction)
: parent(unknown_list, std::move(a_elementary_data)),
porosity_parent(porosity),
internal_friction_(internal_friction),
velocity_solution_(velocity_solution),
velSHalfVhf_(velSHalfVhfr),
pressure_solution_(pressure_solution),
time_manager_(time_manager),
density_(density)
......@@ -73,6 +75,7 @@ namespace HappyHeart
const auto internal_friction = GetInternalFriction();
decltype(auto) velocity_solution = GetVelocitySolution();
decltype(auto) velSHalfVhf = GetVelSHalfVhf();
decltype(auto) pressure_solution = GetPressureSolution();
const auto time_step = GetTimeManager().GetTimeStep();
......@@ -100,6 +103,7 @@ namespace HappyHeart
* infos_at_quad_pt.GetAbsoluteValueJacobianDeterminant();
decltype(auto) current_velocity = velocity_solution.GetValue(quad_pt, current_geom_elt);
decltype(auto) velSHalfVhf_value = velSHalfVhf.GetValue(quad_pt, current_geom_elt);
const auto current_pressure = pressure_solution.GetValue(quad_pt, current_geom_elt);
assert(current_velocity.GetSize() == Ncomponent);
......@@ -118,12 +122,12 @@ namespace HappyHeart
// Equivalent in Freefem:
// int2d(ThF)( weighting(x) * phi^(2.) * [velFr,velFz]' * poroExchange(velFtestr,velFtestz))
// - int2d(ThF)( weighting(x) *phi^(2.) * [velSHalfVhfr,velSHalfVhfz]'*poroExchange(velFtestr,velFtestz) )
if (is_full_darcy == ImplicitStepFluidNS::IsFullDarcy::yes)
current_item += factor
* square_porosity_value
* internal_friction
* current_velocity(component);
* (current_velocity(component) - velSHalfVhf_value(component));
// Equivalent in Freefem:
// int2d(ThF)( weighting(x) * rhoF * tau * phi * (velFr * velFtestr+velFz * velFtestz))
......@@ -142,8 +146,8 @@ namespace HappyHeart
* gradient_felt_phi(node_index, component);
// \todo #820 Still missing:
// - int2d(ThF)( weighting(x) *rhoF*tau* phi* (velFtr*velFtestr+velFtz*velFtestz))
// - int2d(ThF)( weighting(x) *phi^(2.) * [velSHalfVhfr,velSHalfVhfz]'*poroExchange(velFtestr,velFtestz) )
// - int2d(ThF)( weighting(x) *rhoF*tau* phi* (velFtr*velFtestr+velFtz*velFtestz)) NOT SAME AS LINE ABOBE:velFtr and not velFr
}
}
}
......
......@@ -96,6 +96,7 @@ namespace HappyHeart
const ScalarParameter<>& density,
const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>::type& porosity,
const ParameterAtDof<ParameterNS::Type::vector>::type& velocity_solution,
const ParameterAtDof<ParameterNS::Type::vector>::type& velSHalfVhf,
const ParameterAtDof<ParameterNS::Type::scalar>::type& pressure_solution,
double internal_friction);
......@@ -141,6 +142,9 @@ namespace HappyHeart
//! Current velocity (obtained by previous call to Solve()).
const ParameterAtDof<ParameterNS::Type::vector>::type& GetVelocitySolution() const noexcept;
//! \todo #820 Change name...
const ParameterAtDof<ParameterNS::Type::vector>::type& GetVelSHalfVhf() const noexcept;
//! Current pressure (obtained by previous call to Solve()).
const ParameterAtDof<ParameterNS::Type::scalar>::type& GetPressureSolution() const noexcept;
......@@ -162,6 +166,9 @@ namespace HappyHeart
//! Current velocity (obtained by previous call to Solve()).
const ParameterAtDof<ParameterNS::Type::vector>::type& velocity_solution_;
//! \todo #820 Bruno's name currently used!
const ParameterAtDof<ParameterNS::Type::vector>::type& velSHalfVhf_;
//! Current pressure (obtained by previous call to Solve()).
const ParameterAtDof<ParameterNS::Type::scalar>::type& pressure_solution_;
......
......@@ -61,6 +61,13 @@ namespace HappyHeart
}
inline const ParameterAtDof<ParameterNS::Type::vector>::type& Darcy
::GetVelSHalfVhf() const noexcept
{
return velSHalfVhf_;
}
} // namespace LocalVariationalOperatorNS
......
......@@ -100,6 +100,9 @@ namespace HappyHeart
velSHalfVhf_nullified_out_of_robin_interface,
__FILE__, __LINE__);
variable_holder.template DevPrint<DevPhase::none>(velSHalfVhf,
"velSHalfVhfr_at_begin");
variable_holder.template DevPrint<DevPhase::none>(velSHalfVhf_nullified_out_of_robin_interface,
"velSHalfVhfr_non_robin_zero_at_begin");
}
......
Markdown is supported
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