OGS
ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim > Class Template Reference

Detailed Description

template<typename ShapeFunctionDisplacement, typename ShapeFunctionPressure, typename IntegrationMethod, int GlobalDim>
class ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >

Definition at line 35 of file HydroMechanicsLocalAssemblerMatrix.h.

#include <HydroMechanicsLocalAssemblerMatrix.h>

Inheritance diagram for ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >:
[legend]
Collaboration diagram for ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >:
[legend]

Public Member Functions

 HydroMechanicsLocalAssemblerMatrix (HydroMechanicsLocalAssemblerMatrix const &)=delete
 
 HydroMechanicsLocalAssemblerMatrix (HydroMechanicsLocalAssemblerMatrix &&)=delete
 
 HydroMechanicsLocalAssemblerMatrix (MeshLib::Element const &e, std::size_t const n_variables, std::size_t const local_matrix_size, std::vector< unsigned > const &dofIndex_to_localIndex, bool const is_axially_symmetric, unsigned const integration_order, HydroMechanicsProcessData< GlobalDim > &process_data)
 
void preTimestepConcrete (std::vector< double > const &, double const, double const) override
 
- Public Member Functions inherited from ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerInterface
 HydroMechanicsLocalAssemblerInterface (MeshLib::Element const &element, bool const is_axially_symmetric, std::size_t n_local_size, std::vector< unsigned > dofIndex_to_localIndex)
 
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_, const double, const double, std::vector< double > &, std::vector< double > &, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
 
void postTimestepConcrete (Eigen::VectorXd const &local_x_, const double t, double const dt) 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 assembleWithJacobianForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, const double dxdot_dx, const double dx_dx, 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)
 
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 std::vector< double > interpolateNodalValuesToIntegrationPoints (std::vector< double > const &)
 
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...
 

Protected Types

using ShapeMatricesTypeDisplacement = ShapeMatrixPolicyType< ShapeFunctionDisplacement, GlobalDim >
 
using BMatricesType = BMatrixPolicyType< ShapeFunctionDisplacement, GlobalDim >
 
using ShapeMatricesTypePressure = ShapeMatrixPolicyType< ShapeFunctionPressure, GlobalDim >
 
using IntegrationPointDataType = IntegrationPointDataMatrix< BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, GlobalDim, ShapeFunctionDisplacement::NPOINTS >
 
using GlobalDimVector = Eigen::Matrix< double, GlobalDim, 1 >
 

Protected Member Functions

void assembleWithJacobianConcrete (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_dot, Eigen::VectorXd &local_rhs, Eigen::MatrixXd &local_Jac) override
 
void assembleBlockMatricesWithJacobian (double const t, double const dt, Eigen::Ref< const Eigen::VectorXd > const &p, Eigen::Ref< const Eigen::VectorXd > const &p_dot, Eigen::Ref< const Eigen::VectorXd > const &u, Eigen::Ref< const Eigen::VectorXd > const &u_dot, Eigen::Ref< Eigen::VectorXd > rhs_p, Eigen::Ref< Eigen::VectorXd > rhs_u, Eigen::Ref< Eigen::MatrixXd > J_pp, Eigen::Ref< Eigen::MatrixXd > J_pu, Eigen::Ref< Eigen::MatrixXd > J_uu, Eigen::Ref< Eigen::MatrixXd > J_up)
 
void postTimestepConcreteWithVector (double const t, double const dt, Eigen::VectorXd const &local_x) override
 
void postTimestepConcreteWithBlockVectors (double const t, double const dt, Eigen::Ref< const Eigen::VectorXd > const &p, Eigen::Ref< const Eigen::VectorXd > const &u)
 
void setPressureOfInactiveNodes (double const t, Eigen::Ref< Eigen::VectorXd > p)
 
void setPressureDotOfInactiveNodes (Eigen::Ref< Eigen::VectorXd > p_dot)
 

Protected Attributes

HydroMechanicsProcessData< GlobalDim > & _process_data
 
std::vector< IntegrationPointDataType, Eigen::aligned_allocator< IntegrationPointDataType > > _ip_data
 
- Protected Attributes inherited from ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerInterface
MeshLib::Element const & _element
 
bool const _is_axially_symmetric
 

Static Protected 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
 
static const int kelvin_vector_size
 

Member Typedef Documentation

◆ BMatricesType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
using ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::BMatricesType = BMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>
protected

Definition at line 96 of file HydroMechanicsLocalAssemblerMatrix.h.

◆ GlobalDimVector

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
using ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::GlobalDimVector = Eigen::Matrix<double, GlobalDim, 1>
protected

Definition at line 110 of file HydroMechanicsLocalAssemblerMatrix.h.

◆ IntegrationPointDataType

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
using ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::IntegrationPointDataType = IntegrationPointDataMatrix<BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure, GlobalDim, ShapeFunctionDisplacement::NPOINTS>
protected

Definition at line 103 of file HydroMechanicsLocalAssemblerMatrix.h.

◆ ShapeMatricesTypeDisplacement

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
using ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::ShapeMatricesTypeDisplacement = ShapeMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>
protected

Definition at line 94 of file HydroMechanicsLocalAssemblerMatrix.h.

◆ ShapeMatricesTypePressure

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
using ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::ShapeMatricesTypePressure = ShapeMatrixPolicyType<ShapeFunctionPressure, GlobalDim>
protected

Definition at line 100 of file HydroMechanicsLocalAssemblerMatrix.h.

Constructor & Destructor Documentation

◆ HydroMechanicsLocalAssemblerMatrix() [1/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::HydroMechanicsLocalAssemblerMatrix ( HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim > const &  )
delete

◆ HydroMechanicsLocalAssemblerMatrix() [2/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::HydroMechanicsLocalAssemblerMatrix ( HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim > &&  )
delete

◆ HydroMechanicsLocalAssemblerMatrix() [3/3]

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::HydroMechanicsLocalAssemblerMatrix ( MeshLib::Element const &  e,
std::size_t const  n_variables,
std::size_t const  local_matrix_size,
std::vector< unsigned > const &  dofIndex_to_localIndex,
bool const  is_axially_symmetric,
unsigned const  integration_order,
HydroMechanicsProcessData< GlobalDim > &  process_data 
)

Definition at line 32 of file HydroMechanicsLocalAssemblerMatrix-impl.h.

42  e, is_axially_symmetric,
43  (n_variables - 1) * ShapeFunctionDisplacement::NPOINTS * GlobalDim +
44  ShapeFunctionPressure::NPOINTS,
45  dofIndex_to_localIndex),
46  _process_data(process_data)
47 {
48  IntegrationMethod integration_method(integration_order);
49  unsigned const n_integration_points =
50  integration_method.getNumberOfPoints();
51 
52  _ip_data.reserve(n_integration_points);
53 
54  auto const shape_matrices_u =
55  NumLib::initShapeMatrices<ShapeFunctionDisplacement,
57  e, is_axially_symmetric, integration_method);
58 
59  auto const shape_matrices_p =
60  NumLib::initShapeMatrices<ShapeFunctionPressure,
61  ShapeMatricesTypePressure, GlobalDim>(
62  e, is_axially_symmetric, integration_method);
63 
65  _process_data.solid_materials, _process_data.material_ids, e.getID());
66 
68  x_position.setElementID(e.getID());
69  for (unsigned ip = 0; ip < n_integration_points; ip++)
70  {
71  x_position.setIntegrationPoint(ip);
72 
73  _ip_data.emplace_back(solid_material);
74  auto& ip_data = _ip_data[ip];
75  auto const& sm_u = shape_matrices_u[ip];
76  auto const& sm_p = shape_matrices_p[ip];
77  ip_data.integration_weight =
78  sm_u.detJ * sm_u.integralMeasure *
79  integration_method.getWeightedPoint(ip).getWeight();
80  ip_data.darcy_velocity.setZero();
81 
82  ip_data.N_u = sm_u.N;
83  ip_data.dNdx_u = sm_u.dNdx;
84  ip_data.H_u.setZero(GlobalDim, displacement_size);
85  for (int i = 0; i < GlobalDim; ++i)
86  {
87  ip_data.H_u
88  .template block<1, displacement_size / GlobalDim>(
89  i, i * displacement_size / GlobalDim)
90  .noalias() = ip_data.N_u;
91  }
92 
93  ip_data.N_p = sm_p.N;
94  ip_data.dNdx_p = sm_p.dNdx;
95 
96  ip_data.sigma_eff.setZero(kelvin_vector_size);
97  ip_data.eps.setZero(kelvin_vector_size);
98 
99  ip_data.sigma_eff_prev.resize(kelvin_vector_size);
100  ip_data.eps_prev.resize(kelvin_vector_size);
101  ip_data.C.resize(kelvin_vector_size, kelvin_vector_size);
102 
103  auto const initial_effective_stress =
104  _process_data.initial_effective_stress(0, x_position);
105  for (unsigned i = 0; i < kelvin_vector_size; i++)
106  {
107  ip_data.sigma_eff[i] = initial_effective_stress[i];
108  ip_data.sigma_eff_prev[i] = initial_effective_stress[i];
109  }
110  }
111 }
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
HydroMechanicsLocalAssemblerInterface(MeshLib::Element const &element, bool const is_axially_symmetric, std::size_t n_local_size, std::vector< unsigned > dofIndex_to_localIndex)
std::vector< IntegrationPointDataType, Eigen::aligned_allocator< IntegrationPointDataType > > _ip_data
ShapeMatrixPolicyType< ShapeFunctionPressure, GlobalDim > ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionDisplacement, GlobalDim > ShapeMatricesTypeDisplacement
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)
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)

References ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::_ip_data, ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::_process_data, ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::displacement_size, MeshLib::Element::getID(), NumLib::initShapeMatrices(), ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::kelvin_vector_size, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), and ParameterLib::SpatialPosition::setIntegrationPoint().

Member Function Documentation

◆ assembleBlockMatricesWithJacobian()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
void ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::assembleBlockMatricesWithJacobian ( double const  t,
double const  dt,
Eigen::Ref< const Eigen::VectorXd > const &  p,
Eigen::Ref< const Eigen::VectorXd > const &  p_dot,
Eigen::Ref< const Eigen::VectorXd > const &  u,
Eigen::Ref< const Eigen::VectorXd > const &  u_dot,
Eigen::Ref< Eigen::VectorXd >  rhs_p,
Eigen::Ref< Eigen::VectorXd >  rhs_u,
Eigen::Ref< Eigen::MatrixXd >  J_pp,
Eigen::Ref< Eigen::MatrixXd >  J_pu,
Eigen::Ref< Eigen::MatrixXd >  J_uu,
Eigen::Ref< Eigen::MatrixXd >  J_up 
)
protected

Definition at line 158 of file HydroMechanicsLocalAssemblerMatrix-impl.h.

168 {
169  assert(this->_element.getDimension() == GlobalDim);
170 
171  typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
172  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
173  pressure_size);
174 
175  typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
176  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
177  pressure_size);
178 
179  typename ShapeMatricesTypeDisplacement::template MatrixType<
181  Kup = ShapeMatricesTypeDisplacement::template MatrixType<
183  pressure_size);
184 
185  auto const& gravity_vec = _process_data.specific_body_force;
186 
187  MPL::VariableArray variables;
188  MPL::VariableArray variables_prev;
190  x_position.setElementID(_element.getID());
191 
192  unsigned const n_integration_points = _ip_data.size();
193  for (unsigned ip = 0; ip < n_integration_points; ip++)
194  {
195  x_position.setIntegrationPoint(ip);
196 
197  auto& ip_data = _ip_data[ip];
198  auto const& ip_w = ip_data.integration_weight;
199  auto const& N_u = ip_data.N_u;
200  auto const& dNdx_u = ip_data.dNdx_u;
201  auto const& N_p = ip_data.N_p;
202  auto const& dNdx_p = ip_data.dNdx_p;
203  auto const& H_u = ip_data.H_u;
204 
205  auto const x_coord =
206  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
208  _element, N_u);
209  auto const B =
211  ShapeFunctionDisplacement::NPOINTS,
212  typename BMatricesType::BMatrixType>(
213  dNdx_u, N_u, x_coord, _is_axially_symmetric);
214 
215  auto const& eps_prev = ip_data.eps_prev;
216  auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
217  auto& sigma_eff = ip_data.sigma_eff;
218 
219  auto& eps = ip_data.eps;
220  auto& state = ip_data.material_state_variables;
221 
222  auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
223  auto const rho_sr = _process_data.solid_density(t, x_position)[0];
224  auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
225  auto const porosity = _process_data.porosity(t, x_position)[0];
226 
227  double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
228  auto const& identity2 =
230 
231  eps.noalias() = B * u;
232 
233  variables[static_cast<int>(
236 
237  variables_prev[static_cast<int>(MaterialPropertyLib::Variable::stress)]
239  sigma_eff_prev);
240  variables_prev[static_cast<int>(
243  eps_prev);
244  variables_prev[static_cast<int>(
246  .emplace<double>(_process_data.reference_temperature);
247 
248  auto&& solution = _ip_data[ip].solid_material.integrateStress(
249  variables_prev, variables, t, x_position, dt, *state);
250 
251  if (!solution)
252  {
253  OGS_FATAL("Computation of local constitutive relation failed.");
254  }
255 
257  std::tie(sigma_eff, state, C) = std::move(*solution);
258 
259  J_uu.noalias() += B.transpose() * C * B * ip_w;
260 
261  rhs_u.noalias() -= B.transpose() * sigma_eff * ip_w;
262  rhs_u.noalias() -= -H_u.transpose() * rho * gravity_vec * ip_w;
263 
264  //
265  // pressure equation, pressure part and displacement equation, pressure
266  // part
267  //
268  if (!_process_data.deactivate_matrix_in_flow) // Only for hydraulically
269  // active matrix
270  {
271  Kup.noalias() += B.transpose() * alpha * identity2 * N_p * ip_w;
272 
273  double const k_over_mu =
274  _process_data.intrinsic_permeability(t, x_position)[0] /
275  _process_data.fluid_viscosity(t, x_position)[0];
276  double const S = _process_data.specific_storage(t, x_position)[0];
277 
278  auto q = ip_data.darcy_velocity.head(GlobalDim);
279  q.noalias() = -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
280 
281  laplace_p.noalias() +=
282  dNdx_p.transpose() * k_over_mu * dNdx_p * ip_w;
283  storage_p.noalias() += N_p.transpose() * S * N_p * ip_w;
284 
285  rhs_p.noalias() +=
286  dNdx_p.transpose() * rho_fr * k_over_mu * gravity_vec * ip_w;
287  }
288  }
289 
290  // displacement equation, pressure part
291  J_up.noalias() -= Kup;
292 
293  // pressure equation, pressure part.
294  J_pp.noalias() += laplace_p + storage_p / dt;
295 
296  // pressure equation, displacement part.
297  J_pu.noalias() += Kup.transpose() / dt;
298 
299  // pressure equation
300  rhs_p.noalias() -=
301  laplace_p * p + storage_p * p_dot + Kup.transpose() * u_dot;
302 
303  // displacement equation
304  rhs_u.noalias() -= -Kup * p;
305 }
#define OGS_FATAL(...)
Definition: Error.h:26
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:82
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:55
@ rho
density. For some models, rho substitutes p
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Definition: VariableType.h:108
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:48
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:56
static const double q
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

References MaterialPropertyLib::alpha, ProcessLib::LinearBMatrix::computeBMatrix(), NumLib::interpolateXCoordinate(), MaterialPropertyLib::mechanical_strain, OGS_FATAL, MathLib::q, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), MaterialPropertyLib::stress, and MaterialPropertyLib::temperature.

◆ assembleWithJacobianConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
void ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::assembleWithJacobianConcrete ( double const  t,
double const  dt,
Eigen::VectorXd const &  local_x,
Eigen::VectorXd const &  local_x_dot,
Eigen::VectorXd &  local_rhs,
Eigen::MatrixXd &  local_Jac 
)
overrideprotectedvirtual

Implements ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerInterface.

Reimplemented in ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrixNearFracture< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >.

Definition at line 117 of file HydroMechanicsLocalAssemblerMatrix-impl.h.

122 {
123  auto p = const_cast<Eigen::VectorXd&>(local_x).segment(pressure_index,
124  pressure_size);
125  auto p_dot = const_cast<Eigen::VectorXd&>(local_x_dot)
126  .segment(pressure_index, pressure_size);
127 
128  if (_process_data.deactivate_matrix_in_flow)
129  {
132  }
133 
134  auto u = local_x.segment(displacement_index, displacement_size);
135  auto u_dot = local_x_dot.segment(displacement_index, displacement_size);
136 
137  auto rhs_p = local_rhs.template segment<pressure_size>(pressure_index);
138  auto rhs_u =
139  local_rhs.template segment<displacement_size>(displacement_index);
140 
141  auto J_pp = local_Jac.template block<pressure_size, pressure_size>(
143  auto J_pu = local_Jac.template block<pressure_size, displacement_size>(
145  auto J_uu = local_Jac.template block<displacement_size, displacement_size>(
147  auto J_up = local_Jac.template block<displacement_size, pressure_size>(
149 
150  assembleBlockMatricesWithJacobian(t, dt, p, p_dot, u, u_dot, rhs_p, rhs_u,
151  J_pp, J_pu, J_uu, J_up);
152 }
void assembleBlockMatricesWithJacobian(double const t, double const dt, Eigen::Ref< const Eigen::VectorXd > const &p, Eigen::Ref< const Eigen::VectorXd > const &p_dot, Eigen::Ref< const Eigen::VectorXd > const &u, Eigen::Ref< const Eigen::VectorXd > const &u_dot, Eigen::Ref< Eigen::VectorXd > rhs_p, Eigen::Ref< Eigen::VectorXd > rhs_u, Eigen::Ref< Eigen::MatrixXd > J_pp, Eigen::Ref< Eigen::MatrixXd > J_pu, Eigen::Ref< Eigen::MatrixXd > J_uu, Eigen::Ref< Eigen::MatrixXd > J_up)

◆ postTimestepConcreteWithBlockVectors()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
void ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::postTimestepConcreteWithBlockVectors ( double const  t,
double const  dt,
Eigen::Ref< const Eigen::VectorXd > const &  p,
Eigen::Ref< const Eigen::VectorXd > const &  u 
)
protected

Definition at line 329 of file HydroMechanicsLocalAssemblerMatrix-impl.h.

334 {
335  MPL::VariableArray variables;
336  MPL::VariableArray variables_prev;
338  x_position.setElementID(_element.getID());
339 
340  unsigned const n_integration_points = _ip_data.size();
341  for (unsigned ip = 0; ip < n_integration_points; ip++)
342  {
343  x_position.setIntegrationPoint(ip);
344 
345  auto& ip_data = _ip_data[ip];
346 
347  auto const& eps_prev = ip_data.eps_prev;
348  auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
349 
350  auto& eps = ip_data.eps;
351  auto& sigma_eff = ip_data.sigma_eff;
352  auto& state = ip_data.material_state_variables;
353 
354  auto const& N_u = ip_data.N_u;
355  auto const& dNdx_u = ip_data.dNdx_u;
356 
357  auto const x_coord =
358  NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
360  _element, N_u);
361  auto const B =
363  ShapeFunctionDisplacement::NPOINTS,
364  typename BMatricesType::BMatrixType>(
365  dNdx_u, N_u, x_coord, _is_axially_symmetric);
366 
367  eps.noalias() = B * u;
368 
369  variables[static_cast<int>(
372 
373  variables_prev[static_cast<int>(MaterialPropertyLib::Variable::stress)]
375  sigma_eff_prev);
376  variables_prev[static_cast<int>(
379  eps_prev);
380  variables_prev[static_cast<int>(
382  .emplace<double>(_process_data.reference_temperature);
383 
384  auto&& solution = _ip_data[ip].solid_material.integrateStress(
385  variables_prev, variables, t, x_position, dt, *state);
386 
387  if (!solution)
388  {
389  OGS_FATAL("Computation of local constitutive relation failed.");
390  }
391 
393  std::tie(sigma_eff, state, C) = std::move(*solution);
394 
395  if (!_process_data.deactivate_matrix_in_flow) // Only for hydraulically
396  // active matrix
397  {
398  double const k_over_mu =
399  _process_data.intrinsic_permeability(t, x_position)[0] /
400  _process_data.fluid_viscosity(t, x_position)[0];
401  auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
402  auto const& gravity_vec = _process_data.specific_body_force;
403  auto const& dNdx_p = ip_data.dNdx_p;
404 
405  ip_data.darcy_velocity.head(GlobalDim).noalias() =
406  -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
407  }
408  }
409 
410  int n = GlobalDim == 2 ? 4 : 6;
411  Eigen::VectorXd ele_stress = Eigen::VectorXd::Zero(n);
412  Eigen::VectorXd ele_strain = Eigen::VectorXd::Zero(n);
413  GlobalDimVector ele_velocity = GlobalDimVector::Zero();
414 
415  for (auto const& ip_data : _ip_data)
416  {
417  ele_stress += ip_data.sigma_eff;
418  ele_strain += ip_data.eps;
419  ele_velocity += ip_data.darcy_velocity;
420  }
421 
422  ele_stress /= static_cast<double>(n_integration_points);
423  ele_strain /= static_cast<double>(n_integration_points);
424  ele_velocity /= static_cast<double>(n_integration_points);
425 
426  auto const element_id = _element.getID();
427  (*_process_data.mesh_prop_stress_xx)[_element.getID()] = ele_stress[0];
428  (*_process_data.mesh_prop_stress_yy)[_element.getID()] = ele_stress[1];
429  (*_process_data.mesh_prop_stress_zz)[_element.getID()] = ele_stress[2];
430  (*_process_data.mesh_prop_stress_xy)[_element.getID()] = ele_stress[3];
431  if (GlobalDim == 3)
432  {
433  (*_process_data.mesh_prop_stress_yz)[_element.getID()] = ele_stress[4];
434  (*_process_data.mesh_prop_stress_xz)[_element.getID()] = ele_stress[5];
435  }
436 
437  (*_process_data.mesh_prop_strain_xx)[_element.getID()] = ele_strain[0];
438  (*_process_data.mesh_prop_strain_yy)[_element.getID()] = ele_strain[1];
439  (*_process_data.mesh_prop_strain_zz)[_element.getID()] = ele_strain[2];
440  (*_process_data.mesh_prop_strain_xy)[_element.getID()] = ele_strain[3];
441  if (GlobalDim == 3)
442  {
443  (*_process_data.mesh_prop_strain_yz)[_element.getID()] = ele_strain[4];
444  (*_process_data.mesh_prop_strain_xz)[_element.getID()] = ele_strain[5];
445  }
446 
447  for (unsigned i = 0; i < GlobalDim; i++)
448  {
449  (*_process_data.mesh_prop_velocity)[element_id * GlobalDim + i] =
450  ele_velocity[i];
451  }
452 
454  ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
455  GlobalDim>(_element, _is_axially_symmetric, p,
456  *_process_data.mesh_prop_nodal_p);
457 }
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 ProcessLib::LinearBMatrix::computeBMatrix(), NumLib::interpolateToHigherOrderNodes(), NumLib::interpolateXCoordinate(), MaterialPropertyLib::mechanical_strain, OGS_FATAL, ParameterLib::SpatialPosition::setElementID(), ParameterLib::SpatialPosition::setIntegrationPoint(), MaterialPropertyLib::stress, and MaterialPropertyLib::temperature.

◆ postTimestepConcreteWithVector()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
void ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::postTimestepConcreteWithVector ( double const  t,
double const  dt,
Eigen::VectorXd const &  local_x 
)
overrideprotectedvirtual

Implements ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerInterface.

Reimplemented in ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrixNearFracture< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >.

Definition at line 311 of file HydroMechanicsLocalAssemblerMatrix-impl.h.

313 {
314  auto p = const_cast<Eigen::VectorXd&>(local_x).segment(pressure_index,
315  pressure_size);
316  if (_process_data.deactivate_matrix_in_flow)
317  {
319  }
320  auto u = local_x.segment(displacement_index, displacement_size);
321 
323 }
void postTimestepConcreteWithBlockVectors(double const t, double const dt, Eigen::Ref< const Eigen::VectorXd > const &p, Eigen::Ref< const Eigen::VectorXd > const &u)

◆ preTimestepConcrete()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
void ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::preTimestepConcrete ( std::vector< double > const &  ,
double const  ,
double const   
)
inlineoverridevirtual

◆ setPressureDotOfInactiveNodes()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
void ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::setPressureDotOfInactiveNodes ( Eigen::Ref< Eigen::VectorXd >  p_dot)
protected

Definition at line 486 of file HydroMechanicsLocalAssemblerMatrix-impl.h.

487 {
488  for (unsigned i = 0; i < pressure_size; i++)
489  {
490  // only inactive nodes
491  if (_process_data.p_element_status->isActiveNode(_element.getNode(i)))
492  {
493  continue;
494  }
495  p_dot[i] = 0;
496  }
497 }
virtual const Node * getNode(unsigned idx) const =0

◆ setPressureOfInactiveNodes()

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
void ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::setPressureOfInactiveNodes ( double const  t,
Eigen::Ref< Eigen::VectorXd >  p 
)
protected

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

466 {
468  x_position.setElementID(_element.getID());
469  for (unsigned i = 0; i < pressure_size; i++)
470  {
471  // only inactive nodes
472  if (_process_data.p_element_status->isActiveNode(_element.getNode(i)))
473  {
474  continue;
475  }
476  x_position.setNodeID(getNodeIndex(_element, i));
477  auto const p0 = (*_process_data.p0)(t, x_position)[0];
478  p[i] = p0;
479  }
480 }
void setNodeID(std::size_t node_id)
std::size_t getNodeIndex(Element const &element, unsigned const idx)
Definition: Element.cpp:225

References MeshLib::getNodeIndex(), ParameterLib::SpatialPosition::setElementID(), and ParameterLib::SpatialPosition::setNodeID().

Member Data Documentation

◆ _ip_data

◆ _process_data

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
HydroMechanicsProcessData<GlobalDim>& ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::_process_data
protected

◆ displacement_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
const int ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::displacement_index = ShapeFunctionPressure::NPOINTS
staticprotected

Definition at line 120 of file HydroMechanicsLocalAssemblerMatrix.h.

◆ displacement_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
const int ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::displacement_size
staticprotected

◆ kelvin_vector_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
const int ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::kelvin_vector_size
staticprotected

◆ pressure_index

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
const int ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::pressure_index = 0
staticprotected

Definition at line 118 of file HydroMechanicsLocalAssemblerMatrix.h.

◆ pressure_size

template<typename ShapeFunctionDisplacement , typename ShapeFunctionPressure , typename IntegrationMethod , int GlobalDim>
const int ProcessLib::LIE::HydroMechanics::HydroMechanicsLocalAssemblerMatrix< ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod, GlobalDim >::pressure_size = ShapeFunctionPressure::NPOINTS
staticprotected

Definition at line 119 of file HydroMechanicsLocalAssemblerMatrix.h.


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