OGS 6.2.1-97-g73d1aeda3
HydroMechanicsLocalAssemblerFracture-impl.h
Go to the documentation of this file.
1 
10 #pragma once
11 
13 
15 
17 
19 
20 namespace ProcessLib
21 {
22 namespace LIE
23 {
24 namespace HydroMechanics
25 {
26 template <int GlobalDim, typename RotationMatrix>
27 Eigen::Matrix<double, GlobalDim, GlobalDim> createRotatedTensor(
28  RotationMatrix const& R, double const value)
29 {
30  using M = Eigen::Matrix<double, GlobalDim, GlobalDim>;
31  M tensor = M::Zero();
32  tensor.diagonal().head(GlobalDim - 1).setConstant(value);
33  return (R.transpose() * tensor * R).eval();
34 }
35 
36 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
37  typename IntegrationMethod, int GlobalDim>
38 HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
39  ShapeFunctionPressure, IntegrationMethod,
40  GlobalDim>::
41  HydroMechanicsLocalAssemblerFracture(
42  MeshLib::Element const& e,
43  std::size_t const /*local_matrix_size*/,
44  std::vector<unsigned> const& dofIndex_to_localIndex,
45  bool const is_axially_symmetric,
46  unsigned const integration_order,
49  e, is_axially_symmetric,
50  ShapeFunctionDisplacement::NPOINTS * GlobalDim +
51  ShapeFunctionPressure::NPOINTS,
52  dofIndex_to_localIndex),
53  _process_data(process_data)
54 {
55  assert(e.getDimension() == GlobalDim - 1);
56 
57  IntegrationMethod integration_method(integration_order);
58  unsigned const n_integration_points =
59  integration_method.getNumberOfPoints();
60 
61  _ip_data.reserve(n_integration_points);
62 
63  auto const shape_matrices_u =
64  initShapeMatrices<ShapeFunctionDisplacement,
65  ShapeMatricesTypeDisplacement, IntegrationMethod,
66  GlobalDim>(e, is_axially_symmetric,
67  integration_method);
68 
69  auto const shape_matrices_p =
70  initShapeMatrices<ShapeFunctionPressure, ShapeMatricesTypePressure,
71  IntegrationMethod, GlobalDim>(e, is_axially_symmetric,
72  integration_method);
73 
74  auto const& frac_prop = *_process_data.fracture_property;
75 
76  // Get element nodes for aperture0 interpolation from nodes to integration
77  // point. The aperture0 parameter is time-independent.
79  aperture0_node_values = frac_prop.aperture0.getNodalValuesOnElement(
80  e, /*time independent*/ 0);
81 
83  x_position.setElementID(e.getID());
84  for (unsigned ip = 0; ip < n_integration_points; ip++)
85  {
86  x_position.setIntegrationPoint(ip);
87 
88  _ip_data.emplace_back(*_process_data.fracture_model);
89  auto const& sm_u = shape_matrices_u[ip];
90  auto const& sm_p = shape_matrices_p[ip];
91  auto& ip_data = _ip_data[ip];
92  ip_data.integration_weight =
93  sm_u.detJ * sm_u.integralMeasure *
94  integration_method.getWeightedPoint(ip).getWeight();
95 
96  ip_data.H_u.setZero(GlobalDim,
97  ShapeFunctionDisplacement::NPOINTS * GlobalDim);
99  GlobalDim, ShapeFunctionDisplacement::NPOINTS,
101  HMatrixType>(sm_u.N, ip_data.H_u);
102  ip_data.N_p = sm_p.N;
103  ip_data.dNdx_p = sm_p.dNdx;
104 
105  // Initialize current time step values
106  ip_data.w.setZero(GlobalDim);
107  ip_data.sigma_eff.setZero(GlobalDim);
108 
109  // Previous time step values are not initialized and are set later.
110  ip_data.w_prev.resize(GlobalDim);
111  ip_data.sigma_eff_prev.resize(GlobalDim);
112 
113  ip_data.C.resize(GlobalDim, GlobalDim);
114 
115  ip_data.aperture0 = aperture0_node_values.dot(sm_u.N);
116  ip_data.aperture = ip_data.aperture0;
117 
118  ip_data.permeability_state =
119  frac_prop.permeability_model->getNewState();
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 
131 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
132  typename IntegrationMethod, int GlobalDim>
134  ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
135  GlobalDim>::assembleWithJacobianConcrete(double const t,
136  Eigen::VectorXd const& local_x,
137  Eigen::VectorXd const& local_xdot,
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_dot = local_xdot.segment(pressure_index, pressure_size);
143  auto const g = local_x.segment(displacement_index, displacement_size);
144  auto const g_dot =
145  local_xdot.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,
153  auto J_gp = local_J.block(displacement_index, pressure_index,
155  auto J_gg = local_J.block(displacement_index, displacement_index,
157 
158  assembleBlockMatricesWithJacobian(t, p, p_dot, g, g_dot, rhs_p, rhs_g, J_pp,
159  J_pg, J_gg, J_gp);
160 }
161 
162 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
163  typename IntegrationMethod, int GlobalDim>
164 void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
165  ShapeFunctionPressure,
166  IntegrationMethod, GlobalDim>::
168  double const t, Eigen::Ref<const Eigen::VectorXd> const& p,
169  Eigen::Ref<const Eigen::VectorXd> const& p_dot,
170  Eigen::Ref<const Eigen::VectorXd> const& g,
171  Eigen::Ref<const Eigen::VectorXd> const& g_dot,
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  double const& dt = _process_data.dt;
179 
180  // the index of a normal (normal to a fracture plane) component
181  // in a displacement vector
182  auto constexpr index_normal = GlobalDim - 1;
183 
184  typename ShapeMatricesTypePressure::NodalMatrixType laplace_p =
185  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
186  pressure_size);
187 
188  typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
189  ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
190  pressure_size);
191 
192  typename ShapeMatricesTypeDisplacement::template MatrixType<
194  Kgp = ShapeMatricesTypeDisplacement::template MatrixType<
195  displacement_size, pressure_size>::Zero(displacement_size,
196  pressure_size);
197 
198  using GlobalDimMatrix = Eigen::Matrix<double, GlobalDim, GlobalDim>;
199  using GlobalDimVector = Eigen::Matrix<double, GlobalDim, 1>;
200 
201  auto const& gravity_vec = _process_data.specific_body_force;
202 
204  x_position.setElementID(_element.getID());
205 
206  unsigned const n_integration_points = _ip_data.size();
207  for (unsigned ip = 0; ip < n_integration_points; ip++)
208  {
209  x_position.setIntegrationPoint(ip);
210 
211  auto& ip_data = _ip_data[ip];
212  auto const& ip_w = ip_data.integration_weight;
213  auto const& N_p = ip_data.N_p;
214  auto const& dNdx_p = ip_data.dNdx_p;
215  auto const& H_g = ip_data.H_u;
216  auto const& identity2 =
218 
219  auto& mat = ip_data.fracture_material;
220  auto& effective_stress = ip_data.sigma_eff;
221  auto const& effective_stress_prev = ip_data.sigma_eff_prev;
222  auto& w = ip_data.w;
223  auto const& w_prev = ip_data.w_prev;
224  auto& C = ip_data.C;
225  auto& state = *ip_data.material_state_variables;
226  auto& b_m = ip_data.aperture;
227 
228  double const S = frac_prop.specific_storage(t, x_position)[0];
229  double const mu = _process_data.fluid_viscosity(t, x_position)[0];
230  auto const alpha = frac_prop.biot_coefficient(t, x_position)[0];
231  auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
232 
233  // displacement jumps in local coordinates
234  w.noalias() = R * H_g * g;
235 
236  // aperture
237  b_m = ip_data.aperture0 + w[index_normal];
238  if (b_m < 0.0)
239  {
240  DBUG(
241  "Element %d, gp %d: Fracture aperture is %g, but it must be "
242  "non-negative. Setting it to zero.",
243  _element.getID(), ip, b_m);
244  b_m = 0;
245  }
246 
247  auto const initial_effective_stress =
248  _process_data.initial_fracture_effective_stress(0, x_position);
249 
250  Eigen::Map<typename HMatricesType::ForceVectorType const> const stress0(
251  initial_effective_stress.data(), initial_effective_stress.size());
252 
253  // local C, local stress
254  mat.computeConstitutiveRelation(
255  t, x_position, ip_data.aperture0, stress0, w_prev, w,
256  effective_stress_prev, effective_stress, C, state);
257 
258  auto& permeability = ip_data.permeability;
259  permeability = frac_prop.permeability_model->permeability(
260  ip_data.permeability_state.get(), ip_data.aperture0, b_m);
261 
262  GlobalDimMatrix const k =
263  createRotatedTensor<GlobalDim>(R, permeability);
264 
265  // derivative of permeability respect to aperture
266  double const local_dk_db =
267  frac_prop.permeability_model->dpermeability_daperture(
268  ip_data.permeability_state.get(), ip_data.aperture0, b_m);
269  GlobalDimMatrix const dk_db =
270  createRotatedTensor<GlobalDim>(R, local_dk_db);
271 
272  //
273  // displacement equation, displacement jump part
274  //
275  rhs_g.noalias() -=
276  H_g.transpose() * R.transpose() * effective_stress * ip_w;
277  J_gg.noalias() += H_g.transpose() * R.transpose() * C * R * H_g * ip_w;
278 
279  //
280  // displacement equation, pressure part
281  //
282  Kgp.noalias() +=
283  H_g.transpose() * R.transpose() * alpha * identity2 * N_p * ip_w;
284 
285  //
286  // pressure equation, pressure part.
287  //
288  storage_p.noalias() += N_p.transpose() * b_m * S * N_p * ip_w;
289  laplace_p.noalias() +=
290  dNdx_p.transpose() * b_m * k / mu * dNdx_p * ip_w;
291  rhs_p.noalias() +=
292  dNdx_p.transpose() * b_m * k / mu * rho_fr * gravity_vec * ip_w;
293 
294  //
295  // pressure equation, displacement jump part.
296  //
297  GlobalDimVector const grad_head_over_mu =
298  (dNdx_p * p + rho_fr * gravity_vec) / mu;
299  Eigen::Matrix<double, 1, displacement_size> const mT_R_Hg =
300  identity2.transpose() * R * H_g;
301  // velocity in global coordinates
302  ip_data.darcy_velocity.head(GlobalDim).noalias() =
303  -R.transpose() * k * grad_head_over_mu;
304  J_pg.noalias() += N_p.transpose() * S * N_p * p_dot * mT_R_Hg * ip_w;
305  J_pg.noalias() +=
306  dNdx_p.transpose() * k * grad_head_over_mu * mT_R_Hg * ip_w;
307  J_pg.noalias() += dNdx_p.transpose() * b_m * dk_db * grad_head_over_mu *
308  mT_R_Hg * ip_w;
309  }
310 
311  // displacement equation, pressure part
312  J_gp.noalias() -= Kgp;
313 
314  // pressure equation, pressure part.
315  J_pp.noalias() += laplace_p + storage_p / dt;
316 
317  // pressure equation, displacement jump part.
318  J_pg.noalias() += Kgp.transpose() / dt;
319 
320  // pressure equation
321  rhs_p.noalias() -=
322  laplace_p * p + storage_p * p_dot + Kgp.transpose() * g_dot;
323 
324  // displacement equation
325  rhs_g.noalias() -= -Kgp * p;
326 }
327 
328 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
329  typename IntegrationMethod, int GlobalDim>
330 void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
331  ShapeFunctionPressure,
332  IntegrationMethod, GlobalDim>::
334  Eigen::VectorXd const& local_x)
335 {
336  auto const nodal_g = local_x.segment(displacement_index, displacement_size);
337 
338  auto const& frac_prop = *_process_data.fracture_property;
339  auto const& R = frac_prop.R;
340  // the index of a normal (normal to a fracture plane) component
341  // in a displacement vector
342  auto constexpr index_normal = GlobalDim - 1;
343 
345  x_position.setElementID(_element.getID());
346 
347  unsigned const n_integration_points = _ip_data.size();
348  for (unsigned ip = 0; ip < n_integration_points; ip++)
349  {
350  x_position.setIntegrationPoint(ip);
351 
352  auto& ip_data = _ip_data[ip];
353  auto const& H_g = ip_data.H_u;
354  auto& mat = ip_data.fracture_material;
355  auto& effective_stress = ip_data.sigma_eff;
356  auto const& effective_stress_prev = ip_data.sigma_eff_prev;
357  auto& w = ip_data.w;
358  auto const& w_prev = ip_data.w_prev;
359  auto& C = ip_data.C;
360  auto& state = *ip_data.material_state_variables;
361  auto& b_m = ip_data.aperture;
362 
363  // displacement jumps in local coordinates
364  w.noalias() = R * H_g * nodal_g;
365 
366  // aperture
367  b_m = ip_data.aperture0 + w[index_normal];
368  if (b_m < 0.0)
369  {
370  DBUG(
371  "Element %d, gp %d: Fracture aperture is %g, but it is "
372  "expected to be non-negative.",
373  _element.getID(), ip, b_m);
374  }
375 
376  auto const initial_effective_stress =
377  _process_data.initial_fracture_effective_stress(0, x_position);
378 
379  Eigen::Map<typename HMatricesType::ForceVectorType const> const stress0(
380  initial_effective_stress.data(), initial_effective_stress.size());
381 
382  // local C, local stress
383  mat.computeConstitutiveRelation(
384  t, x_position, ip_data.aperture0, stress0, w_prev, w,
385  effective_stress_prev, effective_stress, C, state);
386  }
387 
388  double ele_b = 0;
389  double ele_k = 0;
390  typename HMatricesType::ForceVectorType ele_sigma_eff =
391  HMatricesType::ForceVectorType::Zero(GlobalDim);
392  typename HMatricesType::ForceVectorType ele_w =
393  HMatricesType::ForceVectorType::Zero(GlobalDim);
394  double ele_Fs = -std::numeric_limits<double>::max();
395  Eigen::Vector3d ele_velocity = Eigen::Vector3d::Zero();
396  for (auto const& ip : _ip_data)
397  {
398  ele_b += ip.aperture;
399  ele_k += ip.permeability;
400  ele_w += ip.w;
401  ele_sigma_eff += ip.sigma_eff;
402  ele_Fs = std::max(
403  ele_Fs, ip.material_state_variables->getShearYieldFunctionValue());
404  ele_velocity += ip.darcy_velocity;
405  }
406  ele_b /= static_cast<double>(n_integration_points);
407  ele_k /= static_cast<double>(n_integration_points);
408  ele_w /= static_cast<double>(n_integration_points);
409  ele_sigma_eff /= static_cast<double>(n_integration_points);
410  ele_velocity /= static_cast<double>(n_integration_points);
411  auto const element_id = _element.getID();
412  (*_process_data.mesh_prop_b)[element_id] = ele_b;
413  (*_process_data.mesh_prop_k_f)[element_id] = ele_k;
414  (*_process_data.mesh_prop_w_n)[element_id] = ele_w[index_normal];
415  (*_process_data.mesh_prop_w_s)[element_id] = ele_w[0];
416  (*_process_data.mesh_prop_fracture_stress_normal)[element_id] =
417  ele_sigma_eff[index_normal];
418  (*_process_data.mesh_prop_fracture_stress_shear)[element_id] =
419  ele_sigma_eff[0];
420  (*_process_data.mesh_prop_fracture_shear_failure)[element_id] = ele_Fs;
421 
422  if (GlobalDim == 3)
423  {
424  (*_process_data.mesh_prop_w_s2)[element_id] = ele_w[1];
425  (*_process_data.mesh_prop_fracture_stress_shear2)[element_id] =
426  ele_sigma_eff[1];
427  }
428 
429  for (unsigned i = 0; i < 3; i++)
430  {
431  (*_process_data.mesh_prop_velocity)[element_id * 3 + i] =
432  ele_velocity[i];
433  }
434 }
435 
436 } // namespace HydroMechanics
437 } // namespace LIE
438 } // namespace ProcessLib
void assembleWithJacobianConcrete(double const t, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, Eigen::VectorXd &local_b, Eigen::MatrixXd &local_J) override
void setElementID(std::size_t element_id)
std::vector< IntegrationPointDataType, Eigen::aligned_allocator< IntegrationPointDataType > > _ip_data
ShapeMatrixPolicyType< ShapeFunctionPressure, GlobalDim > ShapeMatricesTypePressure
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool is_axially_symmetric, IntegrationMethod const &integration_method)
virtual unsigned getDimension() const =0
Get dimension of the mesh element.
void setIntegrationPoint(unsigned integration_point)
void computeSecondaryVariableConcreteWithVector(const double t, Eigen::VectorXd const &local_x) override
void assembleBlockMatricesWithJacobian(double const t, Eigen::Ref< const Eigen::VectorXd > const &p, Eigen::Ref< const Eigen::VectorXd > const &p_dot, Eigen::Ref< const Eigen::VectorXd > const &g, Eigen::Ref< const Eigen::VectorXd > const &g_dot, Eigen::Ref< Eigen::VectorXd > rhs_p, Eigen::Ref< Eigen::VectorXd > rhs_g, Eigen::Ref< Eigen::MatrixXd > J_pp, Eigen::Ref< Eigen::MatrixXd > J_pg, Eigen::Ref< Eigen::MatrixXd > J_gg, Eigen::Ref< Eigen::MatrixXd > J_gp)
Eigen::Matrix< double, 3u, 3u, Eigen::RowMajor > RotationMatrix
VectorType< DisplacementDim > ForceVectorType
Definition: HMatrixUtils.h:44
VectorType< ShapeFunction::NPOINTS > NodalVectorType
Eigen::Matrix< double, GlobalDim, GlobalDim > createRotatedTensor(RotationMatrix const &R, double const value)
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
void computeHMatrix(N_Type const &N, HMatrixType &H)
Fills a H-matrix based on given shape function.
Definition: HMatrixUtils.h:52
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:90
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
ShapeMatrixPolicyType< ShapeFunctionDisplacement, GlobalDim > ShapeMatricesTypeDisplacement