Detailed Description

template<typename ShapeFunction, int GlobalDim>
class ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >

Definition at line 36 of file ThermoRichardsFlowFEM.h.

#include <ThermoRichardsFlowFEM.h>

Inheritance diagram for ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >:

[legend]

Collaboration diagram for ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >:

[legend]

Public Types
using	ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>

using	GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType

using	GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType

using	IpData = IntegrationPointData<ShapeMatricesType>

Public Member Functions
	ThermoRichardsFlowLocalAssembler (ThermoRichardsFlowLocalAssembler const &)=delete

	ThermoRichardsFlowLocalAssembler (ThermoRichardsFlowLocalAssembler &&)=delete

	ThermoRichardsFlowLocalAssembler (MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, ThermoRichardsFlowProcessData &process_data)

std::size_t	setIPDataInitialConditions (std::string_view const name, double const *values, int const integration_order) override

void	setInitialConditionsConcrete (Eigen::VectorXd const local_x, double const t, int const process_id) override

void	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	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	initializeConcrete () override

void	postTimestepConcrete (Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override

void	computeSecondaryVariableConcrete (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev) override

Eigen::Map< const Eigen::RowVectorXd >	getShapeMatrix (const unsigned integration_point) const override
	Provides the shape matrix at the given integration point.

std::vector< double > const &	getIntPtDarcyVelocity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

std::vector< double >	getSaturation () const override

std::vector< double > const &	getIntPtSaturation (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

std::vector< double >	getPorosity () const override

std::vector< double > const &	getIntPtPorosity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

std::vector< double > const &	getIntPtDryDensitySolid (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override

Public Member Functions inherited from ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface
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	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)

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

Private Member Functions
unsigned	getNumberOfIntegrationPoints () const override

Private Attributes
ThermoRichardsFlowProcessData &	_process_data

std::vector< IpData, Eigen::aligned_allocator< IpData > >	_ip_data

NumLib::GenericIntegrationMethod const &	_integration_method

MeshLib::Element const &	_element

bool const	_is_axially_symmetric

Static Private Attributes
static const int	temperature_index = 0

static const int	temperature_size = ShapeFunction::NPOINTS

static const int	pressure_index = temperature_size

static const int	pressure_size = ShapeFunction::NPOINTS

Member Typedef Documentation

◆ GlobalDimMatrixType

template<typename ShapeFunction , int GlobalDim>

using ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType

Definition at line 43 of file ThermoRichardsFlowFEM.h.

◆ GlobalDimVectorType

template<typename ShapeFunction , int GlobalDim>

using ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType

Definition at line 44 of file ThermoRichardsFlowFEM.h.

◆ IpData

template<typename ShapeFunction , int GlobalDim>

using ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::IpData = IntegrationPointData<ShapeMatricesType>

Definition at line 46 of file ThermoRichardsFlowFEM.h.

◆ ShapeMatricesType

template<typename ShapeFunction , int GlobalDim>

using ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>

Definition at line 41 of file ThermoRichardsFlowFEM.h.

Constructor & Destructor Documentation

◆ ThermoRichardsFlowLocalAssembler() [1/3]

template<typename ShapeFunction , int GlobalDim>

ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::ThermoRichardsFlowLocalAssembler ( ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim > const & )

delete

◆ ThermoRichardsFlowLocalAssembler() [2/3]

template<typename ShapeFunction , int GlobalDim>

ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::ThermoRichardsFlowLocalAssembler ( ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim > && )

delete

◆ ThermoRichardsFlowLocalAssembler() [3/3]

template<typename ShapeFunction , int GlobalDim>

ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::ThermoRichardsFlowLocalAssembler	(	MeshLib::Element const &	e,
		std::size_t const	,
		NumLib::GenericIntegrationMethod const &	integration_method,
		bool const	is_axially_symmetric,
		ThermoRichardsFlowProcessData &	process_data )

Definition at line 34 of file ThermoRichardsFlowFEM-impl.h.

    : _process_data(process_data),
      _integration_method(integration_method),
      _element(e),
      _is_axially_symmetric(is_axially_symmetric)
{
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
 
    _ip_data.reserve(n_integration_points);
 
    auto const shape_matrices =
        NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType, GlobalDim>(
            e, is_axially_symmetric, _integration_method);
 
    auto const& medium = *_process_data.media_map.getMedium(_element.getID());
 
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        auto const& sm = shape_matrices[ip];
        _ip_data.emplace_back();
        auto& ip_data = _ip_data[ip];
        _ip_data[ip].integration_weight =
            _integration_method.getWeightedPoint(ip).getWeight() *
            sm.integralMeasure * sm.detJ;
 
        ip_data.N = sm.N;
        ip_data.dNdx = sm.dNdx;
 
        ParameterLib::SpatialPosition const x_position{
            std::nullopt, _element.getID(),
            MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
                                                            ShapeMatricesType>(
                _element, sm.N))};
        // Initial porosity. Could be read from integration point data or mesh.
        ip_data.porosity = medium[MPL::porosity].template initialValue<double>(
            x_position,
            std::numeric_limits<double>::quiet_NaN() /* t independent */);
    }
}

Member Function Documentation

◆ assemble()

template<typename ShapeFunction , int GlobalDim>

void ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::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 687 of file ThermoRichardsFlowFEM-impl.h.

{
    auto const local_matrix_dim = pressure_size + temperature_size;
    assert(local_x.size() == local_matrix_dim);
 
    auto const T = Eigen::Map<typename ShapeMatricesType::template VectorType<
        temperature_size> const>(local_x.data() + temperature_index,
                                 temperature_size);
    auto const p_L = Eigen::Map<
        typename ShapeMatricesType::template VectorType<pressure_size> const>(
        local_x.data() + pressure_index, pressure_size);
 
    auto const p_L_prev = Eigen::Map<
        typename ShapeMatricesType::template VectorType<pressure_size> const>(
        local_x_prev.data() + pressure_index, pressure_size);
 
    auto local_K = MathLib::createZeroedMatrix<
        typename ShapeMatricesType::template MatrixType<local_matrix_dim,
                                                        local_matrix_dim>>(
        local_K_data, local_matrix_dim, local_matrix_dim);
 
    auto local_M = MathLib::createZeroedMatrix<
        typename ShapeMatricesType::template MatrixType<local_matrix_dim,
                                                        local_matrix_dim>>(
        local_M_data, local_matrix_dim, local_matrix_dim);
 
    auto local_rhs = MathLib::createZeroedVector<
        typename ShapeMatricesType::template VectorType<local_matrix_dim>>(
        local_rhs_data, local_matrix_dim);
 
    auto const& medium = *_process_data.media_map.getMedium(_element.getID());
    auto const& liquid_phase = medium.phase("AqueousLiquid");
    auto const& solid_phase = medium.phase("Solid");
    MPL::Phase const* gas_phase =
        medium.hasPhase("Gas") ? &medium.phase("Gas") : nullptr;
    MPL::VariableArray variables;
    MPL::VariableArray variables_prev;
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        auto const& w = _ip_data[ip].integration_weight;
 
        auto const& N = _ip_data[ip].N;
        auto const& dNdx = _ip_data[ip].dNdx;
 
        ParameterLib::SpatialPosition const x_position{
            std::nullopt, _element.getID(),
            MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
                                                            ShapeMatricesType>(
                _element, N))};
 
        double T_ip;
        NumLib::shapeFunctionInterpolate(T, N, T_ip);
 
        double p_cap_ip;
        NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
 
        double p_cap_prev_ip;
        NumLib::shapeFunctionInterpolate(-p_L_prev, N, 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;
        variables.temperature = T_ip;
 
        auto& S_L = _ip_data[ip].saturation;
        auto const S_L_prev = _ip_data[ip].saturation_prev;
        auto const alpha =
            medium[MPL::PropertyType::biot_coefficient].template value<double>(
                variables, x_position, t, dt);
 
        auto& solid_elasticity = *_process_data.simplified_elasticity;
        // TODO (buchwaldj)
        // is bulk_modulus good name for bulk modulus of solid skeleton?
        auto const beta_S =
            solid_elasticity.bulkCompressibilityFromYoungsModulus(
                solid_phase, variables, x_position, t, dt);
        auto const beta_SR = (1 - alpha) * beta_S;
        variables.grain_compressibility = beta_SR;
 
        auto const rho_LR =
            liquid_phase[MPL::PropertyType::density].template value<double>(
                variables, x_position, t, dt);
        auto const& b = _process_data.specific_body_force;
 
        double const drho_LR_dp =
            liquid_phase[MPL::PropertyType::density].template dValue<double>(
                variables, MPL::Variable::liquid_phase_pressure, x_position, t,
                dt);
        auto const beta_LR = drho_LR_dp / rho_LR;
 
        S_L = medium[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[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 chi_S_L = S_L;
        auto chi_S_L_prev = S_L_prev;
        if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
        {
            auto const chi = [&medium, x_position, t, dt](double const S_L)
            {
                MPL::VariableArray variables;
                variables.liquid_saturation = S_L;
                return medium[MPL::PropertyType::bishops_effective_stress]
                    .template value<double>(variables, x_position, t, dt);
            };
            chi_S_L = chi(S_L);
            chi_S_L_prev = chi(S_L_prev);
        }
        // TODO (buchwaldj)
        // should solid_grain_pressure or effective_pore_pressure remain?
        // double const p_FR = -chi_S_L * p_cap_ip;
        // variables.solid_grain_pressure = p_FR;
 
        variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
        variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
 
        auto& phi = _ip_data[ip].porosity;
        {  // Porosity update
 
            variables_prev.porosity = _ip_data[ip].porosity_prev;
            phi = medium[MPL::PropertyType::porosity].template value<double>(
                variables, variables_prev, x_position, t, dt);
            variables.porosity = phi;
        }
 
        if (alpha < phi)
        {
            OGS_FATAL(
                "ThermoRichardsFlow: Biot-coefficient {} is smaller than "
                "porosity {} in element/integration point {}/{}.",
                alpha, phi, _element.getID(), ip);
        }
 
        auto const K_pT_thermal_osmosis =
            (solid_phase.hasProperty(
                 MaterialPropertyLib::PropertyType::thermal_osmosis_coefficient)
                 ? MaterialPropertyLib::formEigenTensor<GlobalDim>(
                       solid_phase
                           [MPL::PropertyType::thermal_osmosis_coefficient]
                               .value(variables, x_position, t, dt))
                 : Eigen::MatrixXd::Zero(GlobalDim, GlobalDim));
 
        double const k_rel =
            medium[MPL::PropertyType::relative_permeability]
                .template value<double>(variables, x_position, t, dt);
        auto const mu =
            liquid_phase[MPL::PropertyType::viscosity].template value<double>(
                variables, x_position, t, dt);
 
        auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
            medium[MPL::PropertyType::permeability].value(variables, x_position,
                                                          t, dt));
 
        GlobalDimMatrixType const Ki_over_mu = K_intrinsic / mu;
        GlobalDimMatrixType const rho_Ki_over_mu = rho_LR * Ki_over_mu;
 
        // Consider anisotropic thermal expansion.
        // Read in 3x3 tensor. 2D case also requires expansion coeff. for z-
        // component.
        Eigen::Matrix<double, 3, 3> const
            solid_linear_thermal_expansion_coefficient =
                MaterialPropertyLib::formEigenTensor<3>(
                    solid_phase
                        [MaterialPropertyLib::PropertyType::thermal_expansivity]
                            .value(variables, x_position, t, dt));
 
        auto const rho_SR =
            solid_phase[MPL::PropertyType::density].template value<double>(
                variables, x_position, t, dt);
 
        //
        // pressure equation, pressure part.
        //
        local_K
            .template block<pressure_size, pressure_size>(pressure_index,
                                                          pressure_index)
            .noalias() += dNdx.transpose() * k_rel * rho_Ki_over_mu * dNdx * w;
 
        const double alphaB_minus_phi = alpha - phi;
        double const a0 = alphaB_minus_phi * beta_SR;
        double const specific_storage_a_p =
            S_L * (phi * beta_LR + S_L * a0 +
                   chi_S_L * alpha * alpha *
                       solid_elasticity.storageContribution(
                           solid_phase, variables, x_position, t, dt));
        double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
 
        local_M
            .template block<pressure_size, pressure_size>(pressure_index,
                                                          pressure_index)
            .noalias() += N.transpose() * rho_LR *
                          (specific_storage_a_p -
                           specific_storage_a_S * DeltaS_L_Deltap_cap) *
                          N * w;
 
        local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
            dNdx.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
 
        //
        // pressure equation, temperature part.
        //
        double const fluid_volumetric_thermal_expansion_coefficient =
            MPL::getLiquidThermalExpansivity(liquid_phase, variables, rho_LR,
                                             x_position, t, dt);
        const double eff_thermal_expansion =
            S_L * (alphaB_minus_phi *
                       solid_linear_thermal_expansion_coefficient.trace() +
                   phi * fluid_volumetric_thermal_expansion_coefficient +
                   alpha * solid_elasticity.thermalExpansivityContribution(
                               solid_linear_thermal_expansion_coefficient,
                               solid_phase, variables, x_position, t, dt));
 
        local_K
            .template block<pressure_size, temperature_size>(pressure_index,
                                                             temperature_index)
            .noalias() +=
            dNdx.transpose() * rho_LR * K_pT_thermal_osmosis * dNdx * w;
 
        local_M
            .template block<pressure_size, temperature_size>(pressure_index,
                                                             temperature_index)
            .noalias() -=
            N.transpose() * rho_LR * eff_thermal_expansion * N * w;
 
        //
        // temperature equation.
        //
        {
            auto const specific_heat_capacity_fluid =
                liquid_phase[MaterialPropertyLib::specific_heat_capacity]
                    .template value<double>(variables, x_position, t, dt);
 
            auto const specific_heat_capacity_solid =
                solid_phase
                    [MaterialPropertyLib::PropertyType::specific_heat_capacity]
                        .template value<double>(variables, x_position, t, dt);
 
            local_M
                .template block<temperature_size, temperature_size>(
                    temperature_index, temperature_index)
                .noalias() +=
                w *
                (rho_SR * specific_heat_capacity_solid * (1 - phi) +
                 (S_L * rho_LR * specific_heat_capacity_fluid) * phi) *
                N.transpose() * N;
 
            auto const thermal_conductivity =
                MaterialPropertyLib::formEigenTensor<GlobalDim>(
                    medium[MaterialPropertyLib::PropertyType::
                               thermal_conductivity]
                        .value(variables, x_position, t, dt));
 
            GlobalDimVectorType const velocity_L = GlobalDimVectorType(
                -Ki_over_mu * k_rel * (dNdx * p_L - rho_LR * b) -
                K_pT_thermal_osmosis * dNdx * T);
 
            local_K
                .template block<temperature_size, temperature_size>(
                    temperature_index, temperature_index)
                .noalias() += (dNdx.transpose() * thermal_conductivity * dNdx +
                               N.transpose() * velocity_L.transpose() * dNdx *
                                   rho_LR * specific_heat_capacity_fluid) *
                              w;
            local_K
                .template block<temperature_size, pressure_size>(
                    temperature_index, pressure_index)
                .noalias() +=
                dNdx.transpose() * T_ip * K_pT_thermal_osmosis * dNdx * w;
        }
        if (gas_phase && S_L < 1.0)
        {
            variables.density = rho_LR;
 
            double const rho_wv =
                gas_phase->property(MaterialPropertyLib::density)
                    .template value<double>(variables, x_position, t, dt);
 
            double const drho_wv_dT =
                gas_phase->property(MaterialPropertyLib::density)
                    .template dValue<double>(variables,
                                             MPL::Variable::temperature,
                                             x_position, t, dt);
            double const drho_wv_dp =
                gas_phase->property(MaterialPropertyLib::density)
                    .template dValue<double>(
                        variables, MPL::Variable::liquid_phase_pressure,
                        x_position, t, dt);
            auto const f_Tv =
                gas_phase
                    ->property(
                        MPL::PropertyType::thermal_diffusion_enhancement_factor)
                    .template value<double>(variables, x_position, t, dt);
 
            variables.porosity = phi;
            auto const tortuosity =
                medium.property(MPL::PropertyType::tortuosity)
                    .template value<double>(variables, x_position, t, dt);
            double const D_v =
                phi * (1.0 - S_L) * tortuosity *
                gas_phase->property(MPL::PropertyType::diffusion)
                    .template value<double>(variables, x_position, t, dt);
 
            double const f_Tv_D_Tv = f_Tv * D_v * drho_wv_dT;
            double const D_pv = D_v * drho_wv_dp;
 
            GlobalDimVectorType const grad_T = dNdx * T;
            GlobalDimVectorType const grad_p_cap = -dNdx * p_L;
            GlobalDimVectorType const vapour_flux =
                -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap);
            double const specific_heat_capacity_vapour =
                gas_phase->property(MaterialPropertyLib::specific_heat_capacity)
                    .template value<double>(variables, x_position, t, dt);
 
            local_M
                .template block<temperature_size, temperature_size>(
                    temperature_index, temperature_index)
                .noalias() +=
                w * (rho_wv * specific_heat_capacity_vapour * (1 - S_L) * phi) *
                N.transpose() * N;
 
            local_K
                .template block<temperature_size, temperature_size>(
                    temperature_index, temperature_index)
                .noalias() += N.transpose() * vapour_flux.transpose() * dNdx *
                              specific_heat_capacity_vapour * w;
 
            double const storage_coefficient_by_water_vapor =
                phi * (rho_wv * dS_L_dp_cap + (1 - S_L) * drho_wv_dp);
            local_M
                .template block<pressure_size, pressure_size>(pressure_index,
                                                              pressure_index)
                .noalias() +=
                N.transpose() * storage_coefficient_by_water_vapor * N * w;
 
            double const vapor_expansion_factor = phi * (1 - S_L) * drho_wv_dT;
            local_M
                .template block<pressure_size, temperature_size>(
                    pressure_index, temperature_index)
                .noalias() += N.transpose() * vapor_expansion_factor * N * w;
 
            local_rhs.template segment<pressure_size>(pressure_index)
                .noalias() -= f_Tv_D_Tv * dNdx.transpose() * (dNdx * T) * w;
 
            local_K
                .template block<pressure_size, pressure_size>(pressure_index,
                                                              pressure_index)
                .noalias() += dNdx.transpose() * D_pv * dNdx * w;
 
            //
            // Latent heat term
            //
            if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
            {
                double const factor = phi * (1 - S_L) / rho_LR;
                // The volumetric latent heat of vaporization of liquid water
                double const L0 =
                    gas_phase->property(MPL::PropertyType::specific_latent_heat)
                        .template value<double>(variables, x_position, t, dt) *
                    rho_LR;
 
                double const drho_LR_dT =
                    liquid_phase.property(MPL::PropertyType::density)
                        .template dValue<double>(variables,
                                                 MPL::Variable::temperature,
                                                 x_position, t, dt);
 
                double const rho_wv_over_rho_L = rho_wv / rho_LR;
                local_M
                    .template block<temperature_size, temperature_size>(
                        temperature_index, temperature_index)
                    .noalias() +=
                    factor * L0 *
                    (drho_wv_dT - rho_wv_over_rho_L * drho_LR_dT) *
                    N.transpose() * N * w;
 
                local_M
                    .template block<temperature_size, pressure_size>(
                        temperature_index, pressure_index)
                    .noalias() +=
                    (factor * L0 *
                         (drho_wv_dp - rho_wv_over_rho_L * drho_LR_dp) +
                     L0 * phi * rho_wv_over_rho_L * dS_L_dp_cap) *
                    N.transpose() * N * w;
 
                // temperature equation, temperature part
                local_K
                    .template block<temperature_size, temperature_size>(
                        temperature_index, temperature_index)
                    .noalias() +=
                    L0 * f_Tv_D_Tv * dNdx.transpose() * dNdx * w / rho_LR;
                // temperature equation, pressure part
                local_K
                    .template block<temperature_size, pressure_size>(
                        temperature_index, pressure_index)
                    .noalias() +=
                    L0 * D_pv * dNdx.transpose() * dNdx * w / rho_LR;
            }
        }
    }
 
    if (_process_data.apply_mass_lumping)
    {
        auto Mpp = local_M.template block<pressure_size, pressure_size>(
            pressure_index, pressure_index);
        Mpp = Mpp.colwise().sum().eval().asDiagonal();
    }
}

◆ assembleWithJacobian()

template<typename ShapeFunction , int GlobalDim>

void ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::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 154 of file ThermoRichardsFlowFEM-impl.h.

{
    auto const local_matrix_dim = pressure_size + temperature_size;
    assert(local_x.size() == local_matrix_dim);
 
    auto const T = Eigen::Map<typename ShapeMatricesType::template VectorType<
        temperature_size> const>(local_x.data() + temperature_index,
                                 temperature_size);
    auto const p_L = Eigen::Map<
        typename ShapeMatricesType::template VectorType<pressure_size> const>(
        local_x.data() + pressure_index, pressure_size);
 
    auto const T_prev =
        Eigen::Map<typename ShapeMatricesType::template VectorType<
            temperature_size> const>(local_x_prev.data() + temperature_index,
                                     temperature_size);
    auto const p_L_prev = Eigen::Map<
        typename ShapeMatricesType::template VectorType<pressure_size> const>(
        local_x_prev.data() + pressure_index, pressure_size);
 
    auto local_Jac = MathLib::createZeroedMatrix<
        typename ShapeMatricesType::template MatrixType<local_matrix_dim,
                                                        local_matrix_dim>>(
        local_Jac_data, local_matrix_dim, local_matrix_dim);
 
    auto local_rhs = MathLib::createZeroedVector<
        typename ShapeMatricesType::template VectorType<local_matrix_dim>>(
        local_rhs_data, local_matrix_dim);
 
    typename ShapeMatricesType::NodalMatrixType M_TT =
        ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
                                                 temperature_size);
    typename ShapeMatricesType::NodalMatrixType M_Tp =
        ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
                                                 pressure_size);
    typename ShapeMatricesType::NodalMatrixType K_TT =
        ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
                                                 temperature_size);
    typename ShapeMatricesType::NodalMatrixType K_Tp =
        ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
                                                 pressure_size);
    typename ShapeMatricesType::NodalMatrixType dK_TT_dp =
        ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
                                                 pressure_size);
    typename ShapeMatricesType::NodalMatrixType M_pT =
        ShapeMatricesType::NodalMatrixType::Zero(pressure_size,
                                                 temperature_size);
    typename ShapeMatricesType::NodalMatrixType laplace_p =
        ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
    typename ShapeMatricesType::NodalMatrixType laplace_T =
        ShapeMatricesType::NodalMatrixType::Zero(pressure_size,
                                                 temperature_size);
    typename ShapeMatricesType::NodalMatrixType storage_p_a_p =
        ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
 
    typename ShapeMatricesType::NodalMatrixType storage_p_a_S_Jpp =
        ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
 
    typename ShapeMatricesType::NodalMatrixType storage_p_a_S =
        ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
 
    auto const& medium = *_process_data.media_map.getMedium(_element.getID());
    auto const& liquid_phase = medium.phase("AqueousLiquid");
    auto const& solid_phase = medium.phase("Solid");
    MPL::Phase const* gas_phase =
        medium.hasPhase("Gas") ? &medium.phase("Gas") : nullptr;
    MPL::VariableArray variables;
    MPL::VariableArray variables_prev;
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        auto const& w = _ip_data[ip].integration_weight;
 
        auto const& N = _ip_data[ip].N;
        auto const& dNdx = _ip_data[ip].dNdx;
 
        ParameterLib::SpatialPosition const x_position{
            std::nullopt, _element.getID(),
            MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
                                                            ShapeMatricesType>(
                _element, N))};
 
        double T_ip;
        NumLib::shapeFunctionInterpolate(T, N, T_ip);
 
        double p_cap_ip;
        NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
 
        double p_cap_prev_ip;
        NumLib::shapeFunctionInterpolate(-p_L_prev, N, 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;
        variables.temperature = T_ip;
 
        auto& S_L = _ip_data[ip].saturation;
        auto const S_L_prev = _ip_data[ip].saturation_prev;
        auto const alpha =
            medium[MPL::PropertyType::biot_coefficient].template value<double>(
                variables, x_position, t, dt);
 
        auto& solid_elasticity = *_process_data.simplified_elasticity;
        // TODO (buchwaldj)
        // is bulk_modulus good name for bulk modulus of solid skeleton?
        auto const beta_S =
            solid_elasticity.bulkCompressibilityFromYoungsModulus(
                solid_phase, variables, x_position, t, dt);
        auto const beta_SR = (1 - alpha) * beta_S;
        variables.grain_compressibility = beta_SR;
 
        auto const rho_LR =
            liquid_phase[MPL::PropertyType::density].template value<double>(
                variables, x_position, t, dt);
        variables.density = rho_LR;
        auto const& b = _process_data.specific_body_force;
 
        double const drho_LR_dp =
            liquid_phase[MPL::PropertyType::density].template dValue<double>(
                variables, MPL::Variable::liquid_phase_pressure, x_position, t,
                dt);
        auto const beta_LR = drho_LR_dp / rho_LR;
 
        S_L = medium[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[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 chi_S_L = S_L;
        auto chi_S_L_prev = S_L_prev;
        auto dchi_dS_L = 1.0;
        if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
        {
            auto const chi = [&medium, x_position, t, dt](double const S_L)
            {
                MPL::VariableArray variables;
                variables.liquid_saturation = S_L;
                return medium[MPL::PropertyType::bishops_effective_stress]
                    .template value<double>(variables, x_position, t, dt);
            };
            chi_S_L = chi(S_L);
            chi_S_L_prev = chi(S_L_prev);
 
            dchi_dS_L = medium[MPL::PropertyType::bishops_effective_stress]
                            .template dValue<double>(
                                variables, MPL::Variable::liquid_saturation,
                                x_position, t, dt);
        }
        // TODO (buchwaldj)
        // should solid_grain_pressure or effective_pore_pressure remain?
        // double const p_FR = -chi_S_L * p_cap_ip;
        // variables.solid_grain_pressure = p_FR;
 
        variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
        variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
 
        auto& phi = _ip_data[ip].porosity;
        {  // Porosity update
 
            variables_prev.porosity = _ip_data[ip].porosity_prev;
            phi = medium[MPL::PropertyType::porosity].template value<double>(
                variables, variables_prev, x_position, t, dt);
            variables.porosity = phi;
        }
 
        if (alpha < phi)
        {
            OGS_FATAL(
                "ThermoRichardsFlow: Biot-coefficient {} is smaller than "
                "porosity {} in element/integration point {}/{}.",
                alpha, phi, _element.getID(), ip);
        }
 
        double const k_rel =
            medium[MPL::PropertyType::relative_permeability]
                .template value<double>(variables, x_position, t, dt);
        auto const mu =
            liquid_phase[MPL::PropertyType::viscosity].template value<double>(
                variables, x_position, t, dt);
 
        auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
            medium[MPL::PropertyType::permeability].value(variables, x_position,
                                                          t, dt));
 
        GlobalDimMatrixType const Ki_over_mu = K_intrinsic / mu;
        GlobalDimMatrixType const rho_Ki_over_mu = rho_LR * Ki_over_mu;
 
        auto const K_pT_thermal_osmosis =
            (solid_phase.hasProperty(
                 MaterialPropertyLib::PropertyType::thermal_osmosis_coefficient)
                 ? MaterialPropertyLib::formEigenTensor<GlobalDim>(
                       solid_phase
                           [MPL::PropertyType::thermal_osmosis_coefficient]
                               .value(variables, x_position, t, dt))
                 : Eigen::MatrixXd::Zero(GlobalDim, GlobalDim));
 
        // Consider anisotropic thermal expansion.
        // Read in 3x3 tensor. 2D case also requires expansion coeff. for z-
        // component.
        Eigen::Matrix<double, 3, 3> const
            solid_linear_thermal_expansion_coefficient =
                MaterialPropertyLib::formEigenTensor<3>(
                    solid_phase
                        [MaterialPropertyLib::PropertyType::thermal_expansivity]
                            .value(variables, x_position, t, dt));
 
        auto const rho_SR =
            solid_phase[MPL::PropertyType::density].template value<double>(
                variables, x_position, t, dt);
 
        //
        // pressure equation, pressure part.
        //
        laplace_p.noalias() +=
            dNdx.transpose() * k_rel * rho_Ki_over_mu * dNdx * w;
        laplace_T.noalias() +=
            dNdx.transpose() * rho_LR * K_pT_thermal_osmosis * dNdx * w;
        const double alphaB_minus_phi = alpha - phi;
        double const a0 = alphaB_minus_phi * beta_SR;
        double const specific_storage_a_p =
            S_L * (phi * beta_LR + S_L * a0 +
                   chi_S_L * alpha * alpha *
                       solid_elasticity.storageContribution(
                           solid_phase, variables, x_position, t, dt));
        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 +
                           alpha * alpha *
                               solid_elasticity.storageContribution(
                                   solid_phase, variables, x_position, t, dt) *
                               (chi_S_L + dchi_dS_L * S_L));
        double const dspecific_storage_a_S_dp_cap =
            -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
 
        storage_p_a_p.noalias() +=
            N.transpose() * rho_LR * specific_storage_a_p * N * w;
 
        storage_p_a_S.noalias() -= N.transpose() * rho_LR *
                                   specific_storage_a_S * DeltaS_L_Deltap_cap *
                                   N * w;
 
        local_Jac
            .template block<pressure_size, pressure_size>(pressure_index,
                                                          pressure_index)
            .noalias() += N.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
                          rho_LR * dspecific_storage_a_p_dp_cap * N * w;
 
        storage_p_a_S_Jpp.noalias() -=
            N.transpose() * rho_LR *
            ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
             specific_storage_a_S * dS_L_dp_cap) /
            dt * N * w;
 
        double const dk_rel_dS_L =
            medium[MPL::PropertyType::relative_permeability]
                .template dValue<double>(variables,
                                         MPL::Variable::liquid_saturation,
                                         x_position, t, dt);
        GlobalDimVectorType const grad_p_cap = -dNdx * p_L;
        local_Jac
            .template block<pressure_size, pressure_size>(pressure_index,
                                                          pressure_index)
            .noalias() += dNdx.transpose() * rho_Ki_over_mu * grad_p_cap *
                          dk_rel_dS_L * dS_L_dp_cap * N * w;
 
        local_Jac
            .template block<pressure_size, pressure_size>(pressure_index,
                                                          pressure_index)
            .noalias() += dNdx.transpose() * rho_LR * rho_Ki_over_mu * b *
                          dk_rel_dS_L * dS_L_dp_cap * N * w;
 
        local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
            dNdx.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
 
        //
        // pressure equation, temperature part.
        //
        double const fluid_volumetric_thermal_expansion_coefficient =
            MPL::getLiquidThermalExpansivity(liquid_phase, variables, rho_LR,
                                             x_position, t, dt);
        const double eff_thermal_expansion =
            S_L * (alphaB_minus_phi *
                       solid_linear_thermal_expansion_coefficient.trace() +
                   phi * fluid_volumetric_thermal_expansion_coefficient +
                   alpha * solid_elasticity.thermalExpansivityContribution(
                               solid_linear_thermal_expansion_coefficient,
                               solid_phase, variables, x_position, t, dt));
        M_pT.noalias() -=
            N.transpose() * rho_LR * eff_thermal_expansion * N * w;
 
        //
        // temperature equation.
        //
        {
            auto const specific_heat_capacity_fluid =
                liquid_phase[MaterialPropertyLib::specific_heat_capacity]
                    .template value<double>(variables, x_position, t, dt);
 
            auto const specific_heat_capacity_solid =
                solid_phase
                    [MaterialPropertyLib::PropertyType::specific_heat_capacity]
                        .template value<double>(variables, x_position, t, dt);
 
            M_TT.noalias() +=
                w *
                (rho_SR * specific_heat_capacity_solid * (1 - phi) +
                 (S_L * rho_LR * specific_heat_capacity_fluid) * phi) *
                N.transpose() * N;
 
            auto const thermal_conductivity =
                MaterialPropertyLib::formEigenTensor<GlobalDim>(
                    medium[MaterialPropertyLib::PropertyType::
                               thermal_conductivity]
                        .value(variables, x_position, t, dt));
 
            GlobalDimVectorType const velocity_L = GlobalDimVectorType(
                -Ki_over_mu * k_rel * (dNdx * p_L - rho_LR * b) -
                K_pT_thermal_osmosis * dNdx * T);
 
            K_TT.noalias() += (dNdx.transpose() * thermal_conductivity * dNdx +
                               N.transpose() * velocity_L.transpose() * dNdx *
                                   rho_LR * specific_heat_capacity_fluid) *
                              w;
 
            //
            // temperature equation, pressure part
            //
            K_Tp.noalias() +=
                dNdx.transpose() * T_ip * K_pT_thermal_osmosis * dNdx * w;
            dK_TT_dp.noalias() -= rho_LR * specific_heat_capacity_fluid *
                                  N.transpose() * (dNdx * T).transpose() *
                                  k_rel * Ki_over_mu * dNdx * w;
 
            dK_TT_dp.noalias() -= rho_LR * specific_heat_capacity_fluid *
                                  N.transpose() * velocity_L.dot(dNdx * T) /
                                  k_rel * dk_rel_dS_L * dS_L_dp_cap * N * w;
        }
        if (gas_phase && S_L < 1.0)
        {
            variables.density = rho_LR;
 
            double const rho_wv =
                gas_phase->property(MPL::PropertyType::density)
                    .template value<double>(variables, x_position, t, dt);
 
            double const drho_wv_dT =
                gas_phase->property(MPL::PropertyType::density)
                    .template dValue<double>(variables,
                                             MPL::Variable::temperature,
                                             x_position, t, dt);
            double const drho_wv_dp =
                gas_phase->property(MPL::PropertyType::density)
                    .template dValue<double>(
                        variables, MPL::Variable::liquid_phase_pressure,
                        x_position, t, dt);
            auto const f_Tv =
                gas_phase
                    ->property(
                        MPL::PropertyType::thermal_diffusion_enhancement_factor)
                    .template value<double>(variables, x_position, t, dt);
 
            variables.porosity = phi;
            auto const tortuosity =
                medium.property(MPL::PropertyType::tortuosity)
                    .template value<double>(variables, x_position, t, dt);
            double const D_v =
                phi * (1.0 - S_L) * tortuosity *
                gas_phase->property(MPL::PropertyType::diffusion)
                    .template value<double>(variables, x_position, t, dt);
 
            double const f_Tv_D_Tv = f_Tv * D_v * drho_wv_dT;
            double const D_pv = D_v * drho_wv_dp;
 
            GlobalDimVectorType const grad_T = dNdx * T;
            GlobalDimVectorType const vapour_flux =
                -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap);
            double const specific_heat_capacity_vapour =
                gas_phase->property(MaterialPropertyLib::specific_heat_capacity)
                    .template value<double>(variables, x_position, t, dt);
 
            M_TT.noalias() +=
                w * (rho_wv * specific_heat_capacity_vapour * (1 - S_L) * phi) *
                N.transpose() * N;
 
            K_TT.noalias() += N.transpose() * vapour_flux.transpose() * dNdx *
                              specific_heat_capacity_vapour * w;
 
            double const storage_coefficient_by_water_vapor =
                phi * (rho_wv * dS_L_dp_cap + (1 - S_L) * drho_wv_dp);
 
            storage_p_a_p.noalias() +=
                N.transpose() * storage_coefficient_by_water_vapor * N * w;
 
            double const vapor_expansion_factor = phi * (1 - S_L) * drho_wv_dT;
            M_pT.noalias() += N.transpose() * vapor_expansion_factor * N * w;
 
            local_Jac
                .template block<pressure_size, temperature_size>(
                    pressure_index, temperature_index)
                .noalias() += dNdx.transpose() * f_Tv_D_Tv * dNdx * w;
 
            local_rhs.template segment<pressure_size>(pressure_index)
                .noalias() -= f_Tv_D_Tv * dNdx.transpose() * (dNdx * T) * w;
 
            laplace_p.noalias() += dNdx.transpose() * D_pv * dNdx * w;
 
            //
            // Latent heat term
            //
            if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
            {
                double const factor = phi * (1 - S_L) / rho_LR;
                // The volumetric latent heat of vaporization of liquid water
                double const L0 =
                    gas_phase->property(MPL::PropertyType::specific_latent_heat)
                        .template value<double>(variables, x_position, t, dt) *
                    rho_LR;
 
                double const drho_LR_dT =
                    liquid_phase.property(MPL::PropertyType::density)
                        .template dValue<double>(variables,
                                                 MPL::Variable::temperature,
                                                 x_position, t, dt);
 
                double const rho_wv_over_rho_L = rho_wv / rho_LR;
                M_TT.noalias() +=
                    factor * L0 *
                    (drho_wv_dT - rho_wv_over_rho_L * drho_LR_dT) *
                    N.transpose() * N * w;
 
                M_Tp.noalias() +=
                    (factor * L0 *
                         (drho_wv_dp - rho_wv_over_rho_L * drho_LR_dp) +
                     L0 * phi * rho_wv_over_rho_L * dS_L_dp_cap) *
                    N.transpose() * N * w;
 
                // temperature equation, temperature part
                K_TT.noalias() +=
                    L0 * f_Tv_D_Tv * dNdx.transpose() * dNdx * w / rho_LR;
                // temperature equation, pressure part
                K_Tp.noalias() +=
                    L0 * D_pv * dNdx.transpose() * dNdx * w / rho_LR;
            }
        }
    }
 
    if (_process_data.apply_mass_lumping)
    {
        storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
        storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
        storage_p_a_S_Jpp =
            storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
    }
 
    //
    // -- Jacobian
    //
    // temperature equation.
    local_Jac
        .template block<temperature_size, temperature_size>(temperature_index,
                                                            temperature_index)
        .noalias() += M_TT / dt + K_TT;
    // temperature equation, pressure part
    local_Jac
        .template block<temperature_size, pressure_size>(temperature_index,
                                                         pressure_index)
        .noalias() += K_Tp + dK_TT_dp;
 
    // 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, temperature part (contributed by thermal expansion).
    local_Jac
        .template block<pressure_size, temperature_size>(pressure_index,
                                                         temperature_index)
        .noalias() += M_pT / dt + laplace_T;
 
    //
    // -- Residual
    //
    // temperature equation
    local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
        M_TT * (T - T_prev) / dt + K_TT * T;
    local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
        K_Tp * p_L;
 
    // pressure equation
    local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
        laplace_p * p_L + laplace_T * T +
        (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
        M_pT * (T - T_prev) / dt;
    if (gas_phase)
    {
        if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
        {
            // Jacobian: temperature equation, pressure part
            local_Jac
                .template block<temperature_size, pressure_size>(
                    temperature_index, pressure_index)
                .noalias() += M_Tp / dt;
            // RHS: temperature part
            local_rhs.template segment<temperature_size>(temperature_index)
                .noalias() -= M_Tp * (p_L - p_L_prev) / dt;
        }
    }
}

◆ computeSecondaryVariableConcrete()

template<typename ShapeFunction , int GlobalDim>

void ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::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 1197 of file ThermoRichardsFlowFEM-impl.h.

{
    auto const T =
        local_x.template segment<temperature_size>(temperature_index);
 
    auto const p_L = local_x.template segment<pressure_size>(pressure_index);
 
    auto p_L_prev =
        local_x_prev.template segment<pressure_size>(pressure_index);
 
    auto const& medium = *_process_data.media_map.getMedium(_element.getID());
    auto const& liquid_phase = medium.phase("AqueousLiquid");
    auto const& solid_phase = medium.phase("Solid");
    MPL::VariableArray variables;
    MPL::VariableArray variables_prev;
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
 
    double saturation_avg = 0;
    double porosity_avg = 0;
 
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        auto const& N = _ip_data[ip].N;
 
        ParameterLib::SpatialPosition const x_position{
            std::nullopt, _element.getID(),
            MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
                                                            ShapeMatricesType>(
                _element, N))};
 
        double T_ip;
        NumLib::shapeFunctionInterpolate(T, N, T_ip);
 
        double p_cap_ip;
        NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
 
        double p_cap_prev_ip;
        NumLib::shapeFunctionInterpolate(-p_L_prev, N, 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;
 
        variables.temperature = T_ip;
 
        auto& S_L = _ip_data[ip].saturation;
        auto const S_L_prev = _ip_data[ip].saturation_prev;
        S_L = medium[MPL::PropertyType::saturation].template value<double>(
            variables, x_position, t, dt);
        variables.liquid_saturation = S_L;
        variables_prev.liquid_saturation = S_L_prev;
 
        auto chi_S_L = S_L;
        auto chi_S_L_prev = S_L_prev;
        if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
        {
            auto const chi = [&medium, x_position, t, dt](double const S_L)
            {
                MPL::VariableArray variables;
                variables.liquid_saturation = S_L;
                return medium[MPL::PropertyType::bishops_effective_stress]
                    .template value<double>(variables, x_position, t, dt);
            };
            chi_S_L = chi(S_L);
            chi_S_L_prev = chi(S_L_prev);
        }
        variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
        variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
 
        auto const alpha =
            medium[MPL::PropertyType::biot_coefficient].template value<double>(
                variables, x_position, t, dt);
 
        auto& solid_elasticity = *_process_data.simplified_elasticity;
        auto const beta_S =
            solid_elasticity.bulkCompressibilityFromYoungsModulus(
                solid_phase, variables, x_position, t, dt);
        auto const beta_SR = (1 - alpha) * beta_S;
        variables.grain_compressibility = beta_SR;
 
        auto& phi = _ip_data[ip].porosity;
        {  // Porosity update
            variables_prev.porosity = _ip_data[ip].porosity_prev;
            phi = medium[MPL::PropertyType::porosity].template value<double>(
                variables, variables_prev, x_position, t, dt);
            variables.porosity = phi;
        }
 
        auto const mu =
            liquid_phase[MPL::PropertyType::viscosity].template value<double>(
                variables, x_position, t, dt);
        auto const rho_LR =
            liquid_phase[MPL::PropertyType::density].template value<double>(
                variables, x_position, t, dt);
 
        auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
            medium[MPL::PropertyType::permeability].value(variables, x_position,
                                                          t, dt));
 
        double const k_rel =
            medium[MPL::PropertyType::relative_permeability]
                .template value<double>(variables, x_position, t, dt);
 
        GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
 
        auto const rho_SR =
            solid_phase[MPL::PropertyType::density].template value<double>(
                variables, x_position, t, dt);
        _ip_data[ip].dry_density_solid = (1 - phi) * rho_SR;
 
        auto const& b = _process_data.specific_body_force;
 
        auto const K_pT_thermal_osmosis =
            (solid_phase.hasProperty(
                 MaterialPropertyLib::PropertyType::thermal_osmosis_coefficient)
                 ? MaterialPropertyLib::formEigenTensor<GlobalDim>(
                       solid_phase
                           [MPL::PropertyType::thermal_osmosis_coefficient]
                               .value(variables, x_position, t, dt))
                 : Eigen::MatrixXd::Zero(GlobalDim, GlobalDim));
 
        // Compute the velocity
        auto const& dNdx = _ip_data[ip].dNdx;
        _ip_data[ip].v_darcy.noalias() = -K_over_mu * dNdx * p_L -
                                         K_pT_thermal_osmosis * dNdx * T +
                                         rho_LR * K_over_mu * b;
 
        saturation_avg += S_L;
        porosity_avg += phi;
    }
    saturation_avg /= n_integration_points;
    porosity_avg /= n_integration_points;
 
    (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
    (*_process_data.element_porosity)[_element.getID()] = porosity_avg;
}

◆ getIntPtDarcyVelocity()

template<typename ShapeFunction , int GlobalDim>

std::vector< double > const & ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getIntPtDarcyVelocity	(	const double	t,
		std::vector< GlobalVector * > const &	x,
		std::vector< NumLib::LocalToGlobalIndexMap const * > const &	dof_table,
		std::vector< double > &	cache ) const

overridevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1118 of file ThermoRichardsFlowFEM-impl.h.

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

References MathLib::createZeroedMatrix().

◆ getIntPtDryDensitySolid()

template<typename ShapeFunction , int GlobalDim>

std::vector< double > const & ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getIntPtDryDensitySolid	(	const double	t,
		std::vector< GlobalVector * > const &	x,
		std::vector< NumLib::LocalToGlobalIndexMap const * > const &	dof_table,
		std::vector< double > &	cache ) const

overridevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1185 of file ThermoRichardsFlowFEM-impl.h.

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

References ProcessLib::getIntegrationPointScalarData().

◆ getIntPtPorosity()

template<typename ShapeFunction , int GlobalDim>

std::vector< double > const & ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getIntPtPorosity	(	const double	t,
		std::vector< GlobalVector * > const &	x,
		std::vector< NumLib::LocalToGlobalIndexMap const * > const &	dof_table,
		std::vector< double > &	cache ) const

overridevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1173 of file ThermoRichardsFlowFEM-impl.h.

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

References ProcessLib::getIntegrationPointScalarData().

◆ getIntPtSaturation()

template<typename ShapeFunction , int GlobalDim>

std::vector< double > const & ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getIntPtSaturation	(	const double	t,
		std::vector< GlobalVector * > const &	x,
		std::vector< NumLib::LocalToGlobalIndexMap const * > const &	dof_table,
		std::vector< double > &	cache ) const

overridevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1152 of file ThermoRichardsFlowFEM-impl.h.

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

References ProcessLib::getIntegrationPointScalarData().

◆ getNumberOfIntegrationPoints()

template<typename ShapeFunction , int GlobalDim>

unsigned ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getNumberOfIntegrationPoints ( ) const

overrideprivatevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1343 of file ThermoRichardsFlowFEM-impl.h.

{
    return _integration_method.getNumberOfPoints();
}

◆ getPorosity()

template<typename ShapeFunction , int GlobalDim>

std::vector< double > ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getPorosity ( ) const

overridevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1164 of file ThermoRichardsFlowFEM-impl.h.

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

◆ getSaturation()

template<typename ShapeFunction , int GlobalDim>

std::vector< double > ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getSaturation ( ) const

overridevirtual

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 1143 of file ThermoRichardsFlowFEM-impl.h.

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

◆ getShapeMatrix()

template<typename ShapeFunction , int GlobalDim>

Eigen::Map< const Eigen::RowVectorXd > ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getShapeMatrix ( const unsigned integration_point ) const

inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 113 of file ThermoRichardsFlowFEM.h.

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

References ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_ip_data.

◆ initializeConcrete()

template<typename ShapeFunction , int GlobalDim>

void ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::initializeConcrete ( )

inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 83 of file ThermoRichardsFlowFEM.h.

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

References ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_integration_method, ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_ip_data, and NumLib::GenericIntegrationMethod::getNumberOfPoints().

◆ postTimestepConcrete()

template<typename ShapeFunction , int GlobalDim>

void ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::postTimestepConcrete	(	Eigen::VectorXd const &	,
		Eigen::VectorXd const &	,
		double const	,
		double const	,
		int const	)

inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 95 of file ThermoRichardsFlowFEM.h.

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

References ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_integration_method, ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_ip_data, and NumLib::GenericIntegrationMethod::getNumberOfPoints().

◆ setInitialConditionsConcrete()

template<typename ShapeFunction , int GlobalDim>

void ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::setInitialConditionsConcrete	(	Eigen::VectorXd const	local_x,
		double const	t,
		int const	process_id )

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 111 of file ThermoRichardsFlowFEM-impl.h.

{
    assert(local_x.size() == temperature_size + pressure_size);
 
    auto const p_L = local_x.template segment<pressure_size>(pressure_index);
 
    auto const& medium = *_process_data.media_map.getMedium(_element.getID());
    MPL::VariableArray variables;
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        auto const& N = _ip_data[ip].N;
 
        ParameterLib::SpatialPosition const x_position{
            std::nullopt, _element.getID(),
            MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
                                                            ShapeMatricesType>(
                _element, N))};
 
        double p_cap_ip;
        NumLib::shapeFunctionInterpolate(-p_L, N, 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;
 
        // Note: temperature dependent saturation model is not considered so
        // far.
        _ip_data[ip].saturation_prev =
            medium[MPL::PropertyType::saturation].template value<double>(
                variables, x_position, t,
                std::numeric_limits<double>::quiet_NaN());
    }
}

References MaterialPropertyLib::VariableArray::capillary_pressure, MaterialPropertyLib::VariableArray::gas_phase_pressure, NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, and NumLib::detail::shapeFunctionInterpolate().

◆ setIPDataInitialConditions()

template<typename ShapeFunction , int GlobalDim>

std::size_t ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::setIPDataInitialConditions	(	std::string_view const	name,
		double const *	values,
		int const	integration_order )

overridevirtual

Returns: the number of read integration points.

Implements ProcessLib::ThermoRichardsFlow::LocalAssemblerInterface.

Definition at line 82 of file ThermoRichardsFlowFEM-impl.h.

{
    if (integration_order !=
        static_cast<int>(_integration_method.getIntegrationOrder()))
    {
        OGS_FATAL(
            "Setting integration point initial conditions; The integration "
            "order of the local assembler for element {:d} is different "
            "from the integration order in the initial condition.",
            _element.getID());
    }
 
    if (name == "saturation")
    {
        return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
                                                         &IpData::saturation);
    }
    if (name == "porosity")
    {
        return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
                                                         &IpData::porosity);
    }
    return 0;
}

References OGS_FATAL, and ProcessLib::setIntegrationPointScalarData().

Member Data Documentation

◆ _element

template<typename ShapeFunction , int GlobalDim>

MeshLib::Element const& ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_element

private

Definition at line 157 of file ThermoRichardsFlowFEM.h.

Referenced by ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::ThermoRichardsFlowLocalAssembler().

◆ _integration_method

template<typename ShapeFunction , int GlobalDim>

NumLib::GenericIntegrationMethod const& ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_integration_method

private

Definition at line 156 of file ThermoRichardsFlowFEM.h.

Referenced by ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::ThermoRichardsFlowLocalAssembler(), ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::initializeConcrete(), and ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::postTimestepConcrete().

◆ _ip_data

template<typename ShapeFunction , int GlobalDim>

std::vector<IpData, Eigen::aligned_allocator<IpData> > ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_ip_data

private

Definition at line 154 of file ThermoRichardsFlowFEM.h.

Referenced by ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::ThermoRichardsFlowLocalAssembler(), ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::getShapeMatrix(), ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::initializeConcrete(), and ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::postTimestepConcrete().

◆ _is_axially_symmetric

template<typename ShapeFunction , int GlobalDim>

bool const ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_is_axially_symmetric

private

Definition at line 158 of file ThermoRichardsFlowFEM.h.

◆ _process_data

template<typename ShapeFunction , int GlobalDim>

ThermoRichardsFlowProcessData& ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::_process_data

private

Definition at line 152 of file ThermoRichardsFlowFEM.h.

Referenced by ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::ThermoRichardsFlowLocalAssembler().

◆ pressure_index

template<typename ShapeFunction , int GlobalDim>

const int ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::pressure_index = temperature_size

staticprivate

Definition at line 162 of file ThermoRichardsFlowFEM.h.

◆ pressure_size

template<typename ShapeFunction , int GlobalDim>

const int ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::pressure_size = ShapeFunction::NPOINTS

staticprivate

Definition at line 163 of file ThermoRichardsFlowFEM.h.

◆ temperature_index

template<typename ShapeFunction , int GlobalDim>

const int ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::temperature_index = 0

staticprivate

Definition at line 160 of file ThermoRichardsFlowFEM.h.

◆ temperature_size

template<typename ShapeFunction , int GlobalDim>

const int ProcessLib::ThermoRichardsFlow::ThermoRichardsFlowLocalAssembler< ShapeFunction, GlobalDim >::temperature_size = ShapeFunction::NPOINTS

staticprivate

Definition at line 161 of file ThermoRichardsFlowFEM.h.

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

ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM.h
ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h

Detailed Description

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Member Typedef Documentation

◆ GlobalDimMatrixType

◆ GlobalDimVectorType

◆ IpData

◆ ShapeMatricesType

Constructor & Destructor Documentation

◆ ThermoRichardsFlowLocalAssembler() [1/3]

◆ ThermoRichardsFlowLocalAssembler() [2/3]

◆ ThermoRichardsFlowLocalAssembler() [3/3]

Member Function Documentation

◆ assemble()

◆ assembleWithJacobian()

◆ computeSecondaryVariableConcrete()

◆ getIntPtDarcyVelocity()

◆ getIntPtDryDensitySolid()

◆ getIntPtPorosity()

◆ getIntPtSaturation()

◆ getNumberOfIntegrationPoints()

◆ getPorosity()

◆ getSaturation()

◆ getShapeMatrix()

◆ initializeConcrete()

◆ postTimestepConcrete()

◆ setInitialConditionsConcrete()

◆ setIPDataInitialConditions()

Member Data Documentation

◆ _element

◆ _integration_method

◆ _ip_data

◆ _is_axially_symmetric

◆ _process_data

◆ pressure_index

◆ pressure_size

◆ temperature_index

◆ temperature_size