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 45 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)
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 > &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_b_data, std::vector< double > &local_Jac_data) override
void	initializeConcrete () 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.
int	getNumberOfVectorElementsForDeformation () const override
Public Member Functions inherited from ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >
	LocalAssemblerInterface (MeshLib::Element const &e, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, RichardsMechanicsProcessData< DisplacementDim > &process_data)
std::size_t	setIPDataInitialConditions (std::string_view name, double const *values, int const integration_order)
std::vector< double >	getMaterialStateVariableInternalState (std::function< std::span< double >(typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables &)> const &get_values_span, int const &n_components) const
unsigned	getNumberOfIntegrationPoints () const
int	getMaterialID () const
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const &	getMaterialStateVariablesAt (unsigned integration_point) const
void	postTimestepConcrete (Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override final
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_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_b_data, std::vector< double > &local_Jac_data)

Static Private Member Functions
static void	assembleWithJacobianEvalConstitutiveSetting (double const t, double const dt, ParameterLib::SpatialPosition const &x_position, IpData &ip_data, MPL::VariableArray &variables, MPL::VariableArray &variables_prev, MPL::Medium const *const medium, TemperatureData const T_data, CapillaryPressureData< DisplacementDim > const &p_cap_data, ConstitutiveData< DisplacementDim > &CD, StatefulData< DisplacementDim > &SD, StatefulDataPrev< DisplacementDim > const &SD_prev, std::optional< MicroPorosityParameters > const &micro_porosity_parameters, MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material, ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > &material_state_data)
static constexpr auto	localDOF (auto const &x)

Private Attributes
std::vector< IpData, Eigen::aligned_allocator< IpData > >	ip_data_
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

Additional Inherited Members
Static Public Member Functions inherited from ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >
static auto	getReflectionDataForOutput ()
Protected Attributes inherited from ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >
RichardsMechanicsProcessData< DisplacementDim > &	process_data_
NumLib::GenericIntegrationMethod const &	integration_method_
MeshLib::Element const &	element_
bool const	is_axially_symmetric_
MaterialLib::Solids::MechanicsBase< DisplacementDim > const &	solid_material_
std::vector< StatefulData< DisplacementDim > >	current_states_
std::vector< StatefulDataPrev< DisplacementDim > >	prev_states_
std::vector< ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > >	material_states_
std::vector< OutputData< DisplacementDim > >	output_data_

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 57 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 54 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 68 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 61 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 59 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 49 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 51 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 70 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

References RichardsMechanicsLocalAssembler().

Referenced by RichardsMechanicsLocalAssembler(), RichardsMechanicsLocalAssembler(), and assembleWithJacobianEvalConstitutiveSetting().

RichardsMechanicsLocalAssembler() [2/3]

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

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

delete

References RichardsMechanicsLocalAssembler().

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 122 of file RichardsMechanicsFEM-impl.h.

    : LocalAssemblerInterface<DisplacementDim>{
          e, integration_method, is_axially_symmetric, process_data}
{
    unsigned const n_integration_points =
        this->integration_method_.getNumberOfPoints();
 
    ip_data_.resize(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,
                                                   this->integration_method_);
 
    auto const shape_matrices_p =
        NumLib::initShapeMatrices<ShapeFunctionPressure,
                                  ShapeMatricesTypePressure, DisplacementDim>(
            e, is_axially_symmetric, this->integration_method_);
 
    auto const& medium =
        this->process_data_.media_map.getMedium(this->element_.getID());
 
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        auto& ip_data = ip_data_[ip];
        auto const& sm_u = shape_matrices_u[ip];
        ip_data_[ip].integration_weight =
            this->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;
 
        ParameterLib::SpatialPosition x_position = {
            std::nullopt, this->element_.getID(),
            MathLib::Point3d(
                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
                                               ShapeMatricesTypeDisplacement>(
                    this->element_, ip_data.N_u))};
 
        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.
        auto& porosity =
            std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
                this->current_states_[ip])
                .phi;
        porosity = medium->property(MPL::porosity)
                       .template initialValue<double>(
                           x_position,
                           std::numeric_limits<
                               double>::quiet_NaN() /* t independent */);
 
        auto& transport_porosity =
            std::get<
                ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(
                this->current_states_[ip])
                .phi;
        transport_porosity = porosity;
        if (medium->hasProperty(MPL::PropertyType::transport_porosity))
        {
            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::LocalAssemblerInterface< DisplacementDim >::LocalAssemblerInterface().

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 371 of file RichardsMechanicsFEM-impl.h.

{
    assert(local_x.size() == pressure_size + displacement_size);
 
    auto const [p_L, u] = localDOF(local_x);
    auto const [p_L_prev, u_prev] = localDOF(local_x_prev);
 
    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 =
        this->process_data_.media_map.getMedium(this->element_.getID());
    auto const& liquid_phase =
        medium->phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
    auto const& solid_phase =
        medium->phase(MaterialPropertyLib::PhaseName::Solid);
    MPL::VariableArray variables;
    MPL::VariableArray variables_prev;
 
    ParameterLib::SpatialPosition x_position;
    x_position.setElementID(this->element_.getID());
 
    unsigned const n_integration_points =
        this->integration_method_.getNumberOfPoints();
    for (unsigned ip = 0; ip < n_integration_points; 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;
 
        x_position = {
            std::nullopt, this->element_.getID(),
            MathLib::Point3d(
                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
                                               ShapeMatricesTypeDisplacement>(
                    this->element_, N_u))};
        auto const x_coord = x_position.getCoordinates().value()[0];
 
        auto const B =
            LinearBMatrix::computeBMatrix<DisplacementDim,
                                          ShapeFunctionDisplacement::NPOINTS,
                                          typename BMatricesType::BMatrixType>(
                dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
 
        auto& eps =
            std::get<StrainData<DisplacementDim>>(this->current_states_[ip]);
        eps.eps.noalias() = B * u;
 
        auto& S_L =
            std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
                this->current_states_[ip])
                .S_L;
        auto const S_L_prev =
            std::get<
                PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
                this->prev_states_[ip])
                ->S_L;
 
        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& SD = this->current_states_[ip];
        variables.stress =
            std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                DisplacementDim>>(SD)
                .sigma_eff;
        // Set mechanical strain temporary to compute tangent stiffness.
        variables.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                eps.eps);
        auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
            variables, t, x_position, dt, this->solid_material_,
            *this->material_states_[ip].material_state_variables);
 
        auto const beta_SR = (1 - alpha) / this->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 = this->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(eps.eps);
        variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
 
        auto& phi = std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
                        this->current_states_[ip])
                        .phi;
        {  // Porosity update
            auto const phi_prev = std::get<PrevState<
                ProcessLib::ThermoRichardsMechanics::PorosityData>>(
                                      this->prev_states_[ip])
                                      ->phi;
            variables_prev.porosity = phi_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, this->element_.getID(), ip);
        }
 
        // Swelling and possibly volumetric strain rate update.
        {
            auto& sigma_sw =
                std::get<ProcessLib::ThermoRichardsMechanics::
                             ConstitutiveStress_StrainTemperature::
                                 SwellingDataStateful<DisplacementDim>>(
                    this->current_states_[ip])
                    .sigma_sw;
            auto const& sigma_sw_prev = std::get<PrevState<
                ProcessLib::ThermoRichardsMechanics::
                    ConstitutiveStress_StrainTemperature::SwellingDataStateful<
                        DisplacementDim>>>(this->prev_states_[ip])
                                            ->sigma_sw;
 
            // 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;
 
                variables.volumetric_mechanical_strain =
                    variables.volumetric_strain +
                    identity2.transpose() * C_el.inverse() * sigma_sw;
                variables_prev.volumetric_mechanical_strain =
                    variables_prev.volumetric_strain +
                    identity2.transpose() * C_el.inverse() * sigma_sw_prev;
            }
            else
            {
                variables.volumetric_mechanical_strain =
                    variables.volumetric_strain;
                variables_prev.volumetric_mechanical_strain =
                    variables_prev.volumetric_strain;
            }
 
            if (medium->hasProperty(MPL::PropertyType::transport_porosity))
            {
                auto& transport_porosity =
                    std::get<ProcessLib::ThermoRichardsMechanics::
                                 TransportPorosityData>(
                        this->current_states_[ip])
                        .phi;
                auto const transport_porosity_prev =
                    std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
                                           TransportPorosityData>>(
                        this->prev_states_[ip])
                        ->phi;
                variables_prev.transport_porosity = transport_porosity_prev;
 
                transport_porosity =
                    medium->property(MPL::PropertyType::transport_porosity)
                        .template value<double>(variables, variables_prev,
                                                x_position, t, dt);
                variables.transport_porosity = 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 =
            std::get<ProcessLib::ThermoRichardsMechanics::
                         ConstitutiveStress_StrainTemperature::
                             SwellingDataStateful<DisplacementDim>>(
                this->current_states_[ip])
                .sigma_sw;
        auto const& sigma_eff =
            std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                DisplacementDim>>(this->current_states_[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 =
            this->material_states_[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
        //
        {
            auto& eps_m = std::get<ProcessLib::ConstitutiveRelations::
                                       MechanicalStrainData<DisplacementDim>>(
                              this->current_states_[ip])
                              .eps_m;
            eps_m.noalias() =
                solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
                    ? eps.eps + C_el.inverse() * sigma_sw
                    : eps.eps;
            variables.mechanical_strain.emplace<
                MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                eps_m);
        }
 
        {
            auto& SD = this->current_states_[ip];
            auto const& SD_prev = this->prev_states_[ip];
            auto& sigma_eff =
                std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                    DisplacementDim>>(SD);
            auto const& sigma_eff_prev =
                std::get<PrevState<ProcessLib::ConstitutiveRelations::
                                       EffectiveStressData<DisplacementDim>>>(
                    SD_prev);
            auto const& eps_m =
                std::get<ProcessLib::ConstitutiveRelations::
                             MechanicalStrainData<DisplacementDim>>(SD);
            auto& eps_m_prev =
                std::get<PrevState<ProcessLib::ConstitutiveRelations::
                                       MechanicalStrainData<DisplacementDim>>>(
                    SD_prev);
 
            auto const C = ip_data_[ip].updateConstitutiveRelation(
                variables, t, x_position, dt, temperature, sigma_eff,
                sigma_eff_prev, eps_m, eps_m_prev, this->solid_material_,
                this->material_states_[ip].material_state_variables);
 
            if (this->process_data_.use_numerical_jacobian)
            {
                K.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);
 
        //
        // 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 (this->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::AqueousLiquid, assemble(), MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::capillary_pressure, MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::current_states_, MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, displacement_index, displacement_size, MaterialPropertyLib::VariableArray::effective_pore_pressure, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::element_, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, ParameterLib::SpatialPosition::getCoordinates(), MaterialPropertyLib::VariableArray::grain_compressibility, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), ip_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::is_axially_symmetric_, MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, localDOF(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::material_states_, MaterialPropertyLib::VariableArray::mechanical_strain, N_u_op, OGS_FATAL, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, pressure_size, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::prev_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::process_data_, MaterialPropertyLib::reference_temperature, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, ParameterLib::SpatialPosition::setElementID(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::Solid, MaterialPropertyLib::VariableArray::solid_grain_pressure, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::solid_material_, MaterialPropertyLib::VariableArray::stress, MaterialPropertyLib::swelling_stress_rate, MaterialPropertyLib::VariableArray::temperature, MathLib::KelvinVector::tensorToKelvin(), MaterialPropertyLib::VariableArray::total_stress, MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace(), MaterialPropertyLib::transport_porosity, MaterialPropertyLib::VariableArray::transport_porosity, MaterialPropertyLib::viscosity, MaterialPropertyLib::VariableArray::volumetric_mechanical_strain, and MaterialPropertyLib::VariableArray::volumetric_strain.

Referenced by assemble().

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 > & local_rhs_data, std::vector< double > & local_Jac_data )

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

{
    assert(local_x.size() == pressure_size + displacement_size);
 
    auto const [p_L, u] = localDOF(local_x);
    auto const [p_L_prev, u_prev] = localDOF(local_x_prev);
 
    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 =
        this->process_data_.media_map.getMedium(this->element_.getID());
    auto const& liquid_phase =
        medium->phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
    auto const& solid_phase =
        medium->phase(MaterialPropertyLib::PhaseName::Solid);
    MPL::VariableArray variables;
    MPL::VariableArray variables_prev;
 
    unsigned const n_integration_points =
        this->integration_method_.getNumberOfPoints();
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        ConstitutiveData<DisplacementDim> CD;
        auto& SD = this->current_states_[ip];
        auto const& SD_prev = this->prev_states_[ip];
        [[maybe_unused]] auto models = createConstitutiveModels(
            this->process_data_, this->solid_material_);
 
        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;
 
        ParameterLib::SpatialPosition x_position = {
            std::nullopt, this->element_.getID(),
            MathLib::Point3d(
                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
                                               ShapeMatricesTypeDisplacement>(
                    this->element_, N_u))};
        auto const x_coord = x_position.getCoordinates().value()[0];
 
        auto const B =
            LinearBMatrix::computeBMatrix<DisplacementDim,
                                          ShapeFunctionDisplacement::NPOINTS,
                                          typename BMatricesType::BMatrixType>(
                dNdx_u, N_u, x_coord, this->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;
 
        std::get<StrainData<DisplacementDim>>(SD).eps.noalias() = B * u;
 
        assembleWithJacobianEvalConstitutiveSetting(
            t, dt, x_position, ip_data_[ip], variables, variables_prev, medium,
            TemperatureData{temperature},
            CapillaryPressureData<DisplacementDim>{
                p_cap_ip, p_cap_prev_ip,
                Eigen::Vector<double, DisplacementDim>::Zero()},
            CD, SD, SD_prev, this->process_data_.micro_porosity_parameters,
            this->solid_material_, this->material_states_[ip]);
 
        {
            auto const& C = *std::get<StiffnessTensor<DisplacementDim>>(CD);
            local_Jac
                .template block<displacement_size, displacement_size>(
                    displacement_index, displacement_index)
                .noalias() += B.transpose() * C * B * w;
        }
 
        auto const& b = this->process_data_.specific_body_force;
 
        {
            auto const& sigma_eff =
                std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                    DisplacementDim>>(this->current_states_[ip])
                    .sigma_eff;
            double const rho = *std::get<Density>(CD);
            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
        //
 
        double const alpha =
            *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD);
        double const dS_L_dp_cap =
            std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(
                CD)
                .dS_L_dp_cap;
 
        {
            double const chi_S_L =
                std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
                    .chi_S_L;
            Kup.noalias() +=
                B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
            double const dchi_dS_L =
                std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
                    .dchi_dS_L;
 
            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;
        }
 
        double const phi =
            std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi;
        double const rho_LR = *std::get<LiquidDensity>(CD);
        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.
        //
        double const S_L =
            std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
                this->current_states_[ip])
                .S_L;
        if (this->process_data_.explicit_hm_coupling_in_unsaturated_zone)
        {
            double const chi_S_L_prev = std::get<PrevState<
                ProcessLib::ThermoRichardsMechanics::BishopsData>>(CD)
                                            ->chi_S_L;
            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.
        //
 
        double const k_rel =
            std::get<ProcessLib::ThermoRichardsMechanics::PermeabilityData<
                DisplacementDim>>(CD)
                .k_rel;
        auto const& K_intrinsic =
            std::get<ProcessLib::ThermoRichardsMechanics::PermeabilityData<
                DisplacementDim>>(CD)
                .Ki;
        double const mu =
            *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(
                CD);
 
        GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu;
 
        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 beta_SR =
            std::get<
                ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData>(
                CD)
                .beta_SR;
        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;
 
        double const DeltaS_L_Deltap_cap =
            std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap;
        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;
 
        double const S_L_prev =
            std::get<
                PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
                this->prev_states_[ip])
                ->S_L;
        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 (!this->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 =
                this->process_data_.micro_porosity_parameters
                    ->mass_exchange_coefficient;
            auto const p_L_m =
                *std::get<MicroPressure>(this->current_states_[ip]);
            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)
            {
                auto const p_L_m_prev = **std::get<PrevState<MicroPressure>>(
                    this->prev_states_[ip]);
                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 (this->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::AqueousLiquid, assembleWithJacobianEvalConstitutiveSetting(), MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), ProcessLib::RichardsMechanics::createConstitutiveModels(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::current_states_, MaterialPropertyLib::density, displacement_index, displacement_size, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::element_, MaterialPropertyLib::VariableArray::gas_phase_pressure, ParameterLib::SpatialPosition::getCoordinates(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), ip_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::is_axially_symmetric_, MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::liquid_saturation, localDOF(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::material_states_, N_u_op, pressure_index, pressure_size, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::prev_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::process_data_, MaterialPropertyLib::reference_temperature, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation_micro, NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::Solid, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::solid_material_, MaterialPropertyLib::swelling_stress_rate, MaterialPropertyLib::VariableArray::temperature, and MathLib::KelvinVector::tensorToKelvin().

assembleWithJacobianEvalConstitutiveSetting()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianEvalConstitutiveSetting ( double const t, double const dt, ParameterLib::SpatialPosition const & x_position, IpData & ip_data, MPL::VariableArray & variables, MPL::VariableArray & variables_prev, MPL::Medium const *const medium, TemperatureData const T_data, CapillaryPressureData< DisplacementDim > const & p_cap_data, ConstitutiveData< DisplacementDim > & CD, StatefulData< DisplacementDim > & SD, StatefulDataPrev< DisplacementDim > const & SD_prev, std::optional< MicroPorosityParameters > const & micro_porosity_parameters, MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material, ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > & material_state_data )

staticprivate

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

{
    auto const& liquid_phase =
        medium->phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
    auto const& solid_phase =
        medium->phase(MaterialPropertyLib::PhaseName::Solid);
 
    auto const& identity2 = MathLib::KelvinVector::Invariants<
        MathLib::KelvinVector::kelvin_vector_dimensions(
            DisplacementDim)>::identity2;
 
    double const temperature = T_data();
    double const p_cap_ip = p_cap_data.p_cap;
    double const p_cap_prev_ip = p_cap_data.p_cap_prev;
 
    auto const& eps = std::get<StrainData<DisplacementDim>>(SD);
    auto& S_L =
        std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(SD).S_L;
    auto const S_L_prev =
        std::get<
            PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
            SD_prev)
            ->S_L;
    auto const alpha =
        medium->property(MPL::PropertyType::biot_coefficient)
            .template value<double>(variables, x_position, t, dt);
    *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD) = alpha;
 
    variables.stress =
        std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
            DisplacementDim>>(SD)
            .sigma_eff;
    // Set mechanical strain temporary to compute tangent stiffness.
    variables.mechanical_strain
        .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
            eps.eps);
    auto const C_el = ip_data.computeElasticTangentStiffness(
        variables, t, x_position, dt, solid_material,
        *material_state_data.material_state_variables);
 
    auto const beta_SR =
        (1 - alpha) / solid_material.getBulkModulus(t, x_position, &C_el);
    variables.grain_compressibility = beta_SR;
    std::get<ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData>(CD)
        .beta_SR = beta_SR;
 
    auto const rho_LR =
        liquid_phase.property(MPL::PropertyType::density)
            .template value<double>(variables, x_position, t, dt);
    variables.density = rho_LR;
    *std::get<LiquidDensity>(CD) = rho_LR;
 
    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);
    std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(CD)
        .dS_L_dp_cap = dS_L_dp_cap;
    // 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);
    std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap =
        DeltaS_L_Deltap_cap;
 
    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);
    std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD).chi_S_L =
        chi_S_L;
    double const chi_S_L_prev = chi(S_L_prev);
    std::get<PrevState<ProcessLib::ThermoRichardsMechanics::BishopsData>>(CD)
        ->chi_S_L = chi_S_L_prev;
 
    auto const dchi_dS_L =
        medium->property(MPL::PropertyType::bishops_effective_stress)
            .template dValue<double>(
                variables, MPL::Variable::liquid_saturation, x_position, t, dt);
    std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD).dchi_dS_L =
        dchi_dS_L;
 
    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(eps.eps);
    // TODO (CL) changed that, using eps_prev for the moment, not B * u_prev
    // variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
    variables_prev.volumetric_strain = Invariants::trace(
        std::get<PrevState<StrainData<DisplacementDim>>>(SD_prev)->eps);
 
    auto& phi =
        std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(SD).phi;
    {  // Porosity update
        auto const phi_prev =
            std::get<
                PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData>>(
                SD_prev)
                ->phi;
        variables_prev.porosity = phi_prev;
        phi = medium->property(MPL::PropertyType::porosity)
                  .template value<double>(variables, variables_prev, x_position,
                                          t, dt);
        variables.porosity = phi;
    }
    std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi = phi;
 
    if (alpha < phi)
    {
        auto const eid =
            x_position.getElementID()
                ? static_cast<std::ptrdiff_t>(*x_position.getElementID())
                : static_cast<std::ptrdiff_t>(-1);
        OGS_FATAL(
            "RichardsMechanics: Biot-coefficient {} is smaller than porosity "
            "{} in element {}.",
            alpha, phi, eid);
    }
 
    auto const mu = liquid_phase.property(MPL::PropertyType::viscosity)
                        .template value<double>(variables, x_position, t, dt);
    *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(CD) =
        mu;
 
    {
        // Swelling and possibly volumetric strain rate update.
        auto& sigma_sw =
            std::get<ProcessLib::ThermoRichardsMechanics::
                         ConstitutiveStress_StrainTemperature::
                             SwellingDataStateful<DisplacementDim>>(SD);
        auto const& sigma_sw_prev =
            std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
                                   ConstitutiveStress_StrainTemperature::
                                       SwellingDataStateful<DisplacementDim>>>(
                SD_prev);
        auto const transport_porosity_prev = std::get<PrevState<
            ProcessLib::ThermoRichardsMechanics::TransportPorosityData>>(
            SD_prev);
        auto const phi_prev = std::get<
            PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData>>(
            SD_prev);
        auto& transport_porosity = std::get<
            ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(SD);
        auto& p_L_m = std::get<MicroPressure>(SD);
        auto const p_L_m_prev = std::get<PrevState<MicroPressure>>(SD_prev);
        auto& S_L_m = std::get<MicroSaturation>(SD);
        auto const S_L_m_prev = std::get<PrevState<MicroSaturation>>(SD_prev);
 
        updateSwellingStressAndVolumetricStrain<DisplacementDim>(
            *medium, solid_phase, C_el, rho_LR, mu, micro_porosity_parameters,
            alpha, phi, p_cap_ip, variables, variables_prev, x_position, t, dt,
            sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev,
            transport_porosity, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m);
    }
 
    if (medium->hasProperty(MPL::PropertyType::transport_porosity))
    {
        if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
        {
            auto& transport_porosity =
                std::get<
                    ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(
                    SD)
                    .phi;
            auto const transport_porosity_prev = std::get<PrevState<
                ProcessLib::ThermoRichardsMechanics::TransportPorosityData>>(
                                                     SD_prev)
                                                     ->phi;
            variables_prev.transport_porosity = transport_porosity_prev;
 
            transport_porosity =
                medium->property(MPL::PropertyType::transport_porosity)
                    .template value<double>(variables, variables_prev,
                                            x_position, t, dt);
            variables.transport_porosity = transport_porosity;
        }
    }
    else
    {
        variables.transport_porosity = phi;
    }
 
    // Set mechanical variables for the intrinsic permeability model
    // For stress dependent permeability.
    {
        // TODO mechanical constitutive relation will be evaluated afterwards
        auto const sigma_total =
            (std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                 DisplacementDim>>(SD)
                 .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 =
        material_state_data.material_state_variables
            ->getEquivalentPlasticStrain();
 
    double const k_rel =
        medium->property(MPL::PropertyType::relative_permeability)
            .template value<double>(variables, x_position, t, dt);
 
    auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
        medium->property(MPL::PropertyType::permeability)
            .value(variables, x_position, t, dt));
 
    std::get<
        ProcessLib::ThermoRichardsMechanics::PermeabilityData<DisplacementDim>>(
        CD)
        .k_rel = k_rel;
    std::get<
        ProcessLib::ThermoRichardsMechanics::PermeabilityData<DisplacementDim>>(
        CD)
        .Ki = K_intrinsic;
 
    //
    // displacement equation, displacement part
    //
 
    {
        auto& sigma_sw =
            std::get<ProcessLib::ThermoRichardsMechanics::
                         ConstitutiveStress_StrainTemperature::
                             SwellingDataStateful<DisplacementDim>>(SD)
                .sigma_sw;
 
        auto& eps_m =
            std::get<ProcessLib::ConstitutiveRelations::MechanicalStrainData<
                DisplacementDim>>(SD)
                .eps_m;
        eps_m.noalias() =
            solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
                ? eps.eps + C_el.inverse() * sigma_sw
                : eps.eps;
        variables.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                eps_m);
    }
 
    {
        auto& sigma_eff =
            std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                DisplacementDim>>(SD);
        auto const& sigma_eff_prev =
            std::get<PrevState<ProcessLib::ConstitutiveRelations::
                                   EffectiveStressData<DisplacementDim>>>(
                SD_prev);
        auto const& eps_m =
            std::get<ProcessLib::ConstitutiveRelations::MechanicalStrainData<
                DisplacementDim>>(SD);
        auto& eps_m_prev =
            std::get<PrevState<ProcessLib::ConstitutiveRelations::
                                   MechanicalStrainData<DisplacementDim>>>(
                SD_prev);
 
        auto C = ip_data.updateConstitutiveRelation(
            variables, t, x_position, dt, temperature, sigma_eff,
            sigma_eff_prev, eps_m, eps_m_prev, solid_material,
            material_state_data.material_state_variables);
 
        *std::get<StiffnessTensor<DisplacementDim>>(CD) = std::move(C);
    }
 
    // p_SR
    variables.solid_grain_pressure =
        p_FR - std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                   DisplacementDim>>(SD)
                       .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;
    *std::get<Density>(CD) = rho;
}

References RichardsMechanicsLocalAssembler(), MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::capillary_pressure, MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialLib::Solids::MechanicsBase< DisplacementDim >::getBulkModulus(), ParameterLib::SpatialPosition::getElementID(), MaterialPropertyLib::VariableArray::grain_compressibility, MaterialPropertyLib::Medium::hasProperty(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::liquid_saturation, MaterialPropertyLib::VariableArray::liquid_saturation, ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim >::material_state_variables, MaterialPropertyLib::VariableArray::mechanical_strain, OGS_FATAL, ProcessLib::RichardsMechanics::CapillaryPressureData< DisplacementDim >::p_cap, ProcessLib::RichardsMechanics::CapillaryPressureData< DisplacementDim >::p_cap_prev, MaterialPropertyLib::permeability, MaterialPropertyLib::Medium::phase(), MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, MaterialPropertyLib::Medium::property(), MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, MaterialPropertyLib::saturation_micro, ProcessLib::ConstitutiveRelations::EffectiveStressData< DisplacementDim >::sigma_eff, MaterialPropertyLib::Solid, MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::stress, MaterialPropertyLib::swelling_stress_rate, MaterialPropertyLib::VariableArray::total_stress, MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace(), MaterialPropertyLib::transport_porosity, MaterialPropertyLib::VariableArray::transport_porosity, ProcessLib::RichardsMechanics::updateSwellingStressAndVolumetricStrain(), MaterialPropertyLib::Property::value(), MaterialPropertyLib::viscosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

Referenced by assembleWithJacobian().

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_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_b_data	Right hand side vector of an element.
local_Jac_data	Element Jacobian matrix for the Newton-Raphson method.

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

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

References OGS_FATAL.

Referenced by assembleWithJacobianForStaggeredScheme().

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_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_b_data	Right hand side vector of an element.
local_Jac_data	Element Jacobian matrix for the Newton-Raphson method.

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

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

References OGS_FATAL.

Referenced by assembleWithJacobianForStaggeredScheme().

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_b_data, std::vector< double > & local_Jac_data )

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

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

References assembleWithJacobianForDeformationEquations(), and assembleWithJacobianForPressureEquations().

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 1549 of file RichardsMechanicsFEM-impl.h.

{
    auto const [p_L, u] = localDOF(local_x);
    auto const [p_L_prev, u_prev] = localDOF(local_x_prev);
 
    auto const& identity2 = MathLib::KelvinVector::Invariants<
        MathLib::KelvinVector::kelvin_vector_dimensions(
            DisplacementDim)>::identity2;
 
    auto const& medium =
        this->process_data_.media_map.getMedium(this->element_.getID());
    auto const& liquid_phase =
        medium->phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
    auto const& solid_phase =
        medium->phase(MaterialPropertyLib::PhaseName::Solid);
    MPL::VariableArray variables;
    MPL::VariableArray variables_prev;
 
    unsigned const n_integration_points =
        this->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++)
    {
        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;
 
        ParameterLib::SpatialPosition x_position = {
            std::nullopt, this->element_.getID(),
            MathLib::Point3d(
                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
                                               ShapeMatricesTypeDisplacement>(
                    this->element_, N_u))};
        auto const x_coord = x_position.getCoordinates().value()[0];
 
        auto const B =
            LinearBMatrix::computeBMatrix<DisplacementDim,
                                          ShapeFunctionDisplacement::NPOINTS,
                                          typename BMatricesType::BMatrixType>(
                dNdx_u, N_u, x_coord, this->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 =
            std::get<StrainData<DisplacementDim>>(this->current_states_[ip])
                .eps;
        eps.noalias() = B * u;
        auto& S_L =
            std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
                this->current_states_[ip])
                .S_L;
        auto const S_L_prev =
            std::get<
                PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
                this->prev_states_[ip])
                ->S_L;
        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& SD = this->current_states_[ip];
        variables.stress =
            std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                DisplacementDim>>(SD)
                .sigma_eff;
        // Set mechanical strain temporary to compute tangent stiffness.
        variables.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                eps);
        auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
            variables, t, x_position, dt, this->solid_material_,
            *this->material_states_[ip].material_state_variables);
 
        auto const beta_SR = (1 - alpha) / this->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(eps);
        variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
 
        auto& phi = std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(
                        this->current_states_[ip])
                        .phi;
        {  // Porosity update
            auto const phi_prev = std::get<PrevState<
                ProcessLib::ThermoRichardsMechanics::PorosityData>>(
                                      this->prev_states_[ip])
                                      ->phi;
            variables_prev.porosity = phi_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.
            auto& sigma_sw =
                std::get<ProcessLib::ThermoRichardsMechanics::
                             ConstitutiveStress_StrainTemperature::
                                 SwellingDataStateful<DisplacementDim>>(
                    this->current_states_[ip]);
            auto const& sigma_sw_prev = std::get<
                PrevState<ProcessLib::ThermoRichardsMechanics::
                              ConstitutiveStress_StrainTemperature::
                                  SwellingDataStateful<DisplacementDim>>>(
                this->prev_states_[ip]);
            auto const transport_porosity_prev = std::get<PrevState<
                ProcessLib::ThermoRichardsMechanics::TransportPorosityData>>(
                this->prev_states_[ip]);
            auto const phi_prev = std::get<
                PrevState<ProcessLib::ThermoRichardsMechanics::PorosityData>>(
                this->prev_states_[ip]);
            auto& transport_porosity = std::get<
                ProcessLib::ThermoRichardsMechanics::TransportPorosityData>(
                this->current_states_[ip]);
            auto& p_L_m = std::get<MicroPressure>(this->current_states_[ip]);
            auto const p_L_m_prev =
                std::get<PrevState<MicroPressure>>(this->prev_states_[ip]);
            auto& S_L_m = std::get<MicroSaturation>(this->current_states_[ip]);
            auto const S_L_m_prev =
                std::get<PrevState<MicroSaturation>>(this->prev_states_[ip]);
 
            updateSwellingStressAndVolumetricStrain<DisplacementDim>(
                *medium, solid_phase, C_el, rho_LR, mu,
                this->process_data_.micro_porosity_parameters, alpha, phi,
                p_cap_ip, variables, variables_prev, x_position, t, dt,
                sigma_sw, sigma_sw_prev, transport_porosity_prev, phi_prev,
                transport_porosity, p_L_m_prev, S_L_m_prev, p_L_m, S_L_m);
        }
 
        if (medium->hasProperty(MPL::PropertyType::transport_porosity))
        {
            if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
            {
                auto& transport_porosity =
                    std::get<ProcessLib::ThermoRichardsMechanics::
                                 TransportPorosityData>(
                        this->current_states_[ip])
                        .phi;
                auto const transport_porosity_prev =
                    std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
                                           TransportPorosityData>>(
                        this->prev_states_[ip])
                        ->phi;
 
                variables_prev.transport_porosity = transport_porosity_prev;
 
                transport_porosity =
                    medium->property(MPL::PropertyType::transport_porosity)
                        .template value<double>(variables, variables_prev,
                                                x_position, t, dt);
                variables.transport_porosity = transport_porosity;
            }
        }
        else
        {
            variables.transport_porosity = phi;
        }
 
        auto const& sigma_eff =
            std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                DisplacementDim>>(this->current_states_[ip])
                .sigma_eff;
 
        // Set mechanical variables for the intrinsic permeability model
        // For stress dependent permeability.
        {
            auto const sigma_total =
                (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 =
            this->material_states_[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;
 
        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);
        *std::get<DrySolidDensity>(this->output_data_[ip]) = (1 - phi) * rho_SR;
 
        {
            auto& SD = this->current_states_[ip];
            auto const& sigma_sw =
                std::get<ProcessLib::ThermoRichardsMechanics::
                             ConstitutiveStress_StrainTemperature::
                                 SwellingDataStateful<DisplacementDim>>(SD)
                    .sigma_sw;
            auto& eps_m =
                std::get<ProcessLib::ConstitutiveRelations::
                             MechanicalStrainData<DisplacementDim>>(SD)
                    .eps_m;
            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& SD = this->current_states_[ip];
            auto const& SD_prev = this->prev_states_[ip];
            auto& sigma_eff =
                std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                    DisplacementDim>>(SD);
            auto const& sigma_eff_prev =
                std::get<PrevState<ProcessLib::ConstitutiveRelations::
                                       EffectiveStressData<DisplacementDim>>>(
                    SD_prev);
            auto const& eps_m =
                std::get<ProcessLib::ConstitutiveRelations::
                             MechanicalStrainData<DisplacementDim>>(SD);
            auto const& eps_m_prev =
                std::get<PrevState<ProcessLib::ConstitutiveRelations::
                                       MechanicalStrainData<DisplacementDim>>>(
                    SD_prev);
 
            ip_data_[ip].updateConstitutiveRelation(
                variables, t, x_position, dt, temperature, sigma_eff,
                sigma_eff_prev, eps_m, eps_m_prev, this->solid_material_,
                this->material_states_[ip].material_state_variables);
        }
 
        auto const& b = this->process_data_.specific_body_force;
 
        // Compute the velocity
        auto const& dNdx_p = ip_data_[ip].dNdx_p;
        std::get<
            ProcessLib::ThermoRichardsMechanics::DarcyLawData<DisplacementDim>>(
            this->output_data_[ip])
            ->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;
 
    (*this->process_data_.element_saturation)[this->element_.getID()] =
        saturation_avg;
    (*this->process_data_.element_porosity)[this->element_.getID()] =
        porosity_avg;
 
    Eigen::Map<KV>(
        &(*this->process_data_.element_stresses)[this->element_.getID() *
                                                 KV::RowsAtCompileTime]) =
        MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_avg);
 
    NumLib::interpolateToHigherOrderNodes<
        ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
        DisplacementDim>(this->element_, this->is_axially_symmetric_, p_L,
                         *this->process_data_.pressure_interpolated);
}

References MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::current_states_, MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::element_, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, ParameterLib::SpatialPosition::getCoordinates(), MaterialPropertyLib::VariableArray::grain_compressibility, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), NumLib::interpolateToHigherOrderNodes(), ip_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::is_axially_symmetric_, MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, localDOF(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::material_states_, MaterialPropertyLib::VariableArray::mechanical_strain, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::output_data_, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::prev_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::process_data_, MaterialPropertyLib::reference_temperature, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, MaterialPropertyLib::saturation_micro, NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::Solid, MaterialPropertyLib::VariableArray::solid_grain_pressure, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::solid_material_, MaterialPropertyLib::VariableArray::stress, MaterialPropertyLib::swelling_stress_rate, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace(), MaterialPropertyLib::transport_porosity, MaterialPropertyLib::VariableArray::transport_porosity, ProcessLib::RichardsMechanics::updateSwellingStressAndVolumetricStrain(), MaterialPropertyLib::viscosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

getNumberOfVectorElementsForDeformation()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

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

inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 167 of file RichardsMechanicsFEM.h.

    {
        return displacement_size;
    }

References displacement_size.

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 158 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 secondary_data_.

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 110 of file RichardsMechanicsFEM.h.

    {
        unsigned const n_integration_points =
            this->integration_method_.getNumberOfPoints();
 
        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            auto& SD = this->current_states_[ip];
            auto& ip_data = ip_data_[ip];
 
            ParameterLib::SpatialPosition const x_position{
                std::nullopt, this->element_.getID(),
                MathLib::Point3d(NumLib::interpolateCoordinates<
                                 ShapeFunctionDisplacement,
                                 ShapeMatricesTypeDisplacement>(this->element_,
                                                                ip_data.N_u))};

            if (this->process_data_.initial_stress.value)
            {
                std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                    DisplacementDim>>(SD)
                    .sigma_eff =
                    MathLib::KelvinVector::symmetricTensorToKelvinVector<
                        DisplacementDim>(
                        // The data in process_data_.initial_stress.value can
                        // be total stress or effective stress.
                        (*this->process_data_.initial_stress.value)(
                            std::numeric_limits<
                                double>::quiet_NaN() /* time independent */,
                            x_position));
            }
 
            double const t = 0;  // TODO (naumov) pass t from top
            this->solid_material_.initializeInternalStateVariables(
                t, x_position,
                *this->material_states_[ip].material_state_variables);
 
            this->material_states_[ip].pushBackState();
 
            this->prev_states_[ip] = SD;
        }
    }

References ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::current_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::element_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), ip_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::material_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::prev_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::process_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::solid_material_, and MathLib::KelvinVector::symmetricTensorToKelvinVector().

localDOF()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

constexpr auto ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::localDOF ( auto const & x )

inlinestaticconstexprprivate

Definition at line 238 of file RichardsMechanicsFEM.h.

    {
        return NumLib::localDOF<
            ShapeFunctionPressure,
            NumLib::Vectorial<ShapeFunctionDisplacement, DisplacementDim>>(x);
    }

References NumLib::localDOF().

Referenced by assemble(), assembleWithJacobian(), computeSecondaryVariableConcrete(), and setInitialConditionsConcrete().

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 207 of file RichardsMechanicsFEM-impl.h.

{
    assert(local_x.size() == pressure_size + displacement_size);
 
    auto const [p_L, u] = localDOF(local_x);
 
    constexpr double dt = std::numeric_limits<double>::quiet_NaN();
    auto const& medium =
        this->process_data_.media_map.getMedium(this->element_.getID());
    MPL::VariableArray variables;
 
    auto const& solid_phase =
        medium->phase(MaterialPropertyLib::PhaseName::Solid);
 
    auto const& identity2 = MathLib::KelvinVector::Invariants<
        MathLib::KelvinVector::kelvin_vector_dimensions(
            DisplacementDim)>::identity2;
 
    unsigned const n_integration_points =
        this->integration_method_.getNumberOfPoints();
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        auto const& N_p = ip_data_[ip].N_p;
 
        ParameterLib::SpatialPosition x_position = {
            std::nullopt, this->element_.getID(),
            MathLib::Point3d(
                NumLib::interpolateCoordinates<ShapeFunctionPressure,
                                               ShapeMatricesTypePressure>(
                    this->element_, 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;
 
        {
            auto& p_L_m = std::get<MicroPressure>(this->current_states_[ip]);
            auto& p_L_m_prev =
                std::get<PrevState<MicroPressure>>(this->prev_states_[ip]);
            **p_L_m_prev = -p_cap_ip;
            *p_L_m = -p_cap_ip;
        }
 
        auto const temperature =
            medium->property(MPL::PropertyType::reference_temperature)
                .template value<double>(variables, x_position, t, dt);
        variables.temperature = temperature;
 
        auto& S_L_prev =
            std::get<
                PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
                this->prev_states_[ip])
                ->S_L;
        S_L_prev = medium->property(MPL::PropertyType::saturation)
                       .template value<double>(variables, x_position, t, dt);
 
        if (this->process_data_.initial_stress.isTotalStress())
        {
            auto const alpha_b =
                medium->property(MPL::PropertyType::biot_coefficient)
                    .template value<double>(variables, x_position, t, dt);
 
            variables.liquid_saturation = S_L_prev;
            double const chi_S_L =
                medium->property(MPL::PropertyType::bishops_effective_stress)
                    .template value<double>(variables, x_position, t, dt);
 
            // Initial stresses are total stress, which were assigned to
            // sigma_eff in
            // RichardsMechanicsLocalAssembler::initializeConcrete().
            auto& sigma_eff =
                std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                    DisplacementDim>>(this->current_states_[ip]);
 
            auto& sigma_eff_prev =
                std::get<PrevState<ProcessLib::ConstitutiveRelations::
                                       EffectiveStressData<DisplacementDim>>>(
                    this->prev_states_[ip]);
 
            // Reset sigma_eff to effective stress
            sigma_eff.sigma_eff.noalias() +=
                chi_S_L * alpha_b * (-p_cap_ip) * identity2;
            sigma_eff_prev->sigma_eff = sigma_eff.sigma_eff;
        }
 
        if (medium->hasProperty(MPL::PropertyType::saturation_micro))
        {
            MPL::VariableArray vars;
            vars.capillary_pressure = p_cap_ip;
 
            auto& S_L_m = std::get<MicroSaturation>(this->current_states_[ip]);
            auto& S_L_m_prev =
                std::get<PrevState<MicroSaturation>>(this->prev_states_[ip]);
 
            *S_L_m = medium->property(MPL::PropertyType::saturation_micro)
                         .template value<double>(vars, x_position, t, dt);
            *S_L_m_prev = S_L_m;
        }
 
        // Set eps_m_prev from potentially non-zero eps and sigma_sw from
        // restart.
        auto& SD = this->current_states_[ip];
        variables.stress =
            std::get<ProcessLib::ConstitutiveRelations::EffectiveStressData<
                DisplacementDim>>(SD)
                .sigma_eff;
 
        auto const& N_u = ip_data_[ip].N_u;
        auto const& dNdx_u = ip_data_[ip].dNdx_u;
        auto const x_coord =
            x_position.getCoordinates().value()[0];  // r for axisymetric
        auto const B =
            LinearBMatrix::computeBMatrix<DisplacementDim,
                                          ShapeFunctionDisplacement::NPOINTS,
                                          typename BMatricesType::BMatrixType>(
                dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
        auto& eps =
            std::get<StrainData<DisplacementDim>>(this->current_states_[ip])
                .eps;
        eps.noalias() = B * u;
 
        // Set mechanical strain temporary to compute tangent stiffness.
        variables.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                eps);
 
        // dt = 0 at initialization: there is no time step yet, which yields
        // the elastic tangent. The function-wide dt is NaN to keep
        // initialization and integration strictly separated.
        auto const C_el = ip_data_[ip].computeElasticTangentStiffness(
            variables, t, x_position, 0.0 /*dt*/, this->solid_material_,
            *this->material_states_[ip].material_state_variables);
 
        auto const& sigma_sw =
            std::get<ProcessLib::ThermoRichardsMechanics::
                         ConstitutiveStress_StrainTemperature::
                             SwellingDataStateful<DisplacementDim>>(
                this->current_states_[ip])
                .sigma_sw;
        auto& eps_m_prev =
            std::get<PrevState<ProcessLib::ConstitutiveRelations::
                                   MechanicalStrainData<DisplacementDim>>>(
                this->prev_states_[ip])
                ->eps_m;
 
        eps_m_prev.noalias() =
            solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
                ? eps + C_el.inverse() * sigma_sw
                : eps;
    }
}

References MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::current_states_, displacement_size, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::element_, MaterialPropertyLib::VariableArray::gas_phase_pressure, ParameterLib::SpatialPosition::getCoordinates(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), ip_data_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::is_axially_symmetric_, MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, localDOF(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::material_states_, MaterialPropertyLib::VariableArray::mechanical_strain, pressure_size, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::prev_states_, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::process_data_, MaterialPropertyLib::reference_temperature, MaterialPropertyLib::saturation, MaterialPropertyLib::saturation_micro, NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::Solid, ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::solid_material_, MaterialPropertyLib::VariableArray::stress, MaterialPropertyLib::swelling_stress_rate, and MaterialPropertyLib::VariableArray::temperature.

Member Data Documentation

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 253 of file RichardsMechanicsFEM.h.

Referenced by assemble(), and assembleWithJacobian().

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 254 of file RichardsMechanicsFEM.h.

Referenced by assemble(), assembleWithJacobian(), getNumberOfVectorElementsForDeformation(), and setInitialConditionsConcrete().

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 245 of file RichardsMechanicsFEM.h.

Referenced by assemble(), assembleWithJacobian(), computeSecondaryVariableConcrete(), initializeConcrete(), and setInitialConditionsConcrete().

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 66 of file RichardsMechanicsFEM.h.

N_u_op

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

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

staticconstexpr

Initial value:

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

Definition at line 72 of file RichardsMechanicsFEM.h.

Referenced by assemble(), and assembleWithJacobian().

pressure_index

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

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

staticprivate

Definition at line 251 of file RichardsMechanicsFEM.h.

Referenced by assemble(), and assembleWithJacobian().

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 252 of file RichardsMechanicsFEM.h.

Referenced by assemble(), assembleWithJacobian(), and setInitialConditionsConcrete().

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 249 of file RichardsMechanicsFEM.h.

Referenced by getShapeMatrix().

---

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

• RichardsMechanicsFEM.h
• RichardsMechanicsFEM-impl.h