OGS
HydroMechanicsLocalAssemblerFracture-impl.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <variant>
14
22
23namespace ProcessLib
24{
25namespace LIE
26{
27namespace HydroMechanics
28{
29namespace MPL = MaterialPropertyLib;
30
31template <int DisplacementDim, typename RotationMatrix>
32Eigen::Matrix<double, DisplacementDim, DisplacementDim> createRotatedTensor(
33 RotationMatrix const& R, double const value)
34{
35 using M = Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
36 M tensor = M::Zero();
37 tensor.diagonal().head(DisplacementDim - 1).setConstant(value);
38 return (R.transpose() * tensor * R).eval();
39}
40
41template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
42 int DisplacementDim>
43HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
44 ShapeFunctionPressure, DisplacementDim>::
45 HydroMechanicsLocalAssemblerFracture(
46 MeshLib::Element const& e,
47 std::size_t const /*local_matrix_size*/,
48 std::vector<unsigned> const& dofIndex_to_localIndex,
49 NumLib::GenericIntegrationMethod const& integration_method,
50 bool const is_axially_symmetric,
53 e, is_axially_symmetric, integration_method,
54 ShapeFunctionDisplacement::NPOINTS * DisplacementDim +
55 ShapeFunctionPressure::NPOINTS,
56 dofIndex_to_localIndex),
57 _process_data(process_data)
58{
59 assert(e.getDimension() == DisplacementDim - 1);
60
61 unsigned const n_integration_points =
62 integration_method.getNumberOfPoints();
63
64 _ip_data.reserve(n_integration_points);
65 _secondary_data.N.resize(n_integration_points);
66
67 auto const shape_matrices_u =
68 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
70 DisplacementDim>(e, is_axially_symmetric,
71 integration_method);
72
73 auto const shape_matrices_p =
74 NumLib::initShapeMatrices<ShapeFunctionPressure,
75 ShapeMatricesTypePressure, DisplacementDim>(
76 e, is_axially_symmetric, integration_method);
77
78 auto mat_id = (*_process_data.mesh_prop_materialIDs)[e.getID()];
79 auto frac_id = _process_data.map_materialID_to_fractureID[mat_id];
80 _fracture_property = &_process_data.fracture_properties[frac_id];
81
82 // Get element nodes for aperture0 interpolation from nodes to integration
83 // point. The aperture0 parameter is time-independent.
85 aperture0_node_values =
87 e, /*time independent*/ 0);
88
89 for (unsigned ip = 0; ip < n_integration_points; ip++)
90 {
91 _ip_data.emplace_back(*_process_data.fracture_model);
92 auto const& sm_u = shape_matrices_u[ip];
93 auto const& sm_p = shape_matrices_p[ip];
95 std::nullopt, _element.getID(),
97 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
99 _element, sm_u.N))};
100
101 auto& ip_data = _ip_data[ip];
102 ip_data.integration_weight =
103 sm_u.detJ * sm_u.integralMeasure *
104 integration_method.getWeightedPoint(ip).getWeight();
105
106 ip_data.H_u.setZero(
107 DisplacementDim,
108 ShapeFunctionDisplacement::NPOINTS * DisplacementDim);
110 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
112 HMatrixType>(sm_u.N, ip_data.H_u);
113 ip_data.N_p = sm_p.N;
114 ip_data.dNdx_p = sm_p.dNdx;
115
116 _secondary_data.N[ip] = sm_u.N;
117
118 // Initialize current time step values
119 ip_data.w.setZero(DisplacementDim);
120 ip_data.sigma_eff.setZero(DisplacementDim);
121
122 // Previous time step values are not initialized and are set later.
123 ip_data.w_prev.resize(DisplacementDim);
124 ip_data.sigma_eff_prev.resize(DisplacementDim);
125
126 ip_data.C.resize(DisplacementDim, DisplacementDim);
127
128 ip_data.aperture0 = aperture0_node_values.dot(sm_u.N);
129 ip_data.aperture = ip_data.aperture0;
130
131 auto const initial_effective_stress =
132 _process_data.initial_fracture_effective_stress(0, x_position);
133 for (int i = 0; i < DisplacementDim; i++)
134 {
135 ip_data.sigma_eff[i] = initial_effective_stress[i];
136 ip_data.sigma_eff_prev[i] = initial_effective_stress[i];
137 }
138 }
139}
140
141template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
142 int DisplacementDim>
144 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
145 assembleWithJacobianConcrete(double const t, double const dt,
146 Eigen::VectorXd const& local_x,
147 Eigen::VectorXd const& local_x_prev,
148 Eigen::VectorXd& local_b,
149 Eigen::MatrixXd& local_J)
150{
151 auto const p = local_x.segment(pressure_index, pressure_size);
152 auto const p_prev = local_x_prev.segment(pressure_index, pressure_size);
153 auto const g = local_x.segment(displacement_index, displacement_size);
154 auto const g_prev =
155 local_x_prev.segment(displacement_index, displacement_size);
156
157 auto rhs_p = local_b.segment(pressure_index, pressure_size);
158 auto rhs_g = local_b.segment(displacement_index, displacement_size);
159 auto J_pp = local_J.block(pressure_index, pressure_index, pressure_size,
160 pressure_size);
161 auto J_pg = local_J.block(pressure_index, displacement_index, pressure_size,
162 displacement_size);
163 auto J_gp = local_J.block(displacement_index, pressure_index,
164 displacement_size, pressure_size);
165 auto J_gg = local_J.block(displacement_index, displacement_index,
166 displacement_size, displacement_size);
167
168 assembleBlockMatricesWithJacobian(t, dt, p, p_prev, g, g_prev, rhs_p, rhs_g,
169 J_pp, J_pg, J_gg, J_gp);
170}
171
172template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
173 int DisplacementDim>
175 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
176 assembleBlockMatricesWithJacobian(
177 double const t, double const dt,
178 Eigen::Ref<const Eigen::VectorXd> const& p,
179 Eigen::Ref<const Eigen::VectorXd> const& p_prev,
180 Eigen::Ref<const Eigen::VectorXd> const& g,
181 Eigen::Ref<const Eigen::VectorXd> const& g_prev,
182 Eigen::Ref<Eigen::VectorXd> rhs_p, Eigen::Ref<Eigen::VectorXd> rhs_g,
183 Eigen::Ref<Eigen::MatrixXd> J_pp, Eigen::Ref<Eigen::MatrixXd> J_pg,
184 Eigen::Ref<Eigen::MatrixXd> J_gg, Eigen::Ref<Eigen::MatrixXd> J_gp)
185{
186 auto const& R = _fracture_property->R;
187
188 // the index of a normal (normal to a fracture plane) component
189 // in a displacement vector
190 auto constexpr index_normal = DisplacementDim - 1;
191
193 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
194 pressure_size);
195
197 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
198 pressure_size);
199
200 typename ShapeMatricesTypeDisplacement::template MatrixType<
201 displacement_size, pressure_size>
202 Kgp = ShapeMatricesTypeDisplacement::template MatrixType<
203 displacement_size, pressure_size>::Zero(displacement_size,
204 pressure_size);
205
206 // Projection of the body force vector at the element.
207 Eigen::MatrixXd const global2local_rotation =
208 R.template topLeftCorner<ShapeFunctionPressure::DIM, DisplacementDim>();
209
210 GlobalDimVectorType const gravity_vec = global2local_rotation.transpose() *
211 global2local_rotation *
212 _process_data.specific_body_force;
213
215 x_position.setElementID(_element.getID());
216
217 MPL::VariableArray variables;
218 auto const& medium = _process_data.media_map.getMedium(_element.getID());
219 auto const& liquid_phase = medium->phase("AqueousLiquid");
220 auto const T_ref =
221 medium->property(MPL::PropertyType::reference_temperature)
222 .template value<double>(variables, x_position, t, dt);
223 variables.temperature = T_ref;
224
225 unsigned const n_integration_points = _ip_data.size();
226 for (unsigned ip = 0; ip < n_integration_points; ip++)
227 {
228 auto& ip_data = _ip_data[ip];
229 auto const& ip_w = ip_data.integration_weight;
230 auto const& N_p = ip_data.N_p;
231 auto const& dNdx_p = ip_data.dNdx_p;
232 auto const& H_g = ip_data.H_u;
233 auto const& identity2 =
235
236 x_position = {
237 std::nullopt, _element.getID(),
239 NumLib::interpolateCoordinates<ShapeFunctionPressure,
241 _element, N_p))};
242
243 auto& mat = ip_data.fracture_material;
244 auto& effective_stress = ip_data.sigma_eff;
245 auto const& effective_stress_prev = ip_data.sigma_eff_prev;
246 auto& w = ip_data.w;
247 auto const& w_prev = ip_data.w_prev;
248 auto& C = ip_data.C;
249 auto& state = *ip_data.material_state_variables;
250 auto& b_m = ip_data.aperture;
251
252 auto const rho_fr =
253 liquid_phase.property(MPL::PropertyType::density)
254 .template value<double>(variables, x_position, t, dt);
255 variables.density = rho_fr;
256
257 auto const alpha =
258 medium->property(MPL::PropertyType::biot_coefficient)
259 .template value<double>(variables, x_position, t, dt);
260
261 double const S =
262 medium->property(MPL::PropertyType::storage)
263 .template value<double>(variables, x_position, t, dt);
264
265 auto const mu =
266 liquid_phase.property(MPL::PropertyType::viscosity)
267 .template value<double>(variables, x_position, t, dt);
268
269 // displacement jumps in local coordinates
270 w.noalias() = R * H_g * g;
271
272 // aperture
273 b_m = ip_data.aperture0 + w[index_normal];
274 if (b_m < 0.0)
275 {
276 DBUG(
277 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it must "
278 "be "
279 "non-negative. Setting it to zero.",
280 _element.getID(), ip, b_m);
281 b_m = 0;
282 }
283
284 auto const initial_effective_stress =
285 _process_data.initial_fracture_effective_stress(0, x_position);
286
287 Eigen::Map<typename HMatricesType::ForceVectorType const> const stress0(
288 initial_effective_stress.data(), initial_effective_stress.size());
289
290 // local C, local stress
291 mat.computeConstitutiveRelation(
292 t, x_position, ip_data.aperture0, stress0, w_prev, w,
293 effective_stress_prev, effective_stress, C, state);
294
295 //
296 // displacement equation, displacement jump part
297 //
298 rhs_g.noalias() -=
299 H_g.transpose() * R.transpose() * effective_stress * ip_w;
300 J_gg.noalias() += H_g.transpose() * R.transpose() * C * R * H_g * ip_w;
301
302 //
303 // displacement equation, pressure part
304 //
305 Kgp.noalias() +=
306 H_g.transpose() * R.transpose() * alpha * identity2 * N_p * ip_w;
307
308 //
309 // pressure equation, pressure part.
310 //
311
312 variables.fracture_aperture = b_m;
313 // Assume that the fracture permeability is isotropic
314 auto const permeability =
315 medium->property(MPL::PropertyType::permeability)
316 .value(variables, x_position, t, dt);
317
318 auto& k = ip_data.permeability;
319 k = std::get<double>(permeability);
320 double const k_over_mu = k / mu;
321 storage_p.noalias() += N_p.transpose() * b_m * S * N_p * ip_w;
322 laplace_p.noalias() +=
323 dNdx_p.transpose() * b_m * k_over_mu * dNdx_p * ip_w;
324 rhs_p.noalias() +=
325 dNdx_p.transpose() * b_m * k_over_mu * rho_fr * gravity_vec * ip_w;
326
327 //
328 // pressure equation, displacement jump part.
329 //
330 GlobalDimVectorType const grad_head = dNdx_p * p + rho_fr * gravity_vec;
331 Eigen::Matrix<double, 1, displacement_size> const mT_R_Hg =
332 identity2.transpose() * R * H_g;
333 // velocity in global coordinates
334 ip_data.darcy_velocity = -k_over_mu * grad_head;
335 J_pg.noalias() +=
336 N_p.transpose() * S * N_p * (p - p_prev) / dt * mT_R_Hg * ip_w;
337
338 // derivative of permeability with respect to aperture
339 double const dk_db_over_mu =
340 medium->property(MPL::PropertyType::permeability)
341 .template dValue<double>(variables,
342 MPL::Variable::fracture_aperture,
343 x_position, t, dt) /
344 mu;
345 J_pg.noalias() +=
346 dNdx_p.transpose() * k_over_mu * grad_head * mT_R_Hg * ip_w;
347 J_pg.noalias() += dNdx_p.transpose() * b_m * dk_db_over_mu * grad_head *
348 mT_R_Hg * ip_w;
349 }
350
351 // displacement equation, pressure part
352 J_gp.noalias() -= Kgp;
353
354 // pressure equation, pressure part.
355 J_pp.noalias() += laplace_p + storage_p / dt;
356
357 // pressure equation, displacement jump part.
358 J_pg.noalias() += Kgp.transpose() / dt;
359
360 // pressure equation
361 rhs_p.noalias() -= laplace_p * p + storage_p * (p - p_prev) / dt +
362 Kgp.transpose() * (g - g_prev) / dt;
363
364 // displacement equation
365 rhs_g.noalias() -= -Kgp * p;
366}
367
368template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
369 int DisplacementDim>
371 ShapeFunctionDisplacement, ShapeFunctionPressure,
372 DisplacementDim>::postTimestepConcreteWithVector(const double t,
373 double const /*dt*/,
374 Eigen::VectorXd const&
375 local_x)
376{
377 auto const nodal_g = local_x.segment(displacement_index, displacement_size);
378
379 auto const& R = _fracture_property->R;
380 // the index of a normal (normal to a fracture plane) component
381 // in a displacement vector
382 auto constexpr index_normal = DisplacementDim - 1;
383
385 auto const e_id = _element.getID();
386 x_position.setElementID(e_id);
387
388 unsigned const n_integration_points = _ip_data.size();
389 for (unsigned ip = 0; ip < n_integration_points; ip++)
390 {
391 auto& ip_data = _ip_data[ip];
392 auto const& H_g = ip_data.H_u;
393 auto& mat = ip_data.fracture_material;
394 auto& effective_stress = ip_data.sigma_eff;
395 auto const& effective_stress_prev = ip_data.sigma_eff_prev;
396 auto& w = ip_data.w;
397 auto const& w_prev = ip_data.w_prev;
398 auto& C = ip_data.C;
399 auto& state = *ip_data.material_state_variables;
400 auto& b_m = ip_data.aperture;
401
402 x_position = {
403 std::nullopt, e_id,
405 NumLib::interpolateCoordinates<ShapeFunctionPressure,
407 _element, ip_data.N_p))};
408
409 // displacement jumps in local coordinates
410 w.noalias() = R * H_g * nodal_g;
411
412 // aperture
413 b_m = ip_data.aperture0 + w[index_normal];
414 if (b_m < 0.0)
415 {
416 DBUG(
417 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it is "
418 "expected to be non-negative. Setting it to zero now.",
419 _element.getID(), ip, b_m);
420 b_m = 0;
421 }
422
423 auto const initial_effective_stress =
424 _process_data.initial_fracture_effective_stress(0, x_position);
425
426 Eigen::Map<typename HMatricesType::ForceVectorType const> const stress0(
427 initial_effective_stress.data(), initial_effective_stress.size());
428
429 // local C, local stress
430 mat.computeConstitutiveRelation(
431 t, x_position, ip_data.aperture0, stress0, w_prev, w,
432 effective_stress_prev, effective_stress, C, state);
433 }
434
435 double ele_b = 0;
436 double ele_k = 0;
437 typename HMatricesType::ForceVectorType ele_sigma_eff =
438 HMatricesType::ForceVectorType::Zero(DisplacementDim);
439 typename HMatricesType::ForceVectorType ele_w =
440 HMatricesType::ForceVectorType::Zero(DisplacementDim);
441 GlobalDimVectorType ele_velocity = GlobalDimVectorType::Zero();
442
443 double ele_Fs = -std::numeric_limits<double>::max();
444 for (auto const& ip : _ip_data)
445 {
446 ele_b += ip.aperture;
447 ele_k += ip.permeability;
448 ele_w += ip.w;
449 ele_sigma_eff += ip.sigma_eff;
450 ele_velocity += ip.darcy_velocity;
451 ele_Fs = std::max(
452 ele_Fs, ip.material_state_variables->getShearYieldFunctionValue());
453 }
454 ele_b /= static_cast<double>(n_integration_points);
455 ele_k /= static_cast<double>(n_integration_points);
456 ele_w /= static_cast<double>(n_integration_points);
457 ele_sigma_eff /= static_cast<double>(n_integration_points);
458 ele_velocity /= static_cast<double>(n_integration_points);
459 auto const element_id = _element.getID();
460 (*_process_data.mesh_prop_b)[element_id] = ele_b;
461 (*_process_data.mesh_prop_k_f)[element_id] = ele_k;
462
463 Eigen::Map<GlobalDimVectorType>(
464 &(*_process_data.element_fracture_stresses)[e_id * DisplacementDim]) =
465 ele_sigma_eff;
466
467 Eigen::Map<GlobalDimVectorType>(
468 &(*_process_data.element_fracture_velocities)[e_id * DisplacementDim]) =
469 ele_velocity;
470
471 Eigen::Map<GlobalDimVectorType>(
472 &(*_process_data.element_local_jumps)[e_id * DisplacementDim]) = ele_w;
473
474 (*_process_data.mesh_prop_fracture_shear_failure)[element_id] = ele_Fs;
475}
476
477template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
478 int DisplacementDim>
479std::vector<double> const& HydroMechanicsLocalAssemblerFracture<
480 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
481 getIntPtFractureVelocity(
482 const double /*t*/,
483 std::vector<GlobalVector*> const& /*x*/,
484 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
485 std::vector<double>& cache) const
486{
487 unsigned const n_integration_points = _ip_data.size();
488 cache.clear();
489 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
490 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
491 cache, DisplacementDim, 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
501template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
502 int DisplacementDim>
503std::vector<double> const& HydroMechanicsLocalAssemblerFracture<
504 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
505 getIntPtFractureStress(
506 const double /*t*/,
507 std::vector<GlobalVector*> const& /*x*/,
508 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
509 std::vector<double>& cache) const
510{
511 unsigned const n_integration_points = _ip_data.size();
512 cache.clear();
513 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
514 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
515 cache, DisplacementDim, n_integration_points);
516
517 for (unsigned ip = 0; ip < n_integration_points; ip++)
518 {
519 cache_matrix.col(ip).noalias() = _ip_data[ip].sigma_eff;
520 }
521
522 return cache;
523}
524
525template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
526 int DisplacementDim>
527std::vector<double> const& HydroMechanicsLocalAssemblerFracture<
528 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
529 getIntPtFractureAperture(
530 const double /*t*/,
531 std::vector<GlobalVector*> const& /*x*/,
532 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
533 std::vector<double>& cache) const
534{
536 _ip_data, &IntegrationPointDataType::aperture, cache);
537}
538
539template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
540 int DisplacementDim>
541std::vector<double> const& HydroMechanicsLocalAssemblerFracture<
542 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
543 getIntPtFracturePermeability(
544 const double /*t*/,
545 std::vector<GlobalVector*> const& /*x*/,
546 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
547 std::vector<double>& cache) const
548{
550 _ip_data, &IntegrationPointDataType::permeability, cache);
551}
552
553} // namespace HydroMechanics
554} // namespace LIE
555} // namespace ProcessLib
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
constexpr double getWeight() const
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:91
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
void setElementID(std::size_t element_id)
VectorType< DisplacementDim > ForceVectorType
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
std::vector< IntegrationPointDataType, Eigen::aligned_allocator< IntegrationPointDataType > > _ip_data
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
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)
Eigen::Matrix< double, DisplacementDim, DisplacementDim > createRotatedTensor(RotationMatrix const &R, double const value)
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
void computeHMatrix(N_Type const &N, HMatrixType &H)
Fills a H-matrix based on given shape function.
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
virtual Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > getNodalValuesOnElement(MeshLib::Element const &element, double const t) const
Returns a matrix of values for all nodes of the given element.
Definition Parameter.h:164
ParameterLib::Parameter< double > const & aperture0
Initial aperture.
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N