Commit afc60bd1 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#823 Poromechanics: write version of scalar div vectorial with porosity.

parent 3211e057
......@@ -1250,6 +1250,8 @@
BEABCEE21AE66F5F00817D37 /* libFormulationSolver.a in Frameworks */ = {isa = PBXBuildFile; fileRef = BEABCEBD1AE65A9A00817D37 /* libFormulationSolver.a */; };
BEADAD2A1B2ED2C3001D9111 /* AtQuadraturePoint.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BEADAD271B2ED2C3001D9111 /* AtQuadraturePoint.hpp */; };
BEADAD2B1B2ED2C3001D9111 /* AtQuadraturePoint.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BEADAD281B2ED2C3001D9111 /* AtQuadraturePoint.hxx */; };
BEAE0D201CB2AF0B00CE333D /* ScalarDivVectorial.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BEAE0D1D1CB2AF0B00CE333D /* ScalarDivVectorial.cpp */; };
BEAE0D251CB2AF1600CE333D /* ScalarDivVectorial.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BEAE0D221CB2AF1600CE333D /* ScalarDivVectorial.cpp */; };
BEAE12761C7CA858004FB5A7 /* GradOnGradientBasedElasticityTensor.cpp in Sources */ = {isa = PBXBuildFile; fileRef = BEAE12731C7CA858004FB5A7 /* GradOnGradientBasedElasticityTensor.cpp */; };
BEAE12771C7CA858004FB5A7 /* GradOnGradientBasedElasticityTensor.hpp in Headers */ = {isa = PBXBuildFile; fileRef = BEAE12741C7CA858004FB5A7 /* GradOnGradientBasedElasticityTensor.hpp */; };
BEAE12781C7CA858004FB5A7 /* GradOnGradientBasedElasticityTensor.hxx in Headers */ = {isa = PBXBuildFile; fileRef = BEAE12751C7CA858004FB5A7 /* GradOnGradientBasedElasticityTensor.hxx */; };
......@@ -5058,6 +5060,12 @@
BEABCEBD1AE65A9A00817D37 /* libFormulationSolver.a */ = {isa = PBXFileReference; explicitFileType = archive.ar; includeInIndex = 0; path = libFormulationSolver.a; sourceTree = BUILT_PRODUCTS_DIR; };
BEADAD271B2ED2C3001D9111 /* AtQuadraturePoint.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = AtQuadraturePoint.hpp; sourceTree = "<group>"; };
BEADAD281B2ED2C3001D9111 /* AtQuadraturePoint.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = AtQuadraturePoint.hxx; sourceTree = "<group>"; };
BEAE0D1D1CB2AF0B00CE333D /* ScalarDivVectorial.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = ScalarDivVectorial.cpp; sourceTree = "<group>"; };
BEAE0D1E1CB2AF0B00CE333D /* ScalarDivVectorial.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = ScalarDivVectorial.hpp; sourceTree = "<group>"; };
BEAE0D1F1CB2AF0B00CE333D /* ScalarDivVectorial.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = ScalarDivVectorial.hxx; sourceTree = "<group>"; };
BEAE0D221CB2AF1600CE333D /* ScalarDivVectorial.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = ScalarDivVectorial.cpp; sourceTree = "<group>"; };
BEAE0D231CB2AF1600CE333D /* ScalarDivVectorial.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = ScalarDivVectorial.hpp; sourceTree = "<group>"; };
BEAE0D241CB2AF1600CE333D /* ScalarDivVectorial.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = ScalarDivVectorial.hxx; sourceTree = "<group>"; };
BEAE12731C7CA858004FB5A7 /* GradOnGradientBasedElasticityTensor.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = GradOnGradientBasedElasticityTensor.cpp; sourceTree = "<group>"; };
BEAE12741C7CA858004FB5A7 /* GradOnGradientBasedElasticityTensor.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = GradOnGradientBasedElasticityTensor.hpp; sourceTree = "<group>"; };
BEAE12751C7CA858004FB5A7 /* GradOnGradientBasedElasticityTensor.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = GradOnGradientBasedElasticityTensor.hxx; sourceTree = "<group>"; };
......@@ -7234,6 +7242,9 @@
BE60DBFD1C7F417C007B334C /* Porosity.hxx */,
BE51CDDE1CA4560100B66D2A /* TransientSource.hpp */,
BE51CDDF1CA4560100B66D2A /* TransientSource.hxx */,
BEAE0D1D1CB2AF0B00CE333D /* ScalarDivVectorial.cpp */,
BEAE0D1E1CB2AF0B00CE333D /* ScalarDivVectorial.hpp */,
BEAE0D1F1CB2AF0B00CE333D /* ScalarDivVectorial.hxx */,
);
path = GlobalVariationalOperatorInstances;
sourceTree = "<group>";
......@@ -7255,6 +7266,9 @@
BE60DBFF1C7F4222007B334C /* Porosity.cpp */,
BE60DC001C7F4222007B334C /* Porosity.hpp */,
BE60DC011C7F4222007B334C /* Porosity.hxx */,
BEAE0D221CB2AF1600CE333D /* ScalarDivVectorial.cpp */,
BEAE0D231CB2AF1600CE333D /* ScalarDivVectorial.hpp */,
BEAE0D241CB2AF1600CE333D /* ScalarDivVectorial.hxx */,
BE83AF2C1C733E15002C6FA3 /* HyperelasticLaw */,
);
path = LocalVariationalOperatorInstances;
......@@ -12310,6 +12324,7 @@
BE5842131C6B8C920040C694 /* Mass.cpp in Sources */,
BEAEA8C71C71DFD30054FEF9 /* Porosity.cpp in Sources */,
BE669E1E1C7C9FDE00C521CA /* GradOnGradientBasedElasticityTensor.cpp in Sources */,
BEAE0D201CB2AF0B00CE333D /* ScalarDivVectorial.cpp in Sources */,
BE5842171C6B8C980040C694 /* Mass.cpp in Sources */,
BE60DC021C7F4222007B334C /* Porosity.cpp in Sources */,
BE0B0BE11C6A196C0038A9B9 /* main.cpp in Sources */,
......@@ -12317,6 +12332,7 @@
BE669E1A1C7C9FCE00C521CA /* GradOnGradientBasedElasticityTensor.cpp in Sources */,
BEF183281C6DE3A3008A6F1E /* Ale.cpp in Sources */,
BE60DBFE1C7F417C007B334C /* Porosity.cpp in Sources */,
BEAE0D251CB2AF1600CE333D /* ScalarDivVectorial.cpp in Sources */,
BE83AF3F1C73669D002C6FA3 /* BulkSolid.cpp in Sources */,
);
runOnlyForDeploymentPostprocessing = 0;
//! \file
//
//
// ScalarDivVectorial2.cpp
// HappyHeart
//
// Created by Sebastien Gilles on 30/04/15.
// Copyright (c) 2015 Inria. All rights reserved.
//
#include "ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/ScalarDivVectorial.hpp"
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace GlobalVariationalOperatorNS
{
ScalarDivVectorial::ScalarDivVectorial(const FEltSpace& felt_space,
const Unknown& vectorial_unknown,
const Unknown& scalar_unknown,
unsigned int geom_mesh_region_dimension,
const quadrature_rule_per_topology_type& quadrature_rule_per_topology,
const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>::type& porosity,
DoComputeProcessorWiseLocal2Global do_consider_processor_wise_local_2_global)
: parent(felt_space,
vectorial_unknown,
scalar_unknown,
geom_mesh_region_dimension,
quadrature_rule_per_topology,
AllocateGradientFEltPhi::yes,
do_consider_processor_wise_local_2_global,
porosity
)
{
assert(vectorial_unknown.GetNature() == UnknownNS::Nature::vectorial);
assert(scalar_unknown.GetNature() == UnknownNS::Nature::scalar);
}
const std::string& ScalarDivVectorial::ClassName()
{
static std::string name("ScalarDivVectorial");
return name;
}
} // namespace GlobalVariationalOperatorNS
} // namespace PoromechanicsNS
} // namespace HappyHeart
//! \file
//
//
// ScalarDivVectorial.hpp
// HappyHeart
//
// Created by Sebastien Gilles on 08/09/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_SCALAR_DIV_VECTORIAL_HPP_
# define HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_SCALAR_DIV_VECTORIAL_HPP_
# include "Parameters/ParameterAtDof.hpp"
# include "Operators/GlobalVariationalOperator/GlobalVariationalOperator.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/LocalVariationalOperatorInstances/ScalarDivVectorial.hpp"
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace GlobalVariationalOperatorNS
{
/*!
* \brief Implementation of global Stokes operator.
*
* \todo Improve the comment by writing its mathematical definition!
*/
class ScalarDivVectorial final : public GlobalVariationalOperator
<
ScalarDivVectorial,
LocalVariationalOperatorNS::ScalarDivVectorial
>
{
public:
//! Returns the name of the operator.
static const std::string& ClassName();
//! Convenient alias to pinpoint the GlobalVariationalOperator parent.
using parent = GlobalVariationalOperator
<
ScalarDivVectorial,
LocalVariationalOperatorNS::ScalarDivVectorial
>;
//! Friendship to the parent class so that the CRTP can reach private methods defined below.
friend parent;
//! Unique ptr.
using const_unique_ptr = std::unique_ptr<const ScalarDivVectorial>;
public:
/// \name Special members.
///@{
//! Constructor.
explicit ScalarDivVectorial(const FEltSpace& felt_space,
const Unknown& vectorial_unknown,
const Unknown& scalar_unknown,
unsigned int geom_mesh_region_dimension,
const quadrature_rule_per_topology_type& quadrature_rule_per_topology,
const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>::type& porosity,
DoComputeProcessorWiseLocal2Global do_consider_processor_wise_local_2_global);
//! Destructor.
~ScalarDivVectorial() = default;
//! Move constructor.
ScalarDivVectorial(ScalarDivVectorial&&) = default;
//! Copy constructor.
ScalarDivVectorial(const ScalarDivVectorial&) = delete;
//! Copy affectation.
ScalarDivVectorial& operator=(const ScalarDivVectorial&) = delete;
//! Move affectation.
ScalarDivVectorial& operator=(ScalarDivVectorial&&) = delete;
///@}
/*!
* \brief Assemble into one or several matrices.
*
* \tparam LinearAlgebraTupleT A tuple that may include \a GlobalMatrixWithCoefficient objects.
*
* \param[in] global_matrix_with_coeff_list List of global matrices into which the operator is
* assembled. These matrices 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.
*
*/
template<class LinearAlgebraTupleT>
void Assemble(LinearAlgebraTupleT&& global_matrix_with_coeff_list, const Domain& domain = Domain()) const;
};
} // namespace GlobalVariationalOperatorNS
} // namespace PoromechanicsNS
} // namespace HappyHeart
# include "ModelInstances/UnderDevelopment/Poromechanics/GlobalVariationalOperatorInstances/ScalarDivVectorial.hxx"
#endif // HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_SCALAR_DIV_VECTORIAL_HPP_
//! \file
//
//
// ScalarDivVectorial.hxx
// HappyHeart
//
// Created by Sebastien Gilles on 08/09/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_SCALAR_DIV_VECTORIAL_HXX_
# define HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_SCALAR_DIV_VECTORIAL_HXX_
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace GlobalVariationalOperatorNS
{
template<class LinearAlgebraTupleT>
inline void ScalarDivVectorial
::Assemble(LinearAlgebraTupleT&& linear_algebra_tuple, const Domain& domain) const
{
return parent::AssembleImpl(std::move(linear_algebra_tuple), domain);
}
} // namespace GlobalVariationalOperatorNS
} // namespace PoromechanicsNS
} // namespace HappyHeart
#endif // HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_GLOBAL_VARIATIONAL_OPERATOR_INSTANCES_x_SCALAR_DIV_VECTORIAL_HXX_
//! \file
//
//
// ScalarDivVectorial.cpp
// HappyHeart
//
// Created by Sebastien Gilles on 08/09/14.
// Copyright (c) 2014 Inria. All rights reserved.
//
#include "ThirdParty/Wrappers/Seldon/SeldonFunctions.hpp"
#ifndef NDEBUG
# include "ThirdParty/Wrappers/Seldon/MatrixOperations.hpp"
#endif // NDEBUG
#include "ModelInstances/UnderDevelopment/Poromechanics/LocalVariationalOperatorInstances/ScalarDivVectorial.hpp"
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace LocalVariationalOperatorNS
{
ScalarDivVectorial::ScalarDivVectorial(const ExtendedUnknown::vector_const_shared_ptr& a_unknown_storage,
elementary_data_type&& a_elementary_data,
const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>::type& porosity)
: BilinearLocalVariationalOperator(a_unknown_storage, std::move(a_elementary_data)),
porosity_parent(porosity)
{
#ifndef NDEBUG
{
const auto& unknown_storage = this->GetExtendedUnknownList();
assert(unknown_storage.size() == 2ul && "Namesake GlobalVariationalOperator should provide the correct one.");
assert(!(!unknown_storage.front()));
assert("Namesake GlobalVariationalOperator should provide two unknowns, the first being vectorial"
&& unknown_storage.front()->GetUnknown().GetNature() == UnknownNS::Nature::vectorial);
assert(!(!unknown_storage.back()));
assert("Namesake GlobalVariationalOperator should provide two unknowns, the second being scalar"
&& unknown_storage.back()->GetUnknown().GetNature() == UnknownNS::Nature::scalar);
}
#endif // NDEBUG
}
ScalarDivVectorial::~ScalarDivVectorial() = default;
const std::string& ScalarDivVectorial::ClassName()
{
static std::string name("ScalarDivVectorial");
return name;
}
void ScalarDivVectorial::ComputeEltArray()
{
static_assert(std::is_same<matrix_type, LocalSymmetricMatrix>(),
"Current implementation relies on the symmetry on the matrix (only upper right part "
"of the matrix is filled.");
auto& elementary_data = GetNonCstElementaryData();
auto& matrix_result = elementary_data.GetNonCstMatrixResult();
const auto& scalar_ref_felt =
elementary_data.GetRefFElt(GetNthUnknown(EnumUnderlyingType(UnknownIndex::scalar)));
const auto& vectorial_ref_felt =
elementary_data.GetRefFElt(GetNthUnknown(EnumUnderlyingType(UnknownIndex::vectorial)));
const auto Nnode_velocity = static_cast<int>(vectorial_ref_felt.Nnode());
const auto Ndof_velocity = static_cast<int>(vectorial_ref_felt.Ndof());
const auto Nnode_pressure = static_cast<int>(scalar_ref_felt.Nnode());
const auto Ncomponent = static_cast<int>(vectorial_ref_felt.Ncomponent());
const auto& infos_at_quad_pt_list = elementary_data.GetInformationsAtQuadraturePointList();
const auto& porosity = porosity_parent::GetPorosity();
decltype(auto) geom_elt = elementary_data.GetCurrentGeomElt();
for (const auto& infos_at_quad_pt : infos_at_quad_pt_list)
{
const auto& felt_phi = infos_at_quad_pt.GetFEltPhi();
const auto& gradient_felt_phi = infos_at_quad_pt.GetGradientFEltPhi();
const double factor = infos_at_quad_pt.GetQuadraturePoint().GetWeight()
* infos_at_quad_pt.GetAbsoluteValueJacobianDeterminant()
* porosity.GetValue(infos_at_quad_pt.GetQuadraturePoint(), geom_elt);
for (int node_pressure_index = 0; node_pressure_index < Nnode_pressure; ++node_pressure_index)
{
const int& dof_pressure = node_pressure_index; // alias
const auto pressure_term = felt_phi(Nnode_velocity + node_pressure_index);
for (int node_velocity_index = 0; node_velocity_index < Nnode_velocity; ++node_velocity_index)
{
int dof_velocity_index = node_velocity_index;
for (int component = 0; component < Ncomponent; ++component, dof_velocity_index += Nnode_velocity)
{
matrix_result(dof_velocity_index, Ndof_velocity + dof_pressure)
+= factor * pressure_term * gradient_felt_phi(node_velocity_index, component);
}
}
}
}
}
} // namespace LocalVariationalOperatorNS
} // namespace PoromechanicsNS
} // namespace HappyHeart
//! \file
//
//
// ScalarDivVectorial.hpp
// HappyHeart
//
// Created by Sebastien Gilles on 08/09/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_SCALAR_DIV_VECTORIAL_HPP_
# define HAPPY_HEART_x_MODEL_INSTANCES_x_UNDER_DEVELOPMENT_x_POROMECHANICS_x_LOCAL_VARIATIONAL_OPERATOR_INSTANCES_x_SCALAR_DIV_VECTORIAL_HPP_
# include "Parameters/ParameterAtDof.hpp"
# include "Operators/LocalVariationalOperator/BilinearLocalVariationalOperator.hpp"
# include "ModelInstances/UnderDevelopment/Poromechanics/Crtp/Porosity.hpp"
namespace HappyHeart
{
namespace PoromechanicsNS
{
namespace LocalVariationalOperatorNS
{
/*!
* \brief Implementation of Stokes operator.
*
* \todo Improve the comment by writing its mathematical definition!
*/
class ScalarDivVectorial final
: public BilinearLocalVariationalOperator<LocalSymmetricMatrix>,
private Crtp::Porosity<ScalarDivVectorial>
{
public:
//! Alias to self.
using self = ScalarDivVectorial;
//! Alias to unique pointer.
using unique_ptr = std::unique_ptr<self>;
//! Returns the name of the operator.
static const std::string& ClassName();
//! Alias to crtp parent.
using porosity_parent = Crtp::Porosity<self>;
//! Alias to parent.
using parent = BilinearLocalVariationalOperator<LocalSymmetricMatrix>;
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
* two unknowns (the first one vectorial and the second one scalar).
* \param[in] elementary_data Elementary matrices and vectors that will perform the calculations.
*
* \internal <b><tt>[internal]</tt></b> This constructor must not be called manually: it is involved only in
* GlobalVariationalOperator<DerivedT, LocalVariationalOperatorT>::CreateLocalOperatorList() method.
*/
explicit ScalarDivVectorial(const ExtendedUnknown::vector_const_shared_ptr& unknown_list,
elementary_data_type&& elementary_data,
const ParameterAtDof<ParameterNS::Type::scalar, ParameterNS::TimeDependencyNS::None, 2u>::type& porosity);
//! Destructor.
~ScalarDivVectorial();
//! Copy constructor.
ScalarDivVectorial(const ScalarDivVectorial&) = delete;
//! Move constructor.
ScalarDivVectorial(ScalarDivVectorial&&) = delete;
//! Copy affectation.
ScalarDivVectorial& operator=(const ScalarDivVectorial&) = delete;
//! Move affectation.
ScalarDivVectorial& operator=(ScalarDivVectorial&&) = delete;
///@}
/*!
* \brief Compute the elementary matrix.
*
* Structure of the elementary matrix (with dimension = 3 and respectively (m, n) dofs for velocity
* and pressure):
*
* vx1, vx2, ..., vxm, vy1, vy2, ..., vym, vz1, vz2, ..., vzm, p1, p2, ..., pn
* ( ) // Same ordering for rows.
* ( )
* ( )
* ( )
* ( )
* ( )
* ( )
* ( )
* ( )
* ( )
* ( )
* ( )
* ( )
* ( )
* ( )
* ( )
*
* This operator fills only the mixed terms (v, p) and (p, v); the matrix is symmetric.
*
*/
void ComputeEltArray();
private:
enum class UnknownIndex
{
vectorial = 0,
scalar = 1
};
};