diff --git a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h
index e909848e70bafabeab9d5a83ec2fda13c6958f49..7783b02fc7d93841a2b2ac54e9617f968e7ee2f8 100644
--- a/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h
@@ -19,15 +19,18 @@
#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Porosity.h"
#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Saturation.h"
#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/SolidCompressibilityData.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/SolidMechanics.h"
#include "SaturationSecantDerivative.h"
#include "StiffnessTensor.h"
namespace ProcessLib::RichardsMechanics
{
-// TODO directly declare these type aliases in Traits.h
/// Data whose state must be tracked by the TRM process.
template <int DisplacementDim>
-using StatefulData = std::tuple<>;
+using StatefulData = std::tuple<
+ StrainData<DisplacementDim>,
+ ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>;
template <int DisplacementDim>
using StatefulDataPrev = ProcessLib::ConstitutiveRelations::PrevStateOf<
diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h
index 1845dd7481a6b918efad44a5d62b9505603fa571..6363455b5f7a234607894a9f2a0dd369f349542e 100644
--- a/ProcessLib/RichardsMechanics/IntegrationPointData.h
+++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h
@@ -32,20 +32,14 @@ struct IntegrationPointData final
// Initialize current time step values
static const int kelvin_vector_size =
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
- 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 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;
@@ -80,9 +74,7 @@ struct IntegrationPointData final
void pushBackState()
{
- eps_prev = eps;
eps_m_prev = eps_m;
- sigma_eff_prev = sigma_eff;
sigma_sw_prev = sigma_sw;
saturation_prev = saturation;
saturation_m_prev = saturation_m;
@@ -134,12 +126,16 @@ struct IntegrationPointData final
double const t,
ParameterLib::SpatialPosition const& x_position,
double const dt,
- double const temperature)
+ double const temperature,
+ ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::EffectiveStressData<
+ DisplacementDim>& sigma_eff,
+ PrevState<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::EffectiveStressData<
+ DisplacementDim>> const& sigma_eff_prev)
{
MaterialPropertyLib::VariableArray variable_array_prev;
- variable_array_prev.stress
- .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
- sigma_eff_prev);
+ variable_array_prev.stress = sigma_eff_prev->sigma_eff;
variable_array_prev.mechanical_strain
.emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
eps_m_prev);
@@ -155,7 +151,8 @@ struct IntegrationPointData final
}
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C;
- std::tie(sigma_eff, material_state_variables, C) = std::move(*solution);
+ std::tie(sigma_eff.sigma_eff, material_state_variables, C) =
+ std::move(*solution);
return C;
}
diff --git a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
index 15b57eaafcd6dd6208e47040e50df3043fee3e56..529002d1f651b7f59957c2d24c19c659abbe5ee6 100644
--- a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
+++ b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
@@ -10,9 +10,14 @@
#pragma once
+#include "ConstitutiveRelations/ConstitutiveData.h"
#include "MaterialLib/SolidModels/MechanicsBase.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/Integration/GenericIntegrationMethod.h"
#include "ProcessLib/LocalAssemblerInterface.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/MaterialState.h"
+#include "RichardsMechanicsProcessData.h"
namespace ProcessLib
{
@@ -22,6 +27,41 @@ template <int DisplacementDim>
struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
public NumLib::ExtrapolatableElement
{
+ LocalAssemblerInterface(
+ MeshLib::Element const& e,
+ NumLib::GenericIntegrationMethod const& integration_method,
+ bool const is_axially_symmetric,
+ RichardsMechanicsProcessData<DisplacementDim>& process_data)
+ : process_data_(process_data),
+ integration_method_(integration_method),
+ element_(e),
+ is_axially_symmetric_(is_axially_symmetric)
+ {
+ unsigned const n_integration_points =
+ integration_method_.getNumberOfPoints();
+
+ current_states_.resize(n_integration_points);
+ prev_states_.resize(n_integration_points);
+ output_data_.resize(n_integration_points);
+
+ material_states_.reserve(n_integration_points);
+
+ 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)
+ {
+ material_states_.emplace_back(
+ solid_material.createMaterialStateVariables());
+
+ // Set initial strain field to zero.
+ std::get<StrainData<DisplacementDim>>(current_states_[ip]).eps =
+ KelvinVector<DisplacementDim>::Zero();
+ }
+ }
+
virtual std::size_t setIPDataInitialConditions(
std::string_view const name, double const* values,
int const integration_order) = 0;
@@ -116,7 +156,19 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
virtual typename MaterialLib::Solids::MechanicsBase<
DisplacementDim>::MaterialStateVariables const&
getMaterialStateVariablesAt(unsigned /*integration_point*/) const = 0;
-};
+protected:
+ RichardsMechanicsProcessData<DisplacementDim>& process_data_;
+ NumLib::GenericIntegrationMethod const& integration_method_;
+ MeshLib::Element const& element_;
+ bool const is_axially_symmetric_;
+
+ std::vector<StatefulData<DisplacementDim>> current_states_;
+ std::vector<StatefulDataPrev<DisplacementDim>> prev_states_;
+ std::vector<
+ ProcessLib::ThermoRichardsMechanics::MaterialStateData<DisplacementDim>>
+ material_states_;
+ std::vector<OutputData<DisplacementDim>> output_data_;
+};
} // namespace RichardsMechanics
} // namespace ProcessLib
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 75cac8331605503d0cf72986aa73323c7fe976f8..13728943eda598e0c24875a796e11d38409b1cb2 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -25,6 +25,7 @@
#include "NumLib/Fem/Interpolation.h"
#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
#include "ProcessLib/Utils/TransposeInPlace.h"
+#include "RichardsMechanicsFEM.h"
namespace ProcessLib
{
@@ -122,13 +123,11 @@ RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
NumLib::GenericIntegrationMethod const& integration_method,
bool const is_axially_symmetric,
RichardsMechanicsProcessData<DisplacementDim>& process_data)
- : _process_data(process_data),
- _integration_method(integration_method),
- _element(e),
- _is_axially_symmetric(is_axially_symmetric)
+ : LocalAssemblerInterface<DisplacementDim>{
+ e, integration_method, is_axially_symmetric, process_data}
{
unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
+ this->integration_method_.getNumberOfPoints();
_ip_data.reserve(n_integration_points);
_secondary_data.N_u.resize(n_integration_points);
@@ -137,22 +136,23 @@ RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
NumLib::initShapeMatrices<ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement,
DisplacementDim>(e, is_axially_symmetric,
- _integration_method);
+ this->integration_method_);
auto const shape_matrices_p =
NumLib::initShapeMatrices<ShapeFunctionPressure,
ShapeMatricesTypePressure, DisplacementDim>(
- e, is_axially_symmetric, _integration_method);
+ e, is_axially_symmetric, this->integration_method_);
auto const& solid_material =
MaterialLib::Solids::selectSolidConstitutiveRelation(
- _process_data.solid_materials, _process_data.material_ids,
- e.getID());
+ this->process_data_.solid_materials,
+ this->process_data_.material_ids, e.getID());
- auto const& medium = _process_data.media_map.getMedium(_element.getID());
+ auto const& medium =
+ this->process_data_.media_map.getMedium(this->element_.getID());
ParameterLib::SpatialPosition x_position;
- x_position.setElementID(_element.getID());
+ x_position.setElementID(this->element_.getID());
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
@@ -160,7 +160,7 @@ RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto& ip_data = _ip_data[ip];
auto const& sm_u = shape_matrices_u[ip];
_ip_data[ip].integration_weight =
- _integration_method.getWeightedPoint(ip).getWeight() *
+ this->integration_method_.getWeightedPoint(ip).getWeight() *
sm_u.integralMeasure * sm_u.detJ;
ip_data.N_u = sm_u.N;
@@ -201,27 +201,33 @@ std::size_t RichardsMechanicsLocalAssembler<
int const integration_order)
{
if (integration_order !=
- static_cast<int>(_integration_method.getIntegrationOrder()))
+ static_cast<int>(this->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());
+ this->element_.getID());
}
if (name == "sigma")
{
- if (_process_data.initial_stress != nullptr)
+ if (this->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);
+ this->process_data_.initial_stress->name);
}
return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
- values, _ip_data, &IpData::sigma_eff);
+ values, this->current_states_, [](auto& tuple) -> auto& {
+ return std::get<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>(
+ tuple)
+ .sigma_eff;
+ });
}
if (name == "saturation")
@@ -247,7 +253,9 @@ std::size_t RichardsMechanicsLocalAssembler<
if (name == "epsilon")
{
return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
- values, _ip_data, &IpData::eps);
+ values, this->current_states_, [](auto& tuple) -> auto& {
+ return std::get<StrainData<DisplacementDim>>(tuple).eps;
+ });
}
if (name.starts_with("material_state_variable_"))
{
@@ -288,16 +296,17 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const p_L = local_x.template segment<pressure_size>(pressure_index);
constexpr double dt = std::numeric_limits<double>::quiet_NaN();
- auto const& medium = _process_data.media_map.getMedium(_element.getID());
+ auto const& medium =
+ this->process_data_.media_map.getMedium(this->element_.getID());
MPL::VariableArray variables;
ParameterLib::SpatialPosition x_position;
- x_position.setElementID(_element.getID());
+ x_position.setElementID(this->element_.getID());
auto const& solid_phase = medium->phase("Solid");
unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
+ this->integration_method_.getNumberOfPoints();
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
@@ -338,7 +347,9 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
// restart.
auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
t, x_position, dt, temperature);
- auto& eps = _ip_data[ip].eps;
+ auto& eps =
+ std::get<StrainData<DisplacementDim>>(this->current_states_[ip])
+ .eps;
auto& sigma_sw = _ip_data[ip].sigma_sw;
_ip_data[ip].eps_m_prev.noalias() =
@@ -404,17 +415,18 @@ void RichardsMechanicsLocalAssembler<
MathLib::KelvinVector::kelvin_vector_dimensions(
DisplacementDim)>::identity2;
- auto const& medium = _process_data.media_map.getMedium(_element.getID());
+ auto const& medium =
+ this->process_data_.media_map.getMedium(this->element_.getID());
auto const& liquid_phase = medium->phase("AqueousLiquid");
auto const& solid_phase = medium->phase("Solid");
MPL::VariableArray variables;
MPL::VariableArray variables_prev;
ParameterLib::SpatialPosition x_position;
- x_position.setElementID(_element.getID());
+ x_position.setElementID(this->element_.getID());
unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
+ this->integration_method_.getNumberOfPoints();
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
@@ -429,16 +441,17 @@ void RichardsMechanicsLocalAssembler<
auto const x_coord =
NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement>(
- _element, N_u);
+ this->element_, N_u);
auto const B =
LinearBMatrix::computeBMatrix<DisplacementDim,
ShapeFunctionDisplacement::NPOINTS,
typename BMatricesType::BMatrixType>(
- dNdx_u, N_u, x_coord, _is_axially_symmetric);
+ dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
- auto& eps = _ip_data[ip].eps;
+ auto& eps =
+ std::get<StrainData<DisplacementDim>>(this->current_states_[ip]);
auto& eps_m = _ip_data[ip].eps_m;
- eps.noalias() = B * u;
+ eps.eps.noalias() = B * u;
auto& S_L = _ip_data[ip].saturation;
auto const S_L_prev = _ip_data[ip].saturation_prev;
@@ -476,7 +489,7 @@ void RichardsMechanicsLocalAssembler<
.template value<double>(variables, x_position, t, dt);
variables.density = rho_LR;
- auto const& b = _process_data.specific_body_force;
+ auto const& b = this->process_data_.specific_body_force;
S_L = medium->property(MPL::PropertyType::saturation)
.template value<double>(variables, x_position, t, dt);
@@ -511,7 +524,7 @@ void RichardsMechanicsLocalAssembler<
variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
// Set volumetric strain rate for the general case without swelling.
- variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
+ variables.volumetric_strain = Invariants::trace(eps.eps);
variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
auto& phi = _ip_data[ip].porosity;
@@ -528,7 +541,7 @@ void RichardsMechanicsLocalAssembler<
OGS_FATAL(
"RichardsMechanics: Biot-coefficient {} is smaller than "
"porosity {} in element/integration point {}/{}.",
- alpha, phi, _element.getID(), ip);
+ alpha, phi, this->element_.getID(), ip);
}
// Swelling and possibly volumetric strain rate update.
@@ -583,7 +596,12 @@ void RichardsMechanicsLocalAssembler<
.template value<double>(variables, x_position, t, dt);
auto const& sigma_sw = _ip_data[ip].sigma_sw;
- auto const& sigma_eff = _ip_data[ip].sigma_eff;
+ auto const& sigma_eff =
+ std::get<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>(
+ this->current_states_[ip])
+ .sigma_eff;
// Set mechanical variables for the intrinsic permeability model
// For stress dependent permeability.
@@ -612,14 +630,29 @@ void RichardsMechanicsLocalAssembler<
//
eps_m.noalias() =
solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
- ? eps + C_el.inverse() * sigma_sw
- : eps;
+ ? eps.eps + C_el.inverse() * sigma_sw
+ : eps.eps;
variables.mechanical_strain
.emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
eps_m);
- _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
- temperature);
+ {
+ auto& SD = this->current_states_[ip];
+ auto const& SD_prev = this->prev_states_[ip];
+ auto& sigma_eff =
+ std::get<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>(SD);
+ auto const& sigma_eff_prev = std::get<
+ PrevState<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>>(
+ SD_prev);
+
+ _ip_data[ip].updateConstitutiveRelation(variables, t, x_position,
+ dt, temperature, sigma_eff,
+ sigma_eff_prev);
+ }
// p_SR
variables.solid_grain_pressure =
@@ -677,7 +710,7 @@ void RichardsMechanicsLocalAssembler<
identity2.transpose() * B * w;
}
- if (_process_data.apply_mass_lumping)
+ if (this->process_data_.apply_mass_lumping)
{
auto Mpp = M.template block<pressure_size, pressure_size>(
pressure_index, pressure_index);
@@ -698,7 +731,9 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
MPL::VariableArray& variables, MPL::VariableArray& variables_prev,
MPL::Medium const* const medium, TemperatureData const T_data,
CapillaryPressureData<DisplacementDim> const& p_cap_data,
- ConstitutiveData<DisplacementDim>& CD)
+ ConstitutiveData<DisplacementDim>& CD,
+ StatefulData<DisplacementDim>& SD,
+ StatefulDataPrev<DisplacementDim> const& SD_prev)
{
auto const& liquid_phase = medium->phase("AqueousLiquid");
auto const& solid_phase = medium->phase("Solid");
@@ -711,7 +746,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
double const p_cap_ip = p_cap_data.p_cap;
double const p_cap_prev_ip = p_cap_data.p_cap_prev;
- auto& eps = ip_data.eps;
+ auto const& eps = std::get<StrainData<DisplacementDim>>(SD);
auto& eps_m = ip_data.eps_m;
auto& S_L = ip_data.saturation;
auto const S_L_prev = ip_data.saturation_prev;
@@ -783,10 +818,11 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
// Set volumetric strain rate for the general case without swelling.
- variables.volumetric_strain = Invariants::trace(ip_data.eps);
+ variables.volumetric_strain = Invariants::trace(eps.eps);
// TODO (CL) changed that, using eps_prev for the moment, not B * u_prev
// variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
- variables_prev.volumetric_strain = Invariants::trace(ip_data.eps_prev);
+ variables_prev.volumetric_strain = Invariants::trace(
+ std::get<PrevState<StrainData<DisplacementDim>>>(SD_prev)->eps);
auto& phi = ip_data.porosity;
{ // Porosity update
@@ -804,7 +840,8 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
OGS_FATAL(
"RichardsMechanics: Biot-coefficient {} is smaller than "
"porosity {} in element/integration point {}/{}.",
- alpha, phi, _element.getID(), -1 /*ip*/ /* TODO (CL) re-enable */);
+ alpha, phi, this->element_.getID(),
+ -1 /*ip*/ /* TODO (CL) re-enable */);
}
auto const mu = liquid_phase.property(MPL::PropertyType::viscosity)
@@ -815,7 +852,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
// Swelling and possibly volumetric strain rate update.
updateSwellingStressAndVolumetricStrain<DisplacementDim>(
ip_data, *medium, solid_phase, C_el, rho_LR, mu,
- _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
+ this->process_data_.micro_porosity_parameters, alpha, phi, p_cap_ip,
variables, variables_prev, x_position, t, dt);
if (medium->hasProperty(MPL::PropertyType::transport_porosity))
@@ -839,8 +876,14 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
// Set mechanical variables for the intrinsic permeability model
// For stress dependent permeability.
{
+ // TODO mechanical constitutive relation will be evaluated afterwards
auto const sigma_total =
- (ip_data.sigma_eff + alpha * p_FR * identity2).eval();
+ (std::get<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>(SD)
+ .sigma_eff +
+ alpha * p_FR * identity2)
+ .eval();
// For stress dependent permeability.
variables.total_stress.emplace<SymmetricTensor>(
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_total));
@@ -873,20 +916,38 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
eps_m.noalias() =
solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
- ? eps + C_el.inverse() * sigma_sw
- : eps;
+ ? eps.eps + C_el.inverse() * sigma_sw
+ : eps.eps;
variables.mechanical_strain
.emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
eps_m);
- auto C = ip_data.updateConstitutiveRelation(variables, t, x_position, dt,
- temperature);
-
- *std::get<StiffnessTensor<DisplacementDim>>(CD) = std::move(C);
+ {
+ auto& sigma_eff =
+ std::get<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>(SD);
+ auto const& sigma_eff_prev =
+ std::get<PrevState<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>>(
+ SD_prev);
+
+ auto C = ip_data.updateConstitutiveRelation(variables, t, x_position,
+ dt, temperature, sigma_eff,
+ sigma_eff_prev);
+
+ *std::get<StiffnessTensor<DisplacementDim>>(CD) = std::move(C);
+ }
// p_SR
variables.solid_grain_pressure =
- p_FR - ip_data.sigma_eff.dot(identity2) / (3 * (1 - phi));
+ p_FR -
+ std::get<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::EffectiveStressData<
+ DisplacementDim>>(SD)
+ .sigma_eff.dot(identity2) /
+ (3 * (1 - phi));
auto const rho_SR =
solid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
@@ -972,22 +1033,25 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
pressure_size, displacement_size>::Zero(pressure_size,
displacement_size);
- auto const& medium = _process_data.media_map.getMedium(_element.getID());
+ auto const& medium =
+ this->process_data_.media_map.getMedium(this->element_.getID());
auto const& liquid_phase = medium->phase("AqueousLiquid");
auto const& solid_phase = medium->phase("Solid");
MPL::VariableArray variables;
MPL::VariableArray variables_prev;
ParameterLib::SpatialPosition x_position;
- x_position.setElementID(_element.getID());
+ x_position.setElementID(this->element_.getID());
unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
+ this->integration_method_.getNumberOfPoints();
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
ConstitutiveData<DisplacementDim> CD;
+ auto& SD = this->current_states_[ip];
+ auto const& SD_prev = this->prev_states_[ip];
[[maybe_unused]] auto models = createConstitutiveModels(
- _process_data, _ip_data[ip].solid_material);
+ this->process_data_, _ip_data[ip].solid_material);
x_position.setIntegrationPoint(ip);
auto const& w = _ip_data[ip].integration_weight;
@@ -1001,12 +1065,12 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const x_coord =
NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement>(
- _element, N_u);
+ this->element_, N_u);
auto const B =
LinearBMatrix::computeBMatrix<DisplacementDim,
ShapeFunctionDisplacement::NPOINTS,
typename BMatricesType::BMatrixType>(
- dNdx_u, N_u, x_coord, _is_axially_symmetric);
+ dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
double p_cap_ip;
NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
@@ -1025,8 +1089,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
.template value<double>(variables, x_position, t, dt);
variables.temperature = temperature;
- auto& eps = _ip_data[ip].eps;
- eps.noalias() = B * u;
+ std::get<StrainData<DisplacementDim>>(SD).eps.noalias() = B * u;
assembleWithJacobianEvalConstitutiveSetting(
t, dt, x_position, _ip_data[ip], variables, variables_prev, medium,
@@ -1034,20 +1097,31 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
CapillaryPressureData<DisplacementDim>{
p_cap_ip, p_cap_prev_ip,
Eigen::Vector<double, DisplacementDim>::Zero()},
- CD);
+ CD, SD, SD_prev);
- auto const& C = *std::get<StiffnessTensor<DisplacementDim>>(CD);
- local_Jac
- .template block<displacement_size, displacement_size>(
- displacement_index, displacement_index)
- .noalias() += B.transpose() * C * B * w;
-
- auto const& sigma_eff = _ip_data[ip].sigma_eff;
- double const rho = *std::get<Density>(CD);
- auto const& b = _process_data.specific_body_force;
- local_rhs.template segment<displacement_size>(displacement_index)
- .noalias() -=
- (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
+ {
+ auto const& C = *std::get<StiffnessTensor<DisplacementDim>>(CD);
+ local_Jac
+ .template block<displacement_size, displacement_size>(
+ displacement_index, displacement_index)
+ .noalias() += B.transpose() * C * B * w;
+ }
+
+ auto const& b = this->process_data_.specific_body_force;
+
+ {
+ auto const& sigma_eff =
+ std::get<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>(
+ this->current_states_[ip])
+ .sigma_eff;
+ double const rho = *std::get<Density>(CD);
+ local_rhs.template segment<displacement_size>(displacement_index)
+ .noalias() -= (B.transpose() * sigma_eff -
+ N_u_op(N_u).transpose() * rho * b) *
+ w;
+ }
//
// displacement equation, pressure part
@@ -1055,24 +1129,28 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
double const alpha =
*std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD);
- double const chi_S_L =
- std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
- .chi_S_L;
- Kup.noalias() += B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
-
- double const dchi_dS_L =
- std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
- .dchi_dS_L;
double const dS_L_dp_cap =
std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(
CD)
.dS_L_dp_cap;
- local_Jac
- .template block<displacement_size, pressure_size>(
- displacement_index, pressure_index)
- .noalias() -= B.transpose() * alpha *
- (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
- identity2 * N_p * w;
+
+ {
+ double const chi_S_L =
+ std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
+ .chi_S_L;
+ Kup.noalias() +=
+ B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
+ double const dchi_dS_L =
+ std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
+ .dchi_dS_L;
+
+ local_Jac
+ .template block<displacement_size, pressure_size>(
+ displacement_index, pressure_index)
+ .noalias() -= B.transpose() * alpha *
+ (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
+ identity2 * N_p * w;
+ }
double const phi =
std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi;
@@ -1114,7 +1192,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
// pressure equation, displacement part.
//
double const S_L = _ip_data[ip].saturation;
- if (_process_data.explicit_hm_coupling_in_unsaturated_zone)
+ if (this->process_data_.explicit_hm_coupling_in_unsaturated_zone)
{
double const chi_S_L_prev = std::get<PrevState<
ProcessLib::ThermoRichardsMechanics::BishopsData>>(CD)
@@ -1192,7 +1270,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
specific_storage_a_S * dS_L_dp_cap) /
dt * N_p * w;
- if (!_process_data.explicit_hm_coupling_in_unsaturated_zone)
+ if (!this->process_data_.explicit_hm_coupling_in_unsaturated_zone)
{
local_Jac
.template block<pressure_size, pressure_size>(pressure_index,
@@ -1226,8 +1304,9 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
if (medium->hasProperty(MPL::PropertyType::saturation_micro))
{
- double const alpha_bar = _process_data.micro_porosity_parameters
- ->mass_exchange_coefficient;
+ double const alpha_bar =
+ this->process_data_.micro_porosity_parameters
+ ->mass_exchange_coefficient;
auto const p_L_m = _ip_data[ip].liquid_pressure_m;
local_rhs.template segment<pressure_size>(pressure_index)
.noalias() -=
@@ -1250,7 +1329,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
}
}
- if (_process_data.apply_mass_lumping)
+ if (this->process_data_.apply_mass_lumping)
{
storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
@@ -1306,7 +1385,14 @@ std::vector<double> const& RichardsMechanicsLocalAssembler<
std::vector<double>& cache) const
{
return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
- _ip_data, &IpData::sigma_eff, cache);
+ this->current_states_,
+ [](auto& tuple) -> auto& {
+ return std::get<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>(tuple)
+ .sigma_eff;
+ },
+ cache);
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
@@ -1377,7 +1463,11 @@ std::vector<double> const& RichardsMechanicsLocalAssembler<
std::vector<double>& cache) const
{
return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
- _ip_data, &IpData::eps, cache);
+ this->current_states_,
+ [](auto& tuple) -> auto& {
+ return std::get<StrainData<DisplacementDim>>(tuple).eps;
+ },
+ cache);
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
@@ -1386,9 +1476,9 @@ int RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
ShapeFunctionPressure,
DisplacementDim>::getMaterialID() const
{
- return _process_data.material_ids == nullptr
+ return this->process_data_.material_ids == nullptr
? 0
- : (*_process_data.material_ids)[_element.getID()];
+ : (*this->process_data_.material_ids)[this->element_.getID()];
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
@@ -1417,7 +1507,7 @@ std::vector<double> const& RichardsMechanicsLocalAssembler<
std::vector<double>& cache) const
{
unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
+ this->integration_method_.getNumberOfPoints();
cache.clear();
auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
@@ -1648,17 +1738,18 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
MathLib::KelvinVector::kelvin_vector_dimensions(
DisplacementDim)>::identity2;
- auto const& medium = _process_data.media_map.getMedium(_element.getID());
+ auto const& medium =
+ this->process_data_.media_map.getMedium(this->element_.getID());
auto const& liquid_phase = medium->phase("AqueousLiquid");
auto const& solid_phase = medium->phase("Solid");
MPL::VariableArray variables;
MPL::VariableArray variables_prev;
ParameterLib::SpatialPosition x_position;
- x_position.setElementID(_element.getID());
+ x_position.setElementID(this->element_.getID());
unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
+ this->integration_method_.getNumberOfPoints();
double saturation_avg = 0;
double porosity_avg = 0;
@@ -1676,12 +1767,12 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const x_coord =
NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement>(
- _element, N_u);
+ this->element_, N_u);
auto const B =
LinearBMatrix::computeBMatrix<DisplacementDim,
ShapeFunctionDisplacement::NPOINTS,
typename BMatricesType::BMatrixType>(
- dNdx_u, N_u, x_coord, _is_axially_symmetric);
+ dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
double p_cap_ip;
NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
@@ -1700,7 +1791,9 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
.template value<double>(variables, x_position, t, dt);
variables.temperature = temperature;
- auto& eps = _ip_data[ip].eps;
+ auto& eps =
+ std::get<StrainData<DisplacementDim>>(this->current_states_[ip])
+ .eps;
eps.noalias() = B * u;
auto& eps_m = _ip_data[ip].eps_m;
auto& S_L = _ip_data[ip].saturation;
@@ -1736,7 +1829,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
// Set volumetric strain rate for the general case without swelling.
- variables.volumetric_strain = Invariants::trace(_ip_data[ip].eps);
+ variables.volumetric_strain = Invariants::trace(eps);
variables_prev.volumetric_strain = Invariants::trace(B * u_prev);
auto& phi = _ip_data[ip].porosity;
@@ -1759,7 +1852,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
// Swelling and possibly volumetric strain rate update.
updateSwellingStressAndVolumetricStrain<DisplacementDim>(
_ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu,
- _process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
+ this->process_data_.micro_porosity_parameters, alpha, phi, p_cap_ip,
variables, variables_prev, x_position, t, dt);
if (medium->hasProperty(MPL::PropertyType::transport_porosity))
@@ -1781,12 +1874,18 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables.transport_porosity = phi;
}
+ auto const& sigma_eff =
+ std::get<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>(
+ this->current_states_[ip])
+ .sigma_eff;
+
// Set mechanical variables for the intrinsic permeability model
// For stress dependent permeability.
{
- auto const sigma_total = (_ip_data[ip].sigma_eff +
- alpha * chi_S_L * identity2 * p_cap_ip)
- .eval();
+ auto const sigma_total =
+ (sigma_eff + alpha * chi_S_L * identity2 * p_cap_ip).eval();
// For stress dependent permeability.
variables.total_stress.emplace<SymmetricTensor>(
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
@@ -1806,7 +1905,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
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.solid_grain_pressure =
@@ -1826,10 +1924,25 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
.emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
eps_m);
- _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
- temperature);
+ {
+ auto& SD = this->current_states_[ip];
+ auto const& SD_prev = this->prev_states_[ip];
+ auto& sigma_eff =
+ std::get<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>(SD);
+ auto const& sigma_eff_prev = std::get<
+ PrevState<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>>(
+ SD_prev);
+
+ _ip_data[ip].updateConstitutiveRelation(variables, t, x_position,
+ dt, temperature, sigma_eff,
+ sigma_eff_prev);
+ }
- auto const& b = _process_data.specific_body_force;
+ auto const& b = this->process_data_.specific_body_force;
// Compute the velocity
auto const& dNdx_p = _ip_data[ip].dNdx_p;
@@ -1844,17 +1957,20 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
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;
+ (*this->process_data_.element_saturation)[this->element_.getID()] =
+ saturation_avg;
+ (*this->process_data_.element_porosity)[this->element_.getID()] =
+ porosity_avg;
- Eigen::Map<KV>(&(*_process_data.element_stresses)[_element.getID() *
- KV::RowsAtCompileTime]) =
+ Eigen::Map<KV>(
+ &(*this->process_data_.element_stresses)[this->element_.getID() *
+ KV::RowsAtCompileTime]) =
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_avg);
NumLib::interpolateToHigherOrderNodes<
ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
- DisplacementDim>(_element, _is_axially_symmetric, p_L,
- *_process_data.pressure_interpolated);
+ DisplacementDim>(this->element_, this->is_axially_symmetric_, p_L,
+ *this->process_data_.pressure_interpolated);
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
@@ -1863,7 +1979,7 @@ unsigned RichardsMechanicsLocalAssembler<
ShapeFunctionDisplacement, ShapeFunctionPressure,
DisplacementDim>::getNumberOfIntegrationPoints() const
{
- return _integration_method.getNumberOfPoints();
+ return this->integration_method_.getNumberOfPoints();
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
index 1ad821a245edd5f6110d998c383124f03676e706..16b8210746220b487568525e539f34da4aa24392 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -125,25 +125,29 @@ public:
void initializeConcrete() override
{
unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
+ this->integration_method_.getNumberOfPoints();
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
+ auto& SD = this->current_states_[ip];
auto& ip_data = _ip_data[ip];
ParameterLib::SpatialPosition const x_position{
- std::nullopt, _element.getID(), ip,
- MathLib::Point3d(
- NumLib::interpolateCoordinates<
- ShapeFunctionDisplacement,
- ShapeMatricesTypeDisplacement>(_element, ip_data.N_u))};
+ std::nullopt, this->element_.getID(), ip,
+ MathLib::Point3d(NumLib::interpolateCoordinates<
+ ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(this->element_,
+ ip_data.N_u))};
/// Set initial stress from parameter.
- if (_process_data.initial_stress != nullptr)
+ if (this->process_data_.initial_stress != nullptr)
{
- ip_data.sigma_eff =
+ std::get<ProcessLib::ThermoRichardsMechanics::
+ ConstitutiveStress_StrainTemperature::
+ EffectiveStressData<DisplacementDim>>(SD)
+ .sigma_eff =
MathLib::KelvinVector::symmetricTensorToKelvinVector<
- DisplacementDim>((*_process_data.initial_stress)(
+ DisplacementDim>((*this->process_data_.initial_stress)(
std::numeric_limits<
double>::quiet_NaN() /* time independent */,
x_position));
@@ -154,6 +158,8 @@ public:
t, x_position, *ip_data.material_state_variables);
ip_data.pushBackState();
+
+ this->prev_states_[ip] = SD;
}
}
@@ -163,12 +169,18 @@ public:
int const /*process_id*/) override
{
unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
+ this->integration_method_.getNumberOfPoints();
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
_ip_data[ip].pushBackState();
}
+
+ // TODO move to the local assembler interface
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ this->prev_states_[ip] = this->current_states_[ip];
+ }
}
void computeSecondaryVariableConcrete(
@@ -336,15 +348,12 @@ private:
MPL::VariableArray& variables, MPL::VariableArray& variables_prev,
MPL::Medium const* const medium, TemperatureData const T_data,
CapillaryPressureData<DisplacementDim> const& p_cap_data,
- ConstitutiveData<DisplacementDim>& CD);
-
- RichardsMechanicsProcessData<DisplacementDim>& _process_data;
+ ConstitutiveData<DisplacementDim>& CD,
+ StatefulData<DisplacementDim>& SD,
+ StatefulDataPrev<DisplacementDim> const& SD_prev);
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;