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