Commit 68af1e18 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1022 Make new fluid pressure name more accurate: it wasn't really that so far...

#1022 Make new fluid pressure name more accurate: it wasn't really that so far as fluid density was not integrated there, but added afterwards in operators.
parent ed868a18
......@@ -65,6 +65,7 @@ namespace HappyHeart
//! Constructor.
explicit NewFluidPressure(const HyperelasticLawT& hyperelastic_law,
const ScalarParameter<>& fluid_density,
const cauchy_green_tensor_type& cauchy_green_tensor,
const TimeManager& time_manager);
......
......@@ -35,6 +35,7 @@ namespace HappyHeart
template<class HyperelasticLawT>
NewFluidPressure<HyperelasticLawT>::NewFluidPressure(const HyperelasticLawT& a_hyperelastic_law,
const ScalarParameter<>& fluid_density,
const cauchy_green_tensor_type& cauchy_green_tensor,
const TimeManager& time_manager)
: ::HappyHeart::PoromechanicsNS::Crtp::CauchyGreenAccess<self>(cauchy_green_tensor),
......@@ -78,6 +79,7 @@ namespace HappyHeart
std::make_unique<pressure_parameter_operator_type<time_label_type>>(solid_felt_space,
solid_displacement_unknown,
*pressure_parameter_on_solid_[index],
fluid_density,
hyperelastic_law,
this->GetCauchyGreenTensor());
pressure_parameter_on_fluid_[index] =
......@@ -91,6 +93,7 @@ namespace HappyHeart
std::make_unique<pressure_parameter_operator_type<time_label_type>>(fluid_felt_space,
solid_displacement_unknown,
*pressure_parameter_on_fluid_[index],
fluid_density,
hyperelastic_law,
this->GetCauchyGreenTensor());
}
......@@ -112,6 +115,7 @@ namespace HappyHeart
std::make_unique<pressure_parameter_operator_type<time_label_type>>(solid_felt_space,
solid_displacement_unknown,
*pressure_parameter_on_solid_[index],
fluid_density,
hyperelastic_law,
this->GetCauchyGreenTensor());
......@@ -126,6 +130,7 @@ namespace HappyHeart
std::make_unique<pressure_parameter_operator_type<time_label_type>>(fluid_felt_space,
solid_displacement_unknown,
*pressure_parameter_on_fluid_[index],
fluid_density,
hyperelastic_law,
this->GetCauchyGreenTensor());
}
......
......@@ -85,7 +85,7 @@ namespace HappyHeart
decltype(auto) velocity_solution = GetVelocitySolution();
decltype(auto) velSHalfVhf = GetVelSHalfVhf();
// decltype(auto) pressure_solution = GetPressureSolution();
decltype(auto) pressure_solution = GetPressureSolution();
decltype(auto) velFtr = velFtr_;
const auto time_step = GetTimeManager().GetTimeStep();
......@@ -114,14 +114,14 @@ namespace HappyHeart
decltype(auto) new_velocity = velocity_solution.GetValue(quad_pt, geom_elt);
decltype(auto) velSHalfVhf_value = velSHalfVhf.GetValue(quad_pt, geom_elt);
// const auto new_pressure = pressure_solution.GetValue(quad_pt, geom_elt);
const auto new_pressure = pressure_solution.GetValue(quad_pt, geom_elt);
const auto density_value = density.GetValue(quad_pt, geom_elt);
// std::cout << "PRESSURE IN DARCY -> " << new_pressure << "\t"
// << pressure_at_quad_pt.GetValue(quad_pt, geom_elt) * density_value << std::endl;
const auto pressure_value = pressure_at_quad_pt.GetValue(quad_pt, geom_elt) * density_value;
const auto pressure_value = pressure_at_quad_pt.GetValue(quad_pt, geom_elt);
decltype(auto) velFtr_value = velFtr.GetValue(quad_pt, geom_elt);
assert(new_velocity.GetSize() == Ncomponent);
......
//! \file
//
//
// PressureContribToRhs.hxx
// HappyHeart
//
// Created by Sebastien Gilles on 03/02/15.
// Copyright (c) 2015 Inria. All rights reserved.
//
#ifndef HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_IMPLICIT_STEP_x_IMPLICIT_STEP_FLUID_x_NEWTON_FIXED_POINT_x_LOCAL_VARIATIONAL_OPERATOR_INSTANCES_x_PRESSURE_CONTRIB_TO_RHS_HXX_
# define HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_IMPLICIT_STEP_x_IMPLICIT_STEP_FLUID_x_NEWTON_FIXED_POINT_x_LOCAL_VARIATIONAL_OPERATOR_INSTANCES_x_PRESSURE_CONTRIB_TO_RHS_HXX_
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace ImplicitStepFluidNS
{
namespace NewtonFixedPointNS
{
namespace LocalVariationalOperatorNS
{
template<class HyperelasticLawT>
PressureContribToRhs<HyperelasticLawT>
::PressureContribToRhs(const ExtendedUnknown::vector_const_shared_ptr& unknown_list,
elementary_data_type&& a_elementary_data,
const ScalarParameter<>& fluid_density,
const ParameterAtQuadraturePoint<::HappyHeart::ParameterNS::Type::scalar>& pressure_at_quad_pt)
: LinearLocalVariationalOperator(unknown_list, std::move(a_elementary_data)),
fluid_density_(fluid_density),
pressure_at_quad_pt_(pressure_at_quad_pt)
{
assert(unknown_list.size() == 1ul);
}
template<class HyperelasticLawT>
const std::string& PressureContribToRhs<HyperelasticLawT>::ClassName()
{
static std::string name("PressureContribToRhs");
return name;
}
template<class HyperelasticLawT>
void PressureContribToRhs<HyperelasticLawT>::ComputeEltArray()
{
auto& elementary_data = GetNonCstElementaryData();
const auto& infos_at_quad_pt_list = elementary_data.GetInformationsAtQuadraturePointList();
auto& vector_result = elementary_data.GetNonCstVectorResult();
const auto& geom_elt = elementary_data.GetCurrentGeomElt();
const auto& ref_felt = elementary_data.GetRefFElt(GetNthUnknown());
const auto Nnode = static_cast<std::size_t>(ref_felt.Nnode());
decltype(auto) fluid_density = GetFluidDensity();
decltype(auto) pressure_at_quad_pt = GetPressureAtQuadPt();
for (const auto& infos_at_quad_pt : infos_at_quad_pt_list)
{
const auto& phi = infos_at_quad_pt.GetRefFEltPhi();
const auto& quad_pt = infos_at_quad_pt.GetQuadraturePoint();
const double hyperelastic_contrib = pressure_at_quad_pt.GetValue(quad_pt, geom_elt);
const double factor = - quad_pt.GetWeight()
* infos_at_quad_pt.GetAbsoluteValueJacobianDeterminant()
* hyperelastic_contrib;
for (auto node_index = 0ul; node_index < Nnode; ++node_index)
{
int dof_index = static_cast<int>(node_index);
vector_result(dof_index) += factor * phi(dof_index);
}
}
}
template<class HyperelasticLawT>
inline const ScalarParameter<>& PressureContribToRhs<HyperelasticLawT>::GetFluidDensity() const noexcept
{
return fluid_density_;
}
template<class HyperelasticLawT>
inline const ParameterAtQuadraturePoint<::HappyHeart::ParameterNS::Type::scalar>&
PressureContribToRhs<HyperelasticLawT>::GetPressureAtQuadPt() const noexcept
{
return pressure_at_quad_pt_;
}
} // namespace LocalVariationalOperatorNS
} // namespace NewtonFixedPointNS
} // namespace ImplicitStepFluidNS
} // namespace PoromechanicsNS
} // namespace HappyHeart
#endif // HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_IMPLICIT_STEP_x_IMPLICIT_STEP_FLUID_x_NEWTON_FIXED_POINT_x_LOCAL_VARIATIONAL_OPERATOR_INSTANCES_x_PRESSURE_CONTRIB_TO_RHS_HXX_
......@@ -131,6 +131,7 @@ namespace HappyHeart
{
new_fluid_pressure_data_ =
std::make_unique<DataNS::NewFluidPressure<hyperelastic_law_type>>(GetHyperelasticLaw(),
variable_holder.GetFluidDensity(),
this->GetCauchyGreenTensor(),
parent::GetTimeManager());
......
......@@ -81,6 +81,7 @@ namespace HappyHeart
explicit UpdatePressureAtQuadPt(const ExtendedUnknown::vector_const_shared_ptr& unknown_list,
elementary_data_type&& elementary_data,
ParameterAtQuadraturePoint<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None>& pressure,
const ScalarParameter<>& fluid_density,
const PoroHyperelasticLawT& law,
const cauchy_green_tensor_type& cauchy_green_tensor);
......@@ -115,6 +116,10 @@ namespace HappyHeart
//! Non constant accessor to invariant manager.
InvariantHolder& GetNonCstInvariantHolder() noexcept;
//! Fluid density.
const ScalarParameter<>& GetFluidDensity() const noexcept;
private:
......@@ -123,6 +128,9 @@ namespace HappyHeart
//! Hyperelastic law.
const PoroHyperelasticLawT& hyperelastic_law_;
//! Fluid density.
const ScalarParameter<>& fluid_density_;
};
......
......@@ -31,11 +31,13 @@ namespace HappyHeart
::UpdatePressureAtQuadPt(const ExtendedUnknown::vector_const_shared_ptr& a_unknown_storage,
elementary_data_type&& a_elementary_data,
ParameterAtQuadraturePoint<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None>& pressure,
const ScalarParameter<>& fluid_density,
const PoroHyperelasticLawT& law,
const cauchy_green_tensor_type& cauchy_green_tensor)
: parent(a_unknown_storage, std::move(a_elementary_data), pressure),
Crtp::CauchyGreenAccess<self>(cauchy_green_tensor),
hyperelastic_law_(law)
hyperelastic_law_(law),
fluid_density_(fluid_density)
{
const auto& elementary_data = GetElementaryData();
......@@ -96,6 +98,8 @@ namespace HappyHeart
decltype(auto) difference_fluidmass = fluid_mass_data.GetDifferenceAsParam(TimeLabelT);
decltype(auto) fluid_density = GetFluidDensity();
for (const auto& infos_at_quad_pt : infos_at_quad_pt_list)
{
decltype(auto) quad_pt = infos_at_quad_pt.GetQuadraturePoint();
......@@ -107,6 +111,8 @@ namespace HappyHeart
const auto diff_value = difference_fluidmass.GetValue(quad_pt, geom_elt);
const auto fluid_density_value = fluid_density.GetValue(quad_pt, geom_elt);
if (NumericNS::IsZero(diff_value))
{
// So-called 'analytical' case in Freefem script.
......@@ -134,6 +140,7 @@ namespace HappyHeart
new_value /= diff_value;
}
new_value *= fluid_density_value;
auto functor = [new_value](double& pressure)
{
......@@ -180,7 +187,19 @@ namespace HappyHeart
return const_cast<InvariantHolder&>(GetInvariantHolder());
}
template
<
class PoroHyperelasticLawT,
DataNS::TimeLabel2 TimeLabelT
>
inline const ScalarParameter<>& UpdatePressureAtQuadPt<PoroHyperelasticLawT, TimeLabelT>
::GetFluidDensity() const noexcept
{
return fluid_density_;
}
} // namespace LocalParameterOperatorNS
......
......@@ -99,6 +99,7 @@ namespace HappyHeart
explicit UpdatePressureAtQuadPt(const FEltSpace& felt_space,
const Unknown& unknown,
ParameterAtQuadraturePoint<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None>& pressure_at_quad_pt,
const ScalarParameter<>& fluid_density,
const PoroHyperelasticLawT& hyperelastic_law,
const cauchy_green_tensor_type& cauchy_green_tensor,
const QuadratureRulePerTopology* const quadrature_rule_per_topology = nullptr);
......@@ -124,9 +125,8 @@ namespace HappyHeart
/*!
* \brief Assemble into one or several vectors.
*
* \param[in] displacement_increment Vector that includes data from the previous iteration.
*/
void Update(const GlobalVector& displacement_increment) const;
void Update() const;
//!
const ParameterAtQuadraturePoint<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None>&
......@@ -135,25 +135,6 @@ namespace HappyHeart
private:
/*!
* \brief Computes the additional arguments required by the AssembleWithVariadicArguments() method as
* a tuple.
*
* \param[in] local_operator Local operator in charge of the elementary computation. A mutable work
* variable is actually set in this call.
* \param[in] local_felt_space List of finite elements being considered; all those related to the
* same GeometricElt.
* \param[in] additional_arguments Additional arguments given to PerformElementaryCalculation().
* These arguments might need treatment before being given to ComputeEltArray: for instance if there
* is a GlobalVector that reflects a previous state, ComputeEltArray needs only the dofs that are
* relevant locally.
*/
void
SetComputeEltArrayArguments(const LocalFEltSpace& local_felt_space,
local_operator_type& local_operator,
std::tuple<const GlobalVector&>&& additional_arguments) const;
//! Pressure at quadrature point.
ParameterAtQuadraturePoint<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None>& pressure_at_quad_pt_;
......
......@@ -31,6 +31,7 @@ namespace HappyHeart
::UpdatePressureAtQuadPt(const FEltSpace& felt_space,
const Unknown& unknown,
ParameterAtQuadraturePoint<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None>& pressure_at_quad_pt,
const ScalarParameter<>& fluid_density,
const PoroHyperelasticLawT& hyperelastic_law,
const cauchy_green_tensor_type& cauchy_green_tensor,
const QuadratureRulePerTopology* const quadrature_rule_per_topology)
......@@ -39,14 +40,11 @@ namespace HappyHeart
std::move(quadrature_rule_per_topology),
AllocateGradientFEltPhi::yes,
pressure_at_quad_pt,
fluid_density,
hyperelastic_law,
cauchy_green_tensor),
pressure_at_quad_pt_(pressure_at_quad_pt)
{
felt_space.ComputeLocal2Global(Advanced::GlobalVariationalOperatorNS::DetermineExtendedUnknownList(felt_space,
unknown),
DoComputeProcessorWiseLocal2Global::yes);
}
{ }
template
......@@ -67,28 +65,11 @@ namespace HappyHeart
DataNS::TimeLabel2 TimeLabelT
>
inline void UpdatePressureAtQuadPt<PoroHyperelasticLawT, TimeLabelT>
::Update(const GlobalVector& displacement_increment) const
::Update() const
{
return parent::UpdateImpl(displacement_increment);
return parent::UpdateImpl();
}
template
<
class PoroHyperelasticLawT,
DataNS::TimeLabel2 TimeLabelT
>
inline void UpdatePressureAtQuadPt<PoroHyperelasticLawT, TimeLabelT>
::SetComputeEltArrayArguments(const LocalFEltSpace& local_felt_space,
local_operator_type& local_operator,
std::tuple<const GlobalVector&>&& additional_arguments) const
{
::HappyHeart::GlobalVariationalOperatorNS::ExtractLocalDofValues(local_felt_space,
this->GetNthUnknown(),
std::get<0>(additional_arguments),
local_operator.GetNonCstIncrementLocalDisplacement());
}
template
<
......
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