Commit 0fc88797 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#820 Introduce an operator SolidDeltaResidual that should compute the...

#820 Introduce an operator SolidDeltaResidual that should compute the 'deltaResidualFSIMiddlePoint' of dH.edp. At the moment it's just a trimmed down version of the SecondPiolaKirchhoffStressTensor.
parent 68fcf79e
......@@ -4487,6 +4487,10 @@
BE3C40951B0F7DD900396F68 /* ComputeGradientBasedElasticityTensor.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = ComputeGradientBasedElasticityTensor.cpp; sourceTree = "<group>"; };
BE3C40961B0F7DD900396F68 /* ComputeGradientBasedElasticityTensor.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = ComputeGradientBasedElasticityTensor.hpp; sourceTree = "<group>"; };
BE3C40971B0F7DD900396F68 /* ComputeGradientBasedElasticityTensor.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = ComputeGradientBasedElasticityTensor.hxx; sourceTree = "<group>"; };
BE3CC6CB1D3FB41800E5029B /* SolidDeltaResidual.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = SolidDeltaResidual.hpp; sourceTree = "<group>"; };
BE3CC6CC1D3FB41800E5029B /* SolidDeltaResidual.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = SolidDeltaResidual.hxx; sourceTree = "<group>"; };
BE3CC6CD1D3FB42300E5029B /* SolidDeltaResidual.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = SolidDeltaResidual.hpp; sourceTree = "<group>"; };
BE3CC6CE1D3FB42300E5029B /* SolidDeltaResidual.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = SolidDeltaResidual.hxx; sourceTree = "<group>"; };
BE3CD7BD196D7EE1007ADFFB /* SnesFunctions.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = SnesFunctions.hpp; sourceTree = "<group>"; };
BE3CD7BE196D7EE1007ADFFB /* SnesFunctions.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; lineEnding = 0; path = SnesFunctions.hxx; sourceTree = "<group>"; xcLanguageSpecificationIdentifier = xcode.lang.cpp; };
BE3E66421ACEB61100A3F7E2 /* Partition.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = Partition.cpp; sourceTree = "<group>"; };
......@@ -7951,6 +7955,8 @@
BEDAF9CB1D17F50B00D5AA0D /* LoadOnSolid.cpp */,
BEDAF9CC1D17F50B00D5AA0D /* LoadOnSolid.hpp */,
BEDAF9CD1D17F50B00D5AA0D /* LoadOnSolid.hxx */,
BE3CC6CB1D3FB41800E5029B /* SolidDeltaResidual.hpp */,
BE3CC6CC1D3FB41800E5029B /* SolidDeltaResidual.hxx */,
);
path = GlobalVariationalOperatorInstances;
sourceTree = "<group>";
......@@ -7978,6 +7984,8 @@
BEDAF9CF1D17F51600D5AA0D /* LoadOnSolid.cpp */,
BEDAF9D01D17F51600D5AA0D /* LoadOnSolid.hpp */,
BEDAF9D11D17F51600D5AA0D /* LoadOnSolid.hxx */,
BE3CC6CD1D3FB42300E5029B /* SolidDeltaResidual.hpp */,
BE3CC6CE1D3FB42300E5029B /* SolidDeltaResidual.hxx */,
);
path = LocalVariationalOperatorInstances;
sourceTree = "<group>";
//! \file
//
//
// SolidDeltaResidual.hpp
// HappyHeart
//
// Created by Sebastien Gilles on 26/06/14.
// Copyright (c) 2014 Inria. All rights reserved.
//
#ifndef HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_SOLID_DELTA_RESIDUAL_HPP_
# define HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_SOLID_DELTA_RESIDUAL_HPP_
# include <set>
# include "Core/ParameterType.hpp"
# include "Parameters/Parameter.hpp"
# include "Operators/GlobalVariationalOperator/GlobalVariationalOperator.hpp"
# include "Operators/GlobalVariationalOperator/ExtractLocalDofValues.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/LocalVariationalOperatorInstances/SolidDeltaResidual.hpp"
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace GlobalVariationalOperatorNS
{
/*!
* \brief Equivalent of deltaResidualFSIMiddlePoint in dH.edp.
*
*/
template<class HyperelasticLawT>
class SolidDeltaResidual final
: public GlobalVariationalOperator
<
SolidDeltaResidual<HyperelasticLawT>,
typename LocalVariationalOperatorNS::SolidDeltaResidual<HyperelasticLawT>
>
{
public:
//! Convenient alias.
using self = SolidDeltaResidual<HyperelasticLawT>;
//! Convenient alias.
using LocalVariationalOperator =
typename LocalVariationalOperatorNS::SolidDeltaResidual<HyperelasticLawT>;
//! Returns the name of the operator.
static const std::string& ClassName();
//! Convenient alias to pinpoint the GlobalVariationalOperator parent.
using parent = GlobalVariationalOperator
<
self,
LocalVariationalOperator
>;
//! Friendship to parent class, so this one can access private methods defined below through CRTP.
friend parent;
//! Unique ptr.
using const_unique_ptr = std::unique_ptr<const self>;
public:
/// \name Special members.
///@{
/*!
* \brief Constructor.
*
* \param[in] felt_space Finite element space upon which the operator is defined.
* \param[in] vectorial_unknown Vectorial unknown considered (should be a displacement).
* \param[in] geom_mesh_region_dimension Dimension in the geometric mesh region.
* Must be equal or higher than felt_space_dimension.
* \param[in] solid Object which provides the required material parameters for the solid.
* \param[in] time_manager Time manager need for Viscoelasticity and Active Stress.
* \param[in] input_active_stress_policy Object required only for Active Stress to compute U0 and U1 locally.
* \copydoc doxygen_hide_quadrature_rule_per_topology_nullptr_arg
*/
explicit SolidDeltaResidual(const FEltSpace& felt_space,
const Unknown& vectorial_unknown,
unsigned int geom_mesh_region_dimension,
const Solid& solid,
const TimeManager& time_manager,
const HyperelasticLawT* hyperelastic_law,
QuadratureRulePerTopology::const_unique_ptr&& quadrature_rule_per_topology = nullptr);
//! Destructor.
~SolidDeltaResidual() = default;
//! Move constructor.
SolidDeltaResidual(SolidDeltaResidual&&) = delete;
//! Copy constructor.
SolidDeltaResidual(const SolidDeltaResidual&) = delete;
//! Copy affectation.
SolidDeltaResidual& operator=(const SolidDeltaResidual&) = delete;
//! Move affectation.
SolidDeltaResidual& operator=(SolidDeltaResidual&&) = delete;
///@}
/*!
* \brief Assemble into one or several matrices and/or vectors.
*
* This overload is expected to be used when there are neither viscoelasticity nor active policy
* (hyperelasticity doesn't matter as it requires no additional arguments here).
*
* \tparam LinearAlgebraTupleT A tuple that may include \a GlobalMatrixAndCoefficient and/or
* \a GlobalVectorAndCoefficient objects. Ordering doesn't matter.
*
* \param[in] linear_algebra_tuple List of global matrices and/or vectors into which the operator is
* assembled. These objects are assumed to be already properly allocated.
* \param[in] domain Domain upon which the assembling takes place. Beware: if this domain is not a subset
* of the finite element space, assembling can only occur in a subset of the domain defined in the finite
* element space; if current \a domain is not a subset of finite element space one, assembling will occur
* upon the intersection of both.
* \param[in] displacement_previous_iteration Vector that includes data from the previous iteration. (its nature
* varies depending on the time scheme used).
*
*/
template<class LinearAlgebraTupleT>
void Assemble(LinearAlgebraTupleT&& linear_algebra_tuple,
const GlobalVector& displacement_previous_iteration,
const Domain& domain = Domain()) const;
private:
/*!
* \brief Computes the additional arguments required by the AssembleWithVariadicArguments() method as
* a tuple.
*
* These arguments are data attributes of the LocalVariationalOperator, for the sake of simplicity (a former
* scheme where the arguments appeared directly in ComputeEltArray() prototype was much harder to implement
* and understand, and induces unwanted copies).
*
* \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. Here just one vector for the displacement.
*/
void
SetComputeEltArrayArguments(const LocalFEltSpace& local_felt_space,
typename parent::local_variational_operator_type& local_operator,
std::tuple<const GlobalVector&>&& additional_arguments) const;
};
} // namespace GlobalVariationalOperatorNS
} //namespace PoromechanicsNS
} // namespace HappyHeart
# include "ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/SolidDeltaResidual.hpp"
#endif // HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_SOLID_DELTA_RESIDUAL_HPP_
//! \file
//
//
// SolidDeltaResidual.hxx
// HappyHeart
//
// Created by Sebastien Gilles on 26/06/14.
// Copyright (c) 2014 Inria. All rights reserved.
//
#ifndef HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_SOLID_DELTA_RESIDUAL_HXX_
# define HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_SOLID_DELTA_RESIDUAL_HXX_
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace GlobalVariationalOperatorNS
{
template <class HyperelasticLawT>
SolidDeltaResidual<HyperelasticLawT>
::SolidDeltaResidual(const FEltSpace& felt_space,
const Unknown& vectorial_unknown,
unsigned int geom_mesh_region_dimension,
const Solid& solid,
const TimeManager& time_manager,
const typename HyperelasticityPolicyT::law_type* hyperelastic_law,
QuadratureRulePerTopology::const_unique_ptr&& quadrature_rule_per_topology)
: parent(felt_space,
vectorial_unknown,
geom_mesh_region_dimension,
std::move(quadrature_rule_per_topology),
AllocateGradientFEltPhi::yes,
DoComputeProcessorWiseLocal2Global::yes,
solid,
time_manager,
hyperelastic_law)
{
assert(vectorial_unknown.GetNature() == UnknownNS::Nature::vectorial);
}
template <class HyperelasticLawT>
const std::string& SolidDeltaResidual<HyperelasticLawT>
::ClassName()
{
static std::string name("SolidDeltaResidual");
return name;
}
template <class HyperelasticLawT>
template<class LinearAlgebraTupleT>
inline void
SolidDeltaResidual<HyperelasticLawT>
::Assemble(LinearAlgebraTupleT&& linear_algebra_tuple,
const GlobalVector& displacement_previous_iteration,
const Domain& domain) const
{
return parent::AssembleImpl(std::move(linear_algebra_tuple), domain, displacement_previous_iteration);
}
template <class HyperelasticLawT>
void SolidDeltaResidual<HyperelasticLawT>
::SetComputeEltArrayArguments(const LocalFEltSpace& local_felt_space,
typename parent::local_variational_operator_type& local_operator,
std::tuple<const GlobalVector&>&& additional_arguments) const
{
ExtractLocalDofValues(local_felt_space,
this->GetNthUnknown(),
std::get<0>(additional_arguments),
local_operator.GetNonCstFormerLocalDisplacement());
}
} // namespace GlobalVariationalOperatorNS
} //namespace PoromechanicsNS
} // namespace HappyHeart
#endif // HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_SOLID_DELTA_RESIDUAL_HXX_
......@@ -29,6 +29,7 @@
# include "ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/TransientSource.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/Mass.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/ScalarDivVectorial.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/SolidDeltaResidual.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/GlobalVariationalOperatorInstances/T33.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/GlobalVariationalOperatorInstances/T21.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/ImplicitStepFluid/GlobalVariationalOperatorInstances/Darcy.hpp"
......@@ -748,6 +749,10 @@ namespace HappyHeart
//! Pressure div Velocity operator.
GlobalVariationalOperatorNS::ScalarDivVectorial::const_unique_ptr pressure_div_velocity_ = nullptr;
//! SolidDeltaResidual (deltaResidualFSIMiddlePoint in dH.edp).
typename GlobalVariationalOperatorNS::SolidDeltaResidual<HyperelasticLawT>::const_unique_ptr
solid_delta_residual_operator_ = nullptr;
///@}
//! Interpolator fluid mass on solid mesh -> fluid mass on fluid mesh.
......
//! \file
//
//
// SolidDeltaResidual.hpp
// HappyHeart
//
// Created by Sebastien Gilles on 02/05/14.
// Copyright (c) 2014 Inria. All rights reserved.
//
#ifndef HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_LOCAL_VARIATIONAL_OPERATOR_INSTANCES_x_SOLID_DELTA_RESIDUAL_HPP_
# define HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_LOCAL_VARIATIONAL_OPERATOR_INSTANCES_x_SOLID_DELTA_RESIDUAL_HPP_
# include <memory>
# include <vector>
# include <type_traits>
# include "Utilities/Containers/EnumClass.hpp"
# include "Utilities/LocalLinearAlgebraStorage/LocalMatrixStorage.hpp"
# include "Utilities/LocalLinearAlgebraStorage/LocalVectorStorage.hpp"
# ifndef NDEBUG
# include "ThirdParty/Wrappers/Seldon/MatrixOperations.hpp"
# endif // NDEBUG
# include "Parameters/InitParameter.hpp"
# include "Parameters/Compound/Solid/Solid.hpp"
# include "Operators/LocalVariationalOperator/ElementaryData.hpp"
# include "Operators/LocalVariationalOperator/NonlinearLocalVariationalOperator.hpp"
# include "Operators/LocalVariationalOperator/Private/ExtractGradientBasedBlock.hpp"
# include "Operators/LocalVariationalOperator/Private/DerivativeGreenLagrange.hpp"
# include "Operators/LocalVariationalOperator/Private/GradientDisplacementMatrix.hpp"
# include "Operators/LocalVariationalOperatorInstances/NonlinearForm/SecondPiolaKirchhoffStressTensor/Private/Helper.hpp"
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace LocalVariationalOperatorNS
{
/*!
* \brief Define the operator used in the hyperelastic problem.
*
* \todo Improve the comment by writing its mathematical definition!
*
* \tparam HyperelasticityPolicyT Policy that defines if an hyperelastic contribution is present in the tensor.
* \tparam ViscoelasticityPolicyT Policy that defines if an viscoelastic contribution is present in the tensor.
* \tparam ActiveStressPolicyT Policy that defines if an active stress contribution is present in the tensor.
*
*/
template<class HyperelasticLawT>
class SolidDeltaResidual final
: public LinearLocalVariationalOperator<LocalVector>,
public ::HappyHeart::Crtp::LocalMatrixStorage<SolidDeltaResidual<HyperelasticLawT>, 9ul>,
public ::HappyHeart::Crtp::LocalVectorStorage<SolidDeltaResidual<HyperelasticLawT>, 2ul>
{
public:
//! Alias to self.
using self = SolidDeltaResidual<HyperelasticLawT>;
//! Alias to unique pointer.
using unique_ptr = std::unique_ptr<self>;
//! Returns the name of the operator.
static const std::string& ClassName();
//! Type of the elementary matrix.
using matrix_type = LocalMatrix;
//! Type of the elementary vector.
using vector_type = LocalVector;
//! Alias to variational operator parent.
using parent = NonlinearLocalVariationalOperator
<
LocalMatrix,
LocalVector
>;
//! Alias to the parent that provides LocalMatrixStorage.
using matrix_parent = ::HappyHeart::Crtp::LocalMatrixStorage<self, 9ul>;
//! Alias to the parent that provides LocalVectorStorage.
using vector_parent = ::HappyHeart::Crtp::LocalVectorStorage<self, 2ul>;
//! Rejuvenate alias from parent class.
using elementary_data_type = typename parent::elementary_data_type;
//! Alias to Green Lagrange.
using deriv_green_lagrange_type =
::HappyHeart::LocalVariationalOperatorNS::Private::DerivativeGreenLagrange<::HappyHeart::LocalVariationalOperatorNS::Private::GreenLagrangeOrEta::green_lagrange>;
public:
/// \name Special members.
///@{
/*!
* \brief Constructor.
*
* \param[in] unknown_list List of unknowns considered by the operators. Its type (vector_shared_ptr)
* is due to constraints from genericity; for current operator it is expected to hold exactly
* one vectorial unknown (as long as #7 is not implemented at the very least...)
* \param[in] elementary_data Elementary matrices and vectors that will perform the calculations.
* \param[in] solid Object which provides the required material parameters for the solid.
* \param[in] hyperelastic_law The hyperelastic law if there is an hyperelastic part in the operator,
* or nullptr otherwise. Example of hyperelasticlaw is HappyHeart::HyperelasticLawNS::CiarletGeymonat.
*
* \internal <b><tt>[internal]</tt></b> This constructor must not be called manually: it is involved only in
* GlobalVariationalOperator<DerivedT, LocalVariationalOperatorT>::CreateLocalOperatorList() method.
*/
explicit SolidDeltaResidual(const ExtendedUnknown::vector_const_shared_ptr& unknown_list,
elementary_data_type&& elementary_data,
const Solid& solid,
const TimeManager& time_manager,
const HyperelasticLawT* hyperelastic_law);
//! Destructor.
virtual ~SolidDeltaResidual();
//! Copy constructor.
SolidDeltaResidual(const SolidDeltaResidual&) = delete;
//! Move constructor.
SolidDeltaResidual(SolidDeltaResidual&&) = delete;
//! Copy affectation.
SolidDeltaResidual& operator=(const SolidDeltaResidual&) = delete;
//! Move affectation.
SolidDeltaResidual& operator=(SolidDeltaResidual&&) = delete;
///@}
/*!
* \brief Compute the elementary matrix and/or vector.
*/
void ComputeEltArray();
//! Constant accessor to the former local displacement required by ComputeEltArray().
const std::vector<double>& GetFormerLocalDisplacement() const noexcept;
//! Non constant accessor to the former local displacement required by ComputeEltArray().
std::vector<double>& GetNonCstFormerLocalDisplacement() noexcept;
private:
//! Accessor to the displacement unknown.
Unknown::const_shared_ptr GetVectorialUnknownPtr() const noexcept;
using InformationsAtQuadraturePoint = ::HappyHeart::Private::InformationsAtQuadraturePoint;
/*!
* \brief Part of ComputeEltArray once all the dimension depending operations have been performed.
*/
void ComputeAtQuadraturePoint(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const LocalMatrix& tangent_matrix,
const LocalVector& rhs_part,
elementary_data_type& elementary_data,
OperatorNS::Nature what_must_be_built) const;
//! Compute the derivates of W.
template<unsigned int DimensionT>
void ComputeWDerivates(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& ref_felt,
const LocalMatrix& gradient_displacement,
const CauchyGreenTensor& cauchy_green_tensor,
const LocalMatrix& De,
const LocalMatrix& transposed_De,
LocalVector& dW,
LocalMatrix& d2W);
/*!
* \brief Compute internal data such as invariants, De, Cauchy-Green tensor for current quadrature point.
*/
void PrepareInternalDataForQuadraturePoint(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& ref_felt,
const std::vector<double>& local_displacement,
OperatorNS::Nature what_must_be_built,
const unsigned int mesh_dimension);
//! Helper class to manage computation of Cauchy-Green tensor.
CauchyGreenTensor& GetNonCstCauchyGreenTensor() noexcept;
//! Helper class to manage computation of derivative of Green-Lagrange.
deriv_green_lagrange_type& GetNonCstDerivativeGreenLagrange() noexcept;
/*!
* \brief Manage the ComputeEltArray part of the calculation that is dependant on the dimension considered.
*/
template<unsigned int DimensionT>
void PrepareInternalDataForQuadraturePointForDimension(const InformationsAtQuadraturePoint& infos_at_quad_pt,
const GeometricElt& geom_elt,
const Advanced::RefFEltInLocalOperator& ref_felt,
const std::vector<double>& local_displacement,
OperatorNS::Nature what_must_be_built);
private:
//! Helper class to manage computation of derivative of Green-Lagrange.
typename deriv_green_lagrange_type::unique_ptr deriv_green_lagrange_ = nullptr;
//! Helper class to manage computation of Cauchy-Green tensor.
CauchyGreenTensor::unique_ptr cauchy_green_tensor_ = nullptr;
/*!
* \brief Displacement obtained at previous time iteration expressed at the local level.
*
* \internal <b><tt>[internal]</tt></b> This is a work variable that should be used only within ComputeEltArray.
*/
std::vector<double> former_local_displacement_;
/// \name Useful indexes to fetch the work matrices and vectors.
///@{