OGS
HydroMechanicsLocalAssemblerMatrix-impl.h
Go to the documentation of this file.
1
11#pragma once
12
22
23namespace ProcessLib
24{
25namespace LIE
26{
27namespace HydroMechanics
28{
29template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
30 int GlobalDim>
31HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
32 ShapeFunctionPressure, GlobalDim>::
33 HydroMechanicsLocalAssemblerMatrix(
34 MeshLib::Element const& e,
35 std::size_t const n_variables,
36 std::size_t const /*local_matrix_size*/,
37 std::vector<unsigned> const& dofIndex_to_localIndex,
38 NumLib::GenericIntegrationMethod const& integration_method,
39 bool const is_axially_symmetric,
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 unsigned const n_integration_points =
49 integration_method.getNumberOfPoints();
50
51 _ip_data.reserve(n_integration_points);
52 _secondary_data.N.resize(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
78 ip_data.integration_weight =
79 sm_u.detJ * sm_u.integralMeasure *
80 integration_method.getWeightedPoint(ip).getWeight();
81 ip_data.darcy_velocity.setZero();
82
83 ip_data.N_u = sm_u.N;
84 ip_data.dNdx_u = sm_u.dNdx;
85 ip_data.H_u.setZero(GlobalDim, displacement_size);
86 for (int i = 0; i < GlobalDim; ++i)
87 {
88 ip_data.H_u
89 .template block<1, displacement_size / GlobalDim>(
90 i, i * displacement_size / GlobalDim)
91 .noalias() = ip_data.N_u;
92 }
93
94 ip_data.N_p = sm_p.N;
95 ip_data.dNdx_p = sm_p.dNdx;
96
97 _secondary_data.N[ip] = sm_u.N;
98
99 ip_data.sigma_eff.setZero(kelvin_vector_size);
100 ip_data.eps.setZero(kelvin_vector_size);
101
102 ip_data.sigma_eff_prev.resize(kelvin_vector_size);
103 ip_data.eps_prev.resize(kelvin_vector_size);
104 ip_data.C.resize(kelvin_vector_size, kelvin_vector_size);
105
106 auto const initial_effective_stress =
107 _process_data.initial_effective_stress(0, x_position);
108 for (unsigned i = 0; i < kelvin_vector_size; i++)
109 {
110 ip_data.sigma_eff[i] = initial_effective_stress[i];
111 ip_data.sigma_eff_prev[i] = initial_effective_stress[i];
112 }
113 }
114}
115
116template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
117 int GlobalDim>
118void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
119 ShapeFunctionPressure, GlobalDim>::
120 assembleWithJacobianConcrete(double const t, double const dt,
121 Eigen::VectorXd const& local_x,
122 Eigen::VectorXd const& local_x_prev,
123 Eigen::VectorXd& local_rhs,
124 Eigen::MatrixXd& local_Jac)
125{
126 auto p = const_cast<Eigen::VectorXd&>(local_x).segment(pressure_index,
127 pressure_size);
128 auto p_prev = const_cast<Eigen::VectorXd&>(local_x_prev)
129 .segment(pressure_index, pressure_size);
130
131 if (_process_data.deactivate_matrix_in_flow)
132 {
133 setPressureOfInactiveNodes(t, p);
134 }
135
136 auto u = local_x.segment(displacement_index, displacement_size);
137 auto u_prev = local_x_prev.segment(displacement_index, displacement_size);
138
139 auto rhs_p = local_rhs.template segment<pressure_size>(pressure_index);
140 auto rhs_u =
141 local_rhs.template segment<displacement_size>(displacement_index);
142
143 auto J_pp = local_Jac.template block<pressure_size, pressure_size>(
144 pressure_index, pressure_index);
145 auto J_pu = local_Jac.template block<pressure_size, displacement_size>(
146 pressure_index, displacement_index);
147 auto J_uu = local_Jac.template block<displacement_size, displacement_size>(
148 displacement_index, displacement_index);
149 auto J_up = local_Jac.template block<displacement_size, pressure_size>(
150 displacement_index, pressure_index);
151
152 assembleBlockMatricesWithJacobian(t, dt, p, p_prev, u, u_prev, rhs_p, rhs_u,
153 J_pp, J_pu, J_uu, J_up);
154}
155
156template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
157 int GlobalDim>
158void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
159 ShapeFunctionPressure, GlobalDim>::
160 assembleBlockMatricesWithJacobian(
161 double const t, double const dt,
162 Eigen::Ref<const Eigen::VectorXd> const& p,
163 Eigen::Ref<const Eigen::VectorXd> const& p_prev,
164 Eigen::Ref<const Eigen::VectorXd> const& u,
165 Eigen::Ref<const Eigen::VectorXd> const& u_prev,
166 Eigen::Ref<Eigen::VectorXd> rhs_p, Eigen::Ref<Eigen::VectorXd> rhs_u,
167 Eigen::Ref<Eigen::MatrixXd> J_pp, Eigen::Ref<Eigen::MatrixXd> J_pu,
168 Eigen::Ref<Eigen::MatrixXd> J_uu, Eigen::Ref<Eigen::MatrixXd> J_up)
169{
170 assert(this->_element.getDimension() == GlobalDim);
171
173 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
174 pressure_size);
175
177 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
178 pressure_size);
179
180 typename ShapeMatricesTypeDisplacement::template MatrixType<
181 displacement_size, pressure_size>
182 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
183 displacement_size, pressure_size>::Zero(displacement_size,
184 pressure_size);
185
186 auto const& gravity_vec = _process_data.specific_body_force;
187
188 MPL::VariableArray variables;
189 MPL::VariableArray variables_prev;
191 x_position.setElementID(_element.getID());
192
193 unsigned const n_integration_points = _ip_data.size();
194 for (unsigned ip = 0; ip < n_integration_points; ip++)
195 {
196 x_position.setIntegrationPoint(ip);
197
198 auto& ip_data = _ip_data[ip];
199 auto const& ip_w = ip_data.integration_weight;
200 auto const& N_u = ip_data.N_u;
201 auto const& dNdx_u = ip_data.dNdx_u;
202 auto const& N_p = ip_data.N_p;
203 auto const& dNdx_p = ip_data.dNdx_p;
204 auto const& H_u = ip_data.H_u;
205
206 auto const x_coord =
207 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
209 _element, N_u);
210 auto const B =
212 ShapeFunctionDisplacement::NPOINTS,
214 dNdx_u, N_u, x_coord, _is_axially_symmetric);
215
216 auto const& eps_prev = ip_data.eps_prev;
217 auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
218 auto& sigma_eff = ip_data.sigma_eff;
219
220 auto& eps = ip_data.eps;
221 auto& state = ip_data.material_state_variables;
222
223 auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
224 auto const rho_sr = _process_data.solid_density(t, x_position)[0];
225 auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
226 auto const porosity = _process_data.porosity(t, x_position)[0];
227
228 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
229 auto const& identity2 =
231
232 eps.noalias() = B * u;
233
234 variables.mechanical_strain
236
237 variables_prev.stress
239 sigma_eff_prev);
240 variables_prev.mechanical_strain
242 eps_prev);
243 variables_prev.temperature = _process_data.reference_temperature;
244
245 auto&& solution = _ip_data[ip].solid_material.integrateStress(
246 variables_prev, variables, t, x_position, dt, *state);
247
248 if (!solution)
249 {
250 OGS_FATAL("Computation of local constitutive relation failed.");
251 }
252
254 std::tie(sigma_eff, state, C) = std::move(*solution);
255
256 J_uu.noalias() += B.transpose() * C * B * ip_w;
257
258 rhs_u.noalias() -= B.transpose() * sigma_eff * ip_w;
259 rhs_u.noalias() -= -H_u.transpose() * rho * gravity_vec * ip_w;
260
261 //
262 // pressure equation, pressure part and displacement equation, pressure
263 // part
264 //
265 if (!_process_data.deactivate_matrix_in_flow) // Only for hydraulically
266 // active matrix
267 {
268 Kup.noalias() += B.transpose() * alpha * identity2 * N_p * ip_w;
269
270 double const k_over_mu =
271 _process_data.intrinsic_permeability(t, x_position)[0] /
272 _process_data.fluid_viscosity(t, x_position)[0];
273 double const S = _process_data.specific_storage(t, x_position)[0];
274
275 auto q = ip_data.darcy_velocity.head(GlobalDim);
276 q.noalias() = -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
277
278 laplace_p.noalias() +=
279 dNdx_p.transpose() * k_over_mu * dNdx_p * ip_w;
280 storage_p.noalias() += N_p.transpose() * S * N_p * ip_w;
281
282 rhs_p.noalias() +=
283 dNdx_p.transpose() * rho_fr * k_over_mu * gravity_vec * ip_w;
284 }
285 }
286
287 // displacement equation, pressure part
288 J_up.noalias() -= Kup;
289
290 // pressure equation, pressure part.
291 J_pp.noalias() += laplace_p + storage_p / dt;
292
293 // pressure equation, displacement part.
294 J_pu.noalias() += Kup.transpose() / dt;
295
296 // pressure equation
297 rhs_p.noalias() -= laplace_p * p + storage_p * (p - p_prev) / dt +
298 Kup.transpose() * (u - u_prev) / dt;
299
300 // displacement equation
301 rhs_u.noalias() -= -Kup * p;
302}
303
304template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
305 int GlobalDim>
307 ShapeFunctionDisplacement, ShapeFunctionPressure,
308 GlobalDim>::postTimestepConcreteWithVector(double const t, double const dt,
309 Eigen::VectorXd const& local_x)
310{
311 auto p = const_cast<Eigen::VectorXd&>(local_x).segment(pressure_index,
312 pressure_size);
313 if (_process_data.deactivate_matrix_in_flow)
314 {
315 setPressureOfInactiveNodes(t, p);
316 }
317 auto u = local_x.segment(displacement_index, displacement_size);
318
319 postTimestepConcreteWithBlockVectors(t, dt, p, u);
320}
321
322template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
323 int GlobalDim>
324void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
325 ShapeFunctionPressure, GlobalDim>::
326 postTimestepConcreteWithBlockVectors(
327 double const t, double const dt,
328 Eigen::Ref<const Eigen::VectorXd> const& p,
329 Eigen::Ref<const Eigen::VectorXd> const& u)
330{
331 MPL::VariableArray variables;
332 MPL::VariableArray variables_prev;
334
335 auto const e_id = _element.getID();
336 x_position.setElementID(e_id);
337
339 KV sigma_avg = KV::Zero();
340 GlobalDimVector velocity_avg;
341 velocity_avg.setZero();
342
343 unsigned const n_integration_points = _ip_data.size();
344 for (unsigned ip = 0; ip < n_integration_points; ip++)
345 {
346 x_position.setIntegrationPoint(ip);
347
348 auto& ip_data = _ip_data[ip];
349
350 auto const& eps_prev = ip_data.eps_prev;
351 auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
352
353 auto& eps = ip_data.eps;
354 auto& sigma_eff = ip_data.sigma_eff;
355 auto& state = ip_data.material_state_variables;
356
357 auto const& N_u = ip_data.N_u;
358 auto const& dNdx_u = ip_data.dNdx_u;
359
360 auto const x_coord =
361 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
363 _element, N_u);
364 auto const B =
366 ShapeFunctionDisplacement::NPOINTS,
368 dNdx_u, N_u, x_coord, _is_axially_symmetric);
369
370 eps.noalias() = B * u;
371
372 variables.mechanical_strain
374
375 variables_prev.stress
377 sigma_eff_prev);
378 variables_prev.mechanical_strain
380 eps_prev);
381 variables_prev.temperature = _process_data.reference_temperature;
382
383 auto&& solution = _ip_data[ip].solid_material.integrateStress(
384 variables_prev, variables, t, x_position, dt, *state);
385
386 if (!solution)
387 {
388 OGS_FATAL("Computation of local constitutive relation failed.");
389 }
390
392 std::tie(sigma_eff, state, C) = std::move(*solution);
393
394 sigma_avg += ip_data.sigma_eff;
395
396 if (!_process_data.deactivate_matrix_in_flow) // Only for hydraulically
397 // active matrix
398 {
399 double const k_over_mu =
400 _process_data.intrinsic_permeability(t, x_position)[0] /
401 _process_data.fluid_viscosity(t, x_position)[0];
402 auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
403 auto const& gravity_vec = _process_data.specific_body_force;
404 auto const& dNdx_p = ip_data.dNdx_p;
405
406 ip_data.darcy_velocity.head(GlobalDim).noalias() =
407 -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
408 velocity_avg += ip_data.darcy_velocity.head(GlobalDim);
409 }
410 }
411
412 sigma_avg /= n_integration_points;
413 velocity_avg /= n_integration_points;
414
415 Eigen::Map<KV>(
416 &(*_process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) =
418
419 Eigen::Map<GlobalDimVector>(
420 &(*_process_data.element_velocities)[e_id * GlobalDim]) = velocity_avg;
421
423 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
424 GlobalDim>(_element, _is_axially_symmetric, p,
425 *_process_data.mesh_prop_nodal_p);
426}
427
428template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
429 int GlobalDim>
430void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
431 ShapeFunctionPressure, GlobalDim>::
432 setPressureOfInactiveNodes(double const t, Eigen::Ref<Eigen::VectorXd> p)
433{
435 x_position.setElementID(_element.getID());
436 for (unsigned i = 0; i < pressure_size; i++)
437 {
438 // only inactive nodes
439 if (_process_data.p_element_status->isActiveNode(_element.getNode(i)))
440 {
441 continue;
442 }
443 x_position.setNodeID(getNodeIndex(_element, i));
444 auto const p0 = (*_process_data.p0)(t, x_position)[0];
445 p[i] = p0;
446 }
447}
448
449template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
450 int GlobalDim>
451std::vector<double> const& HydroMechanicsLocalAssemblerMatrix<
452 ShapeFunctionDisplacement, ShapeFunctionPressure, GlobalDim>::
453 getIntPtSigma(
454 const double /*t*/,
455 std::vector<GlobalVector*> const& /*x*/,
456 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
457 std::vector<double>& cache) const
458{
459 return ProcessLib::getIntegrationPointKelvinVectorData<GlobalDim>(
460 _ip_data, &IntegrationPointDataType::sigma_eff, cache);
461}
462template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
463 int GlobalDim>
464std::vector<double> const& HydroMechanicsLocalAssemblerMatrix<
465 ShapeFunctionDisplacement, ShapeFunctionPressure, GlobalDim>::
466 getIntPtEpsilon(
467 const double /*t*/,
468 std::vector<GlobalVector*> const& /*x*/,
469 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
470 std::vector<double>& cache) const
471{
472 return ProcessLib::getIntegrationPointKelvinVectorData<GlobalDim>(
473 _ip_data, &IntegrationPointDataType::eps, cache);
474}
475
476template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
477 int GlobalDim>
478std::vector<double> const& HydroMechanicsLocalAssemblerMatrix<
479 ShapeFunctionDisplacement, ShapeFunctionPressure, GlobalDim>::
480 getIntPtDarcyVelocity(
481 const double /*t*/,
482 std::vector<GlobalVector*> const& /*x*/,
483 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
484 std::vector<double>& cache) const
485{
486 unsigned const n_integration_points = _ip_data.size();
487
488 cache.clear();
489 auto cache_matrix = MathLib::createZeroedMatrix<
490 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
491 cache, GlobalDim, n_integration_points);
492
493 for (unsigned ip = 0; ip < n_integration_points; ip++)
494 {
495 cache_matrix.col(ip).noalias() = _ip_data[ip].darcy_velocity;
496 }
497
498 return cache;
499}
500} // namespace HydroMechanics
501} // namespace LIE
502} // namespace ProcessLib
Definition of the ElementStatus class.
#define OGS_FATAL(...)
Definition Error.h:26
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > stress
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)
void setIntegrationPoint(unsigned integration_point)
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
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
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 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