OGS
HydroMechanicsLocalAssemblerMatrix-impl.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <limits>
14
27
28namespace ProcessLib
29{
30namespace LIE
31{
32namespace HydroMechanics
33{
34namespace MPL = MaterialPropertyLib;
35
36template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
37 int GlobalDim>
38HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
39 ShapeFunctionPressure, GlobalDim>::
40 HydroMechanicsLocalAssemblerMatrix(
41 MeshLib::Element const& e,
42 std::size_t const n_variables,
43 std::size_t const /*local_matrix_size*/,
44 std::vector<unsigned> const& dofIndex_to_localIndex,
45 NumLib::GenericIntegrationMethod const& integration_method,
46 bool const is_axially_symmetric,
49 e, is_axially_symmetric, integration_method,
50 (n_variables - 1) * ShapeFunctionDisplacement::NPOINTS * GlobalDim +
51 ShapeFunctionPressure::NPOINTS,
52 dofIndex_to_localIndex),
53 _process_data(process_data)
54{
55 unsigned const n_integration_points =
56 integration_method.getNumberOfPoints();
57
58 _ip_data.reserve(n_integration_points);
59 _secondary_data.N.resize(n_integration_points);
60
61 auto const shape_matrices_u =
62 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
64 e, is_axially_symmetric, integration_method);
65
66 auto const shape_matrices_p =
67 NumLib::initShapeMatrices<ShapeFunctionPressure,
68 ShapeMatricesTypePressure, GlobalDim>(
69 e, is_axially_symmetric, integration_method);
70
72 _process_data.solid_materials, _process_data.material_ids, e.getID());
73
75 x_position.setElementID(e.getID());
76 for (unsigned ip = 0; ip < n_integration_points; ip++)
77 {
78 _ip_data.emplace_back(solid_material);
79 auto& ip_data = _ip_data[ip];
80 auto const& sm_u = shape_matrices_u[ip];
81 auto const& sm_p = shape_matrices_p[ip];
82
83 ip_data.integration_weight =
84 sm_u.detJ * sm_u.integralMeasure *
85 integration_method.getWeightedPoint(ip).getWeight();
86 ip_data.darcy_velocity.setZero();
87
88 ip_data.N_u = sm_u.N;
89 ip_data.dNdx_u = sm_u.dNdx;
90 ip_data.H_u.setZero(GlobalDim, displacement_size);
91 for (int i = 0; i < GlobalDim; ++i)
92 {
93 ip_data.H_u
94 .template block<1, displacement_size / GlobalDim>(
95 i, i * displacement_size / GlobalDim)
96 .noalias() = ip_data.N_u;
97 }
98
99 ip_data.N_p = sm_p.N;
100 ip_data.dNdx_p = sm_p.dNdx;
101
102 _secondary_data.N[ip] = sm_u.N;
103
104 ip_data.sigma_eff.setZero(kelvin_vector_size);
105 ip_data.eps.setZero(kelvin_vector_size);
106
107 ip_data.sigma_eff_prev.resize(kelvin_vector_size);
108 ip_data.eps_prev.resize(kelvin_vector_size);
109 ip_data.C.resize(kelvin_vector_size, kelvin_vector_size);
110
111 auto const initial_effective_stress =
112 _process_data.initial_effective_stress(0, x_position);
113 for (unsigned i = 0; i < kelvin_vector_size; i++)
114 {
115 ip_data.sigma_eff[i] = initial_effective_stress[i];
116 ip_data.sigma_eff_prev[i] = initial_effective_stress[i];
117 }
118 }
119}
120
121template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
122 int GlobalDim>
123void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
124 ShapeFunctionPressure, GlobalDim>::
125 assembleWithJacobianConcrete(double const t, double const dt,
126 Eigen::VectorXd const& local_x,
127 Eigen::VectorXd const& local_x_prev,
128 Eigen::VectorXd& local_rhs,
129 Eigen::MatrixXd& local_Jac)
130{
131 auto p = const_cast<Eigen::VectorXd&>(local_x).segment(pressure_index,
132 pressure_size);
133 auto const p_prev = local_x_prev.segment(pressure_index, pressure_size);
134
135 if (_process_data.deactivate_matrix_in_flow)
136 {
137 setPressureOfInactiveNodes(t, p);
138 }
139
140 auto u = local_x.segment(displacement_index, displacement_size);
141 auto u_prev = local_x_prev.segment(displacement_index, displacement_size);
142
143 auto rhs_p = local_rhs.template segment<pressure_size>(pressure_index);
144 auto rhs_u =
145 local_rhs.template segment<displacement_size>(displacement_index);
146
147 auto J_pp = local_Jac.template block<pressure_size, pressure_size>(
148 pressure_index, pressure_index);
149 auto J_pu = local_Jac.template block<pressure_size, displacement_size>(
150 pressure_index, displacement_index);
151 auto J_uu = local_Jac.template block<displacement_size, displacement_size>(
152 displacement_index, displacement_index);
153 auto J_up = local_Jac.template block<displacement_size, pressure_size>(
154 displacement_index, pressure_index);
155
156 assembleBlockMatricesWithJacobian(t, dt, p, p_prev, u, u_prev, rhs_p, rhs_u,
157 J_pp, J_pu, J_uu, J_up);
158}
159
160template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
161 int GlobalDim>
162void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
163 ShapeFunctionPressure, GlobalDim>::
164 assembleBlockMatricesWithJacobian(
165 double const t, double const dt,
166 Eigen::Ref<const Eigen::VectorXd> const& p,
167 Eigen::Ref<const Eigen::VectorXd> const& p_prev,
168 Eigen::Ref<const Eigen::VectorXd> const& u,
169 Eigen::Ref<const Eigen::VectorXd> const& u_prev,
170 Eigen::Ref<Eigen::VectorXd> rhs_p, Eigen::Ref<Eigen::VectorXd> rhs_u,
171 Eigen::Ref<Eigen::MatrixXd> J_pp, Eigen::Ref<Eigen::MatrixXd> J_pu,
172 Eigen::Ref<Eigen::MatrixXd> J_uu, Eigen::Ref<Eigen::MatrixXd> J_up)
173{
174 assert(this->_element.getDimension() == GlobalDim);
175
177 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
178 pressure_size);
179
181 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
182 pressure_size);
183
184 typename ShapeMatricesTypeDisplacement::template MatrixType<
185 displacement_size, pressure_size>
186 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
187 displacement_size, pressure_size>::Zero(displacement_size,
188 pressure_size);
189
190 auto const& gravity_vec = _process_data.specific_body_force;
191
192 MPL::VariableArray variables;
193 MPL::VariableArray variables_prev;
195 x_position.setElementID(_element.getID());
196
197 auto const& medium = _process_data.media_map.getMedium(_element.getID());
198 auto const& liquid_phase = medium->phase("AqueousLiquid");
199 auto const& solid_phase = medium->phase("Solid");
200
201 auto const T_ref =
202 medium->property(MPL::PropertyType::reference_temperature)
203 .template value<double>(variables, x_position, t, dt);
204 variables.temperature = T_ref;
205 variables_prev.temperature = T_ref;
206
207 bool const has_storage_property =
208 medium->hasProperty(MPL::PropertyType::storage);
209
210 auto const B_dil_bar = getDilatationalBBarMatrix();
211
212 unsigned const n_integration_points = _ip_data.size();
213 for (unsigned ip = 0; ip < n_integration_points; ip++)
214 {
215 auto& ip_data = _ip_data[ip];
216 auto const& ip_w = ip_data.integration_weight;
217 auto const& N_u = ip_data.N_u;
218 auto const& dNdx_u = ip_data.dNdx_u;
219 auto const& N_p = ip_data.N_p;
220 auto const& dNdx_p = ip_data.dNdx_p;
221 auto const& H_u = ip_data.H_u;
222
223 variables.liquid_phase_pressure = N_p.dot(p);
224
225 auto const x_coord =
226 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
228 _element, N_u);
229 auto const B =
231 GlobalDim, ShapeFunctionDisplacement::NPOINTS, BBarMatrixType,
233 dNdx_u, N_u, B_dil_bar, x_coord, this->_is_axially_symmetric)
234 .eval();
235
236 auto const& eps_prev = ip_data.eps_prev;
237 auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
238 auto& sigma_eff = ip_data.sigma_eff;
239
240 auto& eps = ip_data.eps;
241 auto& state = ip_data.material_state_variables;
242
243 auto const rho_sr =
244 solid_phase.property(MPL::PropertyType::density)
245 .template value<double>(variables, x_position, t, dt);
246 auto const rho_fr =
247 liquid_phase.property(MPL::PropertyType::density)
248 .template value<double>(variables, x_position, t, dt);
249
250 auto const alpha =
251 medium->property(MPL::PropertyType::biot_coefficient)
252 .template value<double>(variables, x_position, t, dt);
253 auto const porosity =
254 medium->property(MPL::PropertyType::porosity)
255 .template value<double>(variables, x_position, t, dt);
256
257 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
258 auto const& identity2 =
260
261 eps.noalias() = B * u;
262
263 variables.mechanical_strain
265
266 variables_prev.stress
268 sigma_eff_prev);
269 variables_prev.mechanical_strain
271 eps_prev);
272 variables_prev.temperature = T_ref;
273
274 auto&& solution = _ip_data[ip].solid_material.integrateStress(
275 variables_prev, variables, t, x_position, dt, *state);
276
277 if (!solution)
278 {
279 OGS_FATAL("Computation of local constitutive relation failed.");
280 }
281
283 std::tie(sigma_eff, state, C) = std::move(*solution);
284
285 J_uu.noalias() += B.transpose() * C * B * ip_w;
286
287 rhs_u.noalias() -= B.transpose() * sigma_eff * ip_w;
288 rhs_u.noalias() -= -H_u.transpose() * rho * gravity_vec * ip_w;
289
290 //
291 // pressure equation, pressure part and displacement equation, pressure
292 // part
293 //
294 if (!_process_data.deactivate_matrix_in_flow) // Only for hydraulically
295 // active matrix
296 {
297 Kup.noalias() += B.transpose() * alpha * identity2 * N_p * ip_w;
298
299 variables.density = rho_fr;
300 auto const mu =
301 liquid_phase.property(MPL::PropertyType::viscosity)
302 .template value<double>(variables, x_position, t, dt);
303
304 auto const k_over_mu =
306 medium->property(MPL::PropertyType::permeability)
307 .value(variables, x_position, t, dt)) /
308 mu)
309 .eval();
310
311 double const S =
312 has_storage_property
313 ? medium->property(MPL::PropertyType::storage)
314 .template value<double>(variables, x_position, t, dt)
315 : porosity *
316 (liquid_phase.property(MPL::PropertyType::density)
317 .template dValue<double>(
318 variables,
319 MPL::Variable::liquid_phase_pressure,
320 x_position, t, dt) /
321 rho_fr) +
322 (alpha - porosity) * (1.0 - alpha) /
323 _ip_data[ip].solid_material.getBulkModulus(
324 t, x_position, &C);
325
326 auto q = ip_data.darcy_velocity.head(GlobalDim);
327 q.noalias() = -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
328
329 laplace_p.noalias() +=
330 dNdx_p.transpose() * k_over_mu * dNdx_p * ip_w;
331 storage_p.noalias() += N_p.transpose() * S * N_p * ip_w;
332
333 rhs_p.noalias() +=
334 dNdx_p.transpose() * rho_fr * k_over_mu * gravity_vec * ip_w;
335 }
336 }
337
338 // displacement equation, pressure part
339 J_up.noalias() -= Kup;
340
341 // pressure equation, pressure part.
342 J_pp.noalias() += laplace_p + storage_p / dt;
343
344 // pressure equation, displacement part.
345 J_pu.noalias() += Kup.transpose() / dt;
346
347 // pressure equation
348 rhs_p.noalias() -= laplace_p * p + storage_p * (p - p_prev) / dt +
349 Kup.transpose() * (u - u_prev) / dt;
350
351 // displacement equation
352 rhs_u.noalias() -= -Kup * p;
353}
354
355template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
356 int GlobalDim>
358 ShapeFunctionDisplacement, ShapeFunctionPressure,
359 GlobalDim>::postTimestepConcreteWithVector(double const t, double const dt,
360 Eigen::VectorXd const& local_x)
361{
362 auto p = const_cast<Eigen::VectorXd&>(local_x).segment(pressure_index,
363 pressure_size);
364 if (_process_data.deactivate_matrix_in_flow)
365 {
366 setPressureOfInactiveNodes(t, p);
367 }
368 auto u = local_x.segment(displacement_index, displacement_size);
369
370 postTimestepConcreteWithBlockVectors(t, dt, p, u);
371}
372
373template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
374 int GlobalDim>
375void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
376 ShapeFunctionPressure, GlobalDim>::
377 postTimestepConcreteWithBlockVectors(
378 double const t, double const dt,
379 Eigen::Ref<const Eigen::VectorXd> const& p,
380 Eigen::Ref<const Eigen::VectorXd> const& u)
381{
382 MPL::VariableArray variables;
383 MPL::VariableArray variables_prev;
385
386 auto const e_id = _element.getID();
387 x_position.setElementID(e_id);
388
390 KV sigma_avg = KV::Zero();
391 GlobalDimVector velocity_avg;
392 velocity_avg.setZero();
393
394 unsigned const n_integration_points = _ip_data.size();
395
396 auto const& medium = _process_data.media_map.getMedium(_element.getID());
397 auto const& liquid_phase = medium->phase("AqueousLiquid");
398 auto const T_ref =
399 medium->property(MPL::PropertyType::reference_temperature)
400 .template value<double>(variables, x_position, t, dt);
401 variables.temperature = T_ref;
402 variables_prev.temperature = T_ref;
403
404 auto const B_dil_bar = getDilatationalBBarMatrix();
405
406 for (unsigned ip = 0; ip < n_integration_points; ip++)
407 {
408 auto& ip_data = _ip_data[ip];
409
410 auto const& eps_prev = ip_data.eps_prev;
411 auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
412
413 auto& eps = ip_data.eps;
414 auto& sigma_eff = ip_data.sigma_eff;
415 auto& state = ip_data.material_state_variables;
416
417 auto const& N_u = ip_data.N_u;
418 auto const& N_p = ip_data.N_p;
419 auto const& dNdx_u = ip_data.dNdx_u;
420
421 variables.liquid_phase_pressure = N_p.dot(p);
422
423 auto const x_coord =
424 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
426 _element, N_u);
427 auto const B =
429 GlobalDim, ShapeFunctionDisplacement::NPOINTS, BBarMatrixType,
431 dNdx_u, N_u, B_dil_bar, x_coord, this->_is_axially_symmetric)
432 .eval();
433
434 eps.noalias() = B * u;
435
436 variables.mechanical_strain
438
439 variables_prev.stress
441 sigma_eff_prev);
442 variables_prev.mechanical_strain
444 eps_prev);
445
446 auto&& solution = _ip_data[ip].solid_material.integrateStress(
447 variables_prev, variables, t, x_position, dt, *state);
448
449 if (!solution)
450 {
451 OGS_FATAL("Computation of local constitutive relation failed.");
452 }
453
455 std::tie(sigma_eff, state, C) = std::move(*solution);
456
457 sigma_avg += ip_data.sigma_eff;
458
459 if (!_process_data.deactivate_matrix_in_flow) // Only for hydraulically
460 // active matrix
461 {
462 auto const rho_fr =
463 liquid_phase.property(MPL::PropertyType::density)
464 .template value<double>(variables, x_position, t, dt);
465 variables.density = rho_fr;
466 auto const mu =
467 liquid_phase.property(MPL::PropertyType::viscosity)
468 .template value<double>(variables, x_position, t, dt);
469
470 auto const k = MPL::formEigenTensor<GlobalDim>(
471 medium->property(MPL::PropertyType::permeability)
472 .value(variables, x_position, t, dt))
473 .eval();
474
475 GlobalDimMatrixType const k_over_mu = k / mu;
476
477 auto const& gravity_vec = _process_data.specific_body_force;
478 auto const& dNdx_p = ip_data.dNdx_p;
479
480 ip_data.darcy_velocity.head(GlobalDim).noalias() =
481 -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
482 velocity_avg.noalias() += ip_data.darcy_velocity.head(GlobalDim);
483 }
484 }
485
486 sigma_avg /= n_integration_points;
487 velocity_avg /= n_integration_points;
488
489 Eigen::Map<KV>(
490 &(*_process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) =
492
493 Eigen::Map<GlobalDimVector>(
494 &(*_process_data.element_velocities)[e_id * GlobalDim]) = velocity_avg;
495
497 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
498 GlobalDim>(_element, _is_axially_symmetric, p,
499 *_process_data.mesh_prop_nodal_p);
500}
501
502template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
503 int GlobalDim>
504void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
505 ShapeFunctionPressure, GlobalDim>::
506 setPressureOfInactiveNodes(double const t, Eigen::Ref<Eigen::VectorXd> p)
507{
509 x_position.setElementID(_element.getID());
510 for (unsigned i = 0; i < pressure_size; i++)
511 {
512 // only inactive nodes
513 if (_process_data.p_element_status->isActiveNode(_element.getNode(i)))
514 {
515 continue;
516 }
517 x_position.setNodeID(getNodeIndex(_element, i));
518 auto const p0 = (*_process_data.p0)(t, x_position)[0];
519 p[i] = p0;
520 }
521}
522
523template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
524 int GlobalDim>
525std::vector<double> const& HydroMechanicsLocalAssemblerMatrix<
526 ShapeFunctionDisplacement, ShapeFunctionPressure, GlobalDim>::
527 getIntPtSigma(
528 const double /*t*/,
529 std::vector<GlobalVector*> const& /*x*/,
530 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
531 std::vector<double>& cache) const
532{
534 _ip_data, &IntegrationPointDataType::sigma_eff, cache);
535}
536template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
537 int GlobalDim>
538std::vector<double> const& HydroMechanicsLocalAssemblerMatrix<
539 ShapeFunctionDisplacement, ShapeFunctionPressure, GlobalDim>::
540 getIntPtEpsilon(
541 const double /*t*/,
542 std::vector<GlobalVector*> const& /*x*/,
543 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
544 std::vector<double>& cache) const
545{
547 _ip_data, &IntegrationPointDataType::eps, cache);
548}
549
550template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
551 int GlobalDim>
552std::vector<double> const& HydroMechanicsLocalAssemblerMatrix<
553 ShapeFunctionDisplacement, ShapeFunctionPressure, GlobalDim>::
554 getIntPtDarcyVelocity(
555 const double /*t*/,
556 std::vector<GlobalVector*> const& /*x*/,
557 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
558 std::vector<double>& cache) const
559{
560 unsigned const n_integration_points = _ip_data.size();
561
562 cache.clear();
563 auto cache_matrix = MathLib::createZeroedMatrix<
564 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
565 cache, GlobalDim, n_integration_points);
566
567 for (unsigned ip = 0; ip < n_integration_points; ip++)
568 {
569 cache_matrix.col(ip).noalias() = _ip_data[ip].darcy_velocity;
570 }
571
572 return cache;
573}
574} // namespace HydroMechanics
575} // namespace LIE
576} // namespace ProcessLib
Definition of the ElementStatus class.
#define OGS_FATAL(...)
Definition Error.h:26
constexpr double getWeight() const
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
void setNodeID(std::size_t node_id)
void setElementID(std::size_t element_id)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
ShapeMatrixPolicyType< ShapeFunctionDisplacement, GlobalDim > ShapeMatricesTypeDisplacement
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
std::vector< IntegrationPointDataType, Eigen::aligned_allocator< IntegrationPointDataType > > _ip_data
ShapeMatrixPolicyType< ShapeFunctionPressure, GlobalDim > ShapeMatricesTypePressure
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
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
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
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 computeBMatrixPossiblyWithBbar(DNDX_Type const &dNdx, N_Type const &N, std::optional< BBarMatrixType > const &B_dil_bar, const double radius, const bool is_axially_symmetric)
Fills a B matrix, or a B bar matrix if required.
std::vector< double > const & getIntegrationPointKelvinVectorData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N