OGS 6.2.1-97-g73d1aeda3
SmallDeformationLocalAssemblerFracture-impl.h
Go to the documentation of this file.
1 
10 #pragma once
11 
13 
14 #include <Eigen/Eigen>
15 
17 
19 
21 
22 namespace ProcessLib
23 {
24 namespace LIE
25 {
26 namespace SmallDeformation
27 {
28 template <typename ShapeFunction, typename IntegrationMethod,
29  int DisplacementDim>
30 SmallDeformationLocalAssemblerFracture<ShapeFunction, IntegrationMethod,
31  DisplacementDim>::
32  SmallDeformationLocalAssemblerFracture(
33  MeshLib::Element const& e,
34  std::size_t const n_variables,
35  std::size_t const /*local_matrix_size*/,
36  std::vector<unsigned> const& dofIndex_to_localIndex,
37  bool const is_axially_symmetric,
38  unsigned const integration_order,
41  n_variables * ShapeFunction::NPOINTS * DisplacementDim,
42  dofIndex_to_localIndex),
43  _process_data(process_data),
44  _integration_method(integration_order),
45  _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
46  IntegrationMethod, DisplacementDim>(
47  e, is_axially_symmetric, _integration_method)),
48  _element(e)
49 {
50  assert(_element.getDimension() == DisplacementDim - 1);
51 
52  unsigned const n_integration_points =
53  _integration_method.getNumberOfPoints();
54 
55  _ip_data.reserve(n_integration_points);
56  _secondary_data.N.resize(n_integration_points);
57 
58  auto mat_id = (*_process_data._mesh_prop_materialIDs)[e.getID()];
59  auto frac_id = _process_data._map_materialID_to_fractureID[mat_id];
60  _fracture_property = &_process_data.fracture_properties[frac_id];
61  for (auto fid : process_data._vec_ele_connected_fractureIDs[e.getID()])
62  {
63  _fracID_to_local.insert({fid, _fracture_props.size()});
64  _fracture_props.push_back(&_process_data.fracture_properties[fid]);
65  }
66 
67  for (auto jid : process_data._vec_ele_connected_junctionIDs[e.getID()])
68  {
69  _junction_props.push_back(&_process_data.junction_properties[jid]);
70  }
71 
73  x_position.setElementID(_element.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 = _shape_matrices[ip];
80  auto& ip_data = _ip_data[ip];
81  ip_data.integration_weight =
82  _integration_method.getWeightedPoint(ip).getWeight() *
83  sm.integralMeasure * sm.detJ;
84  ip_data.h_matrices.setZero(DisplacementDim,
85  ShapeFunction::NPOINTS * DisplacementDim);
86 
87  computeHMatrix<DisplacementDim, ShapeFunction::NPOINTS,
89  HMatrixType>(sm.N, ip_data.h_matrices);
90 
91  // Initialize current time step values
92  ip_data.w.setZero(DisplacementDim);
93  ip_data.sigma.setZero(DisplacementDim);
94 
95  // Previous time step values are not initialized and are set later.
96  ip_data.sigma_prev.resize(DisplacementDim);
97  ip_data.w_prev.resize(DisplacementDim);
98 
99  ip_data.C.resize(DisplacementDim, DisplacementDim);
100 
101  ip_data.aperture0 = _fracture_property->aperture0(0, x_position)[0];
102  ip_data.aperture_prev = ip_data.aperture0;
103 
104  _secondary_data.N[ip] = sm.N;
105  }
106 }
107 
108 template <typename ShapeFunction, typename IntegrationMethod,
109  int DisplacementDim>
111  ShapeFunction, IntegrationMethod,
112  DisplacementDim>::assembleWithJacobian(double const t,
113  Eigen::VectorXd const& local_u,
114  Eigen::VectorXd& local_b,
115  Eigen::MatrixXd& local_J)
116 {
117  auto const N_DOF_PER_VAR = ShapeFunction::NPOINTS * DisplacementDim;
118  auto const n_fractures = _fracture_props.size();
119  auto const n_junctions = _junction_props.size();
120  auto const n_enrich_var = n_fractures + n_junctions;
121 
122  //--------------------------------------------------------------------------------------
123  // prepare sub vectors, matrices for regular displacement (u) and
124  // displacement jumps (g)
125  //
126  // example with two fractures with one intersection:
127  // b = |b(g1)|
128  // |b(g2)|
129  //
130  // J = |J(g1,g1) J(g1,g2)|
131  // |J(g2,g1) J(g2,g2)|
132  //--------------------------------------------------------------------------------------
133 
134  using BlockVectorType =
135  typename Eigen::VectorXd::FixedSegmentReturnType<N_DOF_PER_VAR>::Type;
136  using BlockMatrixType =
137  Eigen::Block<Eigen::MatrixXd, N_DOF_PER_VAR, N_DOF_PER_VAR>;
138 
139  std::vector<BlockVectorType> vec_local_b_g;
140  for (unsigned i = 0; i < n_enrich_var; i++)
141  {
142  vec_local_b_g.push_back(
143  local_b.segment<N_DOF_PER_VAR>(N_DOF_PER_VAR * i));
144  }
145  std::vector<std::vector<BlockMatrixType>> vec_local_J_gg(n_enrich_var);
146  for (unsigned i = 0; i < n_enrich_var; i++)
147  {
148  for (unsigned j = 0; j < n_enrich_var; j++)
149  {
150  auto sub_gg = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
151  N_DOF_PER_VAR * i, N_DOF_PER_VAR * j);
152  vec_local_J_gg[i].push_back(sub_gg);
153  }
154  }
155 
156  std::vector<Eigen::VectorXd> vec_nodal_g;
157  for (unsigned i = 0; i < n_enrich_var; i++)
158  {
159  auto sub = const_cast<Eigen::VectorXd&>(local_u).segment<N_DOF_PER_VAR>(
160  N_DOF_PER_VAR * i);
161  vec_nodal_g.push_back(sub);
162  }
163 
164  //------------------------------------------------
165  // integration
166  //------------------------------------------------
167  // the index of a normal (normal to a fracture plane) component
168  // in a displacement vector
169  int const index_normal = DisplacementDim - 1;
170  auto const& R = _fracture_property->R;
171 
172  unsigned const n_integration_points =
173  _integration_method.getNumberOfPoints();
174 
176  x_position.setElementID(_element.getID());
177 
178  for (unsigned ip = 0; ip < n_integration_points; ip++)
179  {
180  x_position.setIntegrationPoint(ip);
181 
182  auto& ip_data = _ip_data[ip];
183  auto const& integration_weight = ip_data.integration_weight;
184  auto const& H = ip_data.h_matrices;
185  auto& mat = ip_data.fracture_material;
186  auto& sigma = ip_data.sigma;
187  auto const& sigma_prev = ip_data.sigma_prev;
188  auto& w = ip_data.w;
189  auto const& w_prev = ip_data.w_prev;
190  auto& C = ip_data.C;
191  auto& state = *ip_data.material_state_variables;
192  auto& N = _secondary_data.N[ip];
193 
194  Eigen::Vector3d const ip_physical_coords(
195  computePhysicalCoordinates(_element, N).getCoords());
196  std::vector<double> const levelsets(duGlobalEnrichments(
198  _fracID_to_local, ip_physical_coords));
199 
200  // du = du^hat + sum_i(enrich^br_i(x) * [u]_i) + sum_i(enrich^junc_i(x)
201  // * [u]_i)
202  Eigen::VectorXd nodal_gap(N_DOF_PER_VAR);
203  nodal_gap.setZero();
204  for (unsigned i = 0; i < n_enrich_var; i++)
205  {
206  nodal_gap += levelsets[i] * vec_nodal_g[i];
207  }
208 
209  // displacement jumps
210  w.noalias() = R * H * nodal_gap;
211 
212  // total aperture
213  ip_data.aperture = ip_data.aperture0 + w[index_normal];
214 
215  // local C, local stress
216  mat.computeConstitutiveRelation(
217  t, x_position, ip_data.aperture0,
218  Eigen::Matrix<double, DisplacementDim, 1>::Zero(), // TODO (naumov)
219  // Replace with
220  // initial
221  // stress values
222  w_prev, w, sigma_prev, sigma, C, state);
223 
224  // r_[u] += H^T*Stress
225  for (unsigned i = 0; i < n_enrich_var; i++)
226  {
227  vec_local_b_g[i].noalias() -= levelsets[i] * H.transpose() *
228  R.transpose() * sigma *
229  integration_weight;
230  }
231 
232  // J_[u][u] += H^T*C*H
233  for (unsigned i = 0; i < n_enrich_var; i++)
234  {
235  for (unsigned j = 0; j < n_enrich_var; j++)
236  {
237  // J_[u][u] += (levelset * B)^T * C * (levelset * B)
238  vec_local_J_gg[i][j].noalias() +=
239  (levelsets[i] * H.transpose() * R.transpose()) * C *
240  (levelsets[j] * R * H) * integration_weight;
241  }
242  }
243  }
244 }
245 
246 template <typename ShapeFunction, typename IntegrationMethod,
247  int DisplacementDim>
248 void SmallDeformationLocalAssemblerFracture<ShapeFunction, IntegrationMethod,
249  DisplacementDim>::
251  Eigen::VectorXd const& local_u)
252 {
253  auto const N_DOF_PER_VAR = ShapeFunction::NPOINTS * DisplacementDim;
254  auto const n_fractures = _fracture_props.size();
255  auto const n_junctions = _junction_props.size();
256  auto const n_enrich_var = n_fractures + n_junctions;
257 
258  auto const& R = _fracture_property->R;
259 
260  // the index of a normal (normal to a fracture plane) component
261  // in a displacement vector
262  int const index_normal = DisplacementDim - 1;
263 
264  unsigned const n_integration_points =
265  _integration_method.getNumberOfPoints();
266 
268  x_position.setElementID(_element.getID());
269 
270  std::vector<Eigen::VectorXd> vec_nodal_g;
271  for (unsigned i = 0; i < n_enrich_var; i++)
272  {
273  auto sub = const_cast<Eigen::VectorXd&>(local_u).segment<N_DOF_PER_VAR>(
274  N_DOF_PER_VAR * i);
275  vec_nodal_g.push_back(sub);
276  }
277 
278  for (unsigned ip = 0; ip < n_integration_points; ip++)
279  {
280  x_position.setIntegrationPoint(ip);
281 
282  auto& ip_data = _ip_data[ip];
283  auto const& H = ip_data.h_matrices;
284  auto& mat = ip_data.fracture_material;
285  auto& sigma = ip_data.sigma;
286  auto const& sigma_prev = ip_data.sigma_prev;
287  auto& w = ip_data.w;
288  auto const& w_prev = ip_data.w_prev;
289  auto& C = ip_data.C;
290  auto& state = *ip_data.material_state_variables;
291  auto& b_m = ip_data.aperture;
292  auto& N = _secondary_data.N[ip];
293 
294  Eigen::Vector3d const ip_physical_coords(
295  computePhysicalCoordinates(_element, N).getCoords());
296  std::vector<double> const levelsets(duGlobalEnrichments(
298  _fracID_to_local, ip_physical_coords));
299 
300  // du = du^hat + sum_i(enrich^br_i(x) * [u]_i) + sum_i(enrich^junc_i(x)
301  // * [u]_i)
302  Eigen::VectorXd nodal_gap(N_DOF_PER_VAR);
303  nodal_gap.setZero();
304  for (unsigned i = 0; i < n_enrich_var; i++)
305  {
306  nodal_gap += levelsets[i] * vec_nodal_g[i];
307  }
308 
309  // displacement jumps in local coordinates
310  w.noalias() = R * H * nodal_gap;
311 
312  // aperture
313  b_m = ip_data.aperture0 + w[index_normal];
314  if (b_m < 0.0)
315  {
316  OGS_FATAL(
317  "Element %d, gp %d: Fracture aperture is %g, but it must be "
318  "non-negative.",
319  _element.getID(), ip, b_m);
320  }
321 
322  // local C, local stress
323  mat.computeConstitutiveRelation(
324  t, x_position, ip_data.aperture0,
325  Eigen::Matrix<double, DisplacementDim, 1>::Zero(), // TODO (naumov)
326  // Replace with
327  // initial
328  // stress values
329  w_prev, w, sigma_prev, sigma, C, state);
330  }
331 
332  double ele_b = 0;
333  typename HMatricesType::ForceVectorType ele_sigma =
334  HMatricesType::ForceVectorType::Zero(DisplacementDim);
335  typename HMatricesType::ForceVectorType ele_w =
336  HMatricesType::ForceVectorType::Zero(DisplacementDim);
337 
338  for (unsigned ip = 0; ip < n_integration_points; ip++)
339  {
340  auto& ip_data = _ip_data[ip];
341  ele_b += ip_data.aperture;
342  ele_w += ip_data.w;
343  ele_sigma += ip_data.sigma;
344  }
345  ele_b /= n_integration_points;
346  ele_w /= n_integration_points;
347  ele_sigma /= n_integration_points;
348  (*_process_data._mesh_prop_b)[_element.getID()] = ele_b;
349  (*_process_data._mesh_prop_w_n)[_element.getID()] = ele_w[index_normal];
350  (*_process_data._mesh_prop_w_s)[_element.getID()] = ele_w[0];
351  (*_process_data._mesh_prop_fracture_stress_normal)[_element.getID()] =
352  ele_sigma[index_normal];
353  (*_process_data._mesh_prop_fracture_stress_shear)[_element.getID()] =
354  ele_sigma[0];
355 }
356 
357 } // namespace SmallDeformation
358 } // namespace LIE
359 } // namespace ProcessLib
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
Definition: SecondaryData.h:29
Eigen::MatrixXd R
Rotation matrix from global to local coordinates.
void setElementID(std::size_t element_id)
std::vector< ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > _shape_matrices
MathLib::Point3d computePhysicalCoordinates(MeshLib::Element const &e, Eigen::MatrixBase< Derived > const &shape)
Definition: Utils.h:46
void assembleWithJacobian(double const t, Eigen::VectorXd const &local_u, Eigen::VectorXd &local_b, Eigen::MatrixXd &local_J) override
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool is_axially_symmetric, IntegrationMethod const &integration_method)
ParameterLib::Parameter< double > const & aperture0
Initial aperture.
virtual unsigned getDimension() const =0
Get dimension of the mesh element.
void setIntegrationPoint(unsigned integration_point)
VectorType< DisplacementDim > ForceVectorType
Definition: HMatrixUtils.h:44
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
std::vector< double > duGlobalEnrichments(std::size_t this_frac_id, std::vector< FractureProperty *> const &frac_props, std::vector< JunctionProperty *> const &junction_props, std::unordered_map< int, int > const &fracID_to_local, Eigen::Vector3d const &x)
void computeHMatrix(N_Type const &N, HMatrixType &H)
Fills a H-matrix based on given shape function.
Definition: HMatrixUtils.h:52
void computeSecondaryVariableConcreteWithVector(const double t, Eigen::VectorXd const &local_u) override
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
std::vector< IntegrationPointDataFracture< HMatricesType, DisplacementDim >, Eigen::aligned_allocator< IntegrationPointDataFracture< HMatricesType, DisplacementDim > > > _ip_data