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

#820 Compute the solid velocity term that acts as a Dirichlet in the implicit...

#820 Compute the solid velocity term that acts as a Dirichlet in the implicit loop. It is still not exact: I need to restrict the data to the Robin interface.
parent e179ff3c
......@@ -83,12 +83,11 @@ namespace HappyHeart
variable_holder.template DevPrint<DevPhase::none>(variable_holder.GetFormerSolidDisplacement(),
"dispSr_old_at_begin");
//velSHalfVhf.Copy(variable_holder.GetSolidDisplacement(), __FILE__, __LINE__);
variable_holder.GetNonCstSolidOnFluidMesh().ComputeSolidVelocity(velSHalfVhf);
variable_holder.template DevPrint<DevPhase::none>(velSHalfVhf,
"velSHalfVhfr_at_begin");
}
......
......@@ -120,7 +120,7 @@ namespace HappyHeart
// dispSr[]=dispSrk0[]-dispDelta;
solid_displacement.Copy(GetVector<SolidIndex::displacement_k0>(),
__FILE__, __LINE__);
__FILE__, __LINE__);
Wrappers::Petsc::AXPY(-1.,
delta_displacement,
......
......@@ -51,8 +51,8 @@ namespace HappyHeart
const auto& solid_vel_numbering_subset = solid_god_of_dof.GetNumberingSubset(EnumUnderlyingType(NumberingSubsetIndex::solid_velocity));
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_);
work_velocity_on_solid_mesh_ = std::make_unique<GlobalVector>(solid_vel_numbering_subset);
AllocateGlobalVector(solid_god_of_dof, *work_velocity_on_solid_mesh_);
{
const auto& numbering_subset =
......@@ -70,10 +70,10 @@ namespace HappyHeart
half_sum_velocity_ = std::make_unique<GlobalVector>(numbering_subset);
}
auto& half_sum_velocity = GetNonCstHalfSumVelocity();
AllocateGlobalVector(fluid_god_of_dof, half_sum_velocity);
InitInterpolationMatrix(input_parameter_data);
}
......@@ -85,7 +85,7 @@ namespace HappyHeart
decltype(auto) former_disp = GetFormerDisplacementOnSolidMesh();
const auto& solid_displacement_on_solid_mesh = GetSolidDisplacementOnSolidMesh();
auto& half_sum_vel_on_solid_mesh = GetNonCstHalfSumVelocityOnSolidMesh();
auto& half_sum_vel_on_solid_mesh = GetNonCstWorkVelocityOnSolidMesh();
// \todo #530 Use rather a swap when it's sorted out!
half_sum_vel_on_solid_mesh.Copy(former_disp, __FILE__, __LINE__);
......@@ -103,7 +103,7 @@ namespace HappyHeart
// Then interpolate on the fluid mesh.
{
const auto& half_sum_vel_on_solid_mesh = GetHalfSumVelocityOnSolidMesh();
const auto& half_sum_vel_on_solid_mesh = GetWorkVelocityOnSolidMesh();
auto& half_sum_vel = GetNonCstHalfSumVelocity();
Wrappers::Petsc::MatMult(GetMatrixSolidToFluid(),
......@@ -197,6 +197,32 @@ namespace HappyHeart
}
void SolidOnFluidMesh::ComputeSolidVelocity(GlobalVector& solid_velocity)
{
auto& delta_disp_on_solid = GetNonCstWorkVelocityOnSolidMesh();
delta_disp_on_solid.Copy(GetSolidDisplacementOnSolidMesh(), __FILE__, __LINE__);
Wrappers::Petsc::AXPY(-1.,
GetFormerDisplacementOnSolidMesh(),
delta_disp_on_solid,
__FILE__, __LINE__);
const auto time_step = GetTimeManager().GetTimeStep();
assert(time_step > 0.);
delta_disp_on_solid.Scale(1. / time_step, __FILE__, __LINE__);
Wrappers::Petsc::MatMult(GetMatrixSolidToFluid(),
delta_disp_on_solid,
solid_velocity,
__FILE__, __LINE__);
solid_velocity.UpdateGhosts(__FILE__, __LINE__);
}
} // namespace Private
......
......@@ -129,6 +129,14 @@ namespace HappyHeart
///@}
/*!
* \brief Compute and return the velocity deduced from present and former solid displacement.
*
* \param[out] solid_velocity Vector defined on NumberingSubsetIndex::fluid_velocity numbering subset.
*
*/
void ComputeSolidVelocity(GlobalVector& solid_velocity);
private:
......@@ -188,14 +196,14 @@ namespace HappyHeart
* \brief Constant accessor to the an internal vector to store the half sum of former and current solid
* velocity on the solid mesh.
*/
const GlobalVector& GetHalfSumVelocityOnSolidMesh() const noexcept;
const GlobalVector& GetWorkVelocityOnSolidMesh() const noexcept;
/*!
* \brief Non constant accessor to the an internal vector to store the half sum of former and current
* solid velocity on the solid mesh.
*/
GlobalVector& GetNonCstHalfSumVelocityOnSolidMesh() noexcept;
GlobalVector& GetNonCstWorkVelocityOnSolidMesh() noexcept;
///@}
......@@ -217,11 +225,12 @@ namespace HappyHeart
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;
GlobalVector::unique_ptr work_velocity_on_solid_mesh_ = nullptr;
//! Matrix that transforms a P1 displacement on solid mesh into a P1b displacement on fluid mesh.
GlobalMatrix::unique_ptr matrix_solid_to_fluid_ = nullptr;
public:
GlobalMatrix::unique_ptr tmp_p1_to_p1b_ = nullptr;
......
......@@ -81,16 +81,16 @@ namespace HappyHeart
}
inline const GlobalVector& SolidOnFluidMesh::GetHalfSumVelocityOnSolidMesh() const noexcept
inline const GlobalVector& SolidOnFluidMesh::GetWorkVelocityOnSolidMesh() const noexcept
{
assert(!(!half_sum_velocity_on_solid_mesh_));
return *half_sum_velocity_on_solid_mesh_;
assert(!(!work_velocity_on_solid_mesh_));
return *work_velocity_on_solid_mesh_;
}
inline GlobalVector& SolidOnFluidMesh::GetNonCstHalfSumVelocityOnSolidMesh() noexcept
inline GlobalVector& SolidOnFluidMesh::GetNonCstWorkVelocityOnSolidMesh() noexcept
{
return const_cast<GlobalVector&>(GetHalfSumVelocityOnSolidMesh());
return const_cast<GlobalVector&>(GetWorkVelocityOnSolidMesh());
}
......
......@@ -269,14 +269,14 @@ namespace HappyHeart
BoundaryConditionMethod::penalization
>(numbering_subset, numbering_subset);
std::string filename = parent::GetOutputDirectory(numbering_subset) + "/solid_solve_residual_call_0_time_"
+ std::to_string(parent::GetTimeManager().NtimeModified()) + "_it_1_"
+ "it_newton_" + std::to_string(parent::GetSnes().GetSnesIteration(__FILE__, __LINE__)) + ".m";
rhs.View(parent::MpiHappyHeart(),
filename,
__FILE__, __LINE__,
PETSC_VIEWER_ASCII_MATLAB);
// std::string filename = parent::GetOutputDirectory(numbering_subset) + "/solid_solve_residual_call_0_time_"
// + std::to_string(parent::GetTimeManager().NtimeModified()) + "_it_1_"
// + "it_newton_" + std::to_string(parent::GetSnes().GetSnesIteration(__FILE__, __LINE__)) + ".m"; NOT THE SNES INDEXED IN CURRENT COMP FILES
//
// rhs.View(parent::MpiHappyHeart(),
// filename,
// __FILE__, __LINE__,
// PETSC_VIEWER_ASCII_MATLAB);
}
......
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