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 52 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.
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.
virtual std::optional< VectorSegment >	getVectorDeformationSegment () const
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 64 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 61 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 75 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 68 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 66 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 56 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 58 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 77 of file RichardsMechanicsFEM.h.

Constructor & Destructor Documentation

RichardsMechanicsLocalAssembler() [1/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

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

delete

RichardsMechanicsLocalAssembler() [2/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

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

delete

RichardsMechanicsLocalAssembler() [3/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

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

Definition at line 129 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;
    }
}

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 365 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("AqueousLiquid");
    auto const& solid_phase = medium->phase("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 const C_el = ip_data_[ip].computeElasticTangentStiffness(
            t, x_position, dt, temperature, 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);
 
            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);
        }
 
        // 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::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::formEigenTensor< 3 >(), MaterialPropertyLib::VariableArray::gas_phase_pressure, ParameterLib::SpatialPosition::getCoordinates(), MaterialPropertyLib::VariableArray::grain_compressibility, NumLib::interpolateCoordinates(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::VariableArray::mechanical_strain, OGS_FATAL, MaterialPropertyLib::VariableArray::porosity, ParameterLib::SpatialPosition::setElementID(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::temperature, MathLib::KelvinVector::tensorToKelvin(), MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::VariableArray::transport_porosity, MaterialPropertyLib::VariableArray::volumetric_mechanical_strain, and MaterialPropertyLib::VariableArray::volumetric_strain.

assembleWithJacobian()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

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

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1067 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("AqueousLiquid");
    auto const& solid_phase = medium->phase("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::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), ProcessLib::RichardsMechanics::createConstitutiveModels(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::VariableArray::gas_phase_pressure, ParameterLib::SpatialPosition::getCoordinates(), NumLib::interpolateCoordinates(), MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, NumLib::detail::shapeFunctionInterpolate(), 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 762 of file RichardsMechanicsFEM-impl.h.

{
    auto const& liquid_phase = medium->phase("AqueousLiquid");
    auto const& solid_phase = medium->phase("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;
 
    auto const C_el = ip_data.computeElasticTangentStiffness(
        t, x_position, dt, temperature, 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 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::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::Medium::phase(), MaterialPropertyLib::VariableArray::porosity, MaterialPropertyLib::Medium::property(), ProcessLib::ConstitutiveRelations::EffectiveStressData< DisplacementDim >::sigma_eff, MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::VariableArray::transport_porosity, ProcessLib::RichardsMechanics::updateSwellingStressAndVolumetricStrain(), MaterialPropertyLib::Property::value(), and MaterialPropertyLib::VariableArray::volumetric_strain.

assembleWithJacobianForDeformationEquations()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForDeformationEquations ( double const t, double const dt, Eigen::VectorXd const & local_x, Eigen::VectorXd const & local_x_prev, std::vector< double > & local_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 1475 of file RichardsMechanicsFEM-impl.h.

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

References OGS_FATAL.

assembleWithJacobianForPressureEquations()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

void ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForPressureEquations ( double const t, double const dt, Eigen::VectorXd const & local_x, Eigen::VectorXd const & local_x_prev, std::vector< double > & local_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 1461 of file RichardsMechanicsFEM-impl.h.

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

References OGS_FATAL.

assembleWithJacobianForStaggeredScheme()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

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

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1489 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);
}

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 1513 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("AqueousLiquid");
    auto const& solid_phase = medium->phase("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 const C_el = ip_data_[ip].computeElasticTangentStiffness(
            t, x_position, dt, temperature, 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::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MaterialPropertyLib::VariableArray::density, MaterialPropertyLib::VariableArray::effective_pore_pressure, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, ParameterLib::SpatialPosition::getCoordinates(), MaterialPropertyLib::VariableArray::grain_compressibility, NumLib::interpolateCoordinates(), NumLib::interpolateToHigherOrderNodes(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::VariableArray::porosity, NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::solid_grain_pressure, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::VariableArray::transport_porosity, ProcessLib::RichardsMechanics::updateSwellingStressAndVolumetricStrain(), and MaterialPropertyLib::VariableArray::volumetric_strain.

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

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

References ProcessLib::RichardsMechanics::SecondaryData< ShapeMatrixType >::N_u, and ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::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 117 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_, MeshLib::Element::getID(), NumLib::GenericIntegrationMethod::getNumberOfPoints(), ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >::integration_method_, NumLib::interpolateCoordinates(), ProcessLib::RichardsMechanics::RichardsMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::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>

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

inlinestaticconstexprprivate

Definition at line 240 of file RichardsMechanicsFEM.h.

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

References NumLib::localDOF().

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 214 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("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 const C_el = ip_data_[ip].computeElasticTangentStiffness(
            t, x_position, dt, temperature, this->solid_material_,
            *this->material_states_[ip].material_state_variables);
 
        auto const& N_u = ip_data_[ip].N_u;
        auto const& dNdx_u = ip_data_[ip].dNdx_u;
        auto const x_coord =
            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
                                           ShapeMatricesTypeDisplacement>(
                this->element_, N_u);
        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 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::VariableArray::capillary_pressure, ProcessLib::LinearBMatrix::computeBMatrix(), MaterialPropertyLib::VariableArray::gas_phase_pressure, NumLib::interpolateCoordinates(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, NumLib::detail::shapeFunctionInterpolate(), 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 255 of file RichardsMechanicsFEM.h.

displacement_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

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

staticprivate

Initial value:

=
        ShapeFunctionDisplacement::NPOINTS * DisplacementDim

Definition at line 256 of file RichardsMechanicsFEM.h.

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

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

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

pressure_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

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

staticprivate

Definition at line 253 of file RichardsMechanicsFEM.h.

pressure_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>

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

staticprivate

Definition at line 254 of file RichardsMechanicsFEM.h.

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

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

---

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

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