Commit 0858d331 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1022 Darcy: remove unused arguments (and local becomes non template as a result).

parent 3c29670a
......@@ -748,6 +748,8 @@
BE3F2EDA1CD259EB00125E90 /* NonZeroPattern.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BE3F2ED81CD259E100125E90 /* NonZeroPattern.hpp */; };
BE3F2EDC1CD260FE00125E90 /* NonZeroPattern.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BE3F2EDB1CD260FE00125E90 /* NonZeroPattern.hpp */; };
BE3F2EDE1CD2616D00125E90 /* NonZeroPattern.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BE3F2EDD1CD2616D00125E90 /* NonZeroPattern.hxx */; };
BE3FF9791DDB6E6600EEB426 /* Darcy.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE3FF9781DDB6E6600EEB426 /* Darcy.cpp */; };
BE3FF97A1DDB6E6600EEB426 /* Darcy.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE3FF9781DDB6E6600EEB426 /* Darcy.cpp */; };
BE4053EB1AC172580024D5F9 /* FEltSpace.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BE4053E81AC172580024D5F9 /* FEltSpace.hpp */; };
BE4053EC1AC172580024D5F9 /* FEltSpace.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BE4053E91AC172580024D5F9 /* FEltSpace.hxx */; };
BE4053F71AC1729B0024D5F9 /* FEltSpace.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BE4053F41AC1729B0024D5F9 /* FEltSpace.cpp */; };
......@@ -5382,6 +5384,7 @@
BE3F2ED81CD259E100125E90 /* NonZeroPattern.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = NonZeroPattern.hpp; sourceTree = "<group>"; };
BE3F2EDB1CD260FE00125E90 /* NonZeroPattern.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = NonZeroPattern.hpp; sourceTree = "<group>"; };
BE3F2EDD1CD2616D00125E90 /* NonZeroPattern.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = NonZeroPattern.hxx; sourceTree = "<group>"; };
BE3FF9781DDB6E6600EEB426 /* Darcy.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = Darcy.cpp; sourceTree = "<group>"; };
BE4053E81AC172580024D5F9 /* FEltSpace.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = FEltSpace.hpp; sourceTree = "<group>"; };
BE4053E91AC172580024D5F9 /* FEltSpace.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = FEltSpace.hxx; sourceTree = "<group>"; };
BE4053F41AC1729B0024D5F9 /* FEltSpace.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = FEltSpace.cpp; sourceTree = "<group>"; };
......@@ -10226,6 +10229,7 @@
BE71B8E01DB75AE5005DCB59 /* LocalVariationalOperatorInstances */ = {
isa = PBXGroup;
children = (
BE3FF9781DDB6E6600EEB426 /* Darcy.cpp */,
BE71B8E21DB75AE5005DCB59 /* Darcy.hpp */,
BE71B8E31DB75AE5005DCB59 /* Darcy.hxx */,
);
......@@ -15470,6 +15474,7 @@
BE17E7CB1D7EA9E300A1AE6A /* Ale.cpp in Sources */,
BE17E7CC1D7EA9E300A1AE6A /* VariationalFormulation.cpp in Sources */,
BE71B90F1DB75AE5005DCB59 /* DifferentialDarcy.cpp in Sources */,
BE3FF97A1DDB6E6600EEB426 /* Darcy.cpp in Sources */,
BE17E7CE1D7EA9E300A1AE6A /* Porosity.cpp in Sources */,
BE17E7CF1D7EA9E300A1AE6A /* InterpolatorHolder.cpp in Sources */,
BE17E7D01D7EA9E300A1AE6A /* InternalFriction.cpp in Sources */,
......@@ -15823,6 +15828,7 @@
BE6FB22C1D8C3BF300B9F0FB /* InternalFriction.cpp in Sources */,
BE6FB2101D8C3BD800B9F0FB /* SolidVelocity.cpp in Sources */,
BE6FB22B1D8C3BF300B9F0FB /* InitialPorosityValue.cpp in Sources */,
BE3FF9791DDB6E6600EEB426 /* Darcy.cpp in Sources */,
);
runOnlyForDeploymentPostprocessing = 0;
};
......@@ -44,7 +44,7 @@ namespace HappyHeart
: public GlobalVariationalOperator
<
Darcy<HyperelasticLawT>,
LocalVariationalOperatorNS::Darcy<HyperelasticLawT>
LocalVariationalOperatorNS::Darcy
>
{
......@@ -60,7 +60,7 @@ namespace HappyHeart
static const std::string& ClassName();
//! Alias to local operator.
using local_operator_type = LocalVariationalOperatorNS::Darcy<HyperelasticLawT>;
using local_operator_type = LocalVariationalOperatorNS::Darcy;
//! Convenient alias to pinpoint the GlobalVariationalOperator parent.
using parent = GlobalVariationalOperator
......
//! \file
//
//
// Darcy.cpp
// HappyHeart
//
// Created by Sebastien Gilles on 03/02/15.
// Copyright (c) 2015 Inria. All rights reserved.
//
# include "HappyHeart/ModelInstances/UnderDevelopment/Poromechanics/ImplicitStep/ImplicitStepFluid/LocalVariationalOperatorInstances/Darcy.hpp"
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace LocalVariationalOperatorNS
{
Darcy::Darcy(const ExtendedUnknown::vector_const_shared_ptr& unknown_list,
typename parent::elementary_data_type&& a_elementary_data,
const TimeManager& time_manager,
const ScalarParameter<>& density,
const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>& porosity,
const ParameterAtDof<ParameterNS::Type::vector>& velocity_solution,
const ParameterAtDof<ParameterNS::Type::vector>& velSHalfVhfr,
const ParameterAtDof<ParameterNS::Type::vector>& velFtr,
const double internal_friction,
const cauchy_green_tensor_type& cauchy_green_tensor,
const ParameterAtQuadraturePoint<::HappyHeart::ParameterNS::Type::scalar>& pressure_at_quad_pt)
: parent(unknown_list, std::move(a_elementary_data)),
porosity_parent(porosity),
Crtp::CauchyGreenAccess<self>(cauchy_green_tensor),
internal_friction_(internal_friction),
velocity_solution_(velocity_solution),
velSHalfVhf_(velSHalfVhfr),
velFtr_(velFtr),
time_manager_(time_manager),
density_(density),
pressure_at_quad_pt_(pressure_at_quad_pt)
{ }
Darcy::~Darcy() = default;
const std::string& Darcy::ClassName()
{
static std::string name("DarcyWithPorosity");
return name;
}
void Darcy::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 int Ncomponent = static_cast<int>(ref_felt.Ncomponent());
const auto Nnode = static_cast<int>(ref_felt.Nnode());
decltype(auto) porosity = porosity_parent::GetPorosity();
const auto internal_friction = GetInternalFriction();
decltype(auto) velocity_solution = GetVelocitySolution();
decltype(auto) velSHalfVhf = GetVelSHalfVhf();
decltype(auto) velFtr = velFtr_;
const auto time_step = GetTimeManager().GetTimeStep();
assert(!NumericNS::IsZero(time_step));
const auto inv_time_step = 1. / time_step;
decltype(auto) density = GetDensity();
const auto is_full_darcy = GetIsFullDarcy();
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& gradient_felt_phi = infos_at_quad_pt.GetGradientFEltPhi();
const auto& quad_pt = infos_at_quad_pt.GetQuadraturePoint();
const auto porosity_value = porosity.GetValue(quad_pt, geom_elt);
const auto square_porosity_value = NumericNS::Square(porosity_value);
const double quad_pt_factor = quad_pt.GetWeight()
* infos_at_quad_pt.GetAbsoluteValueJacobianDeterminant();
decltype(auto) new_velocity = velocity_solution.GetValue(quad_pt, geom_elt);
decltype(auto) velSHalfVhf_value = velSHalfVhf.GetValue(quad_pt, geom_elt);
const auto density_value = density.GetValue(quad_pt, geom_elt);
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);
for (int node_index = 0; node_index < Nnode; ++node_index)
{
int dof_index = node_index;
const auto factor = quad_pt_factor * phi(node_index);
for (int component = 0; component < Ncomponent; ++component, dof_index += Nnode)
{
auto& new_item = vector_result(dof_index);
// Equivalent in Freefem:
// int2d(ThF)( weighting(x) * phi^(2.) * [velFr,velFz]' * poroExchange(velFtestr,velFtestz))
// - int2d(ThF)( weighting(x) *phi^(2.) * [velSHalfVhfr,velSHalfVhfz]'*poroExchange(velFtestr,velFtestz) )
if (is_full_darcy == ImplicitStepFluidNS::IsFullDarcy::yes)
new_item += factor
* square_porosity_value
* internal_friction
* (new_velocity(component) - velSHalfVhf_value(component));
// Equivalent in Freefem:
// int2d(ThF)( weighting(x) * rhoF * tau * phi * (velFr * velFtestr+velFz * velFtestz))
new_item += factor
* porosity_value
* inv_time_step
* new_velocity(component)
* density_value;
// - int2d(ThF)( weighting(x) *rhoF*tau* phi* (velFtr*velFtestr+velFtz*velFtestz))
if (is_full_darcy == ImplicitStepFluidNS::IsFullDarcy::yes) // \todo #820 Maybe to add in other cases; I'm just in FullDarcy now...
new_item -= factor
* porosity_value
* inv_time_step
* velFtr_value(component)
* density_value;
// Equivalent in Freefem:
// - int2d(ThF)( weighting(x) * pressure * divPhi(velFtestr,velFtestz,phi))
if (is_full_darcy == ImplicitStepFluidNS::IsFullDarcy::yes)
new_item -= quad_pt_factor
* pressure_value
* porosity_value
* gradient_felt_phi(node_index, component);
}
}
}
}
void Darcy::SetIsFullDarcy(ImplicitStepFluidNS::IsFullDarcy value)
{
assert(value != ImplicitStepFluidNS::IsFullDarcy::undetermined
&& "Undetermined is only used to check whether it has been uninitialized; should never be set!");
is_full_darcy_ = value;
}
} // namespace LocalVariationalOperatorNS
} // namespace PoromechanicsNS
} // namespace HappyHeart
......@@ -53,17 +53,17 @@ namespace HappyHeart
* \todo Improve the comment by writing its mathematical definition!
* \todo #559 Time dependency has been momentarily dropped.
*/
template<class HyperelasticLawT>
class Darcy final
: public ::HappyHeart::Advanced::LocalVariationalOperatorNS::LinearLocalVariationalOperator<LocalVector>,
private Crtp::Porosity<Darcy<HyperelasticLawT>>,
public Crtp::CauchyGreenAccess<Darcy<HyperelasticLawT>>
private Crtp::Porosity<Darcy>,
public Crtp::CauchyGreenAccess<Darcy>
{
public:
//! \copydoc doxygen_hide_alias_self
using self = Darcy<HyperelasticLawT>;
using self = Darcy;
//! Alias to unique pointer.
using unique_ptr = std::unique_ptr<self>;
......@@ -107,7 +107,7 @@ namespace HappyHeart
const ParameterAtQuadraturePoint<::HappyHeart::ParameterNS::Type::scalar>& pressure_at_quad_pt);
//! Destructor.
virtual ~Darcy() = default;
virtual ~Darcy();
//! Copy constructor.
Darcy(const Darcy&) = delete;
......
......@@ -24,204 +24,52 @@ namespace HappyHeart
{
template<class HyperelasticLawT>
Darcy<HyperelasticLawT>::Darcy(const ExtendedUnknown::vector_const_shared_ptr& unknown_list,
typename parent::elementary_data_type&& a_elementary_data,
const TimeManager& time_manager,
const ScalarParameter<>& density,
const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>& porosity,
const ParameterAtDof<ParameterNS::Type::vector>& velocity_solution,
const ParameterAtDof<ParameterNS::Type::vector>& velSHalfVhfr,
const ParameterAtDof<ParameterNS::Type::vector>& velFtr,
const double internal_friction,
const cauchy_green_tensor_type& cauchy_green_tensor,
const ParameterAtQuadraturePoint<::HappyHeart::ParameterNS::Type::scalar>& pressure_at_quad_pt)
: parent(unknown_list, std::move(a_elementary_data)),
porosity_parent(porosity),
Crtp::CauchyGreenAccess<self>(cauchy_green_tensor),
internal_friction_(internal_friction),
velocity_solution_(velocity_solution),
velSHalfVhf_(velSHalfVhfr),
velFtr_(velFtr),
time_manager_(time_manager),
density_(density),
pressure_at_quad_pt_(pressure_at_quad_pt)
{ }
template<class HyperelasticLawT>
const std::string& Darcy<HyperelasticLawT>::ClassName()
{
static std::string name("DarcyWithPorosity");
return name;
}
template<class HyperelasticLawT>
void Darcy<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 int Ncomponent = static_cast<int>(ref_felt.Ncomponent());
const auto Nnode = static_cast<int>(ref_felt.Nnode());
decltype(auto) porosity = porosity_parent::GetPorosity();
const auto internal_friction = GetInternalFriction();
decltype(auto) velocity_solution = GetVelocitySolution();
decltype(auto) velSHalfVhf = GetVelSHalfVhf();
decltype(auto) velFtr = velFtr_;
const auto time_step = GetTimeManager().GetTimeStep();
assert(!NumericNS::IsZero(time_step));
const auto inv_time_step = 1. / time_step;
decltype(auto) density = GetDensity();
const auto is_full_darcy = GetIsFullDarcy();
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& gradient_felt_phi = infos_at_quad_pt.GetGradientFEltPhi();
const auto& quad_pt = infos_at_quad_pt.GetQuadraturePoint();
const auto porosity_value = porosity.GetValue(quad_pt, geom_elt);
const auto square_porosity_value = NumericNS::Square(porosity_value);
const double quad_pt_factor = quad_pt.GetWeight()
* infos_at_quad_pt.GetAbsoluteValueJacobianDeterminant();
decltype(auto) new_velocity = velocity_solution.GetValue(quad_pt, geom_elt);
decltype(auto) velSHalfVhf_value = velSHalfVhf.GetValue(quad_pt, geom_elt);
const auto density_value = density.GetValue(quad_pt, geom_elt);
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);
for (int node_index = 0; node_index < Nnode; ++node_index)
{
int dof_index = node_index;
const auto factor = quad_pt_factor * phi(node_index);
for (int component = 0; component < Ncomponent; ++component, dof_index += Nnode)
{
auto& new_item = vector_result(dof_index);
// Equivalent in Freefem:
// int2d(ThF)( weighting(x) * phi^(2.) * [velFr,velFz]' * poroExchange(velFtestr,velFtestz))
// - int2d(ThF)( weighting(x) *phi^(2.) * [velSHalfVhfr,velSHalfVhfz]'*poroExchange(velFtestr,velFtestz) )
if (is_full_darcy == ImplicitStepFluidNS::IsFullDarcy::yes)
new_item += factor
* square_porosity_value
* internal_friction
* (new_velocity(component) - velSHalfVhf_value(component));
// Equivalent in Freefem:
// int2d(ThF)( weighting(x) * rhoF * tau * phi * (velFr * velFtestr+velFz * velFtestz))
new_item += factor
* porosity_value
* inv_time_step
* new_velocity(component)
* density_value;
// - int2d(ThF)( weighting(x) *rhoF*tau* phi* (velFtr*velFtestr+velFtz*velFtestz))
if (is_full_darcy == ImplicitStepFluidNS::IsFullDarcy::yes) // \todo #820 Maybe to add in other cases; I'm just in FullDarcy now...
new_item -= factor
* porosity_value
* inv_time_step
* velFtr_value(component)
* density_value;
// Equivalent in Freefem:
// - int2d(ThF)( weighting(x) * pressure * divPhi(velFtestr,velFtestz,phi))
if (is_full_darcy == ImplicitStepFluidNS::IsFullDarcy::yes)
new_item -= quad_pt_factor
* pressure_value
* porosity_value
* gradient_felt_phi(node_index, component);
}
}
}
}
template<class HyperelasticLawT>
void Darcy<HyperelasticLawT>::SetIsFullDarcy(ImplicitStepFluidNS::IsFullDarcy value)
{
assert(value != ImplicitStepFluidNS::IsFullDarcy::undetermined
&& "Undetermined is only used to check whether it has been uninitialized; should never be set!");
is_full_darcy_ = value;
}
template<class HyperelasticLawT>
inline double Darcy<HyperelasticLawT>::GetInternalFriction() const noexcept
inline double Darcy::GetInternalFriction() const noexcept
{
return internal_friction_;
}
template<class HyperelasticLawT>
inline const ParameterAtDof<ParameterNS::Type::vector>& Darcy<HyperelasticLawT>::GetVelocitySolution() const noexcept
inline const ParameterAtDof<ParameterNS::Type::vector>& Darcy::GetVelocitySolution() const noexcept
{
return velocity_solution_;
}
template<class HyperelasticLawT>
inline const TimeManager& Darcy<HyperelasticLawT>::GetTimeManager() const noexcept
inline const TimeManager& Darcy::GetTimeManager() const noexcept
{
return time_manager_;
}
template<class HyperelasticLawT>
inline const ScalarParameter<>& Darcy<HyperelasticLawT>::GetDensity() const noexcept
inline const ScalarParameter<>& Darcy::GetDensity() const noexcept
{
return density_;
}
template<class HyperelasticLawT>
inline ImplicitStepFluidNS::IsFullDarcy Darcy<HyperelasticLawT>::GetIsFullDarcy() const noexcept
inline ImplicitStepFluidNS::IsFullDarcy Darcy::GetIsFullDarcy() const noexcept
{
assert(is_full_darcy_ != ImplicitStepFluidNS::IsFullDarcy::undetermined);
return is_full_darcy_;
}
template<class HyperelasticLawT>
inline const ParameterAtDof<ParameterNS::Type::vector>& Darcy<HyperelasticLawT>
inline const ParameterAtDof<ParameterNS::Type::vector>& Darcy
::GetVelSHalfVhf() const noexcept
{
return velSHalfVhf_;
}
template<class HyperelasticLawT>
inline const ParameterAtQuadraturePoint<::HappyHeart::ParameterNS::Type::scalar>&
Darcy<HyperelasticLawT>::GetPressureAtQuadPt() const noexcept
Darcy::GetPressureAtQuadPt() const noexcept
{
return pressure_at_quad_pt_;
}
......
......@@ -23,6 +23,7 @@ poromechanics_src = Split('''
./ImplicitStep/ImplicitStepFluid/dH/GlobalVariationalOperatorInstances/HybridVector.cpp
./ImplicitStep/ImplicitStepFluid/dH/LocalVariationalOperatorInstances/DifferentialDarcy.cpp
./ImplicitStep/ImplicitStepFluid/dH/LocalVariationalOperatorInstances/HybridVector.cpp
./ImplicitStep/ImplicitStepFluid/LocalVariationalOperatorInstances/Darcy.cpp
./ImplicitStep/ImplicitStepFluid/NewtonFixedPoint/GlobalVariationalOperatorInstances/T22.cpp
./ImplicitStep/ImplicitStepFluid/NewtonFixedPoint/LocalVariationalOperatorInstances/T22.cpp
./InputParameter/BulkSolid.cpp
......
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