ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > Class Template Reference

Detailed Description

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>
class ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >

Definition at line 159 of file HydroMechanicsFEM.h.

#include <HydroMechanicsFEM.h>

Inheritance diagram for ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:

Collaboration diagram for ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:

Public Types
using	ShapeMatricesTypeDisplacement
using	ShapeMatricesTypePressure
using	Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>
using	SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>

Public Member Functions
	HydroMechanicsLocalAssembler (HydroMechanicsLocalAssembler const &)=delete
	HydroMechanicsLocalAssembler (HydroMechanicsLocalAssembler &&)=delete
	HydroMechanicsLocalAssembler (MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, HydroMechanicsProcessData< DisplacementDim > &process_data)
std::size_t	setIPDataInitialConditions (std::string_view const name, double const *values, int const integration_order) override
	Returns number of read integration points.
void	assemble (double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) 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 (const double 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	postTimestepConcrete (Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, int const process_id) override
void	computeSecondaryVariableConcrete (double const t, double const dt, Eigen::VectorXd const &local_xs, Eigen::VectorXd const &local_x_prev) override
void	postNonLinearSolverConcrete (Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, int const process_id) override
void	setInitialConditionsConcrete (Eigen::VectorXd const local_x, double const t, int const process_id) override
Eigen::Map< const Eigen::RowVectorXd >	getShapeMatrix (const unsigned integration_point) const override
	Provides the shape matrix at the given integration point.
std::vector< double >	getSigma () const override
std::vector< double >	getEpsilon () const override
std::vector< double >	getStrainRateVariable () const override
std::vector< double > const &	getIntPtDarcyVelocity (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double > const &	getIntPtSigma (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
std::vector< double > const &	getIntPtEpsilon (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
int	getNumberOfVectorElementsForDeformation () const override
Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual	~LocalAssemblerInterface ()=default
virtual void	setInitialConditions (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, int const process_id)
virtual void	initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
virtual void	preAssemble (double const, double const, std::vector< double > const &)
virtual void	assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
virtual void	computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
virtual void	preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void	postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
void	postNonLinearSolver (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
virtual Eigen::Vector3d	getFlux (MathLib::Point3d const &, double const, std::vector< double > const &) const
virtual Eigen::Vector3d	getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double > > const &) const
	Fits to staggered scheme.
Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual	~ExtrapolatableElement ()=default

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

Private Types
using	BMatricesType
using	IpData

Private Member Functions
void	assembleWithJacobianForDeformationEquations (const double t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void	assembleWithJacobianForPressureEquations (const double 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)
unsigned	getNumberOfIntegrationPoints () const override
int	getMaterialID () const override
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const &	getMaterialStateVariablesAt (unsigned integration_point) const override

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

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

Member Typedef Documentation

BMatricesType

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

using ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::BMatricesType

private

Initial value:

 
        BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 384 of file HydroMechanicsFEM.h.

Invariants

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

using ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>

Definition at line 172 of file HydroMechanicsFEM.h.

IpData

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

using ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::IpData

private

Initial value:

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

Definition at line 386 of file HydroMechanicsFEM.h.

ShapeMatricesTypeDisplacement

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

using ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypeDisplacement

Initial value:

 
        ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 163 of file HydroMechanicsFEM.h.

ShapeMatricesTypePressure

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

using ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypePressure

Initial value:

 
        ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>

Definition at line 167 of file HydroMechanicsFEM.h.

SymmetricTensor

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

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

Definition at line 174 of file HydroMechanicsFEM.h.

Constructor & Destructor Documentation

HydroMechanicsLocalAssembler() [1/3]

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::HydroMechanicsLocalAssembler ( HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > const & )

delete

References HydroMechanicsLocalAssembler().

Referenced by HydroMechanicsLocalAssembler(), and HydroMechanicsLocalAssembler().

HydroMechanicsLocalAssembler() [2/3]

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::HydroMechanicsLocalAssembler ( HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > && )

delete

References HydroMechanicsLocalAssembler().

HydroMechanicsLocalAssembler() [3/3]

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::HydroMechanicsLocalAssembler	(	MeshLib::Element const &	e,
		std::size_t const	,
		NumLib::GenericIntegrationMethod const &	integration_method,
		bool const	is_axially_symmetric,
		HydroMechanicsProcessData< DisplacementDim > &	process_data )

Definition at line 29 of file HydroMechanicsFEM-impl.h.

    : _process_data(process_data),
      _integration_method(integration_method),
      _element(e),
      _is_axially_symmetric(is_axially_symmetric)
{
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
 
    _ip_data.reserve(n_integration_points);
    _secondary_data.N_u.resize(n_integration_points);
 
    auto const shape_matrices_u =
        NumLib::initShapeMatrices<ShapeFunctionDisplacement,
                                  ShapeMatricesTypeDisplacement,
                                  DisplacementDim>(e, is_axially_symmetric,
                                                   _integration_method);
 
    auto const shape_matrices_p =
        NumLib::initShapeMatrices<ShapeFunctionPressure,
                                  ShapeMatricesTypePressure, DisplacementDim>(
            e, is_axially_symmetric, _integration_method);
 
    auto const& solid_material =
        MaterialLib::Solids::selectSolidConstitutiveRelation(
            _process_data.solid_materials, _process_data.material_ids,
            e.getID());
 
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        _ip_data.emplace_back(solid_material);
        auto& ip_data = _ip_data[ip];
        auto const& sm_u = shape_matrices_u[ip];
        ip_data.integration_weight =
            _integration_method.getWeightedPoint(ip).getWeight() *
            sm_u.integralMeasure * sm_u.detJ;
 
        // Initialize current time step values
        static const int kelvin_vector_size =
            MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
        ip_data.sigma_eff.setZero(kelvin_vector_size);
        ip_data.eps.setZero(kelvin_vector_size);
 
        // Previous time step values are not initialized and are set later.
        ip_data.eps_prev.resize(kelvin_vector_size);
        ip_data.sigma_eff_prev.resize(kelvin_vector_size);
 
        ip_data.N_u = sm_u.N;
        ip_data.dNdx_u = sm_u.dNdx;
 
        ip_data.N_p = shape_matrices_p[ip].N;
        ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
 
        _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
    }
}

References _element, _integration_method, _ip_data, _is_axially_symmetric, _process_data, _secondary_data, MeshLib::Element::getID(), NumLib::initShapeMatrices(), MathLib::KelvinVector::kelvin_vector_dimensions(), and MaterialLib::Solids::selectSolidConstitutiveRelation().

Member Function Documentation

assemble()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assemble ( double const , double const , std::vector< double > const & , std::vector< double > const & , std::vector< double > & , std::vector< double > & , std::vector< double > &  )

inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 196 of file HydroMechanicsFEM.h.

    {
        OGS_FATAL(
            "HydroMechanicsLocalAssembler: assembly without jacobian is not "
            "implemented.");
    }

References OGS_FATAL.

assembleWithJacobian()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< 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 95 of file HydroMechanicsFEM-impl.h.

{
    assert(local_x.size() == pressure_size + displacement_size);
 
    auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
        pressure_size> const>(local_x.data() + pressure_index, pressure_size);
 
    auto u =
        Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
            displacement_size> const>(local_x.data() + displacement_index,
                                      displacement_size);
 
    auto p_prev =
        Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
            pressure_size> const>(local_x_prev.data() + pressure_index,
                                  pressure_size);
    auto u_prev =
        Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
            displacement_size> const>(local_x_prev.data() + displacement_index,
                                      displacement_size);
 
    auto local_Jac = MathLib::createZeroedMatrix<
        typename ShapeMatricesTypeDisplacement::template MatrixType<
            displacement_size + pressure_size,
            displacement_size + pressure_size>>(
        local_Jac_data, displacement_size + pressure_size,
        displacement_size + pressure_size);
 
    auto local_rhs = MathLib::createZeroedVector<
        typename ShapeMatricesTypeDisplacement::template VectorType<
            displacement_size + pressure_size>>(
        local_rhs_data, displacement_size + pressure_size);
 
    typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
                                                         pressure_size);
 
    typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
                                                         pressure_size);
 
    typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
        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);
 
    typename ShapeMatricesTypeDisplacement::template MatrixType<
        pressure_size, displacement_size>
        Kpu_k = ShapeMatricesTypeDisplacement::template MatrixType<
            pressure_size, displacement_size>::Zero(pressure_size,
                                                    displacement_size);
 
    auto const& solid_material =
        MaterialLib::Solids::selectSolidConstitutiveRelation(
            _process_data.solid_materials, _process_data.material_ids,
            _element.getID());
 
    ParameterLib::SpatialPosition x_position;
    x_position.setElementID(_element.getID());
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
 
    auto const& b = _process_data.specific_body_force;
    auto const& medium = _process_data.media_map.getMedium(_element.getID());
    auto const& solid = medium->phase("Solid");
    auto const& fluid = fluidPhase(*medium);
    MPL::VariableArray vars;
    auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
                                                   : vars.liquid_phase_pressure;
 
    auto const T_ref =
        medium->property(MPL::PropertyType::reference_temperature)
            .template value<double>(vars, x_position, t, dt);
    vars.temperature = T_ref;
 
    auto const& identity2 = Invariants::identity2;
 
    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, _element.getID(),
            MathLib::Point3d(
                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
                                               ShapeMatricesTypeDisplacement>(
                    _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, _is_axially_symmetric);
 
        auto& eps = _ip_data[ip].eps;
        eps.noalias() = B * u;
        auto const& sigma_eff = _ip_data[ip].sigma_eff;
 
        double const p_int_pt = N_p.dot(p);
 
        phase_pressure = p_int_pt;
 
        auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
            t, x_position, dt, T_ref);
        auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
 
        auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
                               .template value<double>(vars, x_position, t, dt);
 
        auto const rho_sr =
            solid.property(MPL::PropertyType::density)
                .template value<double>(vars, x_position, t, dt);
        auto const porosity =
            medium->property(MPL::PropertyType::porosity)
                .template value<double>(vars, x_position, t, dt);
 
        auto const [rho_fr, mu] =
            MaterialPropertyLib::getFluidDensityAndViscosity(t, dt, x_position,
                                                             fluid, vars);
 
        auto const beta_p =
            fluid.property(MPL::PropertyType::density)
                .template dValue<double>(vars, _process_data.phase_variable,
                                         x_position, t, dt) /
            rho_fr;
 
        // Set mechanical variables for the intrinsic permeability model
        // For stress dependent permeability.
        {
            auto const sigma_total =
                (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
 
            vars.total_stress.emplace<SymmetricTensor>(
                MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
                    sigma_total));
        }
        // For strain dependent permeability
        vars.volumetric_strain = Invariants::trace(eps);
        vars.equivalent_plastic_strain =
            _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
        vars.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                eps);
 
        auto const K = MPL::formEigenTensor<DisplacementDim>(
            medium->property(MPL::PropertyType::permeability)
                .value(vars, x_position, t, dt));
        auto const dkde = MPL::formEigenTensor<
            MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim)>(
            (*medium)[MPL::PropertyType::permeability].dValue(
                vars, MPL::Variable::mechanical_strain, x_position, t, dt));
 
        auto const K_over_mu = K / mu;
 
        auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
                                                         dt, u, T_ref);
 
        //
        // displacement equation, displacement part
        //
        local_Jac
            .template block<displacement_size, displacement_size>(
                displacement_index, displacement_index)
            .noalias() += B.transpose() * C * B * w;
 
        double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
        local_rhs.template segment<displacement_size>(displacement_index)
            .noalias() -=
            (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
 
        //
        // displacement equation, pressure part
        //
        Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
 
        //
        // pressure equation, pressure part.
        //
        laplace_p.noalias() +=
            rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
 
        storage_p.noalias() +=
            rho_fr * N_p.transpose() * N_p * w *
            ((alpha - porosity) * (1.0 - alpha) / K_S + porosity * beta_p);
 
        // density dependence on pressure evaluated for Darcy-term,
        // for laplace and storage terms this dependence is neglected
        add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
                                      K_over_mu *
                                      (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
 
        local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
            dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
 
        //
        // pressure equation, displacement part.
        //
        Kpu.noalias() +=
            rho_fr * alpha * N_p.transpose() * identity2.transpose() * B * w;
 
        Kpu_k.noalias() +=
            dNdx_p.transpose() *
            MathLib::KelvinVector::liftVectorToKelvin<DisplacementDim>(
                dNdx_p * p - rho_fr * b) *
            dkde * B * rho_fr / mu * w;
    }
    // displacement equation, pressure part
    local_Jac
        .template block<displacement_size, pressure_size>(displacement_index,
                                                          pressure_index)
        .noalias() = -Kup;
 
    if (_process_data.mass_lumping)
    {
        storage_p = storage_p.colwise().sum().eval().asDiagonal();
 
        if constexpr (pressure_size == displacement_size)
        {
            Kpu = Kpu.colwise().sum().eval().asDiagonal();
            Kpu_k = Kpu_k.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 / dt + add_p_derivative;
 
    // pressure equation, displacement part.
    local_Jac
        .template block<pressure_size, displacement_size>(pressure_index,
                                                          displacement_index)
        .noalias() += Kpu / dt + Kpu_k;
 
    // pressure equation
    local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
        laplace_p * p + storage_p * (p - p_prev) / dt + Kpu * (u - u_prev) / dt;
 
    // displacement equation
    local_rhs.template segment<displacement_size>(displacement_index)
        .noalias() += Kup * p;
}

References _element, _integration_method, _ip_data, _is_axially_symmetric, _process_data, MaterialPropertyLib::biot_coefficient, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, displacement_index, displacement_size, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, ParameterLib::SpatialPosition::getCoordinates(), MaterialPropertyLib::getFluidDensityAndViscosity(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::identity2, NumLib::interpolateCoordinates(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MathLib::KelvinVector::liftVectorToKelvin(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::mechanical_strain, MaterialPropertyLib::VariableArray::mechanical_strain, N_u_op, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, pressure_index, pressure_size, MaterialPropertyLib::reference_temperature, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace(), and MaterialPropertyLib::VariableArray::volumetric_strain.

assembleWithJacobianForDeformationEquations()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForDeformationEquations ( const double t, double const dt, Eigen::VectorXd const & local_x, 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 of all coupled processes.
local_b_data	Right hand side vector of an element.
local_Jac_data	Element Jacobian matrix for the Newton-Raphson method.

Definition at line 641 of file HydroMechanicsFEM-impl.h.

{
    auto const p = local_x.template segment<pressure_size>(pressure_index);
    auto const u =
        local_x.template segment<displacement_size>(displacement_index);
 
    auto local_Jac = MathLib::createZeroedMatrix<
        typename ShapeMatricesTypeDisplacement::template MatrixType<
            displacement_size, displacement_size>>(
        local_Jac_data, displacement_size, displacement_size);
 
    auto local_rhs =
        MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
                                        template VectorType<displacement_size>>(
            local_b_data, displacement_size);
 
    ParameterLib::SpatialPosition x_position;
    x_position.setElementID(_element.getID());
 
    auto const& medium = _process_data.media_map.getMedium(_element.getID());
    auto const& solid = medium->phase("Solid");
    auto const& fluid = fluidPhase(*medium);
    MPL::VariableArray vars;
    auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
                                                   : vars.liquid_phase_pressure;
 
    auto const T_ref =
        medium->property(MPL::PropertyType::reference_temperature)
            .template value<double>(vars, x_position, t, dt);
    vars.temperature = T_ref;
 
    int const n_integration_points = _integration_method.getNumberOfPoints();
    for (int 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;
 
        x_position = {
            std::nullopt, _element.getID(),
            MathLib::Point3d(
                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
                                               ShapeMatricesTypeDisplacement>(
                    _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, _is_axially_symmetric);
 
        auto& eps = _ip_data[ip].eps;
        auto const& sigma_eff = _ip_data[ip].sigma_eff;
 
        phase_pressure = N_p.dot(p);
 
        auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
                               .template value<double>(vars, x_position, t, dt);
        auto const rho_sr =
            solid.property(MPL::PropertyType::density)
                .template value<double>(vars, x_position, t, dt);
        auto const porosity =
            medium->property(MPL::PropertyType::porosity)
                .template value<double>(vars, x_position, t, dt);
 
        auto const rho_fr = MaterialPropertyLib::getFluidDensity(
            t, dt, x_position, fluid, vars);
 
        auto const& b = _process_data.specific_body_force;
        auto const& identity2 = MathLib::KelvinVector::Invariants<
            MathLib::KelvinVector::kelvin_vector_dimensions(
                DisplacementDim)>::identity2;
 
        eps.noalias() = B * u;
        vars.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                eps);
 
        auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
                                                         dt, u, T_ref);
 
        local_Jac.noalias() += B.transpose() * C * B * w;
 
        double p_at_xi = 0.;
        NumLib::shapeFunctionInterpolate(p, N_p, p_at_xi);
 
        double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
        local_rhs.noalias() -=
            (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
             N_u_op(N_u).transpose() * rho * b) *
            w;
    }
}

References _element, _integration_method, _ip_data, _is_axially_symmetric, _process_data, MaterialPropertyLib::biot_coefficient, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, displacement_index, displacement_size, MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getFluidDensity(), NumLib::interpolateCoordinates(), MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, N_u_op, MaterialPropertyLib::porosity, pressure_index, MaterialPropertyLib::reference_temperature, ParameterLib::SpatialPosition::setElementID(), NumLib::detail::shapeFunctionInterpolate(), and MaterialPropertyLib::VariableArray::temperature.

Referenced by assembleWithJacobianForStaggeredScheme().

assembleWithJacobianForPressureEquations()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForPressureEquations ( const double 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 of all coupled processes.
local_x_prev	Nodal values of \(x^{t-1}\) of an element of all coupled processes.
local_b_data	Right hand side vector of an element.
local_Jac_data	Element Jacobian matrix for the Newton-Raphson method.

Definition at line 463 of file HydroMechanicsFEM-impl.h.

{
    auto local_rhs =
        MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
                                        template VectorType<pressure_size>>(
            local_b_data, pressure_size);
 
    ParameterLib::SpatialPosition pos;
    pos.setElementID(this->_element.getID());
 
    auto const p = local_x.template segment<pressure_size>(pressure_index);
 
    auto const p_prev =
        local_x_prev.template segment<pressure_size>(pressure_index);
 
    auto local_Jac = MathLib::createZeroedMatrix<
        typename ShapeMatricesTypeDisplacement::template MatrixType<
            pressure_size, pressure_size>>(local_Jac_data, pressure_size,
                                           pressure_size);
 
    typename ShapeMatricesTypePressure::NodalMatrixType storage =
        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
                                                         pressure_size);
 
    typename ShapeMatricesTypePressure::NodalMatrixType laplace =
        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
                                                         pressure_size);
 
    typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
        ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
                                                         pressure_size);
 
    auto const& solid_material =
        MaterialLib::Solids::selectSolidConstitutiveRelation(
            _process_data.solid_materials, _process_data.material_ids,
            _element.getID());
 
    ParameterLib::SpatialPosition x_position;
    x_position.setElementID(_element.getID());
 
    auto const& medium = _process_data.media_map.getMedium(_element.getID());
    auto const& fluid = fluidPhase(*medium);
    MPL::VariableArray vars;
    auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
                                                   : vars.liquid_phase_pressure;
 
    auto const T_ref =
        medium->property(MPL::PropertyType::reference_temperature)
            .template value<double>(vars, x_position, t, dt);
    vars.temperature = T_ref;
 
    auto const& identity2 = Invariants::identity2;
 
    auto const staggered_scheme =
        std::get<Staggered>(_process_data.coupling_scheme);
    auto const fixed_stress_stabilization_parameter =
        staggered_scheme.fixed_stress_stabilization_parameter;
    auto const fixed_stress_over_time_step =
        staggered_scheme.fixed_stress_over_time_step;
 
    int const n_integration_points = _integration_method.getNumberOfPoints();
    for (int ip = 0; ip < n_integration_points; ip++)
    {
        auto const& w = _ip_data[ip].integration_weight;
 
        auto const& N_p = _ip_data[ip].N_p;
        auto const& dNdx_p = _ip_data[ip].dNdx_p;
 
        x_position = {
            std::nullopt, _element.getID(),
            MathLib::Point3d(
                NumLib::interpolateCoordinates<ShapeFunctionPressure,
                                               ShapeMatricesTypePressure>(
                    _element, _ip_data[ip].N_p))};
 
        double const p_int_pt = N_p.dot(p);
 
        phase_pressure = p_int_pt;
 
        auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
            t, x_position, dt, T_ref);
        auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
 
        auto const alpha_b =
            medium->property(MPL::PropertyType::biot_coefficient)
                .template value<double>(vars, x_position, t, dt);
 
        // Set mechanical variables for the intrinsic permeability model
        // For stress dependent permeability.
        auto const sigma_total =
            (_ip_data[ip].sigma_eff - alpha_b * p_int_pt * identity2).eval();
        vars.total_stress.emplace<SymmetricTensor>(
            MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
 
        // For strain dependent permeability
        vars.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
        vars.equivalent_plastic_strain =
            _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
        vars.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                _ip_data[ip].eps);
 
        auto const K = MPL::formEigenTensor<DisplacementDim>(
            medium->property(MPL::PropertyType::permeability)
                .value(vars, x_position, t, dt));
        auto const porosity =
            medium->property(MPL::PropertyType::porosity)
                .template value<double>(vars, x_position, t, dt);
 
        auto const [rho_fr, mu] =
            MaterialPropertyLib::getFluidDensityAndViscosity(t, dt, x_position,
                                                             fluid, vars);
 
        auto const beta_p =
            fluid.property(MPL::PropertyType::density)
                .template dValue<double>(vars, _process_data.phase_variable,
                                         x_position, t, dt) /
            rho_fr;
 
        auto const K_over_mu = K / mu;
 
        laplace.noalias() +=
            rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
 
        // Artificial compressibility from the fixed stress splitting:
        auto const beta_FS =
            fixed_stress_stabilization_parameter * alpha_b * alpha_b / K_S;
 
        storage.noalias() += rho_fr * N_p.transpose() * N_p * w *
                             ((alpha_b - porosity) * (1.0 - alpha_b) / K_S +
                              porosity * beta_p + beta_FS);
 
        auto const& b = _process_data.specific_body_force;
 
        // bodyforce-driven Darcy flow
        local_rhs.noalias() +=
            dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
 
        // density dependence on pressure evaluated for Darcy-term,
        // for laplace and storage terms this dependence is neglected (as is
        // done for monolithic too)
        add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
                                      K_over_mu *
                                      (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
 
        if (!fixed_stress_over_time_step)
        {
            auto const& eps = _ip_data[ip].eps;
            auto const& eps_prev = _ip_data[ip].eps_prev;
            const double eps_v_dot =
                (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
 
            // Constant portion of strain rate term:
            double const strain_rate_b =
                alpha_b * eps_v_dot -
                beta_FS * _ip_data[ip].strain_rate_variable;
 
            local_rhs.noalias() -= strain_rate_b * rho_fr * N_p * w;
        }
        else
        {
            // Constant portion of strain rate term:
            local_rhs.noalias() -=
                alpha_b * _ip_data[ip].strain_rate_variable * rho_fr * N_p * w;
        }
    }
    local_Jac.noalias() = laplace + storage / dt + add_p_derivative;
 
    local_rhs.noalias() -= laplace * p + storage * (p - p_prev) / dt;
}

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::biot_coefficient, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getFluidDensityAndViscosity(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::identity2, NumLib::interpolateCoordinates(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, pressure_index, pressure_size, MaterialPropertyLib::reference_temperature, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace(), and MaterialPropertyLib::VariableArray::volumetric_strain.

Referenced by assembleWithJacobianForStaggeredScheme().

assembleWithJacobianForStaggeredScheme()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForStaggeredScheme ( const double 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 745 of file HydroMechanicsFEM-impl.h.

{
    // For the equations with pressure
    if (process_id == _process_data.hydraulic_process_id)
    {
        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_b_data,
                                                local_Jac_data);
}

References _process_data, assembleWithJacobianForDeformationEquations(), and assembleWithJacobianForPressureEquations().

computeSecondaryVariableConcrete()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::computeSecondaryVariableConcrete ( double const t, double const dt, Eigen::VectorXd const & local_xs, Eigen::VectorXd const & local_x_prev )

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1108 of file HydroMechanicsFEM-impl.h.

{
    auto const p = local_x.template segment<pressure_size>(pressure_index);
 
    NumLib::interpolateToHigherOrderNodes<
        ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
        DisplacementDim>(_element, _is_axially_symmetric, p,
                         *_process_data.pressure_interpolated);
 
    int const elem_id = _element.getID();
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
 
    auto const& medium = _process_data.media_map.getMedium(elem_id);
    MPL::VariableArray vars;
    auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
                                                   : vars.liquid_phase_pressure;
 
    SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize);
    auto sigma_eff_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
        Eigen::Matrix<double, 3, 3>::Zero());
 
    auto const& identity2 = Invariants::identity2;
 
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        auto const& eps = _ip_data[ip].eps;
        sigma_eff_sum += _ip_data[ip].sigma_eff;
 
        ParameterLib::SpatialPosition x_position = {
            std::nullopt, _element.getID(),
            MathLib::Point3d(
                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
                                               ShapeMatricesTypeDisplacement>(
                    _element, _ip_data[ip].N_u))};
 
        auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
                               .template value<double>(vars, x_position, t, dt);
        double const p_int_pt = _ip_data[ip].N_p.dot(p);
 
        phase_pressure = p_int_pt;
 
        // Set mechanical variables for the intrinsic permeability model
        // For stress dependent permeability.
        {
            auto const sigma_total =
                (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
            vars.total_stress.emplace<SymmetricTensor>(
                MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
                    sigma_total));
        }
        // For strain dependent permeability
        vars.volumetric_strain = Invariants::trace(eps);
        vars.equivalent_plastic_strain =
            _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
        vars.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                eps);
 
        k_sum += MPL::getSymmetricTensor<DisplacementDim>(
            medium->property(MPL::PropertyType::permeability)
                .value(vars, x_position, t, dt));
    }
 
    Eigen::Map<Eigen::VectorXd>(
        &(*_process_data.permeability)[elem_id * KelvinVectorSize],
        KelvinVectorSize) = k_sum / n_integration_points;
 
    Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
        MathLib::KelvinVector::kelvinVectorToTensor(sigma_eff_sum) /
        n_integration_points;
 
    Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma_avg);
 
    Eigen::Map<Eigen::Vector3d>(
        &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
        e_s.eigenvalues();
 
    auto eigen_vectors = e_s.eigenvectors();
 
    for (auto i = 0; i < 3; i++)
    {
        Eigen::Map<Eigen::Vector3d>(
            &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
            eigen_vectors.col(i);
    }
}

References _element, _integration_method, _ip_data, _is_axially_symmetric, _process_data, MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getSymmetricTensor(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::identity2, NumLib::interpolateCoordinates(), NumLib::interpolateToHigherOrderNodes(), KelvinVectorSize, MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MathLib::KelvinVector::kelvinVectorToTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::permeability, pressure_index, MathLib::KelvinVector::tensorToKelvin(), MaterialPropertyLib::VariableArray::total_stress, MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace(), and MaterialPropertyLib::VariableArray::volumetric_strain.

getEpsilon()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

std::vector< double > ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getEpsilon ( ) const

overridevirtual

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

Definition at line 1064 of file HydroMechanicsFEM-impl.h.

{
    auto const kelvin_vector_size =
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
 
    std::vector<double> ip_epsilon_values;
    auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
        double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
        ip_epsilon_values, n_integration_points, kelvin_vector_size);
 
    for (unsigned ip = 0; ip < n_integration_points; ++ip)
    {
        auto const& eps = _ip_data[ip].eps;
        cache_mat.row(ip) =
            MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
    }
 
    return ip_epsilon_values;
}

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

Referenced by getEpsilon().

getIntPtDarcyVelocity()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

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

overridevirtual

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

Definition at line 367 of file HydroMechanicsFEM-impl.h.

{
    int const hydraulic_process_id = _process_data.hydraulic_process_id;
    auto const indices =
        NumLib::getIndices(_element.getID(), *dof_table[hydraulic_process_id]);
    assert(!indices.empty());
    auto const local_x = x[hydraulic_process_id]->get(indices);
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
    cache.clear();
    auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
        double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
        cache, DisplacementDim, n_integration_points);
 
    auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
        pressure_size> const>(local_x.data() + pressure_index, pressure_size);
 
    ParameterLib::SpatialPosition x_position;
    x_position.setElementID(_element.getID());
 
    auto const& medium = _process_data.media_map.getMedium(_element.getID());
    auto const& fluid = fluidPhase(*medium);
    MPL::VariableArray vars;
    auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
                                                   : vars.liquid_phase_pressure;
 
    // TODO (naumov) Temporary value not used by current material models. Need
    // extension of secondary variables interface.
    double const dt = std::numeric_limits<double>::quiet_NaN();
    vars.temperature =
        medium->property(MPL::PropertyType::reference_temperature)
            .template value<double>(vars, x_position, t, dt);
 
    auto const& identity2 = Invariants::identity2;
 
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        x_position = {
            std::nullopt, _element.getID(),
            MathLib::Point3d(
                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
                                               ShapeMatricesTypeDisplacement>(
                    _element, _ip_data[ip].N_u))};
 
        double const p_int_pt = _ip_data[ip].N_p.dot(p);
 
        phase_pressure = p_int_pt;
 
        auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
                               .template value<double>(vars, x_position, t, dt);
 
        // Set mechanical variables for the intrinsic permeability model
        // For stress dependent permeability.
        auto const sigma_total =
            (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
        vars.total_stress.emplace<SymmetricTensor>(
            MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
 
        // For strain dependent permeability
        vars.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
        vars.equivalent_plastic_strain =
            _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
        vars.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                _ip_data[ip].eps);
 
        auto const K = MPL::formEigenTensor<DisplacementDim>(
            medium->property(MPL::PropertyType::permeability)
                .value(vars, x_position, t, dt));
 
        auto const [rho_fr, mu] =
            MaterialPropertyLib::getFluidDensityAndViscosity(t, dt, x_position,
                                                             fluid, vars);
 
        auto const K_over_mu = K / mu;
 
        auto const& b = _process_data.specific_body_force;
 
        // Compute the velocity
        auto const& dNdx_p = _ip_data[ip].dNdx_p;
        cache_matrix.col(ip).noalias() =
            -K_over_mu * dNdx_p * p + K_over_mu * rho_fr * b;
    }
 
    return cache;
}

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::biot_coefficient, MathLib::createZeroedMatrix(), MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialPropertyLib::getFluidDensityAndViscosity(), NumLib::getIndices(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::identity2, NumLib::interpolateCoordinates(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::permeability, pressure_index, pressure_size, MaterialPropertyLib::reference_temperature, ParameterLib::SpatialPosition::setElementID(), MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace(), and MaterialPropertyLib::VariableArray::volumetric_strain.

getIntPtEpsilon()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

std::vector< double > const & ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtEpsilon ( const double , std::vector< GlobalVector * > const & , std::vector< NumLib::LocalToGlobalIndexMap const * > const & , std::vector< double > & cache ) const

inlineoverridevirtual

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

Definition at line 307 of file HydroMechanicsFEM.h.

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

References _ip_data, ProcessLib::HydroMechanics::IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >::eps, and ProcessLib::getIntegrationPointKelvinVectorData().

getIntPtSigma()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

std::vector< double > const & ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getIntPtSigma ( const double , std::vector< GlobalVector * > const & , std::vector< NumLib::LocalToGlobalIndexMap const * > const & , std::vector< double > & cache ) const

inlineoverridevirtual

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

Definition at line 297 of file HydroMechanicsFEM.h.

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

References _ip_data, ProcessLib::getIntegrationPointKelvinVectorData(), and ProcessLib::HydroMechanics::IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >::sigma_eff.

getMaterialID()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

int ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getMaterialID ( ) const

overrideprivatevirtual

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

Definition at line 1213 of file HydroMechanicsFEM-impl.h.

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

References _element, _process_data, and getMaterialID().

Referenced by getMaterialID().

getMaterialStateVariablesAt()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

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

overrideprivatevirtual

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

Definition at line 1225 of file HydroMechanicsFEM-impl.h.

{
    return *_ip_data[integration_point].material_state_variables;
}

References _ip_data.

getNumberOfIntegrationPoints()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

unsigned ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getNumberOfIntegrationPoints ( ) const

overrideprivatevirtual

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

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

{
    return _integration_method.getNumberOfPoints();
}

References _integration_method, and getNumberOfIntegrationPoints().

Referenced by getNumberOfIntegrationPoints().

getNumberOfVectorElementsForDeformation()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

int ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getNumberOfVectorElementsForDeformation ( ) const

inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 317 of file HydroMechanicsFEM.h.

    {
        return displacement_size;
    }

References displacement_size.

getShapeMatrix()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

Eigen::Map< const Eigen::RowVectorXd > ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< 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 273 of file HydroMechanicsFEM.h.

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

References _secondary_data.

getSigma()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

std::vector< double > ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getSigma ( ) const

overridevirtual

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

Definition at line 1054 of file HydroMechanicsFEM-impl.h.

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

References _ip_data, ProcessLib::getIntegrationPointKelvinVectorData(), getSigma(), and ProcessLib::HydroMechanics::IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >::sigma_eff.

Referenced by getSigma().

getStrainRateVariable()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

std::vector< double > ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getStrainRateVariable ( ) const

overridevirtual

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

Definition at line 1090 of file HydroMechanicsFEM-impl.h.

{
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
 
    std::vector<double> ip_strain_rate_variables(n_integration_points);
 
    for (unsigned ip = 0; ip < n_integration_points; ++ip)
    {
        ip_strain_rate_variables[ip] = _ip_data[ip].strain_rate_variable;
    }
 
    return ip_strain_rate_variables;
}

References _integration_method, _ip_data, and getStrainRateVariable().

Referenced by getStrainRateVariable().

initializeConcrete()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete ( )

inlineoverridevirtual

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 220 of file HydroMechanicsFEM.h.

    {
        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();
 
        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            auto& ip_data = _ip_data[ip];
 
            ParameterLib::SpatialPosition const x_position{
                std::nullopt, _element.getID(),
                MathLib::Point3d(
                    NumLib::interpolateCoordinates<ShapeFunctionPressure,
                                                   ShapeMatricesTypePressure>(
                        _element, ip_data.N_p))};

            if (_process_data.initial_stress.value)
            {
                ip_data.sigma_eff =
                    MathLib::KelvinVector::symmetricTensorToKelvinVector<
                        DisplacementDim>((*_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
            ip_data.solid_material.initializeInternalStateVariables(
                t, x_position, *ip_data.material_state_variables);
 
            ip_data.pushBackState();
        }
    }

References _element, _integration_method, _ip_data, _process_data, NumLib::interpolateCoordinates(), and MathLib::KelvinVector::symmetricTensorToKelvinVector().

postNonLinearSolverConcrete()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::postNonLinearSolverConcrete ( Eigen::VectorXd const & local_x, Eigen::VectorXd const & local_x_prev, double const t, double const dt, int const process_id )

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 831 of file HydroMechanicsFEM-impl.h.

{
    // Note: local_x and local_x_prev only contain the solutions of current
    // process in the staggered scheme. This has to be changed according to the
    // same two arguments in postTimestepConcrete.
 
    int const n_integration_points = _integration_method.getNumberOfPoints();
 
    auto const staggered_scheme_ptr =
        std::get_if<Staggered>(&_process_data.coupling_scheme);
 
    if (staggered_scheme_ptr &&
        process_id == _process_data.hydraulic_process_id)
    {
        if (!staggered_scheme_ptr->fixed_stress_over_time_step)
        {
            auto const p =
                local_x.template segment<pressure_size>(pressure_index);
            auto const p_prev =
                local_x_prev.template segment<pressure_size>(pressure_index);
 
            for (int ip = 0; ip < n_integration_points; ip++)
            {
                auto& ip_data = _ip_data[ip];
 
                auto const& N_p = ip_data.N_p;
 
                ip_data.strain_rate_variable = N_p.dot(p - p_prev) / dt;
            }
        }
    }
 
    if (!staggered_scheme_ptr ||
        process_id == _process_data.mechanics_related_process_id)
    {
        ParameterLib::SpatialPosition x_position;
        x_position.setElementID(_element.getID());
 
        auto const& medium =
            _process_data.media_map.getMedium(_element.getID());
 
        auto const T_ref =
            medium->property(MPL::PropertyType::reference_temperature)
                .template value<double>(MPL::EmptyVariableArray, x_position, t,
                                        dt);
 
        auto const u =
            local_x.template segment<displacement_size>(displacement_index);
 
        MPL::VariableArray vars;
        vars.temperature = T_ref;
 
        for (int ip = 0; ip < n_integration_points; ip++)
        {
            auto const& N_u = _ip_data[ip].N_u;
            auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
            x_position = {std::nullopt, _element.getID(),
                          MathLib::Point3d(NumLib::interpolateCoordinates<
                                           ShapeFunctionDisplacement,
                                           ShapeMatricesTypeDisplacement>(
                              _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,
                                                     _is_axially_symmetric);
 
            auto& eps = _ip_data[ip].eps;
            eps.noalias() = B * u;
            vars.mechanical_strain.emplace<
                MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(eps);
 
            _ip_data[ip].updateConstitutiveRelation(vars, t, x_position, dt, u,
                                                    T_ref);
        }
    }
}

References _element, _integration_method, _ip_data, _is_axially_symmetric, _process_data, ProcessLib::LinearBMatrix::computeBMatrix(), displacement_index, MaterialPropertyLib::EmptyVariableArray, ParameterLib::SpatialPosition::getCoordinates(), NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::mechanical_strain, pressure_index, MaterialPropertyLib::reference_temperature, ParameterLib::SpatialPosition::setElementID(), and MaterialPropertyLib::VariableArray::temperature.

postTimestepConcrete()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::postTimestepConcrete ( Eigen::VectorXd const & local_x, Eigen::VectorXd const & local_x_prev, double const t, double const dt, int const process_id )

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 919 of file HydroMechanicsFEM-impl.h.

{
    auto const staggered_scheme_ptr =
        std::get_if<Staggered>(&_process_data.coupling_scheme);
 
    if (staggered_scheme_ptr &&
        process_id == _process_data.hydraulic_process_id)
    {
        if (staggered_scheme_ptr->fixed_stress_over_time_step)
        {
            auto const fixed_stress_stabilization_parameter =
                staggered_scheme_ptr->fixed_stress_stabilization_parameter;
 
            auto const p =
                local_x.template segment<pressure_size>(pressure_index);
            auto const p_prev =
                local_x_prev.template segment<pressure_size>(pressure_index);
 
            ParameterLib::SpatialPosition x_position;
            x_position.setElementID(_element.getID());
 
            auto const& solid_material =
                MaterialLib::Solids::selectSolidConstitutiveRelation(
                    _process_data.solid_materials, _process_data.material_ids,
                    _element.getID());
 
            auto const& medium =
                _process_data.media_map.getMedium(_element.getID());
            MPL::VariableArray vars;
 
            auto const T_ref =
                medium->property(MPL::PropertyType::reference_temperature)
                    .template value<double>(vars, x_position, t, dt);
            vars.temperature = T_ref;
 
            int const n_integration_points =
                _integration_method.getNumberOfPoints();
            for (int ip = 0; ip < n_integration_points; ip++)
            {
                auto& ip_data = _ip_data[ip];
 
                auto const& N_p = ip_data.N_p;
 
                x_position = {
                    std::nullopt, _element.getID(),
                    MathLib::Point3d(
                        NumLib::interpolateCoordinates<
                            ShapeFunctionPressure, ShapeMatricesTypePressure>(
                            _element, N_p))};
 
                auto const& eps = ip_data.eps;
                auto const& eps_prev = ip_data.eps_prev;
                const double eps_v_dot =
                    (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
 
                auto const C_el = ip_data.computeElasticTangentStiffness(
                    t, x_position, dt, T_ref);
                auto const K_S =
                    solid_material.getBulkModulus(t, x_position, &C_el);
 
                auto const alpha_b =
                    medium->property(MPL::PropertyType::biot_coefficient)
                        .template value<double>(vars, x_position, t, dt);
 
                ip_data.strain_rate_variable =
                    eps_v_dot - fixed_stress_stabilization_parameter * alpha_b *
                                    N_p.dot(p - p_prev) / dt / K_S;
            }
        }
    }
 
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();
 
    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        _ip_data[ip].pushBackState();
    }
}

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::biot_coefficient, NumLib::interpolateCoordinates(), postTimestepConcrete(), pressure_index, MaterialPropertyLib::reference_temperature, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), MaterialPropertyLib::VariableArray::temperature, and MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace().

Referenced by postTimestepConcrete().

setInitialConditionsConcrete()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

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

overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 769 of file HydroMechanicsFEM-impl.h.

{
    ParameterLib::SpatialPosition x_position;
    x_position.setElementID(_element.getID());
 
    auto const& medium = _process_data.media_map.getMedium(_element.getID());
 
    auto const p = local_x.template segment<pressure_size>(pressure_index);
    auto const u =
        local_x.template segment<displacement_size>(displacement_index);
 
    auto const& identity2 = Invariants::identity2;
    const double dt = 0.0;
 
    MPL::VariableArray vars;
 
    int const n_integration_points = _integration_method.getNumberOfPoints();
    for (int ip = 0; ip < n_integration_points; ip++)
    {
        auto const& N_u = _ip_data[ip].N_u;
        auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
        x_position = {
            std::nullopt, _element.getID(),
            MathLib::Point3d(
                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
                                               ShapeMatricesTypeDisplacement>(
                    _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, _is_axially_symmetric);
 
        auto& eps = _ip_data[ip].eps;
        eps.noalias() = B * u;
        vars.mechanical_strain
            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                eps);
 
        if (_process_data.initial_stress.isTotalStress())
        {
            auto const& N_p = _ip_data[ip].N_p;
            auto const alpha_b =
                medium->property(MPL::PropertyType::biot_coefficient)
                    .template value<double>(vars, x_position, t, dt);
 
            auto& sigma_eff = _ip_data[ip].sigma_eff;
            sigma_eff.noalias() += alpha_b * N_p.dot(p) * identity2;
            _ip_data[ip].sigma_eff_prev.noalias() = sigma_eff;
        }
    }
}

References _element, _integration_method, _ip_data, _is_axially_symmetric, _process_data, MaterialPropertyLib::biot_coefficient, ProcessLib::LinearBMatrix::computeBMatrix(), displacement_index, ParameterLib::SpatialPosition::getCoordinates(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::identity2, NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::mechanical_strain, pressure_index, and ParameterLib::SpatialPosition::setElementID().

setIPDataInitialConditions()

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

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

overridevirtual

Returns number of read integration points.

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

Definition at line 1006 of file HydroMechanicsFEM-impl.h.

{
    if (integration_order !=
        static_cast<int>(_integration_method.getIntegrationOrder()))
    {
        OGS_FATAL(
            "Setting integration point initial conditions; The integration "
            "order of the local assembler for element {:d} is different from "
            "the integration order in the initial condition.",
            _element.getID());
    }
 
    if (name == "sigma")
    {
        if (_process_data.initial_stress.value != nullptr)
        {
            OGS_FATAL(
                "Setting initial conditions for stress from integration "
                "point data and from a parameter '{:s}' is not possible "
                "simultaneously.",
                _process_data.initial_stress.value->name);
        }
 
        return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
            values, _ip_data, &IpData::sigma_eff);
    }
 
    if (name == "epsilon")
    {
        return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
            values, _ip_data, &IpData::eps);
    }
 
    if (name == "strain_rate_variable")
    {
        return ProcessLib::setIntegrationPointScalarData(
            values, _ip_data, &IpData::strain_rate_variable);
    }
 
    return 0;
}

References _element, _integration_method, _ip_data, _process_data, ProcessLib::HydroMechanics::IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >::eps, OGS_FATAL, ProcessLib::setIntegrationPointKelvinVectorData(), ProcessLib::setIntegrationPointScalarData(), ProcessLib::HydroMechanics::IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >::sigma_eff, and ProcessLib::HydroMechanics::IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >::strain_rate_variable.

Member Data Documentation

_element

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

MeshLib::Element const& ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_element

private

Definition at line 393 of file HydroMechanicsFEM.h.

Referenced by HydroMechanicsLocalAssembler(), assembleWithJacobian(), assembleWithJacobianForDeformationEquations(), assembleWithJacobianForPressureEquations(), computeSecondaryVariableConcrete(), getIntPtDarcyVelocity(), getMaterialID(), initializeConcrete(), postNonLinearSolverConcrete(), postTimestepConcrete(), setInitialConditionsConcrete(), and setIPDataInitialConditions().

_integration_method

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

NumLib::GenericIntegrationMethod const& ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_integration_method

private

Definition at line 392 of file HydroMechanicsFEM.h.

Referenced by HydroMechanicsLocalAssembler(), assembleWithJacobian(), assembleWithJacobianForDeformationEquations(), assembleWithJacobianForPressureEquations(), computeSecondaryVariableConcrete(), getEpsilon(), getIntPtDarcyVelocity(), getNumberOfIntegrationPoints(), getStrainRateVariable(), initializeConcrete(), postNonLinearSolverConcrete(), postTimestepConcrete(), setInitialConditionsConcrete(), and setIPDataInitialConditions().

_ip_data

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

std::vector<IpData, Eigen::aligned_allocator<IpData> > ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data

private

Definition at line 390 of file HydroMechanicsFEM.h.

Referenced by HydroMechanicsLocalAssembler(), assembleWithJacobian(), assembleWithJacobianForDeformationEquations(), assembleWithJacobianForPressureEquations(), computeSecondaryVariableConcrete(), getEpsilon(), getIntPtDarcyVelocity(), getIntPtEpsilon(), getIntPtSigma(), getMaterialStateVariablesAt(), getSigma(), getStrainRateVariable(), initializeConcrete(), postNonLinearSolverConcrete(), postTimestepConcrete(), setInitialConditionsConcrete(), and setIPDataInitialConditions().

_is_axially_symmetric

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

bool const ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_is_axially_symmetric

private

Definition at line 394 of file HydroMechanicsFEM.h.

Referenced by HydroMechanicsLocalAssembler(), assembleWithJacobian(), assembleWithJacobianForDeformationEquations(), computeSecondaryVariableConcrete(), postNonLinearSolverConcrete(), and setInitialConditionsConcrete().

_process_data

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

HydroMechanicsProcessData<DisplacementDim>& ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_process_data

private

Definition at line 382 of file HydroMechanicsFEM.h.

Referenced by HydroMechanicsLocalAssembler(), assembleWithJacobian(), assembleWithJacobianForDeformationEquations(), assembleWithJacobianForPressureEquations(), assembleWithJacobianForStaggeredScheme(), computeSecondaryVariableConcrete(), getIntPtDarcyVelocity(), getMaterialID(), initializeConcrete(), postNonLinearSolverConcrete(), postTimestepConcrete(), setInitialConditionsConcrete(), and setIPDataInitialConditions().

_secondary_data

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType> ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data

private

Definition at line 397 of file HydroMechanicsFEM.h.

Referenced by HydroMechanicsLocalAssembler(), and getShapeMatrix().

displacement_index

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

const int ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::displacement_index = ShapeFunctionPressure::NPOINTS

staticprivate

Definition at line 401 of file HydroMechanicsFEM.h.

Referenced by assembleWithJacobian(), assembleWithJacobianForDeformationEquations(), postNonLinearSolverConcrete(), and setInitialConditionsConcrete().

displacement_size

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

const int ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::displacement_size

staticprivate

Initial value:

=
        ShapeFunctionDisplacement::NPOINTS * DisplacementDim

Definition at line 402 of file HydroMechanicsFEM.h.

Referenced by assembleWithJacobian(), assembleWithJacobianForDeformationEquations(), and getNumberOfVectorElementsForDeformation().

KelvinVectorSize

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

int const ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::KelvinVectorSize

static

Initial value:

=
        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim)

Definition at line 170 of file HydroMechanicsFEM.h.

Referenced by computeSecondaryVariableConcrete().

N_u_op

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

auto& ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::N_u_op

staticconstexpr

Initial value:

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

Definition at line 176 of file HydroMechanicsFEM.h.

Referenced by assembleWithJacobian(), and assembleWithJacobianForDeformationEquations().

pressure_index

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

const int ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::pressure_index = 0

staticprivate

Definition at line 399 of file HydroMechanicsFEM.h.

Referenced by assembleWithJacobian(), assembleWithJacobianForDeformationEquations(), assembleWithJacobianForPressureEquations(), computeSecondaryVariableConcrete(), getIntPtDarcyVelocity(), postNonLinearSolverConcrete(), postTimestepConcrete(), and setInitialConditionsConcrete().

pressure_size

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, int DisplacementDim>

const int ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::pressure_size = ShapeFunctionPressure::NPOINTS

staticprivate

Definition at line 400 of file HydroMechanicsFEM.h.

Referenced by assembleWithJacobian(), assembleWithJacobianForPressureEquations(), and getIntPtDarcyVelocity().

---

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

• HydroMechanicsFEM.h
• HydroMechanicsFEM-impl.h