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