Skip to content

GitLab

Commit 70438341 authored by Gautier Bureau's avatar Gautier Bureau Committed by GILLES Sebastien
Browse files

#860 Changes made to be able to implement CardiacMechanics.

parent cb3f11b3
Expand all Hide whitespace changes
Inline Side-by-side
HappyHeart.xcodeproj/xcshareddata/xcschemes/EnsightOutput.xcscheme
View file @ 70438341
......@@ -42,7 +42,7 @@
</AdditionalOptions>
</TestAction>
<LaunchAction
buildConfiguration = "Debug"
buildConfiguration = "Release"
selectedDebuggerIdentifier = ""
selectedLauncherIdentifier = "Xcode.IDEFoundation.Launcher.PosixSpawn"
launchStyle = "0"
......@@ -112,7 +112,7 @@
</CommandLineArgument>
<CommandLineArgument
argument = "--input_parameters /Volumes/Data/${USER}/HappyHeart/Results/AcousticWave/2D_DirBCs/input_parameter_file.lua"
isEnabled = "YES">
isEnabled = "NO">
</CommandLineArgument>
<CommandLineArgument
argument = "--input_parameters /Volumes/Data/${USER}/HappyHeart/Results/AcousticWave/1D/input_parameter_file.lua"
......@@ -122,6 +122,10 @@
argument = "--input_parameters /Volumes/Data/${USER}/HappyHeart/Results/AcousticWave/1D_DirBCs/input_parameter_file.lua"
isEnabled = "NO">
</CommandLineArgument>
<CommandLineArgument
argument = "--input_parameters /Volumes/Data/${USER}/HappyHeart/Results/CardiacMechanics/input_parameter_file.lua"
isEnabled = "YES">
</CommandLineArgument>
</CommandLineArguments>
<AdditionalOptions>
</AdditionalOptions>
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/NonLinearSource/ReactionLaw/Instantiations/CourtemancheRamirezNattel.hpp
View file @ 70438341
......@@ -153,11 +153,11 @@ namespace HappyHeart
private:
std::array<ParameterAtQuadaturePoint<ParameterNS::Type::scalar>::unique_ptr, 30> parameter_list_;
std::array<ParameterAtQuadraturePoint<ParameterNS::Type::scalar>::unique_ptr, 30> parameter_list_;
ParameterAtQuadaturePoint<ParameterNS::Type::scalar>& GetNonCstParameter(parameter_index index) noexcept;
ParameterAtQuadraturePoint<ParameterNS::Type::scalar>& GetNonCstParameter(parameter_index index) noexcept;
const ParameterAtQuadaturePoint<ParameterNS::Type::scalar>& GetParameter(parameter_index index) const noexcept;
const ParameterAtQuadraturePoint<ParameterNS::Type::scalar>& GetParameter(parameter_index index) const noexcept;
private:
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/NonLinearSource/ReactionLaw/Instantiations/CourtemancheRamirezNattel.hxx
View file @ 70438341
......@@ -276,7 +276,7 @@ namespace HappyHeart
inline ReactionLaw<ReactionLawName::CourtemancheRamirezNattel>::ScalarParameterAtQuadPt&
ReactionLaw<ReactionLawName::CourtemancheRamirezNattel>::GetNonCstParameter(parameter_index index) noexcept
{
return const_cast<ParameterAtQuadaturePoint<ParameterNS::Type::scalar>&>(GetParameter(index));
return const_cast<ParameterAtQuadraturePoint<ParameterNS::Type::scalar>&>(GetParameter(index));
}
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/NonLinearSource/ReactionLaw/Instantiations/FitzHughNagumo.hpp
View file @ 70438341
......@@ -125,14 +125,14 @@ namespace HappyHeart
double GetC() const noexcept;
ParameterAtQuadaturePoint<ParameterNS::Type::scalar>& GetNonCstGate() noexcept;
ParameterAtQuadraturePoint<ParameterNS::Type::scalar>& GetNonCstGate() noexcept;
const ParameterAtQuadaturePoint<ParameterNS::Type::scalar>& GetGate() const noexcept;
const ParameterAtQuadraturePoint<ParameterNS::Type::scalar>& GetGate() const noexcept;
private:
//! Gate.
ParameterAtQuadaturePoint<ParameterNS::Type::scalar>::unique_ptr gate_ = nullptr;
ParameterAtQuadraturePoint<ParameterNS::Type::scalar>::unique_ptr gate_ = nullptr;
private:
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/NonLinearSource/ReactionLaw/Instantiations/FitzHughNagumo.hxx
View file @ 70438341
......@@ -36,10 +36,10 @@ namespace HappyHeart
const double initial_condition_gate = ipl::Extract<InitialConditionGate::Value>::Value(input_parameter_data);
gate_ = std::make_unique<ParameterAtQuadaturePoint<ParameterNS::Type::scalar>>("Gate",
gate_ = std::make_unique<ParameterAtQuadraturePoint<ParameterNS::Type::scalar>>("Gate",
mesh,
default_quadrature_rule_set,
initial_condtion_gate,
initial_condition_gate,
this->GetTimeManager());
}
......@@ -81,7 +81,7 @@ namespace HappyHeart
inline ReactionLaw<ReactionLawName::FitzHughNagumo>::ScalarParameterAtQuadPt&
ReactionLaw<ReactionLawName::FitzHughNagumo>::GetNonCstGate() noexcept
{
return const_cast<ParameterAtQuadaturePoint<ParameterNS::Type::scalar>&>(this->GetGate());
return const_cast<ParameterAtQuadraturePoint<ParameterNS::Type::scalar>&>(this->GetGate());
}
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/NonLinearSource/ReactionLaw/Instantiations/MitchellSchaeffer.hpp
View file @ 70438341
......@@ -137,14 +137,14 @@ namespace HappyHeart
double GetUMax() const noexcept;
ParameterAtQuadaturePoint<ParameterNS::Type::scalar>& GetNonCstGate() noexcept;
ParameterAtQuadraturePoint<ParameterNS::Type::scalar>& GetNonCstGate() noexcept;
const ParameterAtQuadaturePoint<ParameterNS::Type::scalar>& GetGate() const noexcept;
const ParameterAtQuadraturePoint<ParameterNS::Type::scalar>& GetGate() const noexcept;
private:
//! Gate.
ParameterAtQuadaturePoint<ParameterNS::Type::scalar>::unique_ptr gate_ = nullptr;
ParameterAtQuadraturePoint<ParameterNS::Type::scalar>::unique_ptr gate_ = nullptr;
private:
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/NonLinearSource/ReactionLaw/Instantiations/MitchellSchaeffer.hxx
View file @ 70438341
......@@ -120,7 +120,7 @@ namespace HappyHeart
inline ReactionLaw<ReactionLawName::MitchellSchaeffer>::ScalarParameterAtQuadPt&
ReactionLaw<ReactionLawName::MitchellSchaeffer>::GetNonCstGate() noexcept
{
return const_cast<ParameterAtQuadaturePoint<ParameterNS::Type::scalar>&>(this->GetGate());
return const_cast<ParameterAtQuadraturePoint<ParameterNS::Type::scalar>&>(this->GetGate());
}
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/SecondPiolaKirchhoffStressTensor.hxx
View file @ 70438341
......@@ -46,7 +46,7 @@ namespace HappyHeart
a_elementary_data.GetGeomEltDimension(),
solid,
time_manager),
ActiveStressPolicyT(a_elementary_data.GetGeomMeshRegionDimension(), a_elementary_data.Nnode(), time_manager, input_active_stress_policy),
ActiveStressPolicyT(a_elementary_data.GetGeomMeshRegionDimension(), a_elementary_data.Nnode(), a_elementary_data.NquadraturePoint(), time_manager, input_active_stress_policy),
parent(unknown_list, std::move(a_elementary_data)),
matrix_parent(),
vector_parent()
......@@ -137,9 +137,25 @@ namespace HappyHeart
const auto mesh_dimension = elementary_data.GetGeomMeshRegionDimension();
const auto& geom_elt = elementary_data.GetCurrentGeomElt();
if (geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000)
{
std::cout << "INDEX GEOM ELT" << std::endl;
std::cout << geom_elt.GetIndex() << std::endl;
}
int quad_index = 0;
for (const auto& infos_at_quad_pt : infos_at_quad_pt_list)
{
++quad_index;
if (geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000)
{
std::cout << std::endl;
std::cout << "QUAD NB " << quad_index << std::endl;
}
tangent_matrix.Zero();
PrepareInternalDataForQuadraturePoint(infos_at_quad_pt,
......@@ -157,6 +173,50 @@ namespace HappyHeart
}
if (geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000)
{
std::cout << std::endl;
std::cout << "Matrix result " << std::endl;
const auto& matrix = elementary_data.GetMatrixResult();
for (int i = 0; i < matrix.GetM(); i++)
{
for (int j = 0; j < matrix.GetN(); j++)
std::cout << std::scientific << std::setprecision(12) << matrix(i, j) << "\t";
std::cout << std::endl;
}
std::cout << std::endl;
std::cout << std::endl;
std::cout << "RHS result " << std::endl;
const auto& rhs = elementary_data.GetVectorResult();
for (int i = 0; i < rhs.GetM(); i++)
{
std::cout << std::scientific << std::setprecision(12) << rhs(i) << "\t";
}
std::cout << std::endl;
}
const auto& active_stress_ptr = static_cast<ActiveStressPolicyT*>(this);
auto& active_stress = *active_stress_ptr;
using dispatcher
= typename SecondPiolaKirchhoffStressTensorNS::Private::CorrectRHSWithActiveSchurComplement<ActiveStressPolicyT>;
dispatcher::Perform(active_stress,
elementary_data.GetNonCstVectorResult());
if (geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000)
{
const auto& rhs = elementary_data.GetVectorResult();
std::cout << std::endl;
std::cout << "RHS after correction " << std::endl;
for (int i = 0; i < rhs.GetM(); i++)
{
std::cout << std::scientific << std::setprecision(12) << rhs(i) << "\t";
}
std::cout << std::endl;
}
}
......@@ -440,6 +500,14 @@ namespace HappyHeart
dW,
d2W);
if (geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000 || geom_elt.GetIndex() == 10000)
{
std::cout << "dW after hyper visco active" << std::endl;
for (int i = 0; i < dW.GetLength(); i++)
std::cout << std::scientific << std::setprecision(12) << dW(i) << "\t";
std::cout << std::endl;
}
# ifndef NDEBUG
Wrappers::Seldon::ThrowIfNanInside(d2W, __FILE__, __LINE__);
# endif // NDEBUG
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/SecondPiolaKirchhoffStressTensor/ActiveStressPolicy/AnalyticalPrestressFiber.hpp
View file @ 70438341
......@@ -107,14 +107,10 @@ namespace HappyHeart
using vector_parent = Crtp::LocalVectorStorage<AnalyticalPrestressFiber<FiberIndexT>, 1ul>;
friend ::HappyHeart::GlobalVariationalOperatorNS::SecondPiolaKirchhoffStressTensorNS::ActiveStressPolicyNS::AnalyticalPrestressFiber<FiberIndexT>;
<<<<<<< dbe51396f04a6e795ab23d77796864c98ebcf01a
using ScalarParameterAtQuadPt = ParameterAtQuadraturePoint<ParameterNS::Type::scalar>;
=======
>>>>>>> #860 ParamterAtQuadPoint added.
public:
public:
/// \name Special members.
///@{
......@@ -168,7 +164,7 @@ namespace HappyHeart
void SetDoUpdateSigmaC(const bool do_update);
//! Constant accessor to sigma_c.
const ParameterAtQuadaturePoint<ParameterNS::Type::scalar>& GetSigmaC() const noexcept;
const ParameterAtQuadraturePoint<ParameterNS::Type::scalar>& GetSigmaC() const noexcept;
private:
......@@ -179,10 +175,10 @@ namespace HappyHeart
const ::HappyHeart::Private::FiberListManager<ParameterNS::Type::vector>::fiber_list_type& GetFibers() const noexcept;
//! Non constant accessor to sigma_c.
ParameterAtQuadaturePoint<ParameterNS::Type::scalar>& GetNonCstSigmaC() noexcept;
ParameterAtQuadraturePoint<ParameterNS::Type::scalar>& GetNonCstSigmaC() noexcept;
//! Give sigma_c to the local operator.
void SetSigmaC(ParameterAtQuadaturePoint<ParameterNS::Type::scalar>* sigma_c);
void SetSigmaC(ParameterAtQuadraturePoint<ParameterNS::Type::scalar>* sigma_c);
private:
......@@ -190,7 +186,7 @@ namespace HappyHeart
const ::HappyHeart::Private::FiberListManager<ParameterNS::Type::vector>::fiber_list_type& fibers_;
// Active stress in the fiber direction.
ParameterAtQuadaturePoint<ParameterNS::Type::scalar>* sigma_c_;
ParameterAtQuadraturePoint<ParameterNS::Type::scalar>* sigma_c_;
private:
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/SecondPiolaKirchhoffStressTensor/ActiveStressPolicy/AnalyticalPrestressFiber.hxx
View file @ 70438341
......@@ -160,7 +160,7 @@ namespace HappyHeart
inline typename AnalyticalPrestressFiber<FiberIndexT>::ScalarParameterAtQuadPt&
AnalyticalPrestressFiber<FiberIndexT>::GetNonCstSigmaC() noexcept
{
return const_cast<ParameterAtQuadaturePoint<ParameterNS::Type::scalar>&>(GetSigmaC());
return const_cast<ParameterAtQuadraturePoint<ParameterNS::Type::scalar>&>(GetSigmaC());
}
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/SecondPiolaKirchhoffStressTensor/ActiveStressPolicy/None.cpp
View file @ 70438341
......@@ -30,11 +30,13 @@ namespace HappyHeart
None::None(const unsigned int mesh_dimension,
const unsigned int Nnode,
const unsigned Nquad_point,
const TimeManager& time_manager,
input_active_stress_policy_type* input_active_stress_policy)
{
static_cast<void>(mesh_dimension);
static_cast<void>(Nnode);
static_cast<void>(Nquad_point);
static_cast<void>(time_manager);
static_cast<void>(input_active_stress_policy);
}
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/SecondPiolaKirchhoffStressTensor/ActiveStressPolicy/None.hpp
View file @ 70438341
......@@ -59,6 +59,7 @@ namespace HappyHeart
//! Constructor.
explicit None(const unsigned int mesh_dimension,
const unsigned int Nnode,
const unsigned Nquad_point,
const TimeManager& time_manager,
input_active_stress_policy_type* input_active_stress_policy);
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/SecondPiolaKirchhoffStressTensor/Private/PartialSpecialization.hpp
View file @ 70438341
......@@ -130,6 +130,36 @@ namespace HappyHeart
};
/*! #860
* \brief Helper struct used to call SecondPiolaKirchhoffStressTensor::ComputeWDerivates() if the
* ActiveStressPolicyT is not ActiveStressPolicyT::None.
*
*/
template<class ActiveStressPolicyT>
struct CorrectRHSWithActiveSchurComplement
{
static void Perform(ActiveStressPolicyT& active_stress,
LocalVector& rhs);
};
/*! #860
* \brief Helper struct used to call SecondPiolaKirchhoffStressTensor::ComputeWDerivates() if the
* ActiveStressPolicyT is ActiveStressPolicyT::None.
*
* \internal <b><tt>[internal]</tt></b> Using this struct allows not to define a full-fledged ActiveStressPolicyT::None policy
* in which all the methods would be asked to do nothing.
*
*/
template <>
struct CorrectRHSWithActiveSchurComplement<ActiveStressPolicyNS::None>
{
static void Perform(ActiveStressPolicyNS::None& active_stress,
LocalVector& rhs);
};
/*!
* \brief Helper struct used to call SecondPiolaKirchhoffStressTensor::ComputeWDerivates() if the
......
Sources/Operators/LocalVariationalOperatorInstances/NonlinearForm/SecondPiolaKirchhoffStressTensor/Private/PartialSpecialization.hxx
View file @ 70438341
......@@ -107,6 +107,24 @@ namespace HappyHeart
}
template <class ActiveStressPolicyT>
void CorrectRHSWithActiveSchurComplement<ActiveStressPolicyT>
::Perform(ActiveStressPolicyT& active_stress,
LocalVector& rhs)
{
active_stress.CorrectRHSWithActiveSchurComplement(rhs);
}
inline void CorrectRHSWithActiveSchurComplement<SecondPiolaKirchhoffStressTensorNS::ActiveStressPolicyNS::None>
::Perform(SecondPiolaKirchhoffStressTensorNS::ActiveStressPolicyNS::None& active_stress,
LocalVector& rhs)
{
static_cast<void>(rhs);
static_cast<void>(active_stress);
}
template <unsigned int DimensionT, class ViscoelasticityPolicyT>
void ComputeWDerivatesViscoelasticity<DimensionT, ViscoelasticityPolicyT>
::Perform(const InformationsAtQuadraturePoint& infos_at_quad_pt,
......
Sources/Operators/ParametersOperator/LocalParametersOperatorInstances/UpdateFiberDeformation.cpp
View file @ 70438341
......@@ -60,6 +60,9 @@ namespace HappyHeart
const unsigned int Nnode = ref_felt.Nnode();
//std::cout << "INDEX GEOM ELT" << std::endl;
//std::cout << geom_elt.GetIndex() << std::endl;
for (const auto& infos_at_quad_pt : infos_at_quad_pt_list)
{
const auto& quad_pt = infos_at_quad_pt.GetQuadraturePoint();
......@@ -70,6 +73,16 @@ namespace HappyHeart
double result_scalar_product = 0.;
//std::cout << "Schur" << std::endl; // le 1/2 n'est pas mis ici dans Heartlab
//for (int i = 0; i < schur.GetLength(); i++)
// std::cout << std::scientific << std::setprecision(12) << schur(i) << "\t";
//std::cout << std::endl;
//std::cout << "dY" << std::endl; // le 1/2 n'est pas mis ici dans Heartlab
//for (int i = 0; i < schur.GetLength(); i++)
// std::cout << std::scientific << std::setprecision(12) << increment_local_displacement[i] << "\t";
//std::cout << std::endl;
for (unsigned int node_index = 0; node_index < Nnode; ++node_index)
{
unsigned int dof_index = node_index;
......@@ -81,6 +94,8 @@ namespace HappyHeart
}
}
// std::cout << "result_scalar_product " << '\t' <<std::scientific << std::setprecision(12) << result_scalar_product << std::endl;
auto functor = [residual, result_scalar_product](double& fiber_deformation)
{
fiber_deformation += - (residual + result_scalar_product);
......@@ -89,6 +104,10 @@ namespace HappyHeart
GetNonCstParameter().UpdateValue(quad_pt,
geom_elt,
functor);
//std::cout << "ec_n+1 " << '\t' <<std::scientific << std::setprecision(12) << GetParameter().GetValue(quad_pt, geom_elt) << std::endl;
//std::cout << std::endl;
}
}
......
Sources/Parameters/Policy/AtQuadraturePoint.hpp
View file @ 70438341
......@@ -172,6 +172,7 @@ namespace HappyHeart
// #860
void Copy(const AtQuadraturePoint& parameter_at_quad_point);
......
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