26 namespace HydroMechanics
28 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
29 typename IntegrationMethod,
int GlobalDim>
30 HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
31 ShapeFunctionPressure, IntegrationMethod,
33 HydroMechanicsLocalAssemblerMatrix(
35 std::size_t
const n_variables,
37 std::vector<unsigned>
const& dofIndex_to_localIndex,
38 bool const is_axially_symmetric,
39 unsigned const integration_order,
42 e, is_axially_symmetric,
43 (n_variables - 1) * ShapeFunctionDisplacement::NPOINTS * GlobalDim +
44 ShapeFunctionPressure::NPOINTS,
45 dofIndex_to_localIndex),
46 _process_data(process_data)
48 IntegrationMethod integration_method(integration_order);
49 unsigned const n_integration_points =
50 integration_method.getNumberOfPoints();
52 _ip_data.reserve(n_integration_points);
54 auto const shape_matrices_u =
57 e, is_axially_symmetric, integration_method);
59 auto const shape_matrices_p =
62 e, is_axially_symmetric, integration_method);
69 for (
unsigned ip = 0; ip < n_integration_points; ip++)
73 _ip_data.emplace_back(solid_material);
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();
83 ip_data.dNdx_u = sm_u.dNdx;
85 for (
int i = 0; i < GlobalDim; ++i)
90 .noalias() = ip_data.N_u;
94 ip_data.dNdx_p = sm_p.dNdx;
103 auto const initial_effective_stress =
107 ip_data.sigma_eff[i] = initial_effective_stress[i];
108 ip_data.sigma_eff_prev[i] = initial_effective_stress[i];
113 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
114 typename IntegrationMethod,
int GlobalDim>
116 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
117 GlobalDim>::assembleWithJacobianConcrete(
double const t,
double const dt,
118 Eigen::VectorXd
const& local_x,
119 Eigen::VectorXd
const& local_x_dot,
120 Eigen::VectorXd& local_rhs,
121 Eigen::MatrixXd& local_Jac)
123 auto p =
const_cast<Eigen::VectorXd&
>(local_x).segment(pressure_index,
125 auto p_dot =
const_cast<Eigen::VectorXd&
>(local_x_dot)
126 .segment(pressure_index, pressure_size);
128 if (_process_data.deactivate_matrix_in_flow)
130 setPressureOfInactiveNodes(t, p);
131 setPressureDotOfInactiveNodes(p_dot);
134 auto u = local_x.segment(displacement_index, displacement_size);
135 auto u_dot = local_x_dot.segment(displacement_index, displacement_size);
137 auto rhs_p = local_rhs.template segment<pressure_size>(pressure_index);
139 local_rhs.template segment<displacement_size>(displacement_index);
141 auto J_pp = local_Jac.template block<pressure_size, pressure_size>(
142 pressure_index, pressure_index);
143 auto J_pu = local_Jac.template block<pressure_size, displacement_size>(
144 pressure_index, displacement_index);
145 auto J_uu = local_Jac.template block<displacement_size, displacement_size>(
146 displacement_index, displacement_index);
147 auto J_up = local_Jac.template block<displacement_size, pressure_size>(
148 displacement_index, pressure_index);
150 assembleBlockMatricesWithJacobian(t, dt, p, p_dot, u, u_dot, rhs_p, rhs_u,
151 J_pp, J_pu, J_uu, J_up);
154 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
155 typename IntegrationMethod,
int GlobalDim>
157 ShapeFunctionPressure,
158 IntegrationMethod, GlobalDim>::
159 assembleBlockMatricesWithJacobian(
160 double const t,
double const dt,
161 Eigen::Ref<const Eigen::VectorXd>
const& p,
162 Eigen::Ref<const Eigen::VectorXd>
const& p_dot,
163 Eigen::Ref<const Eigen::VectorXd>
const& u,
164 Eigen::Ref<const Eigen::VectorXd>
const& u_dot,
165 Eigen::Ref<Eigen::VectorXd> rhs_p, Eigen::Ref<Eigen::VectorXd> rhs_u,
166 Eigen::Ref<Eigen::MatrixXd> J_pp, Eigen::Ref<Eigen::MatrixXd> J_pu,
167 Eigen::Ref<Eigen::MatrixXd> J_uu, Eigen::Ref<Eigen::MatrixXd> J_up)
169 assert(this->_element.getDimension() == GlobalDim);
172 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
176 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
179 typename ShapeMatricesTypeDisplacement::template MatrixType<
180 displacement_size, pressure_size>
181 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
182 displacement_size, pressure_size>::Zero(displacement_size,
185 auto const& gravity_vec = _process_data.specific_body_force;
192 unsigned const n_integration_points = _ip_data.size();
193 for (
unsigned ip = 0; ip < n_integration_points; ip++)
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;
211 ShapeFunctionDisplacement::NPOINTS,
213 dNdx_u, N_u, x_coord, _is_axially_symmetric);
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;
219 auto& eps = ip_data.eps;
220 auto& state = ip_data.material_state_variables;
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];
227 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
228 auto const& identity2 =
231 eps.noalias() = B * u;
233 variables[
static_cast<int>(
240 variables_prev[
static_cast<int>(
244 variables_prev[
static_cast<int>(
246 .emplace<double>(_process_data.reference_temperature);
248 auto&& solution = _ip_data[ip].solid_material.integrateStress(
249 variables_prev, variables, t, x_position, dt, *state);
253 OGS_FATAL(
"Computation of local constitutive relation failed.");
257 std::tie(sigma_eff, state, C) = std::move(*solution);
259 J_uu.noalias() += B.transpose() * C * B * ip_w;
261 rhs_u.noalias() -= B.transpose() * sigma_eff * ip_w;
262 rhs_u.noalias() -= -H_u.transpose() * rho * gravity_vec * ip_w;
268 if (!_process_data.deactivate_matrix_in_flow)
271 Kup.noalias() += B.transpose() *
alpha * identity2 * N_p * ip_w;
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];
278 auto q = ip_data.darcy_velocity.head(GlobalDim);
279 q.noalias() = -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
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;
286 dNdx_p.transpose() * rho_fr * k_over_mu * gravity_vec * ip_w;
291 J_up.noalias() -= Kup;
294 J_pp.noalias() += laplace_p + storage_p / dt;
297 J_pu.noalias() += Kup.transpose() / dt;
301 laplace_p * p + storage_p * p_dot + Kup.transpose() * u_dot;
304 rhs_u.noalias() -= -Kup * p;
307 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
308 typename IntegrationMethod,
int GlobalDim>
310 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
311 GlobalDim>::postTimestepConcreteWithVector(
double const t,
double const dt,
312 Eigen::VectorXd
const& local_x)
314 auto p =
const_cast<Eigen::VectorXd&
>(local_x).segment(pressure_index,
316 if (_process_data.deactivate_matrix_in_flow)
318 setPressureOfInactiveNodes(t, p);
320 auto u = local_x.segment(displacement_index, displacement_size);
322 postTimestepConcreteWithBlockVectors(t, dt, p, u);
325 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
326 typename IntegrationMethod,
int GlobalDim>
328 ShapeFunctionPressure,
329 IntegrationMethod, GlobalDim>::
330 postTimestepConcreteWithBlockVectors(
331 double const t,
double const dt,
332 Eigen::Ref<const Eigen::VectorXd>
const& p,
333 Eigen::Ref<const Eigen::VectorXd>
const& u)
340 unsigned const n_integration_points = _ip_data.size();
341 for (
unsigned ip = 0; ip < n_integration_points; ip++)
345 auto& ip_data = _ip_data[ip];
347 auto const& eps_prev = ip_data.eps_prev;
348 auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
350 auto& eps = ip_data.eps;
351 auto& sigma_eff = ip_data.sigma_eff;
352 auto& state = ip_data.material_state_variables;
354 auto const& N_u = ip_data.N_u;
355 auto const& dNdx_u = ip_data.dNdx_u;
363 ShapeFunctionDisplacement::NPOINTS,
365 dNdx_u, N_u, x_coord, _is_axially_symmetric);
367 eps.noalias() = B * u;
369 variables[
static_cast<int>(
376 variables_prev[
static_cast<int>(
380 variables_prev[
static_cast<int>(
382 .emplace<double>(_process_data.reference_temperature);
384 auto&& solution = _ip_data[ip].solid_material.integrateStress(
385 variables_prev, variables, t, x_position, dt, *state);
389 OGS_FATAL(
"Computation of local constitutive relation failed.");
393 std::tie(sigma_eff, state, C) = std::move(*solution);
395 if (!_process_data.deactivate_matrix_in_flow)
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;
405 ip_data.darcy_velocity.head(GlobalDim).noalias() =
406 -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
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);
415 for (
auto const& ip_data : _ip_data)
417 ele_stress += ip_data.sigma_eff;
418 ele_strain += ip_data.eps;
419 ele_velocity += ip_data.darcy_velocity;
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);
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];
433 (*_process_data.mesh_prop_stress_yz)[_element.getID()] = ele_stress[4];
434 (*_process_data.mesh_prop_stress_xz)[_element.getID()] = ele_stress[5];
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];
443 (*_process_data.mesh_prop_strain_yz)[_element.getID()] = ele_strain[4];
444 (*_process_data.mesh_prop_strain_xz)[_element.getID()] = ele_strain[5];
447 for (
unsigned i = 0; i < GlobalDim; i++)
449 (*_process_data.mesh_prop_velocity)[element_id * GlobalDim + i] =
454 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
455 GlobalDim>(_element, _is_axially_symmetric, p,
456 *_process_data.mesh_prop_nodal_p);
459 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
460 typename IntegrationMethod,
int GlobalDim>
462 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
463 GlobalDim>::setPressureOfInactiveNodes(
double const t,
464 Eigen::Ref<Eigen::VectorXd>
469 for (
unsigned i = 0; i < pressure_size; i++)
472 if (_process_data.p_element_status->isActiveNode(_element.getNode(i)))
477 auto const p0 = (*_process_data.p0)(t, x_position)[0];
482 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
483 typename IntegrationMethod,
int GlobalDim>
485 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
486 GlobalDim>::setPressureDotOfInactiveNodes(Eigen::Ref<Eigen::VectorXd> p_dot)
488 for (
unsigned i = 0; i < pressure_size; i++)
491 if (_process_data.p_element_status->isActiveNode(_element.getNode(i)))
Definition of the ElementStatus class.
virtual std::size_t getID() const final
Returns the ID of the element.
void setNodeID(std::size_t node_id)
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
std::vector< IntegrationPointDataType, Eigen::aligned_allocator< IntegrationPointDataType > > _ip_data
static const int kelvin_vector_size
ShapeMatrixPolicyType< ShapeFunctionPressure, GlobalDim > ShapeMatricesTypePressure
Eigen::Matrix< double, GlobalDim, 1 > GlobalDimVector
HydroMechanicsProcessData< GlobalDim > & _process_data
static const int displacement_size
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::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
std::size_t getNodeIndex(Element const &element, unsigned const idx)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
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)
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.
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType