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

#820 In progress: implementing properly the initial value of the...

#820 In progress: implementing properly the initial value of the 'mixtVariables' vector at the beginning of second implicit loop iteration. I cleant-up some variables duplication: there were especially three variables to hold the same solid displacement... Refactoring to regroup variables will indeed be a huge improvement!
parent e00630dc
......@@ -65,6 +65,9 @@ namespace HappyHeart
corrected_values,
__FILE__, __LINE__);
variable_holder.template DevPrint<DevPhase::newton_fp>(variable_holder.GetVelSHalfVhf(),
"dirichlet_just_before_newton_fp_loop");
corrected_values.UpdateGhosts(__FILE__, __LINE__);
}
......@@ -591,8 +594,6 @@ namespace HappyHeart
this->GetVariableHolder().template DevPrint<DevPhase::newton_fp>(vel_only,
"fluid_vel_before_darcy");
}
GetDarcyOperator().Assemble(std::make_tuple(std::ref(contribution_with_coeff)),
......
......@@ -74,6 +74,22 @@ namespace HappyHeart
double absolute_error = 1.e20;
double relative_error = 1.e20;
{
auto& velSHalfVhf = variable_holder.GetNonCstVelSHalfVhf();
variable_holder.template DevPrint<DevPhase::none>(variable_holder.GetSolidDisplacement(),
"dispSr_at_begin");
variable_holder.template DevPrint<DevPhase::none>(variable_holder.GetFormerSolidDisplacement(),
"dispSr_old_at_begin");
//velSHalfVhf.Copy(variable_holder.GetSolidDisplacement(), __FILE__, __LINE__);
}
while (absolute_error > 1.e-10 && relative_error > 1.e-6 && newton_fp_index < 2u) // \todo #820 Should not remain hardcoded!
......
......@@ -182,8 +182,7 @@ namespace HappyHeart
//! Enum class to tag each GlobalVector related to the solid.
enum class SolidIndex : std::size_t
{
displacement = 0,
displacement_k0,
displacement_k0 = 0,
H0,
delta_displacement,
dH_input_vector,
......
......@@ -40,7 +40,7 @@ namespace HappyHeart
iteration_index = 0u;
auto& solid_displacement = GetNonCstVector<SolidIndex::displacement>();
auto& solid_displacement = variable_holder.GetNonCstSolidDisplacement();
......@@ -67,8 +67,7 @@ namespace HappyHeart
decltype(auto) numbering_subset = solid_varf.GetNumberingSubset();
decltype(auto) sol = solid_varf.GetSystemSolution(numbering_subset);
// \todo #820 Check but probably not that useful: relevant value
GetNonCstVector<SolidIndex::displacement>().Copy(sol, __FILE__, __LINE__);
solid_displacement.Copy(sol, __FILE__, __LINE__);
variable_holder.template DevPrint<DevPhase::none>(sol,
"displacement_after_solid_solve");
......@@ -167,6 +166,10 @@ namespace HappyHeart
{
// \todo #820 Update solid displacement from what has been computed by the dedicated varf!
decltype(auto) variable_holder = this->GetNonCstVariableHolder();
variable_holder.GetNonCstFormerSolidDisplacement().Copy(variable_holder.GetSolidDisplacement(),
__FILE__, __LINE__);
auto& solid_on_fluid_mesh = GetNonCstSolidOnFluidMesh();
solid_on_fluid_mesh.Update();
}
......
......@@ -94,7 +94,6 @@ namespace HappyHeart
SolidVariationalFormulationPolicyT::GetVariationalFormulation();
decltype(auto) solution = solid_varf.GetSystemSolution(solid_varf.GetNumberingSubset());
GetNonCstVectorPtr<SolidIndex::displacement>() = std::make_unique<GlobalVector>(solution);
GetNonCstVectorPtr<SolidIndex::displacement_k0>() = std::make_unique<GlobalVector>(solution);
GetNonCstVectorPtr<SolidIndex::H0>() = std::make_unique<GlobalVector>(solution);
GetNonCstVectorPtr<SolidIndex::delta_displacement>() = std::make_unique<GlobalVector>(solution);
......
......@@ -38,10 +38,12 @@ namespace HappyHeart
SolidOnFluidMesh::SolidOnFluidMesh(const GlobalVector& solid_displacement_on_solid_mesh,
const GlobalVector& former_solid_displacement_on_solid_mesh,
const InputParameterList& input_parameter_data,
const TimeManager& time_manager)
: time_manager_(time_manager),
solid_displacement_on_solid_mesh_(solid_displacement_on_solid_mesh)
solid_displacement_on_solid_mesh_(solid_displacement_on_solid_mesh),
former_solid_displacement_on_solid_mesh_(former_solid_displacement_on_solid_mesh)
{
const auto& god_of_dof_manager = GodOfDofManager::GetInstance();
const auto& solid_god_of_dof = god_of_dof_manager.GetGodOfDof(EnumUnderlyingType(MeshIndex::solid));
......@@ -49,9 +51,6 @@ namespace HappyHeart
const auto& solid_vel_numbering_subset = solid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::solid_velocity));
former_displacement_on_solid_mesh_ = std::make_unique<GlobalVector>(solid_displacement_on_solid_mesh);
former_displacement_on_solid_mesh_->ZeroEntries(__FILE__, __LINE__);
half_sum_velocity_on_solid_mesh_ = std::make_unique<GlobalVector>(solid_vel_numbering_subset);
AllocateGlobalVector(solid_god_of_dof, *half_sum_velocity_on_solid_mesh_);
......@@ -83,7 +82,7 @@ namespace HappyHeart
{
// First update quantities on the solid mesh.
{
auto& former_disp = GetNonCstFormerDisplacementOnSolidMesh();
decltype(auto) former_disp = GetFormerDisplacementOnSolidMesh();
const auto& solid_displacement_on_solid_mesh = GetSolidDisplacementOnSolidMesh();
auto& half_sum_vel_on_solid_mesh = GetNonCstHalfSumVelocityOnSolidMesh();
......@@ -99,10 +98,6 @@ namespace HappyHeart
half_sum_vel_on_solid_mesh.Scale(0.5, __FILE__, __LINE__);
half_sum_vel_on_solid_mesh.UpdateGhosts(__FILE__, __LINE__);
// Update former displacement for next call.
// \todo #530 Use rather a swap when it's sorted out!
former_disp.Copy(solid_displacement_on_solid_mesh, __FILE__, __LINE__);
}
......
......@@ -79,6 +79,7 @@ namespace HappyHeart
* \copydetails doxygen_hide_input_parameter_data_arg
*/
explicit SolidOnFluidMesh(const GlobalVector& solid_displacement_on_solid_mesh,
const GlobalVector& former_solid_displacement_on_solid_mesh,
const InputParameterList& input_parameter_data,
const TimeManager& time_manager);
......@@ -119,7 +120,6 @@ namespace HappyHeart
//! Constant accessor to the computation of vs{n - 1/2} (half sum of current and former displacement divided by time step).
const GlobalVector& GetHalfSumVelocity() const noexcept;
/*!
* \brief Constant accessor to the matrix that transforms a P1 displacement on solid mesh into a P1b
* displacement on fluid mesh.
......@@ -151,9 +151,6 @@ namespace HappyHeart
//! Non constant accessor to the current displacement.
GlobalVector& GetNonCstCurrentDisplacement() noexcept;
//! Non constant accessor to the displacement of previous time iteration.
GlobalVector& GetNonCstFormerDisplacement() noexcept;
//! Constant accessor to the solid displacement as computed on the solid mesh.
const GlobalVector& GetSolidDisplacementOnSolidMesh() const noexcept;
......@@ -169,14 +166,14 @@ namespace HappyHeart
/*!
* \brief Constant accessor to the an internal vector to store the half sum of solid displacement prior
* to the itnerpolation.
* to the interpolation.
*/
const GlobalVector& GetHalfSumDisplacement() const noexcept;
/*!
* \brief Non constant accessor to the an internal vector to store the half sum of solid displacement
* prior to the itnerpolation.
* prior to the interpolation.
*/
GlobalVector& GetNonCstHalfSumDisplacement() noexcept;
......@@ -187,14 +184,6 @@ namespace HappyHeart
*/
const GlobalVector& GetFormerDisplacementOnSolidMesh() const noexcept;
/*!
* \brief Non constant accessor to the an internal vector to store the former solid displacement on the
* solid mesh.
*/
GlobalVector& GetNonCstFormerDisplacementOnSolidMesh() noexcept;
/*!
* \brief Constant accessor to the an internal vector to store the half sum of former and current solid
* velocity on the solid mesh.
......@@ -224,9 +213,9 @@ namespace HappyHeart
//! Solid displacement as computed on the solid mesh.
const GlobalVector& solid_displacement_on_solid_mesh_;
//! An internal vector to store the former solid displacement on the solid mesh.
GlobalVector::unique_ptr former_displacement_on_solid_mesh_ = nullptr;
//! Former solid displacement on the solid mesh.
const GlobalVector& former_solid_displacement_on_solid_mesh_;
//! An internal vector to store the half sum of former and current solid velocity on the solid mesh.
GlobalVector::unique_ptr half_sum_velocity_on_solid_mesh_ = nullptr;
......
......@@ -77,17 +77,10 @@ namespace HappyHeart
inline const GlobalVector& SolidOnFluidMesh::GetFormerDisplacementOnSolidMesh() const noexcept
{
assert(!(!former_displacement_on_solid_mesh_));
return *former_displacement_on_solid_mesh_;
return former_solid_displacement_on_solid_mesh_;
}
inline GlobalVector& SolidOnFluidMesh::GetNonCstFormerDisplacementOnSolidMesh() noexcept
{
return const_cast<GlobalVector&>(GetFormerDisplacementOnSolidMesh());
}
inline const GlobalVector& SolidOnFluidMesh::GetHalfSumVelocityOnSolidMesh() const noexcept
{
assert(!(!half_sum_velocity_on_solid_mesh_));
......
......@@ -89,12 +89,15 @@ namespace HappyHeart
solid_displacement_ = std::make_unique<GlobalVector>(numbering_subset);
::HappyHeart::AllocateGlobalVector(solid_god_of_dof, GetNonCstSolidDisplacement());
former_solid_displacement_ = std::make_unique<GlobalVector>(GetSolidDisplacement());
}
{
solid_on_fluid_mesh_ =
std::make_unique<Private::SolidOnFluidMesh>(GetSolidDisplacement(),
GetFormerSolidDisplacement(),
input_parameter_data,
time_manager);
}
......@@ -107,6 +110,13 @@ namespace HappyHeart
);
}
{
decltype(auto) numbering_subset =
fluid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::fluid_velocity));
velSHalfVhf_ = std::make_unique<GlobalVector>(numbering_subset);
AllocateGlobalVector(fluid_god_of_dof, GetNonCstVelSHalfVhf());
}
}
......@@ -116,7 +126,6 @@ namespace HappyHeart
midpoint_position_ = midpoint_position;
}
} // namespace PoromechanicsNS
......
......@@ -130,6 +130,13 @@ namespace HappyHeart
//! Non constant accessor to the computed solid displacement.
GlobalVector& GetNonCstSolidDisplacement() noexcept;
//! Constant accessor to the former solid displacement (i.e. from previous time step).
const GlobalVector& GetFormerSolidDisplacement() const noexcept;
//! Non constant accessor to the former solid displacement (i.e. from previous time step).
GlobalVector& GetNonCstFormerSolidDisplacement() noexcept;
//! Constant accessor to the object which keeps track on solid displacement and velocity on fluid mesh.
const Private::SolidOnFluidMesh& GetSolidOnFluidMesh() const noexcept;
......@@ -140,6 +147,12 @@ namespace HappyHeart
//! Constant accessor to the material parameters related to the solid.
const Solid& GetSolid() const noexcept;
//! Constant accessor to the odo #820 Rename this (Currently same as Bruno's Freefem)
const GlobalVector& GetVelSHalfVhf() const noexcept;
//! Non constant accessor to the odo #820 Rename this (Currently same as Bruno's Freefem)
GlobalVector& GetNonCstVelSHalfVhf() noexcept;
//! Debug indexes...
unsigned int implicit_loop_index = 0u;
......@@ -149,6 +162,30 @@ namespace HappyHeart
unsigned int dh_index = 0u;
private:
///! \name Vectors defined in fluid velocity numbering subset.
///@{
// \todo #820 At the moment Bruno's name is kept!
GlobalVector::unique_ptr velSHalfVhf_ = nullptr;
///@}
private:
///! \name Vectors defined in solid displacement subset.
///@{
//! Current solid displacement (dispSr in Freefem).
GlobalVector::unique_ptr solid_displacement_ = nullptr;
//! Former solid displacement (dispSOldr in Freefem).
GlobalVector::unique_ptr former_solid_displacement_ = nullptr;
///@}
private:
......@@ -165,9 +202,6 @@ namespace HappyHeart
//! Viscosity of the fluid.
ScalarParameter<>::unique_ptr fluid_viscosity_ = nullptr;
//! Computed solid displacement.
GlobalVector::unique_ptr solid_displacement_ = nullptr;
//! Object which keeps track on solid displacement and velocity on fluid mesh.
Private::SolidOnFluidMesh::unique_ptr solid_on_fluid_mesh_ = nullptr;
......
......@@ -73,6 +73,20 @@ namespace HappyHeart
}
inline const GlobalVector& VariableHolder::GetFormerSolidDisplacement() const noexcept
{
assert(!(!former_solid_displacement_));
return *former_solid_displacement_;
}
inline GlobalVector& VariableHolder::GetNonCstFormerSolidDisplacement() noexcept
{
return const_cast<GlobalVector&>(GetFormerSolidDisplacement());
}
inline const Private::SolidOnFluidMesh& VariableHolder::GetSolidOnFluidMesh() const noexcept
{
assert(!(!solid_on_fluid_mesh_));
......@@ -132,9 +146,21 @@ namespace HappyHeart
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
}
inline const GlobalVector& VariableHolder::GetVelSHalfVhf() const noexcept
{
assert(!(!velSHalfVhf_));
return *velSHalfVhf_;
}
inline GlobalVector& VariableHolder::GetNonCstVelSHalfVhf() noexcept
{
return const_cast<GlobalVector&>(GetVelSHalfVhf());
}
} // namespace PoromechanicsNS
......
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