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