Commit 18f1030b authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#718 FSI_EI: rewrite time update of hyperelasticity to make it closer to the...

#718 FSI_EI: rewrite time update of hyperelasticity to make it closer to the one in Elasticity. Remove one step from the policy to put it directly in the model.
parent 078df49f
......@@ -289,12 +289,12 @@ namespace HappyHeart
}
SolidVariationalFormulationPolicyT::FinalizeStepForFormulation(*solid_displacement_end_of_aitken_,
*solid_displacement_previous_time_iteration_,
output_dir);
{
auto& variational_formulation = SolidVariationalFormulationPolicyT::GetNonCstVariationalFormulation();
variational_formulation.UpdatePreviousTimeIterationData(*solid_displacement_end_of_aitken_,
*solid_displacement_previous_time_iteration_,
output_dir);
}
{
const auto& implicit_step_formulation = GetImplicitStepFluidVariationalFormulation();
......
......@@ -62,18 +62,7 @@ namespace HappyHeart
}
void Elasticity::FinalizeStepForFormulation(const GlobalVector& newest_displacement,
const GlobalVector& former_displacement,
const std::string& output_dir)
{
auto& variational_formulation = GetNonCstVariationalFormulation();
variational_formulation.UpdatePreviousTimeIterationData(newest_displacement,
former_displacement,
output_dir);
}
} //namespace Policy
......
......@@ -97,17 +97,7 @@ namespace HappyHeart
//! Call the appropriate solver to deal with the equation.
void Solve();
/*!
* \brief Action to do in the variational formulation when Model::FinalizeStep() is called
*
* \param[in] newest_displacement Displacement that has just been computed at the end of the implicit
* loop.
* \param[in] former_displacement Actually the \a newest_displacement obtained one time iteration ago.
*/
void FinalizeStepForFormulation(const GlobalVector& newest_displacement,
const GlobalVector& former_displacement,
const std::string& output_dir);
protected:
//! Underlying variational formulation.
......
......@@ -98,29 +98,6 @@ namespace HappyHeart
}
void VariationalFormulationElasticity
::AllocateMatricesAndVectors()
{
const auto& numbering_subset = GetNumberingSubset();
parent::AllocateGlobalMatrix(numbering_subset, numbering_subset);
parent::AllocateGlobalVector(numbering_subset);
const auto& system_matrix = parent::GetSystemMatrix(numbering_subset, numbering_subset);
const auto& system_rhs = parent::GetSystemRhs(numbering_subset);
mass_ = std::make_unique<GlobalMatrix>(system_matrix);
stiffness_ = std::make_unique<GlobalMatrix>(system_matrix);
matrix_displacement = std::make_unique<GlobalMatrix>(system_matrix);
matrix_velocity_ = std::make_unique<GlobalMatrix>(system_matrix);
vector_velocity_previous_time_iteration_ = std::make_unique<GlobalVector>(system_rhs);
fixed_contribution_to_rhs_ = std::make_unique<GlobalVector>(system_rhs);
}
void VariationalFormulationElasticity::ComputeMatrices()
{
const auto& numbering_subset = GetNumberingSubset();
......@@ -253,6 +230,27 @@ namespace HappyHeart
}
void VariationalFormulationElasticity
::AllocateMatricesAndVectors()
{
const auto& numbering_subset = GetNumberingSubset();
parent::AllocateGlobalMatrix(numbering_subset, numbering_subset);
parent::AllocateGlobalVector(numbering_subset);
const auto& system_matrix = parent::GetSystemMatrix(numbering_subset, numbering_subset);
const auto& system_rhs = parent::GetSystemRhs(numbering_subset);
mass_ = std::make_unique<GlobalMatrix>(system_matrix);
stiffness_ = std::make_unique<GlobalMatrix>(system_matrix);
matrix_displacement = std::make_unique<GlobalMatrix>(system_matrix);
matrix_velocity_ = std::make_unique<GlobalMatrix>(system_matrix);
vector_velocity_previous_time_iteration_ = std::make_unique<GlobalVector>(system_rhs);
fixed_contribution_to_rhs_ = std::make_unique<GlobalVector>(system_rhs);
}
} //namespace FSI_EINS
......
......@@ -104,17 +104,7 @@ namespace HappyHeart
//! Call the appropriate solver to deal with the equation.
void Solve();
/*!
* \brief Action to do in the variational formulation when Model::FinalizeStep() is called
*
* \param[in] newest_displacement Displacement that has just been computed at the end of the implicit
* loop.
* \param[in] former_displacement Actually the \a newest_displacement obtained one time iteration ago.
*/
void FinalizeStepForFormulation(const GlobalVector& newest_displacement,
const GlobalVector& former_displacement,
const std::string& output_dir);
protected:
//! Underlying variational formulation.
......
......@@ -108,22 +108,7 @@ namespace HappyHeart
}
template
<
class LawPolicyT,
HyperelasticityNS::TimeScheme TimeSchemeT
>
void Hyperelasticity<LawPolicyT, TimeSchemeT>::FinalizeStepForFormulation(const GlobalVector& newest_displacement,
const GlobalVector& former_displacement,
const std::string& output_dir)
{
// auto& variational_formulation = GetNonCstVariationalFormulation();
// variational_formulation.UpdatePreviousTimeIterationData(newest_displacement,
// former_displacement,
// output_dir);
}
} //namespace Policy
......
......@@ -142,10 +142,10 @@ namespace HappyHeart
///@}
//! Update for next time step. (not called after each dynamic iteration).
void UpdateForNextTimeStep();
void UpdatePreviousTimeIterationData(const GlobalVector& newest_displacement,
const GlobalVector& former_displacement,
const std::string& output_dir);
/*!
......@@ -223,9 +223,7 @@ namespace HappyHeart
//! Prepare vectors and matrices before assembling.
void PrepareMatricesAndVectorsBeforeAssembling();
//! Update current displacement. Already called in UpdateForNextTimeStep().
void UpdateDisplacementBetweenTimeStep();
private:
......
......@@ -115,32 +115,38 @@ namespace HappyHeart
HyperelasticityNS::TimeScheme TimeSchemeT
>
void VariationalFormulationHyperElasticity<LawPolicyT, TimeSchemeT>
::UpdateForNextTimeStep()
::UpdatePreviousTimeIterationData(const GlobalVector& newest_displacement,
const GlobalVector& former_displacement,
const std::string& output_dir)
{
auto& vm = this->GetNonCstVectorsAndMatrices();
auto& current_displacement = vm.GetNonCstCurrentDisplacement();
auto& current_velocity = vm.GetNonCstCurrentVelocity();
const auto& numbering_subset = this->GetNumberingSubset();
const auto& system_solution = this->GetSystemSolution(numbering_subset);
// First update velocities.
{
current_velocity.Scale(-1., __FILE__, __LINE__);
std::lock_guard<std::mutex> lock(this->GetMutex());
auto& diff_displ = this->GetNonCstHelperGlobalVector<0>();
{
diff_displ.Copy(system_solution, __FILE__, __LINE__);
Wrappers::Petsc::AXPY(-1., current_displacement, diff_displ, __FILE__, __LINE__);
diff_displ.Scale(2. / this->GetTimeManager().GetTimeStep(), __FILE__, __LINE__);
}
Wrappers::Petsc::AXPY(1., diff_displ, current_velocity, __FILE__, __LINE__);
}
// Update first the velocity: it needs the not-yet-updated previous displacement.
auto& velocity_prev_time_it = vm.GetNonCstCurrentVelocity();
velocity_prev_time_it.Scale(-1., __FILE__, __LINE__);
const double factor = 2. / parent::GetTimeManager().GetTimeStep();
Wrappers::Petsc::AXPY(factor,
newest_displacement,
velocity_prev_time_it,
__FILE__, __LINE__);
Wrappers::Petsc::AXPY(-factor,
former_displacement,
velocity_prev_time_it,
__FILE__, __LINE__);
// Then update displacement.
UpdateDisplacementBetweenTimeStep();
vm.GetNonCstCurrentDisplacement().Copy(newest_displacement, __FILE__, __LINE__);
newest_displacement.template Print<MpiScale::processor_wise>(this->MpiHappyHeart(),
output_dir + "solid_displacement_prev_it.hhdata",
__FILE__, __LINE__);
HyperelasticityNS::Private::TimeSchemeDependantUpdateForNextTimeStep(vm);
}
......@@ -192,39 +198,6 @@ namespace HappyHeart
}
template
<
class LawPolicyT,
HyperelasticityNS::TimeScheme TimeSchemeT
>
void VariationalFormulationHyperElasticity<LawPolicyT, TimeSchemeT>
::AllocateMatricesAndVectors()
{
const auto& numbering_subset = this->GetNumberingSubset();
parent::AllocateGlobalMatrix(numbering_subset, numbering_subset);
parent::AllocateGlobalVector(numbering_subset);
const auto& system_solution = this->GetSystemSolution(numbering_subset);
const auto& system_matrix = this->GetSystemMatrix(numbering_subset, numbering_subset);
{
std::lock_guard<std::mutex> lock(this->GetMutex());
for (auto& item_ptr : helper_global_vector_list_)
{
assert(!item_ptr);
item_ptr = std::make_unique<GlobalVector>(system_solution);
}
}
GetNonCstVectorsAndMatrices().AllocateVectorsAndMatrices(system_matrix,
system_solution);
}
template
<
class LawPolicyT,
......@@ -430,20 +403,6 @@ namespace HappyHeart
}
template
<
class LawPolicyT,
HyperelasticityNS::TimeScheme TimeSchemeT
>
inline void VariationalFormulationHyperElasticity<LawPolicyT, TimeSchemeT>
::UpdateDisplacementBetweenTimeStep()
{
auto& vm = this->GetNonCstVectorsAndMatrices();
vm.GetNonCstCurrentDisplacement().Copy(this->GetSystemSolution(this->GetNumberingSubset()),
__FILE__, __LINE__);
}
template
<
class LawPolicyT,
......@@ -627,7 +586,37 @@ namespace HappyHeart
}
template
<
class LawPolicyT,
HyperelasticityNS::TimeScheme TimeSchemeT
>
void VariationalFormulationHyperElasticity<LawPolicyT, TimeSchemeT>
::AllocateMatricesAndVectors()
{
const auto& numbering_subset = this->GetNumberingSubset();
parent::AllocateGlobalMatrix(numbering_subset, numbering_subset);
parent::AllocateGlobalVector(numbering_subset);
const auto& system_solution = this->GetSystemSolution(numbering_subset);
const auto& system_matrix = this->GetSystemMatrix(numbering_subset, numbering_subset);
{
std::lock_guard<std::mutex> lock(this->GetMutex());
for (auto& item_ptr : helper_global_vector_list_)
{
assert(!item_ptr);
item_ptr = std::make_unique<GlobalVector>(system_solution);
}
}
GetNonCstVectorsAndMatrices().AllocateVectorsAndMatrices(system_matrix,
system_solution);
}
} //namespace Policy
......
......@@ -49,7 +49,7 @@ int main(int argc, char ** argv)
mpi.Barrier();
}
#define HYPERELASTICITY_IN_FSI
//#define HYPERELASTICITY_IN_FSI
#ifdef HYPERELASTICITY_IN_FSI
......
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