ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > Class Template Reference

Detailed Description

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
class ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >

Definition at line 49 of file RichardsMechanicsFEM.h.

#include <RichardsMechanicsFEM.h>

Inheritance diagram for ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:

Collaboration diagram for ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:

Public Types
using	ShapeMatricesTypeDisplacement
using	ShapeMatricesTypePressure
using	GlobalDimMatrixType
using	BMatricesType
using	KelvinVectorType = typename BMatricesType::KelvinVectorType
using	IpData
using	Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>
using	SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>

Public Member Functions
	RichardsMechanicsLocalAssembler (RichardsMechanicsLocalAssembler const &)=delete
	RichardsMechanicsLocalAssembler (RichardsMechanicsLocalAssembler &&)=delete
	RichardsMechanicsLocalAssembler (MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, RichardsMechanicsProcessData< DisplacementDim > &process_data)
std::size_t	setIPDataInitialConditions (std::string_view const name, double const *values, int const integration_order) override
void	setInitialConditionsConcrete (Eigen::VectorXd const local_x, double const t, int const process_id) override
void	assemble (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_rhs_data) override
void	assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &, std::vector< double > &, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
void	assembleWithJacobianForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
void	initializeConcrete () override
void	postTimestepConcrete (Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override
void	computeSecondaryVariableConcrete (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev) override
Eigen::Map< const Eigen::RowVectorXd >	getShapeMatrix (const unsigned integration_point) const override
	Provides the shape matrix at the given integration point.
std::vector< double >	getSigma () const override
std::vector< double > const &	getIntPtDarcyVelocity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double >	getSaturation () const override
std::vector< double > const &	getIntPtSaturation (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double >	getMicroSaturation () const override
std::vector< double > const &	getIntPtMicroSaturation (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double >	getMicroPressure () const override
std::vector< double > const &	getIntPtMicroPressure (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double >	getPorosity () const override
std::vector< double > const &	getIntPtPorosity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double >	getTransportPorosity () const override
std::vector< double > const &	getIntPtTransportPorosity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double > const &	getIntPtSigma (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double >	getSwellingStress () const override
std::vector< double > const &	getIntPtSwellingStress (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double >	getEpsilon () const override
std::vector< double > const &	getIntPtEpsilon (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
int	getMaterialID () const override
std::vector< double >	getMaterialStateVariableInternalState (std::function< std::span< double >(typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables &)> const &get_values_span, int const &n_components) const override
std::vector< double > const &	getIntPtDryDensitySolid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual	~LocalAssemblerInterface ()=default
virtual void	setInitialConditions (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, int const process_id)
virtual void	initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
virtual void	preAssemble (double const, double const, std::vector< double > const &)
virtual void	assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
virtual void	computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
virtual void	preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void	postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
void	postNonLinearSolver (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
virtual Eigen::Vector3d	getFlux (MathLib::Point3d const &, double const, std::vector< double > const &) const
virtual Eigen::Vector3d	getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double > > const &) const
	Fits to staggered scheme.
Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual	~ExtrapolatableElement ()=default

Static Public Attributes
static int const	KelvinVectorSize
static constexpr auto &	N_u_op

Private Member Functions
void	assembleWithJacobianForDeformationEquations (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void	assembleWithJacobianForPressureEquations (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
unsigned	getNumberOfIntegrationPoints () const override
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const &	getMaterialStateVariablesAt (unsigned integration_point) const override

Private Attributes
RichardsMechanicsProcessData< DisplacementDim > &	_process_data
std::vector< IpData, Eigen::aligned_allocator< IpData > >	_ip_data
NumLib::GenericIntegrationMethod const &	_integration_method
MeshLib::Element const &	_element
bool const	_is_axially_symmetric
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType >	_secondary_data

Static Private Attributes
static const int	pressure_index = 0
static const int	pressure_size = ShapeFunctionPressure::NPOINTS
static const int	displacement_index = ShapeFunctionPressure::NPOINTS
static const int	displacement_size

Member Typedef Documentation

BMatricesType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::BMatricesType

Initial value:

 
        BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 61 of file RichardsMechanicsFEM.h.

GlobalDimMatrixType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::GlobalDimMatrixType

Initial value:

 
        typename ShapeMatricesTypePressure::GlobalDimMatrixType

Definition at line 58 of file RichardsMechanicsFEM.h.

Invariants

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>

Definition at line 72 of file RichardsMechanicsFEM.h.

IpData

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::IpData

Initial value:

 
        IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement,
                             ShapeMatricesTypePressure, DisplacementDim,
                             ShapeFunctionDisplacement::NPOINTS>

Definition at line 65 of file RichardsMechanicsFEM.h.

KelvinVectorType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::KelvinVectorType = typename BMatricesType::KelvinVectorType

Definition at line 63 of file RichardsMechanicsFEM.h.

ShapeMatricesTypeDisplacement

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypeDisplacement

Initial value:

 
        ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 53 of file RichardsMechanicsFEM.h.

ShapeMatricesTypePressure

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypePressure

Initial value:

 
        ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>

Definition at line 55 of file RichardsMechanicsFEM.h.

SymmetricTensor

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

using ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>

Definition at line 74 of file RichardsMechanicsFEM.h.

Constructor & Destructor Documentation

RichardsMechanicsLocalAssembler() [1/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::RichardsMechanicsLocalAssembler ( RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > const & )

delete

RichardsMechanicsLocalAssembler() [2/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::RichardsMechanicsLocalAssembler ( RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > && )

delete

RichardsMechanicsLocalAssembler() [3/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::RichardsMechanicsLocalAssembler	(	MeshLib::Element const &	e,
		std::size_t const	,
		NumLib::GenericIntegrationMethod const &	integration_method,
		bool const	is_axially_symmetric,
		RichardsMechanicsProcessData< DisplacementDim > &	process_data )

Definition at line 117 of file RichardsMechanicsFEM-impl.h.

    : _process_data(process_data),
      _integration_method(integration_method),
      _element(e),
      _is_axially_symmetric(is_axially_symmetric)
{
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
 
    _ip_data.reserve(n_integration_points);
    _secondary_data.N_u.resize(n_integration_points);
 
    auto const shape_matrices_u =
        NumLib::initShapeMatrices<ShapeFunctionDisplacement,
                                  ShapeMatricesTypeDisplacement,
                                  DisplacementDim>(e, is_axially_symmetric,
                                                   _integration_method);
 
    auto const shape_matrices_p =
        NumLib::initShapeMatrices<ShapeFunctionPressure,
                                  ShapeMatricesTypePressure, DisplacementDim>(
            e, is_axially_symmetric, _integration_method);
 
    auto const& solid_material =
        MaterialLib::Solids::selectSolidConstitutiveRelation(
            _process_data.solid_materials, _process_data.material_ids,
            e.getID());
 
    auto const& medium = _process_data.media_map.getMedium(_element.getID());
 
    ParameterLib::SpatialPosition x_position;
    x_position.setElementID(_element.getID());
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        x_position.setIntegrationPoint(ip);
        _ip_data.emplace_back(solid_material);
        auto& ip_data = _ip_data[ip];
        auto const& sm_u = shape_matrices_u[ip];
        _ip_data[ip].integration_weight =
            _integration_method.getWeightedPoint(ip).getWeight() *
            sm_u.integralMeasure * sm_u.detJ;
 
        ip_data.N_u = sm_u.N;
        ip_data.dNdx_u = sm_u.dNdx;
 
        ip_data.N_p = shape_matrices_p[ip].N;
        ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
 
        // Initial porosity. Could be read from integration point data or mesh.
        ip_data.porosity =
            medium->property(MPL::porosity)
                .template initialValue<double>(
                    x_position,
                    std::numeric_limits<
                        double>::quiet_NaN() /* t independent */);
 
        ip_data.transport_porosity = ip_data.porosity;
        if (medium->hasProperty(MPL::PropertyType::transport_porosity))
        {
            ip_data.transport_porosity =
                medium->property(MPL::transport_porosity)
                    .template initialValue<double>(
                        x_position,
                        std::numeric_limits<
                            double>::quiet_NaN() /* t independent */);
        }
 
        _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
    }
}

References ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_element, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_integration_method, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_process_data, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data, MeshLib::Element::getID(), NumLib::GenericIntegrationMethod::getNumberOfPoints(), MathLib::WeightedPoint::getWeight(), NumLib::GenericIntegrationMethod::getWeightedPoint(), NumLib::initShapeMatrices(), ProcessLib::HydroMechanics::SecondaryData< ShapeMatrixType >::N_u, MaterialPropertyLib::porosity, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), and MaterialPropertyLib::transport_porosity.

Member Function Documentation

assemble()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assemble ( double const t, double const dt, std::vector< double > const & local_x, std::vector< double > const & local_x_prev, std::vector< double > & local_M_data, std::vector< double > & local_K_data, std::vector< double > & local_rhs_data )

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 354 of file RichardsMechanicsFEM-impl.h.

{
    assert(local_x.size() == pressure_size + displacement_size);
 
    auto p_L =
        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
            pressure_size> const>(local_x.data() + pressure_index,
                                  pressure_size);
 
    auto u =
        Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
            displacement_size> const>(local_x.data() + displacement_index,
                                      displacement_size);
 
    auto p_L_prev =
        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
            pressure_size> const>(local_x_prev.data() + pressure_index,
                                  pressure_size);
 
    auto u_prev =
        Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
            displacement_size> const>(local_x_prev.data() + displacement_index,
                                      displacement_size);
 
    auto K = MathLib::createZeroedMatrix<
        typename ShapeMatricesTypeDisplacement::template MatrixType<
            displacement_size + pressure_size,
            displacement_size + pressure_size>>(
        local_K_data, displacement_size + pressure_size,
        displacement_size + pressure_size);
 
    auto M = MathLib::createZeroedMatrix<
        typename ShapeMatricesTypeDisplacement::template MatrixType<
            displacement_size + pressure_size,
            displacement_size + pressure_size>>(
        local_M_data, displacement_size + pressure_size,
        displacement_size + pressure_size);
 
    auto rhs = MathLib::createZeroedVector<
        typename ShapeMatricesTypeDisplacement::template VectorType<
            displacement_size + pressure_size>>(
        local_rhs_data, displacement_size + pressure_size);
 
    auto const& identity2 = MathLib::KelvinVector::Invariants<
        MathLib::KelvinVector::kelvin_vector_dimensions(
            DisplacementDim)>::identity2;
 
    auto const& medium = _process_data.media_map.getMedium(_element.getID());
    auto const& liquid_phase = medium->phase("AqueousLiquid");
    auto const& solid_phase = medium->phase("Solid");
    MPL::VariableArray variables;
    MPL::VariableArray variables_prev;
 
    ParameterLib::SpatialPosition x_position;
    x_position.setElementID(_element.getID());
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        x_position.setIntegrationPoint(ip);
        auto const& w = _ip_data[ip].integration_weight;
 
        auto const& N_u = _ip_data[ip].N_u;
        auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
        auto const& N_p = _ip_data[ip].N_p;
        auto const& dNdx_p = _ip_data[ip].dNdx_p;
 
        auto const x_coord =
            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
                                           ShapeMatricesTypeDisplacement>(
                _element, N_u);
        auto const B =
            LinearBMatrix::computeBMatrix<DisplacementDim,
                                          ShapeFunctionDisplacement::NPOINTS,
                                          typename BMatricesType::BMatrixType>(
                dNdx_u, N_u, x_coord, _is_axially_symmetric);
 
        auto& eps = _ip_data[ip].eps;
        auto& eps_m = _ip_data[ip].eps_m;
        eps.noalias() = B * u;
 
        auto& S_L = _ip_data[ip].saturation;
        auto const S_L_prev = _ip_data[ip].saturation_prev;
 
        double p_cap_ip;
        NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
 
        double p_cap_prev_ip;
        NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
 
        variables.capillary_pressure = p_cap_ip;
        variables.liquid_phase_pressure = -p_cap_ip;
        // setting pG to 1 atm
        // TODO : rewrite equations s.t. p_L = pG-p_cap
        variables.gas_phase_pressure = 1.0e5;
 
        auto const temperature =
            medium->property(MPL::PropertyType::reference_temperature)
                .template value<double>(variables, x_position, t, dt);
        variables.temperature = temperature;
 
        auto const alpha =
            medium->property(MPL::PropertyType::biot_coefficient)
                .template value<double>(variables, x_position, t, dt);
        auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
            t, x_position, dt, temperature);
 
        auto const beta_SR =
            (1 - alpha) /
            _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
        variables.grain_compressibility = beta_SR;
 
        auto const rho_LR =
            liquid_phase.property(MPL::PropertyType::density)
                .template value<double>(variables, x_position, t, dt);
        variables.density = rho_LR;
 
        auto const& b = _process_data.specific_body_force;
 
        S_L = medium->property(MPL::PropertyType::saturation)
                  .template value<double>(variables, x_position, t, dt);
        variables.liquid_saturation = S_L;
        variables_prev.liquid_saturation = S_L_prev;
 
        // tangent derivative for Jacobian
        double const dS_L_dp_cap =
            medium->property(MPL::PropertyType::saturation)
                .template dValue<double>(variables,
                                         MPL::Variable::capillary_pressure,
                                         x_position, t, dt);
        // secant derivative from time discretization for storage
        // use tangent, if secant is not available
        double const DeltaS_L_Deltap_cap =
            (p_cap_ip == p_cap_prev_ip)
                ? dS_L_dp_cap
                : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
 
        auto const chi = [medium, x_position, t, dt](double const S_L)
        {
            MPL::VariableArray vs;
            vs.liquid_saturation = S_L;
            return medium->property(MPL::PropertyType::bishops_effective_stress)
                .template value<double>(vs, x_position, t, dt);
        };
        double const chi_S_L = chi(S_L);
        double const chi_S_L_prev = chi(S_L_prev);
 
        double const p_FR = -chi_S_L * p_cap_ip;
        variables.effective_pore_pressure = p_FR;
        variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
 
        // Set volumetric strain rate for the general case without swelling.
        variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
        variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
 
        auto& phi = _ip_data[ip].porosity;
        {  // Porosity update
            variables_prev.porosity = _ip_data[ip].porosity_prev;
            phi = medium->property(MPL::PropertyType::porosity)
                      .template value<double>(variables, variables_prev,
                                              x_position, t, dt);
            variables.porosity = phi;
        }
 
        if (alpha < phi)
        {
            OGS_FATAL(
                "RichardsMechanics: Biot-coefficient {} is smaller than "
                "porosity {} in element/integration point {}/{}.",
                alpha, phi, _element.getID(), ip);
        }
 
        // Swelling and possibly volumetric strain rate update.
        {
            auto& sigma_sw = _ip_data[ip].sigma_sw;
            auto const& sigma_sw_prev = _ip_data[ip].sigma_sw_prev;
 
            // If there is swelling, compute it. Update volumetric strain rate,
            // s.t. it corresponds to the mechanical part only.
            sigma_sw = sigma_sw_prev;
            if (solid_phase.hasProperty(
                    MPL::PropertyType::swelling_stress_rate))
            {
                auto const sigma_sw_dot =
                    MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
                        MPL::formEigenTensor<3>(
                            solid_phase[MPL::PropertyType::swelling_stress_rate]
                                .value(variables, variables_prev, x_position, t,
                                       dt)));
                sigma_sw += sigma_sw_dot * dt;
 
                // !!! Misusing volumetric strain for mechanical volumetric
                // strain just to update the transport porosity !!!
                variables.volumetric_strain +=
                    identity2.transpose() * C_el.inverse() * sigma_sw;
                variables_prev.volumetric_strain +=
                    identity2.transpose() * C_el.inverse() * sigma_sw_prev;
            }
 
            if (medium->hasProperty(MPL::PropertyType::transport_porosity))
            {
                variables_prev.transport_porosity =
                    _ip_data[ip].transport_porosity_prev;
 
                _ip_data[ip].transport_porosity =
                    medium->property(MPL::PropertyType::transport_porosity)
                        .template value<double>(variables, variables_prev,
                                                x_position, t, dt);
                variables.transport_porosity = _ip_data[ip].transport_porosity;
            }
            else
            {
                variables.transport_porosity = phi;
            }
        }
 
        double const k_rel =
            medium->property(MPL::PropertyType::relative_permeability)
                .template value<double>(variables, x_position, t, dt);
        auto const mu =
            liquid_phase.property(MPL::PropertyType::viscosity)
                .template value<double>(variables, x_position, t, dt);
 
        auto const& sigma_sw = _ip_data[ip].sigma_sw;
        auto const& sigma_eff = _ip_data[ip].sigma_eff;
 
        // Set mechanical variables for the intrinsic permeability model
        // For stress dependent permeability.
        {
            auto const sigma_total =
                (sigma_eff - alpha * p_FR * identity2).eval();
 
            // For stress dependent permeability.
            variables.total_stress.emplace<SymmetricTensor>(
                MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
                    sigma_total));
        }
 
        variables.equivalent_plastic_strain =
            _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
 
        auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
            medium->property(MPL::PropertyType::permeability)
                .value(variables, x_position, t, dt));
 
        GlobalDimMatrixType const rho_K_over_mu =
            K_intrinsic * rho_LR * k_rel / mu;
 
        //
        // displacement equation, displacement part
        //
        eps_m.noalias() =
            solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
                ? eps + C_el.inverse() * sigma_sw
                : eps;
        variables.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                eps_m);
 
        _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
                                                temperature);
 
        // p_SR
        variables.solid_grain_pressure =
            p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
        auto const rho_SR =
            solid_phase.property(MPL::PropertyType::density)
                .template value<double>(variables, x_position, t, dt);
 
        //
        // displacement equation, displacement part
        //
        double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
        rhs.template segment<displacement_size>(displacement_index).noalias() -=
            (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
 
        //
        // pressure equation, pressure part.
        //
        auto const beta_LR =
            1 / rho_LR *
            liquid_phase.property(MPL::PropertyType::density)
                .template dValue<double>(variables,
                                         MPL::Variable::liquid_phase_pressure,
                                         x_position, t, dt);
 
        double const a0 = S_L * (alpha - phi) * beta_SR;
        // Volumetric average specific storage of the solid and fluid phases.
        double const specific_storage =
            DeltaS_L_Deltap_cap * (p_cap_ip * a0 - phi) +
            S_L * (phi * beta_LR + a0);
        M.template block<pressure_size, pressure_size>(pressure_index,
                                                       pressure_index)
            .noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
 
        K.template block<pressure_size, pressure_size>(pressure_index,
                                                       pressure_index)
            .noalias() += dNdx_p.transpose() * rho_K_over_mu * dNdx_p * w;
 
        rhs.template segment<pressure_size>(pressure_index).noalias() +=
            dNdx_p.transpose() * rho_LR * rho_K_over_mu * b * w;
 
        //
        // displacement equation, pressure part
        //
        K.template block<displacement_size, pressure_size>(displacement_index,
                                                           pressure_index)
            .noalias() -= B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
 
        //
        // pressure equation, displacement part.
        //
        M.template block<pressure_size, displacement_size>(pressure_index,
                                                           displacement_index)
            .noalias() += N_p.transpose() * S_L * rho_LR * alpha *
                          identity2.transpose() * B * w;
    }
 
    if (_process_data.apply_mass_lumping)
    {
        auto Mpp = M.template block<pressure_size, pressure_size>(
            pressure_index, pressure_index);
        Mpp = Mpp.colwise().sum().eval().asDiagonal();
    }
}

References MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor< 3 >(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::VariableArray::grain_compressibility, NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::VariableArray::mechanical_strain, OGS_FATAL, MaterialPropertyLib::VariableArray::porosity, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::VariableArray::transport_porosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

assembleWithJacobian()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobian ( double const t, double const dt, std::vector< double > const & local_x, std::vector< double > const & local_x_prev, std::vector< double > & , std::vector< double > & , std::vector< double > & local_rhs_data, std::vector< double > & local_Jac_data )

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 690 of file RichardsMechanicsFEM-impl.h.

{
    assert(local_x.size() == pressure_size + displacement_size);
 
    auto p_L =
        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
            pressure_size> const>(local_x.data() + pressure_index,
                                  pressure_size);
 
    auto u =
        Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
            displacement_size> const>(local_x.data() + displacement_index,
                                      displacement_size);
 
    auto p_L_prev =
        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
            pressure_size> const>(local_x_prev.data() + pressure_index,
                                  pressure_size);
    auto u_prev =
        Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
            displacement_size> const>(local_x_prev.data() + displacement_index,
                                      displacement_size);
 
    auto local_Jac = MathLib::createZeroedMatrix<
        typename ShapeMatricesTypeDisplacement::template MatrixType<
            displacement_size + pressure_size,
            displacement_size + pressure_size>>(
        local_Jac_data, displacement_size + pressure_size,
        displacement_size + pressure_size);
 
    auto local_rhs = MathLib::createZeroedVector<
        typename ShapeMatricesTypeDisplacement::template VectorType<
            displacement_size + pressure_size>>(
        local_rhs_data, displacement_size + pressure_size);
 
    auto const& identity2 = MathLib::KelvinVector::Invariants<
        MathLib::KelvinVector::kelvin_vector_dimensions(
            DisplacementDim)>::identity2;
 
    typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
                                                         pressure_size);
 
    typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_p =
        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
                                                         pressure_size);
 
    typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S_Jpp =
        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
                                                         pressure_size);
 
    typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S =
        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
                                                         pressure_size);
 
    typename ShapeMatricesTypeDisplacement::template MatrixType<
        displacement_size, pressure_size>
        Kup = ShapeMatricesTypeDisplacement::template MatrixType<
            displacement_size, pressure_size>::Zero(displacement_size,
                                                    pressure_size);
 
    typename ShapeMatricesTypeDisplacement::template MatrixType<
        pressure_size, displacement_size>
        Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
            pressure_size, displacement_size>::Zero(pressure_size,
                                                    displacement_size);
 
    auto const& medium = _process_data.media_map.getMedium(_element.getID());
    auto const& liquid_phase = medium->phase("AqueousLiquid");
    auto const& solid_phase = medium->phase("Solid");
    MPL::VariableArray variables;
    MPL::VariableArray variables_prev;
 
    ParameterLib::SpatialPosition x_position;
    x_position.setElementID(_element.getID());
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        x_position.setIntegrationPoint(ip);
        auto const& w = _ip_data[ip].integration_weight;
 
        auto const& N_u = _ip_data[ip].N_u;
        auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
        auto const& N_p = _ip_data[ip].N_p;
        auto const& dNdx_p = _ip_data[ip].dNdx_p;
 
        auto const x_coord =
            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
                                           ShapeMatricesTypeDisplacement>(
                _element, N_u);
        auto const B =
            LinearBMatrix::computeBMatrix<DisplacementDim,
                                          ShapeFunctionDisplacement::NPOINTS,
                                          typename BMatricesType::BMatrixType>(
                dNdx_u, N_u, x_coord, _is_axially_symmetric);
 
        double p_cap_ip;
        NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
 
        double p_cap_prev_ip;
        NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
 
        variables.capillary_pressure = p_cap_ip;
        variables.liquid_phase_pressure = -p_cap_ip;
        // setting pG to 1 atm
        // TODO : rewrite equations s.t. p_L = pG-p_cap
        variables.gas_phase_pressure = 1.0e5;
 
        auto const temperature =
            medium->property(MPL::PropertyType::reference_temperature)
                .template value<double>(variables, x_position, t, dt);
        variables.temperature = temperature;
 
        auto& eps = _ip_data[ip].eps;
        auto& eps_m = _ip_data[ip].eps_m;
        eps.noalias() = B * u;
        auto const& sigma_eff = _ip_data[ip].sigma_eff;
        auto& S_L = _ip_data[ip].saturation;
        auto const S_L_prev = _ip_data[ip].saturation_prev;
        auto const alpha =
            medium->property(MPL::PropertyType::biot_coefficient)
                .template value<double>(variables, x_position, t, dt);
 
        auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
            t, x_position, dt, temperature);
 
        auto const beta_SR =
            (1 - alpha) /
            _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
        variables.grain_compressibility = beta_SR;
 
        auto const rho_LR =
            liquid_phase.property(MPL::PropertyType::density)
                .template value<double>(variables, x_position, t, dt);
        variables.density = rho_LR;
 
        auto const& b = _process_data.specific_body_force;
 
        S_L = medium->property(MPL::PropertyType::saturation)
                  .template value<double>(variables, x_position, t, dt);
        variables.liquid_saturation = S_L;
        variables_prev.liquid_saturation = S_L_prev;
 
        // tangent derivative for Jacobian
        double const dS_L_dp_cap =
            medium->property(MPL::PropertyType::saturation)
                .template dValue<double>(variables,
                                         MPL::Variable::capillary_pressure,
                                         x_position, t, dt);
        // secant derivative from time discretization for storage
        // use tangent, if secant is not available
        double const DeltaS_L_Deltap_cap =
            (p_cap_ip == p_cap_prev_ip)
                ? dS_L_dp_cap
                : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
 
        auto const chi = [medium, x_position, t, dt](double const S_L)
        {
            MPL::VariableArray vs;
            vs.liquid_saturation = S_L;
            return medium->property(MPL::PropertyType::bishops_effective_stress)
                .template value<double>(vs, x_position, t, dt);
        };
        double const chi_S_L = chi(S_L);
        double const chi_S_L_prev = chi(S_L_prev);
 
        double const p_FR = -chi_S_L * p_cap_ip;
        variables.effective_pore_pressure = p_FR;
        variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
 
        // Set volumetric strain rate for the general case without swelling.
        variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
        variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
 
        auto& phi = _ip_data[ip].porosity;
        {  // Porosity update
 
            variables_prev.porosity = _ip_data[ip].porosity_prev;
            phi = medium->property(MPL::PropertyType::porosity)
                      .template value<double>(variables, variables_prev,
                                              x_position, t, dt);
            variables.porosity = phi;
        }
 
        if (alpha < phi)
        {
            OGS_FATAL(
                "RichardsMechanics: Biot-coefficient {} is smaller than "
                "porosity {} in element/integration point {}/{}.",
                alpha, phi, _element.getID(), ip);
        }
 
        auto const mu =
            liquid_phase.property(MPL::PropertyType::viscosity)
                .template value<double>(variables, x_position, t, dt);
 
        // Swelling and possibly volumetric strain rate update.
        updateSwellingStressAndVolumetricStrain<DisplacementDim>(
            _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu,
            _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
            variables, variables_prev, x_position, t, dt);
 
        if (medium->hasProperty(MPL::PropertyType::transport_porosity))
        {
            if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
            {
                variables_prev.transport_porosity =
                    _ip_data[ip].transport_porosity_prev;
 
                _ip_data[ip].transport_porosity =
                    medium->property(MPL::PropertyType::transport_porosity)
                        .template value<double>(variables, variables_prev,
                                                x_position, t, dt);
                variables.transport_porosity = _ip_data[ip].transport_porosity;
            }
        }
        else
        {
            variables.transport_porosity = phi;
        }
 
        double const k_rel =
            medium->property(MPL::PropertyType::relative_permeability)
                .template value<double>(variables, x_position, t, dt);
 
        // Set mechanical variables for the intrinsic permeability model
        // For stress dependent permeability.
        {
            auto const sigma_total =
                (_ip_data[ip].sigma_eff + alpha * p_FR * identity2).eval();
            // For stress dependent permeability.
            variables.total_stress.emplace<SymmetricTensor>(
                MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
                    sigma_total));
        }
 
        variables.equivalent_plastic_strain =
            _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
 
        auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
            medium->property(MPL::PropertyType::permeability)
                .value(variables, x_position, t, dt));
 
        GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu;
 
        //
        // displacement equation, displacement part
        //
        auto& sigma_sw = _ip_data[ip].sigma_sw;
 
        eps_m.noalias() =
            solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
                ? eps + C_el.inverse() * sigma_sw
                : eps;
        variables.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                eps_m);
 
        auto C = _ip_data[ip].updateConstitutiveRelation(
            variables, t, x_position, dt, temperature);
 
        local_Jac
            .template block<displacement_size, displacement_size>(
                displacement_index, displacement_index)
            .noalias() += B.transpose() * C * B * w;
 
        // p_SR
        variables.solid_grain_pressure =
            p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
        auto const rho_SR =
            solid_phase.property(MPL::PropertyType::density)
                .template value<double>(variables, x_position, t, dt);
 
        double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
        local_rhs.template segment<displacement_size>(displacement_index)
            .noalias() -=
            (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
 
        //
        // displacement equation, pressure part
        //
        Kup.noalias() += B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
 
        auto const dchi_dS_L =
            medium->property(MPL::PropertyType::bishops_effective_stress)
                .template dValue<double>(variables,
                                         MPL::Variable::liquid_saturation,
                                         x_position, t, dt);
        local_Jac
            .template block<displacement_size, pressure_size>(
                displacement_index, pressure_index)
            .noalias() -= B.transpose() * alpha *
                          (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
                          identity2 * N_p * w;
 
        local_Jac
            .template block<displacement_size, pressure_size>(
                displacement_index, pressure_index)
            .noalias() +=
            N_u_op(N_u).transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
 
        // For the swelling stress with double structure model the corresponding
        // Jacobian u-p entry would be required, but it does not improve
        // convergence and sometimes worsens it:
        // if (medium->hasProperty(MPL::PropertyType::saturation_micro))
        // {
        //     -B.transpose() *
        //         dsigma_sw_dS_L_m* dS_L_m_dp_cap_m*(p_L_m - p_L_m_prev) /
        //         (p_cap_ip - p_cap_prev_ip) * N_p* w;
        // }
        if (!medium->hasProperty(MPL::PropertyType::saturation_micro) &&
            solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
        {
            using DimMatrix = Eigen::Matrix<double, 3, 3>;
            auto const dsigma_sw_dS_L =
                MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
                    solid_phase
                        .property(MPL::PropertyType::swelling_stress_rate)
                        .template dValue<DimMatrix>(
                            variables, variables_prev,
                            MPL::Variable::liquid_saturation, x_position, t,
                            dt));
            local_Jac
                .template block<displacement_size, pressure_size>(
                    displacement_index, pressure_index)
                .noalias() +=
                B.transpose() * dsigma_sw_dS_L * dS_L_dp_cap * N_p * w;
        }
        //
        // pressure equation, displacement part.
        //
        if (_process_data.explicit_hm_coupling_in_unsaturated_zone)
        {
            Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR * alpha *
                             identity2.transpose() * B * w;
        }
        else
        {
            Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
                             identity2.transpose() * B * w;
        }
 
        //
        // pressure equation, pressure part.
        //
        laplace_p.noalias() +=
            dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
 
        auto const beta_LR =
            1 / rho_LR *
            liquid_phase.property(MPL::PropertyType::density)
                .template dValue<double>(variables,
                                         MPL::Variable::liquid_phase_pressure,
                                         x_position, t, dt);
 
        double const a0 = (alpha - phi) * beta_SR;
        double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0);
        double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
 
        double const dspecific_storage_a_p_dp_cap =
            dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0);
        double const dspecific_storage_a_S_dp_cap =
            -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
 
        storage_p_a_p.noalias() +=
            N_p.transpose() * rho_LR * specific_storage_a_p * N_p * w;
 
        storage_p_a_S.noalias() -= N_p.transpose() * rho_LR *
                                   specific_storage_a_S * DeltaS_L_Deltap_cap *
                                   N_p * w;
 
        local_Jac
            .template block<pressure_size, pressure_size>(pressure_index,
                                                          pressure_index)
            .noalias() += N_p.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
                          rho_LR * dspecific_storage_a_p_dp_cap * N_p * w;
 
        storage_p_a_S_Jpp.noalias() -=
            N_p.transpose() * rho_LR *
            ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
             specific_storage_a_S * dS_L_dp_cap) /
            dt * N_p * w;
 
        if (!_process_data.explicit_hm_coupling_in_unsaturated_zone)
        {
            local_Jac
                .template block<pressure_size, pressure_size>(pressure_index,
                                                              pressure_index)
                .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
                              identity2.transpose() * B * (u - u_prev) / dt *
                              N_p * w;
        }
 
        double const dk_rel_dS_l =
            medium->property(MPL::PropertyType::relative_permeability)
                .template dValue<double>(variables,
                                         MPL::Variable::liquid_saturation,
                                         x_position, t, dt);
        typename ShapeMatricesTypeDisplacement::GlobalDimVectorType const
            grad_p_cap = -dNdx_p * p_L;
        local_Jac
            .template block<pressure_size, pressure_size>(pressure_index,
                                                          pressure_index)
            .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
                          dk_rel_dS_l * dS_L_dp_cap * N_p * w;
 
        local_Jac
            .template block<pressure_size, pressure_size>(pressure_index,
                                                          pressure_index)
            .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
                          dk_rel_dS_l * dS_L_dp_cap * N_p * w;
 
        local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
            dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
 
        if (medium->hasProperty(MPL::PropertyType::saturation_micro))
        {
            double const alpha_bar = _process_data.micro_porosity_parameters
                                         ->mass_exchange_coefficient;
            auto const p_L_m = _ip_data[ip].liquid_pressure_m;
            local_rhs.template segment<pressure_size>(pressure_index)
                .noalias() -=
                N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w;
 
            local_Jac
                .template block<pressure_size, pressure_size>(pressure_index,
                                                              pressure_index)
                .noalias() += N_p.transpose() * alpha_bar / mu * N_p * w;
            if (p_cap_ip != p_cap_prev_ip)
            {
                double const p_L_m_prev = _ip_data[ip].liquid_pressure_m_prev;
                local_Jac
                    .template block<pressure_size, pressure_size>(
                        pressure_index, pressure_index)
                    .noalias() += N_p.transpose() * alpha_bar / mu *
                                  (p_L_m - p_L_m_prev) /
                                  (p_cap_ip - p_cap_prev_ip) * N_p * w;
            }
        }
    }
 
    if (_process_data.apply_mass_lumping)
    {
        storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
        storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
        storage_p_a_S_Jpp =
            storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
    }
 
    // pressure equation, pressure part.
    local_Jac
        .template block<pressure_size, pressure_size>(pressure_index,
                                                      pressure_index)
        .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
 
    // pressure equation, displacement part.
    local_Jac
        .template block<pressure_size, displacement_size>(pressure_index,
                                                          displacement_index)
        .noalias() = Kpu / dt;
 
    // pressure equation
    local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
        laplace_p * p_L +
        (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
        Kpu * (u - u_prev) / dt;
 
    // displacement equation
    local_rhs.template segment<displacement_size>(displacement_index)
        .noalias() += Kup * p_L;
}

References MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::VariableArray::grain_compressibility, NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::VariableArray::mechanical_strain, OGS_FATAL, MaterialPropertyLib::VariableArray::porosity, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::VariableArray::transport_porosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

assembleWithJacobianForDeformationEquations()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForDeformationEquations ( double const t, double const dt, Eigen::VectorXd const & local_x, Eigen::VectorXd const & local_x_prev, std::vector< double > & local_M_data, std::vector< double > & local_K_data, std::vector< double > & local_b_data, std::vector< double > & local_Jac_data )

private

Assemble local matrices and vectors arise from the linearized discretized weak form of the residual of the momentum balance equation,

\[ \nabla (\sigma - \alpha_b p \mathrm{I}) = f \]

where \( \sigma\) is the effective stress tensor, \(p\) is the pore pressure, \(\alpha_b\) is the Biot constant, \(\mathrm{I}\) is the identity tensor, and \(f\) is the body force.

Parameters

t	Time
dt	Time increment
local_x	Nodal values of \(x\) of an element.
local_x_prev	Nodal values of \(x_{prev}\) of an element.
local_M_data	Mass matrix of an element, which takes the form of \( \int N^T N\mathrm{d}\Omega\). Not used.
local_K_data	Laplacian matrix of an element, which takes the form of \( \int (\nabla N)^T K \nabla N\mathrm{d}\Omega\). Not used.
local_b_data	Right hand side vector of an element.
local_Jac_data	Element Jacobian matrix for the Newton-Raphson method.

Definition at line 1481 of file RichardsMechanicsFEM-impl.h.

{
    OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
}

References OGS_FATAL.

assembleWithJacobianForPressureEquations()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForPressureEquations ( double const t, double const dt, Eigen::VectorXd const & local_x, Eigen::VectorXd const & local_x_prev, std::vector< double > & local_M_data, std::vector< double > & local_K_data, std::vector< double > & local_b_data, std::vector< double > & local_Jac_data )

private

Assemble local matrices and vectors arise from the linearized discretized weak form of the residual of the mass balance equation of single phase flow,

\[ \alpha \cdot{p} - \nabla (K (\nabla p + \rho g \nabla z) + \alpha_b \nabla \cdot \dot{u} = Q \]

where \( alpha\) is a coefficient may depend on storage or the fluid density change, \( \rho\) is the fluid density, \(g\) is the gravitational acceleration, \(z\) is the vertical coordinate, \(u\) is the displacement, and \(Q\) is the source/sink term.

Parameters

t	Time
dt	Time increment
local_x	Nodal values of \(x\) of an element.
local_x_prev	Nodal values of \(x_{prev}\) of an element.
local_M_data	Mass matrix of an element, which takes the form of \( \int N^T N\mathrm{d}\Omega\). Not used.
local_K_data	Laplacian matrix of an element, which takes the form of \( \int (\nabla N)^T K \nabla N\mathrm{d}\Omega\). Not used.
local_b_data	Right hand side vector of an element.
local_Jac_data	Element Jacobian matrix for the Newton-Raphson method.

Definition at line 1465 of file RichardsMechanicsFEM-impl.h.

{
    OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
}

References OGS_FATAL.

assembleWithJacobianForStaggeredScheme()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForStaggeredScheme ( double const t, double const dt, Eigen::VectorXd const & local_x, Eigen::VectorXd const & local_x_prev, int const process_id, std::vector< double > & local_M_data, std::vector< double > & local_K_data, std::vector< double > & local_b_data, std::vector< double > & local_Jac_data )

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1497 of file RichardsMechanicsFEM-impl.h.

{
    // For the equations with pressure
    if (process_id == 0)
    {
        assembleWithJacobianForPressureEquations(t, dt, local_x, local_x_prev,
                                                 local_M_data, local_K_data,
                                                 local_b_data, local_Jac_data);
        return;
    }
 
    // For the equations with deformation
    assembleWithJacobianForDeformationEquations(t, dt, local_x, local_x_prev,
                                                local_M_data, local_K_data,
                                                local_b_data, local_Jac_data);
}

computeSecondaryVariableConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::computeSecondaryVariableConcrete ( double const t, double const dt, Eigen::VectorXd const & local_x, Eigen::VectorXd const & local_x_prev )

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1522 of file RichardsMechanicsFEM-impl.h.

{
    auto p_L = local_x.template segment<pressure_size>(pressure_index);
    auto u = local_x.template segment<displacement_size>(displacement_index);
 
    auto p_L_prev =
        local_x_prev.template segment<pressure_size>(pressure_index);
    auto u_prev =
        local_x_prev.template segment<displacement_size>(displacement_index);
 
    auto const& identity2 = MathLib::KelvinVector::Invariants<
        MathLib::KelvinVector::kelvin_vector_dimensions(
            DisplacementDim)>::identity2;
 
    auto const& medium = _process_data.media_map.getMedium(_element.getID());
    auto const& liquid_phase = medium->phase("AqueousLiquid");
    auto const& solid_phase = medium->phase("Solid");
    MPL::VariableArray variables;
    MPL::VariableArray variables_prev;
 
    ParameterLib::SpatialPosition x_position;
    x_position.setElementID(_element.getID());
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
 
    double saturation_avg = 0;
    double porosity_avg = 0;
 
    using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
    KV sigma_avg = KV::Zero();
 
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        x_position.setIntegrationPoint(ip);
        auto const& N_p = _ip_data[ip].N_p;
        auto const& N_u = _ip_data[ip].N_u;
        auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
        auto const x_coord =
            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
                                           ShapeMatricesTypeDisplacement>(
                _element, N_u);
        auto const B =
            LinearBMatrix::computeBMatrix<DisplacementDim,
                                          ShapeFunctionDisplacement::NPOINTS,
                                          typename BMatricesType::BMatrixType>(
                dNdx_u, N_u, x_coord, _is_axially_symmetric);
 
        double p_cap_ip;
        NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
 
        double p_cap_prev_ip;
        NumLib::shapeFunctionInterpolate(-p_L_prev, N_p, p_cap_prev_ip);
 
        variables.capillary_pressure = p_cap_ip;
        variables.liquid_phase_pressure = -p_cap_ip;
        // setting pG to 1 atm
        // TODO : rewrite equations s.t. p_L = pG-p_cap
        variables.gas_phase_pressure = 1.0e5;
 
        auto const temperature =
            medium->property(MPL::PropertyType::reference_temperature)
                .template value<double>(variables, x_position, t, dt);
        variables.temperature = temperature;
 
        auto& eps = _ip_data[ip].eps;
        eps.noalias() = B * u;
        auto& eps_m = _ip_data[ip].eps_m;
        auto& S_L = _ip_data[ip].saturation;
        auto const S_L_prev = _ip_data[ip].saturation_prev;
        S_L = medium->property(MPL::PropertyType::saturation)
                  .template value<double>(variables, x_position, t, dt);
        variables.liquid_saturation = S_L;
        variables_prev.liquid_saturation = S_L_prev;
 
        auto const chi = [medium, x_position, t, dt](double const S_L)
        {
            MPL::VariableArray vs;
            vs.liquid_saturation = S_L;
            return medium->property(MPL::PropertyType::bishops_effective_stress)
                .template value<double>(vs, x_position, t, dt);
        };
        double const chi_S_L = chi(S_L);
        double const chi_S_L_prev = chi(S_L_prev);
 
        auto const alpha =
            medium->property(MPL::PropertyType::biot_coefficient)
                .template value<double>(variables, x_position, t, dt);
 
        auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
            t, x_position, dt, temperature);
 
        auto const beta_SR =
            (1 - alpha) /
            _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
        variables.grain_compressibility = beta_SR;
 
        variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
        variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
 
        // Set volumetric strain rate for the general case without swelling.
        variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
        variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
 
        auto& phi = _ip_data[ip].porosity;
        {  // Porosity update
            variables_prev.porosity = _ip_data[ip].porosity_prev;
            phi = medium->property(MPL::PropertyType::porosity)
                      .template value<double>(variables, variables_prev,
                                              x_position, t, dt);
            variables.porosity = phi;
        }
 
        auto const rho_LR =
            liquid_phase.property(MPL::PropertyType::density)
                .template value<double>(variables, x_position, t, dt);
        variables.density = rho_LR;
        auto const mu =
            liquid_phase.property(MPL::PropertyType::viscosity)
                .template value<double>(variables, x_position, t, dt);
 
        // Swelling and possibly volumetric strain rate update.
        updateSwellingStressAndVolumetricStrain<DisplacementDim>(
            _ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu,
            _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
            variables, variables_prev, x_position, t, dt);
 
        if (medium->hasProperty(MPL::PropertyType::transport_porosity))
        {
            if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
            {
                variables_prev.transport_porosity =
                    _ip_data[ip].transport_porosity_prev;
 
                _ip_data[ip].transport_porosity =
                    medium->property(MPL::PropertyType::transport_porosity)
                        .template value<double>(variables, variables_prev,
                                                x_position, t, dt);
                variables.transport_porosity = _ip_data[ip].transport_porosity;
            }
        }
        else
        {
            variables.transport_porosity = phi;
        }
 
        // Set mechanical variables for the intrinsic permeability model
        // For stress dependent permeability.
        {
            auto const sigma_total = (_ip_data[ip].sigma_eff +
                                      alpha * chi_S_L * identity2 * p_cap_ip)
                                         .eval();
            // For stress dependent permeability.
            variables.total_stress.emplace<SymmetricTensor>(
                MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
                    sigma_total));
        }
 
        variables.equivalent_plastic_strain =
            _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
 
        auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
            medium->property(MPL::PropertyType::permeability)
                .value(variables, x_position, t, dt));
 
        double const k_rel =
            medium->property(MPL::PropertyType::relative_permeability)
                .template value<double>(variables, x_position, t, dt);
 
        GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
 
        auto const& sigma_eff = _ip_data[ip].sigma_eff;
        double const p_FR = -chi_S_L * p_cap_ip;
        // p_SR
        variables.solid_grain_pressure =
            p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
        auto const rho_SR =
            solid_phase.property(MPL::PropertyType::density)
                .template value<double>(variables, x_position, t, dt);
        _ip_data[ip].dry_density_solid = (1 - phi) * rho_SR;
 
        auto& sigma_sw = _ip_data[ip].sigma_sw;
 
        eps_m.noalias() =
            solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
                ? eps + C_el.inverse() * sigma_sw
                : eps;
        variables.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                eps_m);
 
        _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
                                                temperature);
 
        auto const& b = _process_data.specific_body_force;
 
        // Compute the velocity
        auto const& dNdx_p = _ip_data[ip].dNdx_p;
        _ip_data[ip].v_darcy.noalias() =
            -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
 
        saturation_avg += S_L;
        porosity_avg += phi;
        sigma_avg += sigma_eff;
    }
    saturation_avg /= n_integration_points;
    porosity_avg /= n_integration_points;
    sigma_avg /= n_integration_points;
 
    (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
    (*_process_data.element_porosity)[_element.getID()] = porosity_avg;
 
    Eigen::Map<KV>(&(*_process_data.element_stresses)[_element.getID() *
                                                      KV::RowsAtCompileTime]) =
        MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_avg);
 
    NumLib::interpolateToHigherOrderNodes<
        ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
        DisplacementDim>(_element, _is_axially_symmetric, p_L,
                         *_process_data.pressure_interpolated);
}

References MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::VariableArray::grain_compressibility, NumLib::interpolateToHigherOrderNodes(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::VariableArray::porosity, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::VariableArray::transport_porosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

getEpsilon()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getEpsilon ( ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1247 of file RichardsMechanicsFEM-impl.h.

{
    constexpr int kelvin_vector_size =
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
 
    return transposeInPlace<kelvin_vector_size>(
        [this](std::vector<double>& values)
        { return getIntPtEpsilon(0, {}, {}, values); });
}

References MathLib::KelvinVector::kelvin_vector_dimensions().

getIntPtDarcyVelocity()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtDarcyVelocity ( const double t, std::vector< GlobalVector * > const & x, std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table, std::vector< double > & cache ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1300 of file RichardsMechanicsFEM-impl.h.

{
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
 
    cache.clear();
    auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
        double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
        cache, DisplacementDim, n_integration_points);
 
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        cache_matrix.col(ip).noalias() = _ip_data[ip].v_darcy;
    }
 
    return cache;
}

References MathLib::createZeroedMatrix().

getIntPtDryDensitySolid()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtDryDensitySolid ( const double t, std::vector< GlobalVector * > const & x, std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table, std::vector< double > & cache ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1451 of file RichardsMechanicsFEM-impl.h.

{
    return ProcessLib::getIntegrationPointScalarData(
        _ip_data, &IpData::dry_density_solid, cache);
}

References ProcessLib::getIntegrationPointScalarData().

getIntPtEpsilon()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEpsilon ( const double t, std::vector< GlobalVector * > const & x, std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table, std::vector< double > & cache ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1260 of file RichardsMechanicsFEM-impl.h.

{
    return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
        _ip_data, &IpData::eps, cache);
}

getIntPtMicroPressure()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtMicroPressure ( const double t, std::vector< GlobalVector * > const & x, std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table, std::vector< double > & cache ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1387 of file RichardsMechanicsFEM-impl.h.

{
    return ProcessLib::getIntegrationPointScalarData(
        _ip_data, &IpData::liquid_pressure_m, cache);
}

References ProcessLib::getIntegrationPointScalarData().

getIntPtMicroSaturation()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtMicroSaturation ( const double t, std::vector< GlobalVector * > const & x, std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table, std::vector< double > & cache ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1362 of file RichardsMechanicsFEM-impl.h.

{
    return ProcessLib::getIntegrationPointScalarData(
        _ip_data, &IpData::saturation_m, cache);
}

References ProcessLib::getIntegrationPointScalarData().

getIntPtPorosity()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtPorosity ( const double t, std::vector< GlobalVector * > const & x, std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table, std::vector< double > & cache ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1412 of file RichardsMechanicsFEM-impl.h.

{
    return ProcessLib::getIntegrationPointScalarData(_ip_data,
                                                     &IpData::porosity, cache);
}

References ProcessLib::getIntegrationPointScalarData().

getIntPtSaturation()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtSaturation ( const double t, std::vector< GlobalVector * > const & x, std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table, std::vector< double > & cache ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1337 of file RichardsMechanicsFEM-impl.h.

{
    return ProcessLib::getIntegrationPointScalarData(
        _ip_data, &IpData::saturation, cache);
}

References ProcessLib::getIntegrationPointScalarData().

getIntPtSigma()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtSigma ( const double t, std::vector< GlobalVector * > const & x, std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table, std::vector< double > & cache ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1189 of file RichardsMechanicsFEM-impl.h.

{
    return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
        _ip_data, &IpData::sigma_eff, cache);
}

getIntPtSwellingStress()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtSwellingStress ( const double t, std::vector< GlobalVector * > const & x, std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table, std::vector< double > & cache ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1217 of file RichardsMechanicsFEM-impl.h.

{
    constexpr int kelvin_vector_size =
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
    auto const n_integration_points = _ip_data.size();
 
    cache.clear();
    auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
        double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
        cache, kelvin_vector_size, n_integration_points);
 
    for (unsigned ip = 0; ip < n_integration_points; ++ip)
    {
        auto const& sigma_sw = _ip_data[ip].sigma_sw;
        cache_mat.col(ip) =
            MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_sw);
    }
 
    return cache;
}

References MathLib::createZeroedMatrix(), MathLib::KelvinVector::kelvin_vector_dimensions(), and MathLib::KelvinVector::kelvinVectorToSymmetricTensor().

getIntPtTransportPorosity()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtTransportPorosity ( const double t, std::vector< GlobalVector * > const & x, std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table, std::vector< double > & cache ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1437 of file RichardsMechanicsFEM-impl.h.

{
    return ProcessLib::getIntegrationPointScalarData(
        _ip_data, &IpData::transport_porosity, cache);
}

References ProcessLib::getIntegrationPointScalarData().

getMaterialID()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

int ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMaterialID ( ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1275 of file RichardsMechanicsFEM-impl.h.

{
    return _process_data.material_ids == nullptr
               ? 0
               : (*_process_data.material_ids)[_element.getID()];
}

getMaterialStateVariableInternalState()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMaterialStateVariableInternalState ( std::function< std::span< double >(typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables &)> const & get_values_span, int const & n_components ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1285 of file RichardsMechanicsFEM-impl.h.

{
    return ProcessLib::getIntegrationPointDataMaterialStateVariables(
        _ip_data, &IpData::material_state_variables, get_values_span,
        n_components);
}

References ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::getIntegrationPointDataMaterialStateVariables(), and ProcessLib::RichardsMechanics::IntegrationPointData< BMatricesType, ShapeMatrixTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, NPoints >::material_state_variables.

getMaterialStateVariablesAt()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMaterialStateVariablesAt ( unsigned integration_point ) const

overrideprivatevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1762 of file RichardsMechanicsFEM-impl.h.

{
    return *_ip_data[integration_point].material_state_variables;
}

getMicroPressure()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMicroPressure ( ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1377 of file RichardsMechanicsFEM-impl.h.

{
    std::vector<double> result;
    getIntPtMicroPressure(0, {}, {}, result);
    return result;
}

getMicroSaturation()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMicroSaturation ( ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1352 of file RichardsMechanicsFEM-impl.h.

{
    std::vector<double> result;
    getIntPtMicroSaturation(0, {}, {}, result);
    return result;
}

getNumberOfIntegrationPoints()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

unsigned ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getNumberOfIntegrationPoints ( ) const

overrideprivatevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1752 of file RichardsMechanicsFEM-impl.h.

{
    return _integration_method.getNumberOfPoints();
}

getPorosity()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getPorosity ( ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1402 of file RichardsMechanicsFEM-impl.h.

{
    std::vector<double> result;
    getIntPtPorosity(0, {}, {}, result);
    return result;
}

getSaturation()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getSaturation ( ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1327 of file RichardsMechanicsFEM-impl.h.

{
    std::vector<double> result;
    getIntPtSaturation(0, {}, {}, result);
    return result;
}

getShapeMatrix()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

Eigen::Map< const Eigen::RowVectorXd > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getShapeMatrix ( const unsigned integration_point ) const

inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 176 of file RichardsMechanicsFEM.h.

    {
        auto const& N_u = _secondary_data.N_u[integration_point];
 
        // assumes N is stored contiguously in memory
        return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
    }

References ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data, and ProcessLib::HydroMechanics::SecondaryData< ShapeMatrixType >::N_u.

getSigma()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getSigma ( ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1176 of file RichardsMechanicsFEM-impl.h.

{
    constexpr int kelvin_vector_size =
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
 
    return transposeInPlace<kelvin_vector_size>(
        [this](std::vector<double>& values)
        { return getIntPtSigma(0, {}, {}, values); });
}

References MathLib::KelvinVector::kelvin_vector_dimensions().

getSwellingStress()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getSwellingStress ( ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1204 of file RichardsMechanicsFEM-impl.h.

{
    constexpr int kelvin_vector_size =
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
 
    return transposeInPlace<kelvin_vector_size>(
        [this](std::vector<double>& values)
        { return getIntPtSwellingStress(0, {}, {}, values); });
}

References MathLib::KelvinVector::kelvin_vector_dimensions().

getTransportPorosity()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector< double > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getTransportPorosity ( ) const

overridevirtual

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 1427 of file RichardsMechanicsFEM-impl.h.

{
    std::vector<double> result;
    getIntPtTransportPorosity(0, {}, {}, result);
    return result;
}

initializeConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete ( )

inlineoverridevirtual

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 123 of file RichardsMechanicsFEM.h.

    {
        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();
 
        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            auto& ip_data = _ip_data[ip];
 
            ParameterLib::SpatialPosition const x_position{
                std::nullopt, _element.getID(), ip,
                MathLib::Point3d(
                    NumLib::interpolateCoordinates<
                        ShapeFunctionDisplacement,
                        ShapeMatricesTypeDisplacement>(_element, ip_data.N_u))};
 
            if (_process_data.initial_stress != nullptr)
            {
                ip_data.sigma_eff =
                    MathLib::KelvinVector::symmetricTensorToKelvinVector<
                        DisplacementDim>((*_process_data.initial_stress)(
                        std::numeric_limits<
                            double>::quiet_NaN() /* time independent */,
                        x_position));
            }
 
            double const t = 0;  // TODO (naumov) pass t from top
            ip_data.solid_material.initializeInternalStateVariables(
                t, x_position, *ip_data.material_state_variables);
 
            ip_data.pushBackState();
        }
    }

References ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_element, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_integration_method, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_process_data, MeshLib::Element::getID(), NumLib::GenericIntegrationMethod::getNumberOfPoints(), NumLib::interpolateCoordinates(), and MathLib::KelvinVector::symmetricTensorToKelvinVector().

postTimestepConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::postTimestepConcrete ( Eigen::VectorXd const & , Eigen::VectorXd const & , double const , double const , int const  )

inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 158 of file RichardsMechanicsFEM.h.

    {
        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();
 
        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            _ip_data[ip].pushBackState();
        }
    }

References ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_integration_method, ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, and NumLib::GenericIntegrationMethod::getNumberOfPoints().

setInitialConditionsConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setInitialConditionsConcrete ( Eigen::VectorXd const local_x, double const t, int const process_id )

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 280 of file RichardsMechanicsFEM-impl.h.

{
    assert(local_x.size() == pressure_size + displacement_size);
 
    auto const p_L = local_x.template segment<pressure_size>(pressure_index);
 
    constexpr double dt = std::numeric_limits<double>::quiet_NaN();
    auto const& medium = _process_data.media_map.getMedium(_element.getID());
    MPL::VariableArray variables;
 
    ParameterLib::SpatialPosition x_position;
    x_position.setElementID(_element.getID());
 
    auto const& solid_phase = medium->phase("Solid");
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        x_position.setIntegrationPoint(ip);
 
        auto const& N_p = _ip_data[ip].N_p;
 
        double p_cap_ip;
        NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
 
        variables.capillary_pressure = p_cap_ip;
        variables.liquid_phase_pressure = -p_cap_ip;
        // setting pG to 1 atm
        // TODO : rewrite equations s.t. p_L = pG-p_cap
        variables.gas_phase_pressure = 1.0e5;
        _ip_data[ip].liquid_pressure_m_prev = -p_cap_ip;
        _ip_data[ip].liquid_pressure_m = -p_cap_ip;
 
        auto const temperature =
            medium->property(MPL::PropertyType::reference_temperature)
                .template value<double>(variables, x_position, t, dt);
        variables.temperature = temperature;
 
        _ip_data[ip].saturation_prev =
            medium->property(MPL::PropertyType::saturation)
                .template value<double>(variables, x_position, t, dt);
 
        if (medium->hasProperty(MPL::PropertyType::saturation_micro))
        {
            MPL::VariableArray vars;
            vars.capillary_pressure = p_cap_ip;
            double const S_L_m =
                medium->property(MPL::PropertyType::saturation_micro)
                    .template value<double>(vars, x_position, t, dt);
            _ip_data[ip].saturation_m_prev = S_L_m;
        }
 
        // Set eps_m_prev from potentially non-zero eps and sigma_sw from
        // restart.
        auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
            t, x_position, dt, temperature);
        auto& eps = _ip_data[ip].eps;
        auto& sigma_sw = _ip_data[ip].sigma_sw;
 
        _ip_data[ip].eps_m_prev.noalias() =
            solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
                ? eps + C_el.inverse() * sigma_sw
                : eps;
    }
}

References MaterialPropertyLib::VariableArray::capillary_pressure, MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), NumLib::detail::shapeFunctionInterpolate(), and MaterialPropertyLib::VariableArray::temperature.

setIPDataInitialConditions()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::size_t ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setIPDataInitialConditions ( std::string_view const name, double const * values, int const integration_order )

overridevirtual

Returns

the number of read integration points.

Implements ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >.

Definition at line 198 of file RichardsMechanicsFEM-impl.h.

{
    if (integration_order !=
        static_cast<int>(_integration_method.getIntegrationOrder()))
    {
        OGS_FATAL(
            "Setting integration point initial conditions; The integration "
            "order of the local assembler for element {:d} is different "
            "from the integration order in the initial condition.",
            _element.getID());
    }
 
    if (name == "sigma")
    {
        if (_process_data.initial_stress != nullptr)
        {
            OGS_FATAL(
                "Setting initial conditions for stress from integration "
                "point data and from a parameter '{:s}' is not possible "
                "simultaneously.",
                _process_data.initial_stress->name);
        }
        return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
            values, _ip_data, &IpData::sigma_eff);
    }
 
    if (name == "saturation")
    {
        return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
                                                         &IpData::saturation);
    }
    if (name == "porosity")
    {
        return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
                                                         &IpData::porosity);
    }
    if (name == "transport_porosity")
    {
        return ProcessLib::setIntegrationPointScalarData(
            values, _ip_data, &IpData::transport_porosity);
    }
    if (name == "swelling_stress")
    {
        return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
            values, _ip_data, &IpData::sigma_sw);
    }
    if (name == "epsilon")
    {
        return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
            values, _ip_data, &IpData::eps);
    }
    if (name.starts_with("material_state_variable_"))
    {
        name.remove_prefix(24);
 
        // Using first ip data for solid material. TODO (naumov) move solid
        // material into element, store only material state in IPs.
        auto const& internal_variables =
            _ip_data[0].solid_material.getInternalVariables();
        if (auto const iv = std::find_if(
                begin(internal_variables), end(internal_variables),
                [&name](auto const& iv) { return iv.name == name; });
            iv != end(internal_variables))
        {
            DBUG("Setting material state variable '{:s}'", name);
            return ProcessLib::setIntegrationPointDataMaterialStateVariables(
                values, _ip_data, &IpData::material_state_variables,
                iv->reference);
        }
 
        ERR("Could not find variable {:s} in solid material model's internal "
            "variables.",
            name);
    }
    return 0;
}

References DBUG(), ERR(), OGS_FATAL, ProcessLib::setIntegrationPointDataMaterialStateVariables(), and ProcessLib::setIntegrationPointScalarData().

Member Data Documentation

_element

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

MeshLib::Element const& ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_element

private

Definition at line 336 of file RichardsMechanicsFEM.h.

Referenced by ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::RichardsMechanicsLocalAssembler(), and ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete().

_integration_method

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

NumLib::GenericIntegrationMethod const& ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_integration_method

private

Definition at line 335 of file RichardsMechanicsFEM.h.

Referenced by ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::RichardsMechanicsLocalAssembler(), ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete(), and ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::postTimestepConcrete().

_ip_data

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

std::vector<IpData, Eigen::aligned_allocator<IpData> > ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data

private

Definition at line 333 of file RichardsMechanicsFEM.h.

Referenced by ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::RichardsMechanicsLocalAssembler(), ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMaterialStateVariableInternalState(), ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete(), and ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::postTimestepConcrete().

_is_axially_symmetric

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

bool const ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_is_axially_symmetric

private

Definition at line 337 of file RichardsMechanicsFEM.h.

_process_data

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

RichardsMechanicsProcessData<DisplacementDim>& ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_process_data

private

Definition at line 331 of file RichardsMechanicsFEM.h.

Referenced by ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::RichardsMechanicsLocalAssembler(), and ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete().

_secondary_data

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType> ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data

private

Definition at line 340 of file RichardsMechanicsFEM.h.

Referenced by ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::RichardsMechanicsLocalAssembler(), and ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getShapeMatrix().

displacement_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

const int ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::displacement_index = ShapeFunctionPressure::NPOINTS

staticprivate

Definition at line 344 of file RichardsMechanicsFEM.h.

displacement_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

const int ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::displacement_size

staticprivate

Initial value:

=
        ShapeFunctionDisplacement::NPOINTS * DisplacementDim

Definition at line 345 of file RichardsMechanicsFEM.h.

KelvinVectorSize

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

int const ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::KelvinVectorSize

static

Initial value:

=
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim)

Definition at line 70 of file RichardsMechanicsFEM.h.

N_u_op

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

constexpr auto& ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::N_u_op

staticconstexpr

Initial value:

= MathLib::eigenBlockMatrixView<
        DisplacementDim,
        typename ShapeMatricesTypeDisplacement::NodalRowVectorType>

Definition at line 76 of file RichardsMechanicsFEM.h.

pressure_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

const int ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::pressure_index = 0

staticprivate

Definition at line 342 of file RichardsMechanicsFEM.h.

pressure_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

const int ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::pressure_size = ShapeFunctionPressure::NPOINTS

staticprivate

Definition at line 343 of file RichardsMechanicsFEM.h.

---

The documentation for this class was generated from the following files:

• ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
• ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h