From d1feb884905648cedca93156192a74dfdece597a Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Fri, 16 Oct 2020 15:18:06 +0200
Subject: [PATCH] [Process/TRM] Added local assembler
---
.../IntegrationPointData.h | 182 +++
.../ThermoRichardsMechanicsFEM-impl.h | 1443 +++++++++++++++++
.../ThermoRichardsMechanicsFEM.h | 258 +++
3 files changed, 1883 insertions(+)
create mode 100644 ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h
create mode 100644 ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
create mode 100644 ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h
diff --git a/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h b/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h
new file mode 100644
index 00000000000..caad48f6f7c
--- /dev/null
+++ b/ProcessLib/ThermoRichardsMechanics/IntegrationPointData.h
@@ -0,0 +1,182 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <limits>
+#include <memory>
+
+#include "MaterialLib/MPL/VariableType.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MathLib/KelvinVector.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsMechanics
+{
+template <typename BMatricesType, typename ShapeMatrixTypeDisplacement,
+ typename ShapeMatricesTypePressure, int DisplacementDim, int NPoints>
+struct IntegrationPointData final
+{
+ explicit IntegrationPointData(
+ MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
+ solid_material)
+ : solid_material(solid_material),
+ material_state_variables(
+ solid_material.createMaterialStateVariables())
+ {
+ // Initialize current time step values
+ static const int kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value;
+ sigma_eff.setZero(kelvin_vector_size);
+ sigma_sw.setZero(kelvin_vector_size);
+ eps.setZero(kelvin_vector_size);
+ eps_m.setZero(kelvin_vector_size);
+
+ // Previous time step values are not initialized and are set later.
+ eps_prev.resize(kelvin_vector_size);
+ eps_m_prev.resize(kelvin_vector_size);
+ sigma_eff_prev.resize(kelvin_vector_size);
+ }
+
+ typename ShapeMatrixTypeDisplacement::template MatrixType<
+ DisplacementDim, NPoints * DisplacementDim>
+ N_u_op;
+ typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
+ typename BMatricesType::KelvinVectorType sigma_sw, sigma_sw_prev;
+ typename BMatricesType::KelvinVectorType eps, eps_prev;
+ typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;
+
+ typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
+ typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
+
+ typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
+ typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
+
+ typename ShapeMatricesTypePressure::GlobalDimVectorType v_darcy;
+
+ double saturation = std::numeric_limits<double>::quiet_NaN();
+ double saturation_prev = std::numeric_limits<double>::quiet_NaN();
+ double porosity = std::numeric_limits<double>::quiet_NaN();
+ double porosity_prev = std::numeric_limits<double>::quiet_NaN();
+ double transport_porosity = std::numeric_limits<double>::quiet_NaN();
+ double transport_porosity_prev = std::numeric_limits<double>::quiet_NaN();
+ double dry_density_solid = std::numeric_limits<double>::quiet_NaN();
+ double dry_density_pellet_saturated =
+ std::numeric_limits<double>::quiet_NaN();
+ double dry_density_pellet_unsaturated =
+ std::numeric_limits<double>::quiet_NaN();
+
+ MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
+ std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
+ DisplacementDim>::MaterialStateVariables>
+ material_state_variables;
+
+ double integration_weight = std::numeric_limits<double>::quiet_NaN();
+
+ void pushBackState()
+ {
+ eps_prev = eps;
+ eps_m_prev = eps_m;
+ sigma_eff_prev = sigma_eff;
+ sigma_sw_prev = sigma_sw;
+ saturation_prev = saturation;
+ porosity_prev = porosity;
+ transport_porosity_prev = transport_porosity;
+ material_state_variables->pushBackState();
+ }
+
+ MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>
+ computeElasticTangentStiffness(
+ double const t,
+ ParameterLib::SpatialPosition const& x_position,
+ double const dt,
+ double const temperature)
+ {
+ namespace MPL = MaterialPropertyLib;
+
+ MPL::VariableArray variable_array;
+ MPL::VariableArray variable_array_prev;
+
+ auto const null_state = solid_material.createMaterialStateVariables();
+
+ using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
+
+ variable_array[static_cast<int>(MPL::Variable::stress)].emplace<KV>(
+ KV::Zero());
+ variable_array[static_cast<int>(MPL::Variable::mechanical_strain)]
+ .emplace<KV>(KV::Zero());
+ variable_array[static_cast<int>(MPL::Variable::temperature)]
+ .emplace<double>(temperature);
+
+ variable_array_prev[static_cast<int>(MPL::Variable::stress)]
+ .emplace<KV>(KV::Zero());
+ variable_array_prev[static_cast<int>(MPL::Variable::mechanical_strain)]
+ .emplace<KV>(KV::Zero());
+ variable_array_prev[static_cast<int>(MPL::Variable::temperature)]
+ .emplace<double>(temperature);
+
+ auto&& solution =
+ solid_material.integrateStress(variable_array_prev, variable_array,
+ t, x_position, dt, *null_state);
+
+ if (!solution)
+ {
+ OGS_FATAL("Computation of elastic tangent stiffness failed.");
+ }
+
+ MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C =
+ std::move(std::get<2>(*solution));
+
+ return C;
+ }
+
+ typename BMatricesType::KelvinMatrixType updateConstitutiveRelation(
+ MaterialPropertyLib::VariableArray const& variable_array,
+ double const t,
+ ParameterLib::SpatialPosition const& x_position,
+ double const dt,
+ double const temperature)
+ {
+ MaterialPropertyLib::VariableArray variable_array_prev;
+ variable_array_prev[static_cast<int>(
+ MaterialPropertyLib::Variable::stress)]
+ .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
+ sigma_eff_prev);
+ variable_array_prev[static_cast<int>(MaterialPropertyLib::Variable::
+ mechanical_strain)]
+ .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
+ eps_m_prev);
+ variable_array_prev[static_cast<int>(
+ MaterialPropertyLib::Variable::temperature)]
+ .emplace<double>(temperature);
+
+ auto&& solution = solid_material.integrateStress(
+ variable_array_prev, variable_array, t, x_position, dt,
+ *material_state_variables);
+
+ if (!solution)
+ {
+ OGS_FATAL("Computation of local constitutive relation failed.");
+ }
+
+ MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C;
+ std::tie(sigma_eff, material_state_variables, C) = std::move(*solution);
+
+ return C;
+ }
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+} // namespace ThermoRichardsMechanics
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
new file mode 100644
index 00000000000..2e6df6bcf8f
--- /dev/null
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
@@ -0,0 +1,1443 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file
+ * Created on November 29, 2017, 2:03 PM
+ */
+
+#pragma once
+
+#include <cassert>
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
+#include "MaterialLib/MPL/Utils/GetLiquidThermalExpansivity.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
+#include "MathLib/KelvinVector.h"
+#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsMechanics
+{
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ ThermoRichardsMechanicsLocalAssembler(
+ MeshLib::Element const& e,
+ std::size_t const /*local_matrix_size*/,
+ bool const is_axially_symmetric,
+ unsigned const integration_order,
+ ThermoRichardsMechanicsProcessData<DisplacementDim>& process_data)
+ : _process_data(process_data),
+ _integration_method(integration_order),
+ _element(e),
+ _is_axially_symmetric(is_axially_symmetric)
+{
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ _ip_data.reserve(n_integration_points);
+ _secondary_data.N_u.resize(n_integration_points);
+
+ auto const shape_matrices_u =
+ NumLib::initShapeMatrices<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement,
+ DisplacementDim>(e, is_axially_symmetric,
+ _integration_method);
+
+ auto const shape_matrices_p =
+ NumLib::initShapeMatrices<ShapeFunctionPressure,
+ ShapeMatricesTypePressure, DisplacementDim>(
+ e, is_axially_symmetric, _integration_method);
+
+ auto const& solid_material =
+ MaterialLib::Solids::selectSolidConstitutiveRelation(
+ _process_data.solid_materials, _process_data.material_ids,
+ e.getID());
+
+ auto const& medium = _process_data.media_map->getMedium(_element.getID());
+
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ _ip_data.emplace_back(solid_material);
+ auto& ip_data = _ip_data[ip];
+ auto const& sm_u = shape_matrices_u[ip];
+ _ip_data[ip].integration_weight =
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ sm_u.integralMeasure * sm_u.detJ;
+
+ ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
+ DisplacementDim, displacement_size>::Zero(DisplacementDim,
+ displacement_size);
+ for (int i = 0; i < DisplacementDim; ++i)
+ {
+ ip_data.N_u_op
+ .template block<1, displacement_size / DisplacementDim>(
+ i, i * displacement_size / DisplacementDim)
+ .noalias() = sm_u.N;
+ }
+
+ ip_data.N_u = sm_u.N;
+ ip_data.dNdx_u = sm_u.dNdx;
+
+ ip_data.N_p = shape_matrices_p[ip].N;
+ ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
+
+ // Initial porosity. Could be read from intergration point data or mesh.
+ ip_data.porosity =
+ medium->property(MPL::porosity)
+ .template initialValue<double>(
+ x_position,
+ std::numeric_limits<
+ double>::quiet_NaN() /* t independent */);
+
+ ip_data.transport_porosity = ip_data.porosity;
+ if (medium->hasProperty(MPL::PropertyType::transport_porosity))
+ {
+ ip_data.transport_porosity =
+ medium->property(MPL::transport_porosity)
+ .template initialValue<double>(
+ x_position,
+ std::numeric_limits<
+ double>::quiet_NaN() /* t independent */);
+ }
+
+ _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
+ }
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::size_t ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::setIPDataInitialConditions(std::string const& name,
+ double const* values,
+ int const integration_order)
+{
+ 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_ip")
+ {
+ if (_process_data.initial_stress != nullptr)
+ {
+ OGS_FATAL(
+ "Setting initial conditions for stress from integration "
+ "point data and from a parameter '{:s}' is not possible "
+ "simultaneously.",
+ _process_data.initial_stress->name);
+ }
+ return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
+ values, _ip_data, &IpData::sigma_eff);
+ }
+
+ if (name == "saturation_ip")
+ {
+ return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
+ &IpData::saturation);
+ }
+ if (name == "porosity_ip")
+ {
+ return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
+ &IpData::porosity);
+ }
+ if (name == "transport_porosity_ip")
+ {
+ return ProcessLib::setIntegrationPointScalarData(
+ values, _ip_data, &IpData::transport_porosity);
+ }
+ if (name == "swelling_stress_ip")
+ {
+ return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
+ values, _ip_data, &IpData::sigma_sw);
+ }
+ if (name == "epsilon_ip")
+ {
+ return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
+ values, _ip_data, &IpData::eps);
+ }
+ return 0;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure,
+ IntegrationMethod, DisplacementDim>::
+ setInitialConditionsConcrete(std::vector<double> const& local_x,
+ double const t,
+ bool const /*use_monolithic_scheme*/,
+ int const /*process_id*/)
+{
+ assert(local_x.size() ==
+ temperature_size + pressure_size + displacement_size);
+
+ auto p_L =
+ Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+ pressure_size> const>(local_x.data() + pressure_index,
+ pressure_size);
+
+ auto const& medium = _process_data.media_map->getMedium(_element.getID());
+ MPL::VariableArray variables;
+
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+
+ auto const& N_p = _ip_data[ip].N_p;
+
+ double p_cap_ip;
+ NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
+
+ variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ p_cap_ip;
+ variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
+
+ // Note: temperature dependent saturation model is not considered so
+ // far.
+ _ip_data[ip].saturation_prev =
+ medium->property(MPL::PropertyType::saturation)
+ .template value<double>(
+ variables, x_position, t,
+ std::numeric_limits<double>::quiet_NaN());
+ }
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure,
+ IntegrationMethod, DisplacementDim>::
+ assembleWithJacobian(double const t, double const dt,
+ std::vector<double> const& local_x,
+ std::vector<double> const& local_xdot,
+ const double /*dxdot_dx*/, const double /*dx_dx*/,
+ std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& local_rhs_data,
+ std::vector<double>& local_Jac_data)
+{
+ auto const local_matrix_dim =
+ displacement_size + pressure_size + temperature_size;
+ assert(local_x.size() == local_matrix_dim);
+
+ auto const T =
+ Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+ temperature_size> const>(local_x.data() + temperature_index,
+ temperature_size);
+ auto const p_L =
+ Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+ pressure_size> const>(local_x.data() + pressure_index,
+ pressure_size);
+
+ auto const u =
+ Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
+ displacement_size> const>(local_x.data() + displacement_index,
+ displacement_size);
+
+ auto const T_dot =
+ Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+ temperature_size> const>(local_xdot.data() + temperature_index,
+ temperature_size);
+ auto const p_L_dot =
+ Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+ pressure_size> const>(local_xdot.data() + pressure_index,
+ pressure_size);
+ auto const u_dot =
+ Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
+ displacement_size> const>(local_xdot.data() + displacement_index,
+ displacement_size);
+
+ auto local_Jac = MathLib::createZeroedMatrix<
+ typename ShapeMatricesTypeDisplacement::template MatrixType<
+ local_matrix_dim, local_matrix_dim>>(
+ local_Jac_data, local_matrix_dim, local_matrix_dim);
+
+ auto local_rhs =
+ MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
+ template VectorType<local_matrix_dim>>(
+ local_rhs_data, local_matrix_dim);
+
+ auto const& identity2 = MathLib::KelvinVector::Invariants<
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value>::
+ identity2;
+
+ typename ShapeMatricesTypePressure::NodalMatrixType M_TT =
+ ShapeMatricesTypePressure::NodalMatrixType::Zero(temperature_size,
+ temperature_size);
+ typename ShapeMatricesTypePressure::NodalMatrixType K_TT =
+ ShapeMatricesTypePressure::NodalMatrixType::Zero(temperature_size,
+ temperature_size);
+ typename ShapeMatricesTypePressure::NodalMatrixType K_Tp =
+ ShapeMatricesTypePressure::NodalMatrixType::Zero(temperature_size,
+ pressure_size);
+ typename ShapeMatricesTypePressure::NodalMatrixType M_pT =
+ ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+ temperature_size);
+ typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
+ ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+ pressure_size);
+
+ typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_p =
+ ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+ pressure_size);
+
+ typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S_Jpp =
+ ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+ pressure_size);
+
+ typename ShapeMatricesTypePressure::NodalMatrixType storage_p_a_S =
+ ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+ pressure_size);
+
+ typename ShapeMatricesTypeDisplacement::template MatrixType<
+ displacement_size, pressure_size>
+ Kup = ShapeMatricesTypeDisplacement::template MatrixType<
+ displacement_size, pressure_size>::Zero(displacement_size,
+ pressure_size);
+
+ typename ShapeMatricesTypeDisplacement::template MatrixType<
+ pressure_size, displacement_size>
+ Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
+ pressure_size, displacement_size>::Zero(pressure_size,
+ displacement_size);
+
+ auto const& medium = _process_data.media_map->getMedium(_element.getID());
+ auto const& liquid_phase = medium->phase("AqueousLiquid");
+ auto const& solid_phase = medium->phase("Solid");
+ MPL::VariableArray variables;
+ MPL::VariableArray variables_prev;
+
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& w = _ip_data[ip].integration_weight;
+
+ auto const& N_u_op = _ip_data[ip].N_u_op;
+
+ auto const& N_u = _ip_data[ip].N_u;
+ auto const& dNdx_u = _ip_data[ip].dNdx_u;
+
+ auto const& N_p = _ip_data[ip].N_p;
+ auto const& dNdx_p = _ip_data[ip].dNdx_p;
+
+ auto const x_coord =
+ NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(
+ _element, N_u);
+ auto const B =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunctionDisplacement::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx_u, N_u, x_coord, _is_axially_symmetric);
+
+ double T_ip;
+ NumLib::shapeFunctionInterpolate(T, N_p, T_ip);
+ double T_dot_ip;
+ NumLib::shapeFunctionInterpolate(T_dot, N_p, T_dot_ip);
+
+ double p_cap_ip;
+ NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
+
+ double p_cap_dot_ip;
+ NumLib::shapeFunctionInterpolate(-p_L_dot, N_p, p_cap_dot_ip);
+
+ variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ p_cap_ip;
+ variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
+ variables[static_cast<int>(MPL::Variable::temperature)] = T_ip;
+
+ auto& eps = _ip_data[ip].eps;
+ auto const& sigma_eff = _ip_data[ip].sigma_eff;
+ auto& S_L = _ip_data[ip].saturation;
+ auto const S_L_prev = _ip_data[ip].saturation_prev;
+ auto const alpha =
+ medium->property(MPL::PropertyType::biot_coefficient)
+ .template value<double>(variables, x_position, t, dt);
+
+ auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
+ t, x_position, dt, T_ip);
+
+ auto const beta_SR =
+ (1 - alpha) /
+ _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
+ variables[static_cast<int>(MPL::Variable::grain_compressibility)] =
+ beta_SR;
+
+ auto const K_LR =
+ liquid_phase.property(MPL::PropertyType::bulk_modulus)
+ .template value<double>(variables, x_position, t, dt);
+
+ auto const rho_LR =
+ liquid_phase.property(MPL::PropertyType::density)
+ .template value<double>(variables, x_position, t, dt);
+ auto const& b = _process_data.specific_body_force;
+
+ S_L = medium->property(MPL::PropertyType::saturation)
+ .template value<double>(variables, x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
+ variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L_prev;
+
+ double const dS_L_dp_cap =
+ medium->property(MPL::PropertyType::saturation)
+ .template dValue<double>(variables,
+ MPL::Variable::capillary_pressure,
+ x_position, t, dt);
+
+ auto const chi = [medium, x_position, t, dt](double const S_L) {
+ MPL::VariableArray variables;
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
+ return medium->property(MPL::PropertyType::bishops_effective_stress)
+ .template value<double>(variables, x_position, t, dt);
+ };
+ double const chi_S_L = chi(S_L);
+ double const chi_S_L_prev = chi(S_L_prev);
+
+ variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L * p_cap_ip;
+ variables_prev[static_cast<int>(
+ MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
+
+ // Set volumetric strain rate for the general case without swelling.
+ variables[static_cast<int>(MPL::Variable::volumetric_strain)]
+ .emplace<double>(identity2.transpose() * B * u);
+ variables_prev[static_cast<int>(MPL::Variable::volumetric_strain)]
+ .emplace<double>(identity2.transpose() * B * (u - u_dot * dt));
+
+ auto& phi = _ip_data[ip].porosity;
+ { // Porosity update
+
+ variables_prev[static_cast<int>(MPL::Variable::porosity)] =
+ _ip_data[ip].porosity_prev;
+ phi = medium->property(MPL::PropertyType::porosity)
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::porosity)] = phi;
+ }
+
+ if (alpha < phi)
+ {
+ OGS_FATAL(
+ "ThermoRichardsMechanics: Biot-coefficient {} is smaller than "
+ "porosity {} in element/integration point {}/{}.",
+ alpha, phi, _element.getID(), ip);
+ }
+
+ // Swelling and possibly volumetric strain rate update.
+ auto& sigma_sw = _ip_data[ip].sigma_sw;
+ {
+ auto const& sigma_sw_prev = _ip_data[ip].sigma_sw_prev;
+
+ // If there is swelling, compute it. Update volumetric strain rate,
+ // s.t. it corresponds to the mechanical part only.
+ sigma_sw = sigma_sw_prev;
+ if (solid_phase.hasProperty(
+ MPL::PropertyType::swelling_stress_rate))
+ {
+ using DimMatrix = Eigen::Matrix<double, 3, 3>;
+ auto const sigma_sw_dot =
+ MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
+ solid_phase
+ .property(MPL::PropertyType::swelling_stress_rate)
+ .template value<DimMatrix>(
+ variables, variables_prev, x_position, t, dt));
+ sigma_sw += sigma_sw_dot * dt;
+
+ // !!! Misusing volumetric strain for mechanical volumetric
+ // strain just to update the transport porosity !!!
+ std::get<double>(variables[static_cast<int>(
+ MPL::Variable::volumetric_strain)]) +=
+ identity2.transpose() * C_el.inverse() * sigma_sw;
+ std::get<double>(variables_prev[static_cast<int>(
+ MPL::Variable::volumetric_strain)]) +=
+ identity2.transpose() * C_el.inverse() * sigma_sw_prev;
+ }
+
+ if (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
+ {
+ variables_prev[static_cast<int>(
+ MPL::Variable::transport_porosity)] =
+ _ip_data[ip].transport_porosity_prev;
+
+ _ip_data[ip].transport_porosity =
+ solid_phase.property(MPL::PropertyType::transport_porosity)
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::transport_porosity)] =
+ _ip_data[ip].transport_porosity;
+ }
+ else
+ {
+ variables[static_cast<int>(MPL::Variable::transport_porosity)] =
+ phi;
+ }
+ }
+
+ double const k_rel =
+ medium->property(MPL::PropertyType::relative_permeability)
+ .template value<double>(variables, x_position, t, dt);
+ auto const mu =
+ liquid_phase.property(MPL::PropertyType::viscosity)
+ .template value<double>(variables, x_position, t, dt);
+
+ // For stress dependent permeability.
+ variables[static_cast<int>(MPL::Variable::stress)]
+ .emplace<SymmetricTensor>(
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
+ (_ip_data[ip].sigma_eff +
+ alpha * chi_S_L * identity2 * p_cap_ip)
+ .eval()));
+ auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
+ medium->property(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;
+
+ //
+ // displacement equation, displacement part
+ //
+ eps.noalias() = B * u;
+
+ // 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
+ .property(
+ MaterialPropertyLib::PropertyType::thermal_expansivity)
+ .value(variables, x_position, t, dt));
+
+ MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
+ dthermal_strain =
+ MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
+ solid_linear_thermal_expansion_coefficient) *
+ T_dot_ip * dt;
+
+ variables[static_cast<int>(MPL::Variable::mechanical_strain)]
+ .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
+ eps - dthermal_strain);
+
+ auto C = _ip_data[ip].updateConstitutiveRelation(variables, t,
+ x_position, dt, T_ip);
+
+ local_Jac
+ .template block<displacement_size, displacement_size>(
+ displacement_index, displacement_index)
+ .noalias() += B.transpose() * C * B * w;
+
+ double const p_FR = -chi_S_L * p_cap_ip;
+ // p_SR
+ variables[static_cast<int>(MPL::Variable::solid_grain_pressure)] =
+ p_FR - (sigma_eff + sigma_sw).dot(identity2) / (3 * (1 - phi));
+ auto const rho_SR =
+ solid_phase.property(MPL::PropertyType::density)
+ .template value<double>(variables, x_position, t, dt);
+
+ double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
+ local_rhs.template segment<displacement_size>(displacement_index)
+ .noalias() -= (B.transpose() * (sigma_eff + sigma_sw) -
+ N_u_op.transpose() * rho * b) *
+ w;
+
+ //
+ // displacement equation, pressure part
+ //
+ Kup.noalias() += B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
+
+ auto const dchi_dS_L =
+ medium->property(MPL::PropertyType::bishops_effective_stress)
+ .template dValue<double>(variables,
+ MPL::Variable::liquid_saturation,
+ x_position, t, dt);
+ local_Jac
+ .template block<displacement_size, pressure_size>(
+ displacement_index, pressure_index)
+ .noalias() -= B.transpose() * alpha *
+ (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
+ identity2 * N_p * w;
+
+ local_Jac
+ .template block<displacement_size, pressure_size>(
+ displacement_index, pressure_index)
+ .noalias() +=
+ N_u_op.transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
+
+ if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
+ {
+ using DimMatrix = Eigen::Matrix<double, 3, 3>;
+ auto const dsigma_sw_dS_L =
+ MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
+ solid_phase
+ .property(MPL::PropertyType::swelling_stress_rate)
+ .template dValue<DimMatrix>(
+ variables, variables_prev,
+ MPL::Variable::liquid_saturation, x_position, t,
+ dt));
+ local_Jac
+ .template block<displacement_size, pressure_size>(
+ displacement_index, pressure_index)
+ .noalias() +=
+ B.transpose() * dsigma_sw_dS_L * dS_L_dp_cap * N_p * w;
+ }
+ //
+ // pressure equation, displacement part.
+ //
+ Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
+ identity2.transpose() * B * w;
+
+ //
+ // pressure equation, pressure part.
+ //
+ laplace_p.noalias() +=
+ dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
+
+ const double alphaB_minu_phi = alpha - phi;
+ double const a0 = alphaB_minu_phi * beta_SR;
+ double const specific_storage_a_p = S_L * (phi / K_LR + S_L * a0);
+ double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
+
+ double const dspecific_storage_a_p_dp_cap =
+ dS_L_dp_cap * (phi / K_LR + 2 * S_L * a0);
+ double const dspecific_storage_a_S_dp_cap =
+ -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
+
+ storage_p_a_p.noalias() +=
+ N_p.transpose() * rho_LR * specific_storage_a_p * N_p * w;
+ if (p_cap_dot_ip != 0) // prevent division by zero.
+ {
+ storage_p_a_S.noalias() -= N_p.transpose() * rho_LR *
+ specific_storage_a_S * (S_L - S_L_prev) /
+ (dt * p_cap_dot_ip) * N_p * w;
+ }
+
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += N_p.transpose() * p_cap_dot_ip * rho_LR *
+ dspecific_storage_a_p_dp_cap * N_p * w;
+
+ storage_p_a_S_Jpp.noalias() -=
+ N_p.transpose() * rho_LR *
+ ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
+ specific_storage_a_S * dS_L_dp_cap) /
+ dt * N_p * w;
+
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
+ identity2.transpose() * B * u_dot * N_p * w;
+
+ double const dk_rel_dS_l =
+ medium->property(MPL::PropertyType::relative_permeability)
+ .template dValue<double>(variables,
+ MPL::Variable::liquid_saturation,
+ x_position, t, dt);
+ GlobalDimVectorType const grad_p_cap = -dNdx_p * p_L;
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
+ dk_rel_dS_l * dS_L_dp_cap * N_p * w;
+
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
+ dk_rel_dS_l * dS_L_dp_cap * N_p * w;
+
+ local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
+ dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
+
+ //
+ // pressure equation, temperature part.
+ //
+ // Note: d (rho_l * beta _T)/dp * dotT
+ // Add the thermal expansion term is the point is in the fully
+ // saturated zone.
+ if (p_cap_ip <= 0.0) // p_l>0.0
+ {
+ double const fluid_volumetric_thermal_expansion_coefficient =
+ MPL::getLiquidThermalExpansivity(liquid_phase, variables,
+ rho_LR, x_position, t, dt);
+ const double eff_thermal_expansion =
+ alphaB_minu_phi *
+ solid_linear_thermal_expansion_coefficient.trace() +
+ phi * fluid_volumetric_thermal_expansion_coefficient;
+
+ M_pT.noalias() -=
+ N_p.transpose() * rho_LR * eff_thermal_expansion * N_p * w;
+ }
+
+ //
+ // temperature equation.
+ //
+ {
+ auto const specific_heat_capacity_fluid =
+ liquid_phase
+ .property(MaterialPropertyLib::specific_heat_capacity)
+ .template value<double>(variables, x_position, t, dt);
+
+ auto const specific_heat_capacity_solid =
+ solid_phase
+ .property(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_p.transpose() * N_p;
+
+ auto const thermal_conductivity_solid =
+ solid_phase
+ .property(
+ MaterialPropertyLib::PropertyType::thermal_conductivity)
+ .value(variables, x_position, t, dt);
+
+ auto const thermal_conductivity_fluid =
+ liquid_phase
+ .property(
+ MaterialPropertyLib::PropertyType::thermal_conductivity)
+ .template value<double>(variables, x_position, t, dt) *
+ S_L;
+
+ GlobalDimMatrixType const thermal_conductivity =
+ MaterialPropertyLib::formEffectiveThermalConductivity<
+ DisplacementDim>(thermal_conductivity_solid,
+ thermal_conductivity_fluid, phi);
+
+ GlobalDimVectorType const velocity_l =
+ GlobalDimVectorType(-Ki_over_mu * (dNdx_p * p_L - rho_LR * b));
+
+ K_TT.noalias() +=
+ (dNdx_p.transpose() * thermal_conductivity * dNdx_p +
+ N_p.transpose() * velocity_l.transpose() * dNdx_p * rho_LR *
+ specific_heat_capacity_fluid) *
+ w;
+
+ //
+ // temperature equation, pressure part
+ //
+ K_Tp.noalias() -= rho_LR * specific_heat_capacity_fluid *
+ N_p.transpose() * (dNdx_p * T).transpose() *
+ Ki_over_mu * dNdx_p * w;
+ }
+ }
+
+ if (_process_data.apply_mass_lumping)
+ {
+ storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
+ storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
+ storage_p_a_S_Jpp =
+ storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
+ }
+
+ //
+ // -- 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;
+
+ // 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;
+
+ // pressure equation, displacement part.
+ local_Jac
+ .template block<pressure_size, displacement_size>(pressure_index,
+ displacement_index)
+ .noalias() = Kpu / dt;
+
+ //
+ // -- Residual
+ //
+ // temperature equation
+ local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
+ M_TT * T_dot + K_TT * T;
+
+ // pressure equation
+ local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
+ laplace_p * p_L + (storage_p_a_p + storage_p_a_S) * p_L_dot +
+ Kpu * u_dot + M_pT * T_dot;
+
+ // displacement equation
+ local_rhs.template segment<displacement_size>(displacement_index)
+ .noalias() += Kup * p_L;
+}
+
+template <int Components, typename StoreValuesFunction>
+std::vector<double> transposeInPlace(
+ StoreValuesFunction const& store_values_function)
+{
+ std::vector<double> result;
+ store_values_function(result);
+ // Transpose. For time being Eigen's transposeInPlace doesn't work for
+ // non-square mapped matrices.
+ MathLib::toMatrix<
+ Eigen::Matrix<double, Eigen::Dynamic, Components, Eigen::RowMajor>>(
+ result, result.size() / Components, Components) =
+ MathLib::toMatrix<
+ Eigen::Matrix<double, Components, Eigen::Dynamic, Eigen::RowMajor>>(
+ result, Components, result.size() / Components)
+ .transpose()
+ .eval();
+
+ return result;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::getSigma() const
+{
+ constexpr int kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+
+ return transposeInPlace<kelvin_vector_size>(
+ [this](std::vector<double>& values) {
+ return getIntPtSigma(0, {}, {}, values);
+ });
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ getIntPtSigma(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
+ _ip_data, &IpData::sigma_eff, cache);
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::getSwellingStress() const
+{
+ constexpr int kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+
+ return transposeInPlace<kelvin_vector_size>(
+ [this](std::vector<double>& values) {
+ return getIntPtSwellingStress(0, {}, {}, values);
+ });
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ getIntPtSwellingStress(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ constexpr int kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+ auto const n_integration_points = _ip_data.size();
+
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, kelvin_vector_size, n_integration_points);
+
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ auto const& sigma_sw = _ip_data[ip].sigma_sw;
+ cache_mat.col(ip) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_sw);
+ }
+
+ return cache;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::getEpsilon() const
+{
+ constexpr int kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+
+ return transposeInPlace<kelvin_vector_size>(
+ [this](std::vector<double>& values) {
+ return getIntPtEpsilon(0, {}, {}, values);
+ });
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ getIntPtEpsilon(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
+ _ip_data, &IpData::eps, cache);
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ getIntPtDarcyVelocity(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ cache.clear();
+ auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, DisplacementDim, n_integration_points);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ cache_matrix.col(ip).noalias() = _ip_data[ip].v_darcy;
+ }
+
+ return cache;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::getSaturation() const
+{
+ std::vector<double> result;
+ getIntPtSaturation(0, {}, {}, result);
+ return result;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ getIntPtSaturation(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ return ProcessLib::getIntegrationPointScalarData(
+ _ip_data, &IpData::saturation, cache);
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::getPorosity() const
+{
+ std::vector<double> result;
+ getIntPtPorosity(0, {}, {}, result);
+ return result;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ getIntPtPorosity(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ return ProcessLib::getIntegrationPointScalarData(_ip_data,
+ &IpData::porosity, cache);
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::getTransportPorosity() const
+{
+ std::vector<double> result;
+ getIntPtTransportPorosity(0, {}, {}, result);
+ return result;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ getIntPtTransportPorosity(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ return ProcessLib::getIntegrationPointScalarData(
+ _ip_data, &IpData::transport_porosity, cache);
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ getIntPtDryDensitySolid(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ return ProcessLib::getIntegrationPointScalarData(
+ _ip_data, &IpData::dry_density_solid, cache);
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure,
+ IntegrationMethod, DisplacementDim>::
+ postNonLinearSolverConcrete(std::vector<double> const& local_x,
+ std::vector<double> const& local_xdot,
+ double const t, double const dt,
+ bool const use_monolithic_scheme,
+ int const /*process_id*/)
+{
+ const int displacement_offset =
+ use_monolithic_scheme ? displacement_index : 0;
+
+ auto const T =
+ Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+ temperature_size> const>(local_x.data() + temperature_index,
+ temperature_size);
+ auto const T_dot =
+ Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+ temperature_size> const>(local_xdot.data() + temperature_index,
+ temperature_size);
+
+ auto const u =
+ Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
+ displacement_size> const>(local_x.data() + displacement_offset,
+ displacement_size);
+ auto const& medium = _process_data.media_map->getMedium(_element.getID());
+ MPL::VariableArray variables;
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+ auto const& solid_phase = medium->phase("Solid");
+
+ int const n_integration_points = _integration_method.getNumberOfPoints();
+ for (int ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& N_p = _ip_data[ip].N_p;
+ auto const& N_u = _ip_data[ip].N_u;
+ auto const& dNdx_u = _ip_data[ip].dNdx_u;
+ double T_ip;
+ NumLib::shapeFunctionInterpolate(T, N_p, T_ip);
+ double T_dot_ip;
+ NumLib::shapeFunctionInterpolate(T_dot, N_p, T_dot_ip);
+ variables[static_cast<int>(MPL::Variable::temperature)] = T_ip;
+
+ auto const x_coord =
+ NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(
+ _element, N_u);
+ auto const B =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunctionDisplacement::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx_u, N_u, x_coord, _is_axially_symmetric);
+
+ auto& eps = _ip_data[ip].eps;
+ eps.noalias() = B * u;
+ // 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
+ .property(
+ MaterialPropertyLib::PropertyType::thermal_expansivity)
+ .value(variables, x_position, t, dt));
+
+ MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
+ dthermal_strain =
+ MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
+ solid_linear_thermal_expansion_coefficient) *
+ T_dot_ip * dt;
+
+ variables[static_cast<int>(MPL::Variable::mechanical_strain)]
+ .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
+ eps - dthermal_strain);
+
+ _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
+ T_ip);
+ }
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure,
+ IntegrationMethod, DisplacementDim>::
+ computeSecondaryVariableConcrete(double const t, double const dt,
+ Eigen::VectorXd const& local_x,
+ Eigen::VectorXd const& local_x_dot)
+{
+ auto const T =
+ local_x.template segment<temperature_size>(temperature_index);
+ auto const p_L = local_x.template segment<pressure_size>(pressure_index);
+ auto const u =
+ local_x.template segment<displacement_size>(displacement_index);
+
+ auto const T_dot =
+ local_x_dot.template segment<temperature_size>(temperature_index);
+ auto const p_L_dot =
+ local_x_dot.template segment<pressure_size>(pressure_index);
+ auto const u_dot =
+ local_x_dot.template segment<displacement_size>(displacement_index);
+
+ auto const& identity2 = MathLib::KelvinVector::Invariants<
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value>::
+ identity2;
+
+ auto const& medium = _process_data.media_map->getMedium(_element.getID());
+ auto const& liquid_phase = medium->phase("AqueousLiquid");
+ auto const& solid_phase = medium->phase("Solid");
+ MPL::VariableArray variables;
+ MPL::VariableArray variables_prev;
+
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ double saturation_avg = 0;
+ double porosity_avg = 0;
+
+ using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
+ KV sigma_avg = KV::Zero();
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& N_p = _ip_data[ip].N_p;
+ auto const& N_u = _ip_data[ip].N_u;
+ auto const& dNdx_u = _ip_data[ip].dNdx_u;
+
+ auto const x_coord =
+ NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(
+ _element, N_u);
+ auto const B =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunctionDisplacement::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx_u, N_u, x_coord, _is_axially_symmetric);
+
+ double T_ip;
+ NumLib::shapeFunctionInterpolate(T, N_p, T_ip);
+ double T_dot_ip;
+ NumLib::shapeFunctionInterpolate(T_dot, N_p, T_dot_ip);
+
+ double p_cap_ip;
+ NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
+
+ double p_cap_dot_ip;
+ NumLib::shapeFunctionInterpolate(-p_L_dot, N_p, p_cap_dot_ip);
+
+ variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ p_cap_ip;
+ variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
+
+ variables[static_cast<int>(MPL::Variable::temperature)] = T_ip;
+
+ auto& eps = _ip_data[ip].eps;
+ auto& S_L = _ip_data[ip].saturation;
+ auto const S_L_prev = _ip_data[ip].saturation_prev;
+ S_L = medium->property(MPL::PropertyType::saturation)
+ .template value<double>(variables, x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
+ variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L_prev;
+
+ auto const chi = [medium, x_position, t, dt](double const S_L) {
+ MPL::VariableArray variables;
+ variables.fill(std::numeric_limits<double>::quiet_NaN());
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
+ return medium->property(MPL::PropertyType::bishops_effective_stress)
+ .template value<double>(variables, x_position, t, dt);
+ };
+ double const chi_S_L = chi(S_L);
+ double const chi_S_L_prev = chi(S_L_prev);
+
+ auto const alpha =
+ medium->property(MPL::PropertyType::biot_coefficient)
+ .template value<double>(variables, x_position, t, dt);
+
+ auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
+ t, x_position, dt, T_ip);
+
+ auto const beta_SR =
+ (1 - alpha) /
+ _ip_data[ip].solid_material.getBulkModulus(t, x_position, &C_el);
+ variables[static_cast<int>(MPL::Variable::grain_compressibility)] =
+ beta_SR;
+
+ variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L * p_cap_ip;
+ variables_prev[static_cast<int>(
+ MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
+
+ // Set volumetric strain rate for the general case without swelling.
+ variables[static_cast<int>(MPL::Variable::volumetric_strain)]
+ .emplace<double>(identity2.transpose() * B * u);
+ variables_prev[static_cast<int>(MPL::Variable::volumetric_strain)]
+ .emplace<double>(identity2.transpose() * B * (u - u_dot * dt));
+
+ auto& phi = _ip_data[ip].porosity;
+ { // Porosity update
+ variables_prev[static_cast<int>(MPL::Variable::porosity)] =
+ _ip_data[ip].porosity_prev;
+ phi = medium->property(MPL::PropertyType::porosity)
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::porosity)] = phi;
+ }
+
+ // Swelling and possibly volumetric strain rate update.
+ auto& sigma_sw = _ip_data[ip].sigma_sw;
+ {
+ auto const& sigma_sw_prev = _ip_data[ip].sigma_sw_prev;
+
+ // If there is swelling, compute it. Update volumetric strain rate,
+ // s.t. it corresponds to the mechanical part only.
+ sigma_sw = sigma_sw_prev;
+ if (solid_phase.hasProperty(
+ MPL::PropertyType::swelling_stress_rate))
+ {
+ using DimMatrix = Eigen::Matrix<double, 3, 3>;
+ auto const sigma_sw_dot =
+ MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
+ solid_phase
+ .property(MPL::PropertyType::swelling_stress_rate)
+ .template value<DimMatrix>(
+ variables, variables_prev, x_position, t, dt));
+ sigma_sw += sigma_sw_dot * dt;
+
+ // !!! Misusing volumetric strain for mechanical volumetric
+ // strain just to update the transport porosity !!!
+ std::get<double>(variables[static_cast<int>(
+ MPL::Variable::volumetric_strain)]) +=
+ identity2.transpose() * C_el.inverse() * sigma_sw;
+ std::get<double>(variables_prev[static_cast<int>(
+ MPL::Variable::volumetric_strain)]) +=
+ identity2.transpose() * C_el.inverse() * sigma_sw_prev;
+ }
+
+ if (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
+ {
+ variables_prev[static_cast<int>(
+ MPL::Variable::transport_porosity)] =
+ _ip_data[ip].transport_porosity_prev;
+
+ _ip_data[ip].transport_porosity =
+ solid_phase.property(MPL::PropertyType::transport_porosity)
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::transport_porosity)] =
+ _ip_data[ip].transport_porosity;
+ }
+ else
+ {
+ variables[static_cast<int>(MPL::Variable::transport_porosity)] =
+ phi;
+ }
+ }
+ auto const mu =
+ liquid_phase.property(MPL::PropertyType::viscosity)
+ .template value<double>(variables, x_position, t, dt);
+ auto const rho_LR =
+ liquid_phase.property(MPL::PropertyType::density)
+ .template value<double>(variables, x_position, t, dt);
+
+ // For stress dependent permeability.
+ variables[static_cast<int>(MPL::Variable::stress)]
+ .emplace<SymmetricTensor>(
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
+ (_ip_data[ip].sigma_eff +
+ alpha * chi_S_L * identity2 * p_cap_ip)
+ .eval()));
+ auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
+ medium->property(MPL::PropertyType::permeability)
+ .value(variables, x_position, t, dt));
+
+ double const k_rel =
+ medium->property(MPL::PropertyType::relative_permeability)
+ .template value<double>(variables, x_position, t, dt);
+
+ GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
+
+ auto const& sigma_eff = _ip_data[ip].sigma_eff;
+ double const p_FR = -chi_S_L * p_cap_ip;
+ // p_SR
+ variables[static_cast<int>(MPL::Variable::solid_grain_pressure)] =
+ p_FR - (sigma_eff + sigma_sw).dot(identity2) / (3 * (1 - phi));
+ auto const rho_SR =
+ solid_phase.property(MPL::PropertyType::density)
+ .template value<double>(variables, x_position, t, dt);
+ _ip_data[ip].dry_density_solid = (1 - phi) * rho_SR;
+
+ eps.noalias() = B * u;
+
+ // 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
+ .property(
+ MaterialPropertyLib::PropertyType::thermal_expansivity)
+ .value(variables, x_position, t, dt));
+
+ MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
+ dthermal_strain =
+ MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
+ solid_linear_thermal_expansion_coefficient) *
+ T_dot_ip * dt;
+
+ variables[static_cast<int>(MPL::Variable::mechanical_strain)]
+ .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
+ eps - dthermal_strain);
+
+ _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
+ T_ip);
+
+ auto const& b = _process_data.specific_body_force;
+
+ // Compute the velocity
+ auto const& dNdx_p = _ip_data[ip].dNdx_p;
+ _ip_data[ip].v_darcy.noalias() =
+ -K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
+
+ saturation_avg += S_L;
+ porosity_avg += phi;
+ sigma_avg += sigma_eff;
+ }
+ saturation_avg /= n_integration_points;
+ porosity_avg /= n_integration_points;
+ sigma_avg /= n_integration_points;
+
+ (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
+ (*_process_data.element_porosity)[_element.getID()] = porosity_avg;
+
+ Eigen::Map<KV>(&(*_process_data.element_stresses)[_element.getID() *
+ KV::RowsAtCompileTime]) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_avg);
+
+ NumLib::interpolateToHigherOrderNodes<
+ ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
+ DisplacementDim>(_element, _is_axially_symmetric, p_L,
+ *_process_data.pressure_interpolated);
+ NumLib::interpolateToHigherOrderNodes<
+ ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
+ DisplacementDim>(_element, _is_axially_symmetric, T,
+ *_process_data.temperature_interpolated);
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+unsigned ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::getNumberOfIntegrationPoints() const
+{
+ return _integration_method.getNumberOfPoints();
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+typename MaterialLib::Solids::MechanicsBase<
+ DisplacementDim>::MaterialStateVariables const&
+ThermoRichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::getMaterialStateVariablesAt(unsigned integration_point)
+ const
+{
+ return *_ip_data[integration_point].material_state_variables;
+}
+} // namespace ThermoRichardsMechanics
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h
new file mode 100644
index 00000000000..a0f1cfddc24
--- /dev/null
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h
@@ -0,0 +1,258 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include <vector>
+
+#include "IntegrationPointData.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MathLib/KelvinVector.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ParameterLib/Parameter.h"
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/RichardsMechanics/LocalAssemblerInterface.h"
+#include "ThermoRichardsMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsMechanics
+{
+namespace MPL = MaterialPropertyLib;
+
+/// Used for the extrapolation of the integration point values. It is ordered
+/// (and stored) by integration points.
+template <typename ShapeMatrixType>
+struct SecondaryData
+{
+ std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N_u;
+};
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+class ThermoRichardsMechanicsLocalAssembler
+ : public RichardsMechanics::LocalAssemblerInterface<DisplacementDim>
+{
+public:
+ using ShapeMatricesTypeDisplacement =
+ ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
+ // Note: temperature variable uses the same shape functions as that are used
+ // by pressure variable.
+ using ShapeMatricesTypePressure =
+ ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>;
+
+ using GlobalDimMatrixType =
+ typename ShapeMatricesTypePressure::GlobalDimMatrixType;
+ using GlobalDimVectorType =
+ typename ShapeMatricesTypePressure::GlobalDimVectorType;
+
+ using BMatricesType =
+ BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
+ using KelvinVectorType = typename BMatricesType::KelvinVectorType;
+
+ using IpData =
+ IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement,
+ ShapeMatricesTypePressure, DisplacementDim,
+ ShapeFunctionDisplacement::NPOINTS>;
+
+ static int const KelvinVectorSize =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+ using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+
+ using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
+
+ ThermoRichardsMechanicsLocalAssembler(
+ ThermoRichardsMechanicsLocalAssembler const&) = delete;
+ ThermoRichardsMechanicsLocalAssembler(
+ ThermoRichardsMechanicsLocalAssembler&&) = delete;
+
+ ThermoRichardsMechanicsLocalAssembler(
+ MeshLib::Element const& e,
+ std::size_t const /*local_matrix_size*/,
+ bool const is_axially_symmetric,
+ unsigned const integration_order,
+ ThermoRichardsMechanicsProcessData<DisplacementDim>& process_data);
+
+ /// \return the number of read integration points.
+ std::size_t setIPDataInitialConditions(
+ std::string const& name,
+ double const* values,
+ int const integration_order) override;
+
+ void setInitialConditionsConcrete(std::vector<double> const& local_x,
+ double const t,
+ bool const use_monolithic_scheme,
+ int const process_id) override;
+
+ void assembleWithJacobian(double const t, double const dt,
+ std::vector<double> const& local_x,
+ std::vector<double> const& local_xdot,
+ const double /*dxdot_dx*/, const double /*dx_dx*/,
+ std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& local_rhs_data,
+ std::vector<double>& local_Jac_data) override;
+
+ void initializeConcrete() override
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto& ip_data = _ip_data[ip];
+
+ /// Set initial stress from parameter.
+ if (_process_data.initial_stress != nullptr)
+ {
+ ParameterLib::SpatialPosition const x_position{
+ boost::none, _element.getID(), ip,
+ MathLib::Point3d(NumLib::interpolateCoordinates<
+ ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(
+ _element, ip_data.N_u))};
+
+ ip_data.sigma_eff =
+ MathLib::KelvinVector::symmetricTensorToKelvinVector<
+ DisplacementDim>((*_process_data.initial_stress)(
+ std::numeric_limits<
+ double>::quiet_NaN() /* time independent */,
+ x_position));
+ }
+
+ ip_data.pushBackState();
+ }
+ }
+
+ void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
+ double const /*t*/,
+ double const /*dt*/) override
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ _ip_data[ip].pushBackState();
+ }
+ }
+
+ void computeSecondaryVariableConcrete(
+ double const t, double const dt, Eigen::VectorXd const& local_x,
+ Eigen::VectorXd const& local_x_dot) override;
+
+ void postNonLinearSolverConcrete(std::vector<double> const& local_x,
+ std::vector<double> const& local_xdot,
+ double const t, double const dt,
+ bool const use_monolithic_scheme,
+ int const process_id) override;
+
+ Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+ const unsigned integration_point) const override
+ {
+ 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());
+ }
+
+ std::vector<double> getSigma() const override;
+
+ std::vector<double> const& getIntPtDarcyVelocity(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const override;
+
+ std::vector<double> getSaturation() const override;
+ std::vector<double> const& getIntPtSaturation(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const override;
+
+ std::vector<double> getPorosity() const override;
+ std::vector<double> const& getIntPtPorosity(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const override;
+
+ std::vector<double> getTransportPorosity() const override;
+ std::vector<double> const& getIntPtTransportPorosity(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const override;
+
+ std::vector<double> const& getIntPtSigma(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const override;
+
+ std::vector<double> getSwellingStress() const override;
+ std::vector<double> const& getIntPtSwellingStress(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const override;
+
+ std::vector<double> getEpsilon() const override;
+ std::vector<double> const& getIntPtEpsilon(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const override;
+
+ 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;
+
+private:
+ unsigned getNumberOfIntegrationPoints() const override;
+
+ typename MaterialLib::Solids::MechanicsBase<
+ DisplacementDim>::MaterialStateVariables const&
+ getMaterialStateVariablesAt(unsigned integration_point) const override;
+
+private:
+ ThermoRichardsMechanicsProcessData<DisplacementDim>& _process_data;
+
+ std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
+
+ IntegrationMethod _integration_method;
+ MeshLib::Element const& _element;
+ bool const _is_axially_symmetric;
+ SecondaryData<
+ typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
+ _secondary_data;
+
+ static const int temperature_index = 0;
+ static const int temperature_size = ShapeFunctionPressure::NPOINTS;
+ static const int pressure_index = temperature_size;
+ static const int pressure_size = ShapeFunctionPressure::NPOINTS;
+ static const int displacement_index = 2 * ShapeFunctionPressure::NPOINTS;
+ static const int displacement_size =
+ ShapeFunctionDisplacement::NPOINTS * DisplacementDim;
+};
+
+} // namespace ThermoRichardsMechanics
+} // namespace ProcessLib
+
+#include "ThermoRichardsMechanicsFEM-impl.h"
--
GitLab