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

Detailed Description

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

Definition at line 166 of file HydroMechanicsFEM.h.

#include <HydroMechanicsFEM.h>

Inheritance diagram for ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:
[legend]
Collaboration diagram for ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >:
[legend]

Public Types

using ShapeMatricesTypeDisplacement = ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim >
 
using ShapeMatricesTypePressure = ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim >
 
using Invariants = MathLib::KelvinVector::Invariants< KelvinVectorSize >
 
using SymmetricTensor = Eigen::Matrix< double, KelvinVectorSize, 1 >
 

Public Member Functions

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

Static Public Attributes

static int const KelvinVectorSize
 

Private Types

using BMatricesType = BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim >
 
using IpData = IntegrationPointData< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS >
 

Private Member Functions

void assembleWithJacobianForDeformationEquations (const double t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
 
void assembleWithJacobianForPressureEquations (const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
 
unsigned getNumberOfIntegrationPoints () const override
 
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt (unsigned integration_point) const override
 

Private Attributes

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

Static Private Attributes

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

Member Typedef Documentation

◆ BMatricesType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::BMatricesType = BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>
private

Definition at line 387 of file HydroMechanicsFEM.h.

◆ Invariants

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

Definition at line 179 of file HydroMechanicsFEM.h.

◆ IpData

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::IpData = IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, DisplacementDim, ShapeFunctionDisplacement::NPOINTS>
private

Definition at line 389 of file HydroMechanicsFEM.h.

◆ ShapeMatricesTypeDisplacement

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypeDisplacement = ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>

Definition at line 170 of file HydroMechanicsFEM.h.

◆ ShapeMatricesTypePressure

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
using ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::ShapeMatricesTypePressure = ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>

Definition at line 174 of file HydroMechanicsFEM.h.

◆ SymmetricTensor

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

Definition at line 181 of file HydroMechanicsFEM.h.

Constructor & Destructor Documentation

◆ HydroMechanicsLocalAssembler() [1/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::HydroMechanicsLocalAssembler ( HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > const &  )
delete

◆ HydroMechanicsLocalAssembler() [2/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::HydroMechanicsLocalAssembler ( HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim > &&  )
delete

◆ HydroMechanicsLocalAssembler() [3/3]

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

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

42  : _process_data(process_data),
43  _integration_method(integration_method),
44  _element(e),
45  _is_axially_symmetric(is_axially_symmetric)
46 {
47  unsigned const n_integration_points =
49 
50  _ip_data.reserve(n_integration_points);
51  _secondary_data.N_u.resize(n_integration_points);
52 
53  auto const shape_matrices_u =
54  NumLib::initShapeMatrices<ShapeFunctionDisplacement,
56  DisplacementDim>(e, is_axially_symmetric,
58 
59  auto const shape_matrices_p =
60  NumLib::initShapeMatrices<ShapeFunctionPressure,
61  ShapeMatricesTypePressure, DisplacementDim>(
62  e, is_axially_symmetric, _integration_method);
63 
64  auto const& solid_material =
66  _process_data.solid_materials, _process_data.material_ids,
67  e.getID());
68 
69  for (unsigned ip = 0; ip < n_integration_points; ip++)
70  {
71  _ip_data.emplace_back(solid_material);
72  auto& ip_data = _ip_data[ip];
73  auto const& sm_u = shape_matrices_u[ip];
74  ip_data.integration_weight =
76  sm_u.integralMeasure * sm_u.detJ;
77 
78  // Initialize current time step values
79  static const int kelvin_vector_size =
81  ip_data.sigma_eff.setZero(kelvin_vector_size);
82  ip_data.eps.setZero(kelvin_vector_size);
83 
84  // Previous time step values are not initialized and are set later.
85  ip_data.eps_prev.resize(kelvin_vector_size);
86  ip_data.sigma_eff_prev.resize(kelvin_vector_size);
87 
88  ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
89  DisplacementDim, displacement_size>::Zero(DisplacementDim,
91  for (int i = 0; i < DisplacementDim; ++i)
92  {
93  ip_data.N_u_op
94  .template block<1, displacement_size / DisplacementDim>(
95  i, i * displacement_size / DisplacementDim)
96  .noalias() = sm_u.N;
97  }
98 
99  ip_data.N_u = sm_u.N;
100  ip_data.dNdx_u = sm_u.dNdx;
101 
102  ip_data.N_p = shape_matrices_p[ip].N;
103  ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
104 
105  _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
106  }
107 }
double getWeight() const
Definition: WeightedPoint.h:80
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
NumLib::GenericIntegrationMethod const & _integration_method
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
HydroMechanicsProcessData< DisplacementDim > & _process_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
MechanicsBase< DisplacementDim > & selectSolidConstitutiveRelation(std::map< int, std::unique_ptr< MechanicsBase< DisplacementDim >>> const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:24
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u

References ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_integration_method, ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_process_data, ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data, ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::displacement_size, MeshLib::Element::getID(), NumLib::GenericIntegrationMethod::getNumberOfPoints(), MathLib::WeightedPoint::getWeight(), NumLib::GenericIntegrationMethod::getWeightedPoint(), NumLib::initShapeMatrices(), MathLib::KelvinVector::kelvin_vector_dimensions(), ProcessLib::HydroMechanics::SecondaryData< ShapeMatrixType >::N_u, and MaterialLib::Solids::selectSolidConstitutiveRelation().

Member Function Documentation

◆ assemble()

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

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 199 of file HydroMechanicsFEM.h.

205  {
206  OGS_FATAL(
207  "HydroMechanicsLocalAssembler: assembly without jacobian is not "
208  "implemented.");
209  }
#define OGS_FATAL(...)
Definition: Error.h:26

References OGS_FATAL.

◆ assembleWithJacobian()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobian ( double const  t,
double const  dt,
std::vector< double > const &  local_x,
std::vector< double > const &  local_xdot,
std::vector< double > &  ,
std::vector< double > &  ,
std::vector< double > &  local_rhs_data,
std::vector< double > &  local_Jac_data 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

120 {
121  assert(local_x.size() == pressure_size + displacement_size);
122 
123  auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
124  pressure_size> const>(local_x.data() + pressure_index, pressure_size);
125 
126  auto u =
127  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
128  displacement_size> const>(local_x.data() + displacement_index,
130 
131  auto p_dot =
132  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
133  pressure_size> const>(local_xdot.data() + pressure_index,
134  pressure_size);
135  auto u_dot =
136  Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
137  displacement_size> const>(local_xdot.data() + displacement_index,
139 
140  auto local_Jac = MathLib::createZeroedMatrix<
141  typename ShapeMatricesTypeDisplacement::template MatrixType<
144  local_Jac_data, displacement_size + pressure_size,
146 
147  auto local_rhs = MathLib::createZeroedVector<
148  typename ShapeMatricesTypeDisplacement::template VectorType<
150  local_rhs_data, displacement_size + pressure_size);
151 
152  typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
153  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
154  pressure_size);
155 
156  typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
157  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
158  pressure_size);
159 
160  typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
161  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
162  pressure_size);
163 
164  typename ShapeMatricesTypeDisplacement::template MatrixType<
166  Kup = ShapeMatricesTypeDisplacement::template MatrixType<
168  pressure_size);
169 
170  typename ShapeMatricesTypeDisplacement::template MatrixType<
172  Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
175 
176  typename ShapeMatricesTypeDisplacement::template MatrixType<
178  Kpu_k = ShapeMatricesTypeDisplacement::template MatrixType<
181 
182  auto const& solid_material =
184  _process_data.solid_materials, _process_data.material_ids,
185  _element.getID());
186 
188  x_position.setElementID(_element.getID());
189 
190  unsigned const n_integration_points =
192 
193  auto const& b = _process_data.specific_body_force;
194  auto const& medium = _process_data.media_map->getMedium(_element.getID());
195  auto const& solid = medium->phase("Solid");
196  auto const& fluid = fluidPhase(*medium);
197  MPL::VariableArray vars;
198 
199  auto const T_ref =
201  .template value<double>(vars, x_position, t, dt);
202  vars.temperature = T_ref;
203 
204  auto const& identity2 = Invariants::identity2;
205 
206  for (unsigned ip = 0; ip < n_integration_points; ip++)
207  {
208  x_position.setIntegrationPoint(ip);
209  auto const& w = _ip_data[ip].integration_weight;
210 
211  auto const& N_u_op = _ip_data[ip].N_u_op;
212 
213  auto const& N_u = _ip_data[ip].N_u;
214  auto const& dNdx_u = _ip_data[ip].dNdx_u;
215 
216  auto const& N_p = _ip_data[ip].N_p;
217  auto const& dNdx_p = _ip_data[ip].dNdx_p;
218 
219  auto const x_coord =
220  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
222  _element, N_u);
223  auto const B =
224  LinearBMatrix::computeBMatrix<DisplacementDim,
225  ShapeFunctionDisplacement::NPOINTS,
226  typename BMatricesType::BMatrixType>(
227  dNdx_u, N_u, x_coord, _is_axially_symmetric);
228 
229  auto& eps = _ip_data[ip].eps;
230  eps.noalias() = B * u;
231  auto const& sigma_eff = _ip_data[ip].sigma_eff;
232 
233  double const p_int_pt = N_p.dot(p);
234  vars.phase_pressure = p_int_pt;
235 
236  auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
237  t, x_position, dt, T_ref);
238  auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
239 
240  auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
241  .template value<double>(vars, x_position, t, dt);
242 
243  auto const rho_sr =
244  solid.property(MPL::PropertyType::density)
245  .template value<double>(vars, x_position, t, dt);
246  auto const porosity =
247  medium->property(MPL::PropertyType::porosity)
248  .template value<double>(vars, x_position, t, dt);
249 
250  auto const mu = fluid.property(MPL::PropertyType::viscosity)
251  .template value<double>(vars, x_position, t, dt);
252 
253  // Quick workaround: If fluid density is described as ideal gas, then
254  // the molar mass must be passed to the MPL::IdealGasLaw via the
255  // variable_array and the fluid must have the property
256  // MPL::PropertyType::molar_mass. For other density models (e.g.
257  // Constant), it is not mandatory to specify the molar mass.
258  if (fluid.hasProperty(MPL::PropertyType::molar_mass))
259  {
260  vars.molar_mass =
261  fluid.property(MPL::PropertyType::molar_mass)
262  .template value<double>(vars, x_position, t, dt);
263  }
264  auto const rho_fr =
265  fluid.property(MPL::PropertyType::density)
266  .template value<double>(vars, x_position, t, dt);
267  auto const beta_p =
268  fluid.property(MPL::PropertyType::density)
269  .template dValue<double>(vars, MPL::Variable::phase_pressure,
270  x_position, t, dt) /
271  rho_fr;
272 
273  // Set mechanical variables for the intrinsic permeability model
274  // For stress dependent permeability.
275  {
276  auto const sigma_total =
277  (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
278 
279  vars.total_stress.emplace<SymmetricTensor>(
281  sigma_total));
282  }
283  // For strain dependent permeability
286  _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
287  vars.mechanical_strain
289  eps);
290 
291  auto const K = MPL::formEigenTensor<DisplacementDim>(
292  medium->property(MPL::PropertyType::permeability)
293  .value(vars, x_position, t, dt));
294  auto const dkde = MPL::formEigenTensor<
296  (*medium)[MPL::PropertyType::permeability].dValue(
297  vars, MPL::Variable::mechanical_strain, x_position, t, dt));
298 
299  auto const K_over_mu = K / mu;
300 
301  auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
302  dt, u, T_ref);
303 
304  //
305  // displacement equation, displacement part
306  //
307  local_Jac
308  .template block<displacement_size, displacement_size>(
310  .noalias() += B.transpose() * C * B * w;
311 
312  double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
313  local_rhs.template segment<displacement_size>(displacement_index)
314  .noalias() -=
315  (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
316 
317  //
318  // displacement equation, pressure part
319  //
320  Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
321 
322  //
323  // pressure equation, pressure part.
324  //
325  laplace_p.noalias() +=
326  rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
327 
328  storage_p.noalias() +=
329  rho_fr * N_p.transpose() * N_p * w *
330  ((alpha - porosity) * (1.0 - alpha) / K_S + porosity * beta_p);
331 
332  // density dependence on pressure evaluated for Darcy-term,
333  // for laplace and storage terms this dependence is neglected
334  add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
335  K_over_mu *
336  (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
337 
338  local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
339  dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
340 
341  //
342  // pressure equation, displacement part.
343  //
344  Kpu.noalias() +=
345  rho_fr * alpha * N_p.transpose() * identity2.transpose() * B * w;
346 
347  Kpu_k.noalias() +=
348  dNdx_p.transpose() *
349  MathLib::KelvinVector::liftVectorToKelvin<DisplacementDim>(
350  dNdx_p * p - rho_fr * b) *
351  dkde * B * rho_fr / mu * w;
352  }
353  // displacement equation, pressure part
354  local_Jac
355  .template block<displacement_size, pressure_size>(displacement_index,
357  .noalias() = -Kup;
358 
359  if (_process_data.mass_lumping)
360  {
361  storage_p = storage_p.colwise().sum().eval().asDiagonal();
362 
363  if constexpr (pressure_size == displacement_size)
364  {
365  Kpu = Kpu.colwise().sum().eval().asDiagonal();
366  Kpu_k = Kpu_k.colwise().sum().eval().asDiagonal();
367  }
368  }
369 
370  // pressure equation, pressure part.
371  local_Jac
372  .template block<pressure_size, pressure_size>(pressure_index,
374  .noalias() += laplace_p + storage_p / dt + add_p_derivative;
375 
376  // pressure equation, displacement part.
377  local_Jac
378  .template block<pressure_size, displacement_size>(pressure_index,
380  .noalias() += Kpu / dt + Kpu_k;
381 
382  // pressure equation
383  local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
384  laplace_p * p + storage_p * p_dot + Kpu * u_dot;
385 
386  // displacement equation
387  local_rhs.template segment<displacement_size>(displacement_index)
388  .noalias() += Kup * p;
389 }
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
Definition: VariableType.h:165
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > total_stress
Definition: VariableType.h:184
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:89
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:55
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
@ rho
density. For some models, rho substitutes p
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Phase const & fluidPhase(Medium const &medium)
Returns a gas or aqueous liquid phase of the given medium.
Definition: Medium.cpp:84
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:32
static const double u
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
static const double t
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition: LinearBMatrix.h:42
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Definition: KelvinVector.h:77
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.

References MaterialPropertyLib::alpha, MaterialPropertyLib::biot_coefficient, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::fluidPhase(), MaterialPropertyLib::formEigenTensor(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::mechanical_strain, MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::molar_mass, MaterialPropertyLib::VariableArray::molar_mass, MathLib::p, MaterialPropertyLib::permeability, MaterialPropertyLib::phase_pressure, MaterialPropertyLib::VariableArray::phase_pressure, MaterialPropertyLib::porosity, MaterialPropertyLib::reference_temperature, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), MathLib::t, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MathLib::u, MaterialPropertyLib::viscosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ assembleWithJacobianForDeformationEquations()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForDeformationEquations ( const double  t,
double const  dt,
Eigen::VectorXd const &  local_x,
std::vector< double > &  local_b_data,
std::vector< double > &  local_Jac_data 
)
private

Assemble local matrices and vectors arise from the linearized discretized weak form of the residual of the momentum balance equation,

\[ \nabla (\sigma - \alpha_b p \mathrm{I}) = f \]

where \( \sigma\) is the effective stress tensor, \(p\) is the pore pressure, \(\alpha_b\) is the Biot constant, \(\mathrm{I}\) is the identity tensor, and \(f\) is the body force.

Parameters
tTime
dtTime increment
local_xNodal values of \(x\) of an element of all coupled processes.
local_b_dataRight hand side vector of an element.
local_Jac_dataElement Jacobian matrix for the Newton-Raphson method.

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

673 {
674  auto const p = local_x.template segment<pressure_size>(pressure_index);
675  auto const u =
676  local_x.template segment<displacement_size>(displacement_index);
677 
678  auto local_Jac = MathLib::createZeroedMatrix<
679  typename ShapeMatricesTypeDisplacement::template MatrixType<
681  local_Jac_data, displacement_size, displacement_size);
682 
683  auto local_rhs =
684  MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
685  template VectorType<displacement_size>>(
686  local_b_data, displacement_size);
687 
689  x_position.setElementID(_element.getID());
690 
691  auto const& medium = _process_data.media_map->getMedium(_element.getID());
692  auto const& solid = medium->phase("Solid");
693  auto const& fluid = fluidPhase(*medium);
694  MPL::VariableArray vars;
695 
696  auto const T_ref =
698  .template value<double>(vars, x_position, t, dt);
699  vars.temperature = T_ref;
700 
701  int const n_integration_points = _integration_method.getNumberOfPoints();
702  for (int ip = 0; ip < n_integration_points; ip++)
703  {
704  x_position.setIntegrationPoint(ip);
705  auto const& w = _ip_data[ip].integration_weight;
706 
707  auto const& N_u_op = _ip_data[ip].N_u_op;
708 
709  auto const& N_u = _ip_data[ip].N_u;
710  auto const& dNdx_u = _ip_data[ip].dNdx_u;
711 
712  auto const& N_p = _ip_data[ip].N_p;
713 
714  auto const x_coord =
715  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
717  _element, N_u);
718  auto const B =
719  LinearBMatrix::computeBMatrix<DisplacementDim,
720  ShapeFunctionDisplacement::NPOINTS,
721  typename BMatricesType::BMatrixType>(
722  dNdx_u, N_u, x_coord, _is_axially_symmetric);
723 
724  auto& eps = _ip_data[ip].eps;
725  auto const& sigma_eff = _ip_data[ip].sigma_eff;
726 
727  vars.phase_pressure = N_p.dot(p);
728 
729  auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
730  .template value<double>(vars, x_position, t, dt);
731  auto const rho_sr =
732  solid.property(MPL::PropertyType::density)
733  .template value<double>(vars, x_position, t, dt);
734  auto const porosity =
735  medium->property(MPL::PropertyType::porosity)
736  .template value<double>(vars, x_position, t, dt);
737 
738  // Quick workaround: If fluid density is described as ideal gas, then
739  // the molar mass must be passed to the MPL::IdealGasLaw via the
740  // variable_array and the fluid must have the property
741  // MPL::PropertyType::molar_mass. For other density models (e.g.
742  // Constant), it is not mandatory to specify the molar mass.
743  if (fluid.hasProperty(MPL::PropertyType::molar_mass))
744  {
745  vars.molar_mass =
746  fluid.property(MPL::PropertyType::molar_mass)
747  .template value<double>(vars, x_position, t, dt);
748  }
749  auto const rho_fr =
750  fluid.property(MPL::PropertyType::density)
751  .template value<double>(vars, x_position, t, dt);
752 
753  auto const& b = _process_data.specific_body_force;
754  auto const& identity2 = MathLib::KelvinVector::Invariants<
756  DisplacementDim)>::identity2;
757 
758  eps.noalias() = B * u;
759  vars.mechanical_strain
761  eps);
762 
763  auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
764  dt, u, T_ref);
765 
766  local_Jac.noalias() += B.transpose() * C * B * w;
767 
768  double p_at_xi = 0.;
769  NumLib::shapeFunctionInterpolate(p, N_p, p_at_xi);
770 
771  double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
772  local_rhs.noalias() -=
773  (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
774  N_u_op.transpose() * rho * b) *
775  w;
776  }
777 }
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
Definition: Interpolation.h:80

References MaterialPropertyLib::alpha, MaterialPropertyLib::biot_coefficient, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, MaterialPropertyLib::fluidPhase(), NumLib::interpolateXCoordinate(), MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::molar_mass, MaterialPropertyLib::VariableArray::molar_mass, MathLib::p, MaterialPropertyLib::VariableArray::phase_pressure, MaterialPropertyLib::porosity, MaterialPropertyLib::reference_temperature, ParameterLib::SpatialPosition::setElementID(), NumLib::shapeFunctionInterpolate(), MathLib::t, MaterialPropertyLib::VariableArray::temperature, and MathLib::u.

◆ assembleWithJacobianForPressureEquations()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForPressureEquations ( const double  t,
double const  dt,
Eigen::VectorXd const &  local_x,
Eigen::VectorXd const &  local_xdot,
std::vector< double > &  local_b_data,
std::vector< double > &  local_Jac_data 
)
private

Assemble local matrices and vectors arise from the linearized discretized weak form of the residual of the mass balance equation of single phase flow,

\[ \alpha \cdot{p} - \nabla (K (\nabla p + \rho g \nabla z) + \alpha_b \nabla \cdot \dot{u} = Q \]

where \( alpha\) is a coefficient may depend on storage or the fluid density change, \( \rho\) is the fluid density, \(g\) is the gravitational acceleration, \(z\) is the vertical coordinate, \(u\) is the displacement, and \(Q\) is the source/sink term.

Parameters
tTime
dtTime increment
local_xNodal values of \(x\) of an element of all coupled processes.
local_xdotNodal values of \(\dot{x}\) of an element of all coupled processes.
local_b_dataRight hand side vector of an element.
local_Jac_dataElement Jacobian matrix for the Newton-Raphson method.

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

500 {
501  auto local_rhs =
502  MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
503  template VectorType<pressure_size>>(
504  local_b_data, pressure_size);
505 
507  pos.setElementID(this->_element.getID());
508 
509  auto const p = local_x.template segment<pressure_size>(pressure_index);
510 
511  auto const p_dot =
512  local_xdot.template segment<pressure_size>(pressure_index);
513 
514  auto local_Jac = MathLib::createZeroedMatrix<
515  typename ShapeMatricesTypeDisplacement::template MatrixType<
516  pressure_size, pressure_size>>(local_Jac_data, pressure_size,
517  pressure_size);
518 
520  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
521  pressure_size);
522 
524  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
525  pressure_size);
526 
528  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
529  pressure_size);
530 
531  typename ShapeMatricesTypePressure::NodalMatrixType add_p_derivative =
532  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
533  pressure_size);
534 
535  auto const& solid_material =
537  _process_data.solid_materials, _process_data.material_ids,
538  _element.getID());
539 
541  x_position.setElementID(_element.getID());
542 
543  auto const& medium = _process_data.media_map->getMedium(_element.getID());
544  auto const& fluid = fluidPhase(*medium);
545  MPL::VariableArray vars;
546 
547  auto const T_ref =
549  .template value<double>(vars, x_position, t, dt);
550  vars.temperature = T_ref;
551 
552  auto const& identity2 = Invariants::identity2;
553 
554  auto const fixed_stress_stabilization_parameter =
555  _process_data.fixed_stress_stabilization_parameter;
556 
557  int const n_integration_points = _integration_method.getNumberOfPoints();
558  for (int ip = 0; ip < n_integration_points; ip++)
559  {
560  x_position.setIntegrationPoint(ip);
561  auto const& w = _ip_data[ip].integration_weight;
562 
563  auto const& N_p = _ip_data[ip].N_p;
564  auto const& dNdx_p = _ip_data[ip].dNdx_p;
565 
566  double const p_int_pt = N_p.dot(p);
567  vars.phase_pressure = p_int_pt;
568 
569  auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
570  t, x_position, dt, T_ref);
571  auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
572 
573  auto const alpha_b =
574  medium->property(MPL::PropertyType::biot_coefficient)
575  .template value<double>(vars, x_position, t, dt);
576 
577  // Set mechanical variables for the intrinsic permeability model
578  // For stress dependent permeability.
579  auto const sigma_total =
580  (_ip_data[ip].sigma_eff - alpha_b * p_int_pt * identity2).eval();
581  vars.total_stress.emplace<SymmetricTensor>(
583 
584  // For strain dependent permeability
587  _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
588  vars.mechanical_strain
590  _ip_data[ip].eps);
591 
592  auto const K = MPL::formEigenTensor<DisplacementDim>(
593  medium->property(MPL::PropertyType::permeability)
594  .value(vars, x_position, t, dt));
595  auto const porosity =
596  medium->property(MPL::PropertyType::porosity)
597  .template value<double>(vars, x_position, t, dt);
598 
599  auto const mu = fluid.property(MPL::PropertyType::viscosity)
600  .template value<double>(vars, x_position, t, dt);
601  // Quick workaround: If fluid density is described as ideal gas, then
602  // the molar mass must be passed to the MPL::IdealGasLaw via the
603  // variable_array and the fluid must have the property
604  // MPL::PropertyType::molar_mass. For other density models (e.g.
605  // Constant), it is not mandatory to specify the molar mass.
606  if (fluid.hasProperty(MPL::PropertyType::molar_mass))
607  {
608  vars.molar_mass =
609  fluid.property(MPL::PropertyType::molar_mass)
610  .template value<double>(vars, x_position, t, dt);
611  }
612  auto const rho_fr =
613  fluid.property(MPL::PropertyType::density)
614  .template value<double>(vars, x_position, t, dt);
615  auto const beta_p =
616  fluid.property(MPL::PropertyType::density)
617  .template dValue<double>(vars, MPL::Variable::phase_pressure,
618  x_position, t, dt) /
619  rho_fr;
620 
621  auto const K_over_mu = K / mu;
622 
623  // artificial compressibility to account for coupling
624  auto const beta_FS =
625  fixed_stress_stabilization_parameter * alpha_b * alpha_b / K_S;
626 
627  laplace.noalias() +=
628  rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
629 
630  storage.noalias() +=
631  rho_fr * N_p.transpose() * N_p * w *
632  ((alpha_b - porosity) * (1.0 - alpha_b) / K_S + porosity * beta_p);
633 
634  auto const& b = _process_data.specific_body_force;
635 
636  // bodyforce-driven Darcy flow
637  local_rhs.noalias() +=
638  dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
639 
640  // density dependence on pressure evaluated for Darcy-term,
641  // for laplace and storage terms this dependence is neglected (as is
642  // done for monolithic too)
643  add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
644  K_over_mu *
645  (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
646 
647  auto& eps = _ip_data[ip].eps;
648  auto& eps_prev = _ip_data[ip].eps_prev;
649  const double eps_v_dot =
650  (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
651 
652  const double p_c = _ip_data[ip].coupling_pressure;
653 
654  // pressure-dependent part of fixed-stress coupling
655  coupling.noalias() += rho_fr * N_p.transpose() * N_p * w * beta_FS / dt;
656 
657  // constant part of fixed-stress coupling
658  local_rhs.noalias() -=
659  (alpha_b * eps_v_dot - beta_FS * p_c / dt) * rho_fr * N_p * w;
660  }
661  local_Jac.noalias() = laplace + storage / dt + add_p_derivative + coupling;
662 
663  local_rhs.noalias() -= laplace * p + storage * p_dot + coupling * p;
664 }

References MaterialPropertyLib::biot_coefficient, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::fluidPhase(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::molar_mass, MaterialPropertyLib::VariableArray::molar_mass, MathLib::p, MaterialPropertyLib::permeability, MaterialPropertyLib::phase_pressure, MaterialPropertyLib::VariableArray::phase_pressure, MaterialPropertyLib::porosity, MaterialPropertyLib::reference_temperature, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), MaterialPropertyLib::storage, MathLib::t, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::viscosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ assembleWithJacobianForStaggeredScheme()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::assembleWithJacobianForStaggeredScheme ( const double  t,
double const  dt,
Eigen::VectorXd const &  local_x,
Eigen::VectorXd const &  local_xdot,
int const  process_id,
std::vector< double > &  local_M_data,
std::vector< double > &  local_K_data,
std::vector< double > &  local_b_data,
std::vector< double > &  local_Jac_data 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

789 {
790  // For the equations with pressure
791  if (process_id == _process_data.hydraulic_process_id)
792  {
793  assembleWithJacobianForPressureEquations(t, dt, local_x, local_xdot,
794  local_b_data, local_Jac_data);
795  return;
796  }
797 
798  // For the equations with deformation
799  assembleWithJacobianForDeformationEquations(t, dt, local_x, local_b_data,
800  local_Jac_data);
801 }
void assembleWithJacobianForPressureEquations(const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleWithJacobianForDeformationEquations(const double t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)

References MathLib::t.

◆ computeSecondaryVariableConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::computeSecondaryVariableConcrete ( double const  t,
double const  dt,
Eigen::VectorXd const &  local_xs,
Eigen::VectorXd const &  local_x_dot 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

1023 {
1024  auto const p = local_x.template segment<pressure_size>(pressure_index);
1025 
1027  ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1028  DisplacementDim>(_element, _is_axially_symmetric, p,
1029  *_process_data.pressure_interpolated);
1030 
1031  int const elem_id = _element.getID();
1032  ParameterLib::SpatialPosition x_position;
1033  x_position.setElementID(elem_id);
1034  unsigned const n_integration_points =
1036 
1037  auto const& medium = _process_data.media_map->getMedium(elem_id);
1038  MPL::VariableArray vars;
1039 
1040  SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize);
1041  auto sigma_eff_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
1042  Eigen::Matrix<double, 3, 3>::Zero());
1043 
1044  auto const& identity2 = Invariants::identity2;
1045 
1046  for (unsigned ip = 0; ip < n_integration_points; ip++)
1047  {
1048  x_position.setIntegrationPoint(ip);
1049 
1050  auto const& eps = _ip_data[ip].eps;
1051  sigma_eff_sum += _ip_data[ip].sigma_eff;
1052 
1053  auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
1054  .template value<double>(vars, x_position, t, dt);
1055  double const p_int_pt = _ip_data[ip].N_p.dot(p);
1056  vars.phase_pressure = p_int_pt;
1057 
1058  // Set mechanical variables for the intrinsic permeability model
1059  // For stress dependent permeability.
1060  {
1061  auto const sigma_total =
1062  (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
1063  vars.total_stress.emplace<SymmetricTensor>(
1065  sigma_total));
1066  }
1067  // For strain dependent permeability
1070  _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1071  vars.mechanical_strain
1073  eps);
1074 
1075  k_sum += MPL::getSymmetricTensor<DisplacementDim>(
1076  medium->property(MPL::PropertyType::permeability)
1077  .value(vars, x_position, t, dt));
1078  }
1079 
1080  Eigen::Map<Eigen::VectorXd>(
1081  &(*_process_data.permeability)[elem_id * KelvinVectorSize],
1082  KelvinVectorSize) = k_sum / n_integration_points;
1083 
1084  Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
1086  n_integration_points;
1087 
1088  Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma_avg);
1089 
1090  Eigen::Map<Eigen::Vector3d>(
1091  &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
1092  e_s.eigenvalues();
1093 
1094  auto eigen_vectors = e_s.eigenvectors();
1095 
1096  for (auto i = 0; i < 3; i++)
1097  {
1098  Eigen::Map<Eigen::Vector3d>(
1099  &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
1100  eigen_vectors.col(i);
1101  }
1102 }
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)

References MaterialPropertyLib::alpha, MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, NumLib::interpolateToHigherOrderNodes(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MathLib::KelvinVector::kelvinVectorToTensor(), MaterialPropertyLib::VariableArray::mechanical_strain, MathLib::p, MaterialPropertyLib::permeability, MaterialPropertyLib::VariableArray::phase_pressure, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), MathLib::t, MaterialPropertyLib::VariableArray::total_stress, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ getEpsilon()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getEpsilon
overridevirtual

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

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

995 {
996  auto const kelvin_vector_size =
998  unsigned const n_integration_points =
1000 
1001  std::vector<double> ip_epsilon_values;
1002  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
1003  double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
1004  ip_epsilon_values, n_integration_points, kelvin_vector_size);
1005 
1006  for (unsigned ip = 0; ip < n_integration_points; ++ip)
1007  {
1008  auto const& eps = _ip_data[ip].eps;
1009  cache_mat.row(ip) =
1011  }
1012 
1013  return ip_epsilon_values;
1014 }

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

◆ getIntPtDarcyVelocity()

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

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

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

400 {
401  int const hydraulic_process_id = _process_data.hydraulic_process_id;
402  auto const indices =
403  NumLib::getIndices(_element.getID(), *dof_table[hydraulic_process_id]);
404  assert(!indices.empty());
405  auto const local_x = x[hydraulic_process_id]->get(indices);
406 
407  unsigned const n_integration_points =
409  cache.clear();
410  auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
411  double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
412  cache, DisplacementDim, n_integration_points);
413 
414  auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
415  pressure_size> const>(local_x.data() + pressure_index, pressure_size);
416 
418  x_position.setElementID(_element.getID());
419 
420  auto const& medium = _process_data.media_map->getMedium(_element.getID());
421  auto const& fluid = fluidPhase(*medium);
422  MPL::VariableArray vars;
423 
424  // TODO (naumov) Temporary value not used by current material models. Need
425  // extension of secondary variables interface.
426  double const dt = std::numeric_limits<double>::quiet_NaN();
427  vars.temperature =
429  .template value<double>(vars, x_position, t, dt);
430 
431  auto const& identity2 = Invariants::identity2;
432 
433  for (unsigned ip = 0; ip < n_integration_points; ip++)
434  {
435  x_position.setIntegrationPoint(ip);
436 
437  double const p_int_pt = _ip_data[ip].N_p.dot(p);
438  vars.phase_pressure = p_int_pt;
439 
440  auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
441  .template value<double>(vars, x_position, t, dt);
442 
443  // Set mechanical variables for the intrinsic permeability model
444  // For stress dependent permeability.
445  auto const sigma_total =
446  (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
447  vars.total_stress.emplace<SymmetricTensor>(
449 
450  // For strain dependent permeability
453  _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
454  vars.mechanical_strain
456  _ip_data[ip].eps);
457 
458  auto const K = MPL::formEigenTensor<DisplacementDim>(
459  medium->property(MPL::PropertyType::permeability)
460  .value(vars, x_position, t, dt));
461 
462  auto const mu = fluid.property(MPL::PropertyType::viscosity)
463  .template value<double>(vars, x_position, t, dt);
464  // Quick workaround: If fluid density is described as ideal gas, then
465  // the molar mass must be passed to the MPL::IdealGasLaw via the
466  // variable_array and the fluid must have the property
467  // MPL::PropertyType::molar_mass. For other density models (e.g.
468  // Constant), it is not mandatory to specify the molar mass.
469  if (fluid.hasProperty(MPL::PropertyType::molar_mass))
470  {
471  vars.molar_mass =
472  fluid.property(MPL::PropertyType::molar_mass)
473  .template value<double>(vars, x_position, t, dt);
474  }
475  auto const rho_fr =
476  fluid.property(MPL::PropertyType::density)
477  .template value<double>(vars, x_position, t, dt);
478 
479  auto const K_over_mu = K / mu;
480 
481  auto const& b = _process_data.specific_body_force;
482 
483  // Compute the velocity
484  auto const& dNdx_p = _ip_data[ip].dNdx_p;
485  cache_matrix.col(ip).noalias() =
486  -K_over_mu * dNdx_p * p + K_over_mu * rho_fr * b;
487  }
488 
489  return cache;
490 }
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)

References MaterialPropertyLib::alpha, MaterialPropertyLib::biot_coefficient, MathLib::createZeroedMatrix(), MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::equivalent_plastic_strain, MaterialPropertyLib::fluidPhase(), NumLib::getIndices(), MathLib::KelvinVector::kelvinVectorToSymmetricTensor(), MaterialPropertyLib::VariableArray::mechanical_strain, MaterialPropertyLib::molar_mass, MaterialPropertyLib::VariableArray::molar_mass, MathLib::p, MaterialPropertyLib::permeability, MaterialPropertyLib::VariableArray::phase_pressure, MaterialPropertyLib::reference_temperature, ParameterLib::SpatialPosition::setElementID(), MathLib::t, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::VariableArray::total_stress, MaterialPropertyLib::viscosity, and MaterialPropertyLib::VariableArray::volumetric_strain.

◆ getIntPtEpsilon()

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

◆ getIntPtSigma()

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

◆ getMaterialStateVariablesAt()

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

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

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

1120 {
1121  return *_ip_data[integration_point].material_state_variables;
1122 }

◆ getNumberOfIntegrationPoints()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
unsigned ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getNumberOfIntegrationPoints
overrideprivatevirtual

◆ getShapeMatrix()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
Eigen::Map<const Eigen::RowVectorXd> ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getShapeMatrix ( const unsigned  integration_point) const
inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 285 of file HydroMechanicsFEM.h.

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

References ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_secondary_data, and ProcessLib::HydroMechanics::SecondaryData< ShapeMatrixType >::N_u.

◆ getSigma()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::vector< double > ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::getSigma
overridevirtual

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

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

985 {
986  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
988 }

◆ initializeConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::initializeConcrete ( )
inlineoverridevirtual

Set initial stress from parameter.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 226 of file HydroMechanicsFEM.h.

227  {
228  unsigned const n_integration_points =
230 
231  for (unsigned ip = 0; ip < n_integration_points; ip++)
232  {
233  auto& ip_data = _ip_data[ip];
234 
236  if (_process_data.initial_stress != nullptr)
237  {
238  ParameterLib::SpatialPosition const x_position{
239  std::nullopt, _element.getID(), ip,
241  ShapeFunctionDisplacement,
243  _element, ip_data.N_u))};
244 
245  ip_data.sigma_eff =
247  DisplacementDim>((*_process_data.initial_stress)(
248  std::numeric_limits<
249  double>::quiet_NaN() /* time independent */,
250  x_position));
251  }
252 
253  ip_data.pushBackState();
254  }
255  }
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
Definition: KelvinVector.h:171
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)

References ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_element, ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_integration_method, ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_ip_data, ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_process_data, MeshLib::Element::getID(), NumLib::GenericIntegrationMethod::getNumberOfPoints(), NumLib::interpolateCoordinates(), and MathLib::KelvinVector::symmetricTensorToKelvinVector().

◆ postNonLinearSolverConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, 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 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

873 {
874  int const n_integration_points = _integration_method.getNumberOfPoints();
875 
876  if (use_monolithic_scheme || process_id == 0)
877  {
878  auto p =
879  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
880  pressure_size> const>(local_x.data(), pressure_size);
881 
882  for (int ip = 0; ip < n_integration_points; ip++)
883  {
884  auto const& N_p = _ip_data[ip].N_p;
885  _ip_data[ip].coupling_pressure = N_p.dot(p);
886  }
887  }
888 
889  if (use_monolithic_scheme || process_id == 1)
890  {
892  x_position.setElementID(_element.getID());
893 
894  auto const& medium =
895  _process_data.media_map->getMedium(_element.getID());
896 
897  auto const T_ref =
899  .template value<double>(MPL::VariableArray(), x_position, t,
900  dt);
901 
902  const int displacement_offset =
903  use_monolithic_scheme ? displacement_index : 0;
904 
905  auto u = Eigen::Map<typename ShapeMatricesTypeDisplacement::
906  template VectorType<displacement_size> const>(
907  local_x.data() + displacement_offset, displacement_size);
908 
909  MPL::VariableArray vars;
910  vars.temperature = T_ref;
911 
912  for (int ip = 0; ip < n_integration_points; ip++)
913  {
914  x_position.setIntegrationPoint(ip);
915  auto const& N_u = _ip_data[ip].N_u;
916  auto const& dNdx_u = _ip_data[ip].dNdx_u;
917 
918  auto const x_coord =
919  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
921  _element, N_u);
922  auto const B = LinearBMatrix::computeBMatrix<
923  DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
924  typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
926 
927  auto& eps = _ip_data[ip].eps;
928  eps.noalias() = B * u;
929  vars.mechanical_strain.emplace<
931 
932  _ip_data[ip].updateConstitutiveRelation(vars, t, x_position, dt, u,
933  T_ref);
934  }
935  }
936 }

References ProcessLib::LinearBMatrix::computeBMatrix(), NumLib::interpolateXCoordinate(), MathLib::p, MaterialPropertyLib::reference_temperature, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), MathLib::t, MaterialPropertyLib::VariableArray::temperature, and MathLib::u.

◆ postTimestepConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::postTimestepConcrete ( Eigen::VectorXd const &  ,
double const  ,
double const   
)
inlineoverridevirtual

◆ setInitialConditionsConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
void ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setInitialConditionsConcrete ( std::vector< double > const &  local_x,
double const  t,
bool const  use_monolithic_scheme,
int const  process_id 
)
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

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

811 {
812  int const n_integration_points = _integration_method.getNumberOfPoints();
813 
814  if (use_monolithic_scheme || process_id == 0)
815  {
816  auto p =
817  Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
818  pressure_size> const>(local_x.data(), pressure_size);
819 
820  for (int ip = 0; ip < n_integration_points; ip++)
821  {
822  auto const& N_p = _ip_data[ip].N_p;
823  _ip_data[ip].coupling_pressure = N_p.dot(p);
824  }
825  }
826 
827  if (use_monolithic_scheme || process_id == 1)
828  {
830  x_position.setElementID(_element.getID());
831 
832  const int displacement_offset =
833  use_monolithic_scheme ? displacement_index : 0;
834 
835  auto u = Eigen::Map<typename ShapeMatricesTypeDisplacement::
836  template VectorType<displacement_size> const>(
837  local_x.data() + displacement_offset, displacement_size);
838 
839  MPL::VariableArray vars;
840 
841  for (int ip = 0; ip < n_integration_points; ip++)
842  {
843  x_position.setIntegrationPoint(ip);
844  auto const& N_u = _ip_data[ip].N_u;
845  auto const& dNdx_u = _ip_data[ip].dNdx_u;
846 
847  auto const x_coord =
848  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
850  _element, N_u);
851  auto const B = LinearBMatrix::computeBMatrix<
852  DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
853  typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
855 
856  auto& eps = _ip_data[ip].eps;
857  eps.noalias() = B * u;
858  vars.mechanical_strain.emplace<
860  }
861  }
862 }

References ProcessLib::LinearBMatrix::computeBMatrix(), NumLib::interpolateXCoordinate(), MathLib::p, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), and MathLib::u.

◆ setIPDataInitialConditions()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
std::size_t ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::setIPDataInitialConditions ( std::string const &  name,
double const *  values,
int const  integration_order 
)
overridevirtual

Returns number of read integration points.

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

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

945 {
946  if (integration_order !=
947  static_cast<int>(_integration_method.getIntegrationOrder()))
948  {
949  OGS_FATAL(
950  "Setting integration point initial conditions; The integration "
951  "order of the local assembler for element {:d} is different from "
952  "the integration order in the initial condition.",
953  _element.getID());
954  }
955 
956  if (name == "sigma_ip")
957  {
958  if (_process_data.initial_stress != nullptr)
959  {
960  OGS_FATAL(
961  "Setting initial conditions for stress from integration "
962  "point data and from a parameter '{:s}' is not possible "
963  "simultaneously.",
964  _process_data.initial_stress->name);
965  }
966 
967  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
968  values, _ip_data, &IpData::sigma_eff);
969  }
970 
971  if (name == "epsilon_ip")
972  {
973  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
974  values, _ip_data, &IpData::eps);
975  }
976 
977  return 0;
978 }

References MaterialPropertyLib::name, and OGS_FATAL.

Member Data Documentation

◆ _element

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
MeshLib::Element const& ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_element
private

◆ _integration_method

◆ _ip_data

◆ _is_axially_symmetric

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
bool const ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::_is_axially_symmetric
private

Definition at line 397 of file HydroMechanicsFEM.h.

◆ _process_data

◆ _secondary_data

◆ displacement_index

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

Definition at line 404 of file HydroMechanicsFEM.h.

◆ displacement_size

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

◆ KelvinVectorSize

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , int DisplacementDim>
int const ProcessLib::HydroMechanics::HydroMechanicsLocalAssembler< ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim >::KelvinVectorSize
static
Initial value:

Definition at line 177 of file HydroMechanicsFEM.h.

◆ pressure_index

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

Definition at line 402 of file HydroMechanicsFEM.h.

◆ pressure_size

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

Definition at line 403 of file HydroMechanicsFEM.h.


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