OGS 6.1.0-1699-ge946d4c5f
SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
Go to the documentation of this file.
1 
10 #pragma once
11 
13 
14 #include <valarray>
15 #include <vector>
16 
17 #include <Eigen/Eigen>
18 
20 
21 #include "MathLib/KelvinVector.h"
23 #include "MathLib/Point3d.h"
24 
25 #include "MeshLib/Node.h"
26 
29 
33 
36 
38 #include "SecondaryData.h"
40 
41 namespace ProcessLib
42 {
43 namespace LIE
44 {
45 namespace SmallDeformation
46 {
47 template <typename ShapeFunction,
48  typename IntegrationMethod,
49  int DisplacementDim>
50 SmallDeformationLocalAssemblerMatrixNearFracture<ShapeFunction,
51  IntegrationMethod,
52  DisplacementDim>::
53  SmallDeformationLocalAssemblerMatrixNearFracture(
54  MeshLib::Element const& e,
55  std::size_t const n_variables,
56  std::size_t const /*local_matrix_size*/,
57  std::vector<unsigned> const& dofIndex_to_localIndex,
58  bool const is_axially_symmetric,
59  unsigned const integration_order,
62  n_variables * ShapeFunction::NPOINTS * DisplacementDim,
63  dofIndex_to_localIndex),
64  _process_data(process_data),
65  _integration_method(integration_order),
66  _element(e),
67  _is_axially_symmetric(is_axially_symmetric)
68 {
69  std::vector<
71  Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
72  shape_matrices = initShapeMatrices<ShapeFunction,
74  IntegrationMethod,
75  DisplacementDim>(
76  e, is_axially_symmetric, _integration_method);
77 
78  unsigned const n_integration_points =
79  _integration_method.getNumberOfPoints();
80 
81  _ip_data.reserve(n_integration_points);
82  _secondary_data.N.resize(n_integration_points);
83 
85  _process_data.solid_materials, _process_data.material_ids, e.getID());
86 
87  for (unsigned ip = 0; ip < n_integration_points; ip++)
88  {
89  _ip_data.emplace_back(solid_material);
90  auto& ip_data = _ip_data[ip];
91  auto const& sm = shape_matrices[ip];
92  ip_data.N = sm.N;
93  ip_data.dNdx = sm.dNdx;
94  ip_data.integration_weight =
95  _integration_method.getWeightedPoint(ip).getWeight() *
96  sm.integralMeasure * sm.detJ;
97 
98  // Initialize current time step values
99  static const int kelvin_vector_size =
101  DisplacementDim>::value;
102  ip_data._sigma.setZero(kelvin_vector_size);
103  ip_data._eps.setZero(kelvin_vector_size);
104 
105  // Previous time step values are not initialized and are set later.
106  ip_data._sigma_prev.resize(kelvin_vector_size);
107  ip_data._eps_prev.resize(kelvin_vector_size);
108 
109  ip_data._C.resize(kelvin_vector_size, kelvin_vector_size);
110 
111  _secondary_data.N[ip] = sm.N;
112  }
113 
114  for (auto fid : process_data._vec_ele_connected_fractureIDs[e.getID()])
115  {
116  _fracID_to_local.insert({fid, _fracture_props.size()});
117  _fracture_props.push_back(&_process_data.fracture_properties[fid]);
118  }
119 
120  for (auto jid : process_data._vec_ele_connected_junctionIDs[e.getID()])
121  {
122  _junction_props.push_back(&_process_data.junction_properties[jid]);
123  }
124 }
125 
126 template <typename ShapeFunction,
127  typename IntegrationMethod,
128  int DisplacementDim>
130  ShapeFunction,
131  IntegrationMethod,
132  DisplacementDim>::assembleWithJacobian(double const t,
133  Eigen::VectorXd const& local_u,
134  Eigen::VectorXd& local_b,
135  Eigen::MatrixXd& local_J)
136 {
137  assert(_element.getDimension() == DisplacementDim);
138 
139  auto const N_DOF_PER_VAR = ShapeFunction::NPOINTS * DisplacementDim;
140  auto const n_fractures = _fracture_props.size();
141  auto const n_junctions = _junction_props.size();
142  auto const n_enrich_var = n_fractures + n_junctions;
143 
144  using BlockVectorType =
145  typename Eigen::VectorXd::FixedSegmentReturnType<N_DOF_PER_VAR>::Type;
146  using BlockMatrixType =
147  Eigen::Block<Eigen::MatrixXd, N_DOF_PER_VAR, N_DOF_PER_VAR>;
148 
149  //--------------------------------------------------------------------------------------
150  // prepare sub vectors, matrices for regular displacement (u) and
151  // displacement jumps (g)
152  //
153  // example with two fractures with one intersection:
154  // |b(u)|
155  // b = |b(g1)|
156  // |b(g2)|
157  // |b(j1)|
158  //
159  // |J(u,u) J(u,g1) J(u,g2) J(u,j1) |
160  // J = |J(g1,u) J(g1,g1) J(g1,g2) J(g1,j1)|
161  // |J(g2,u) J(g2,g1) J(g2,g2) J(g2,j1)|
162  // |J(j1,u) J(j1,g1) J(j1,g2) J(j1,j1)|
163  //--------------------------------------------------------------------------------------
164  auto local_b_u = local_b.segment<N_DOF_PER_VAR>(0);
165  std::vector<BlockVectorType> vec_local_b_g;
166  for (unsigned i = 0; i < n_enrich_var; i++)
167  {
168  vec_local_b_g.push_back(
169  local_b.segment<N_DOF_PER_VAR>(N_DOF_PER_VAR * (i + 1)));
170  }
171 
172  auto local_J_uu = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(0, 0);
173  std::vector<BlockMatrixType> vec_local_J_ug;
174  std::vector<BlockMatrixType> vec_local_J_gu;
175  std::vector<std::vector<BlockMatrixType>> vec_local_J_gg(n_enrich_var);
176  for (unsigned i = 0; i < n_enrich_var; i++)
177  {
178  auto sub_ug = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
179  0, N_DOF_PER_VAR * (i + 1));
180  vec_local_J_ug.push_back(sub_ug);
181 
182  auto sub_gu = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
183  N_DOF_PER_VAR * (i + 1), 0);
184  vec_local_J_gu.push_back(sub_gu);
185 
186  for (unsigned j = 0; j < n_enrich_var; j++)
187  {
188  auto sub_gg = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
189  N_DOF_PER_VAR * (i + 1), N_DOF_PER_VAR * (j + 1));
190  vec_local_J_gg[i].push_back(sub_gg);
191  }
192  }
193 
194  auto const nodal_u = local_u.segment<N_DOF_PER_VAR>(0);
195  std::vector<BlockVectorType> vec_nodal_g;
196  for (unsigned i = 0; i < n_enrich_var; i++)
197  {
198  auto sub = const_cast<Eigen::VectorXd&>(local_u).segment<N_DOF_PER_VAR>(
199  N_DOF_PER_VAR * (i + 1));
200  vec_nodal_g.push_back(sub);
201  }
202 
203  //------------------------------------------------
204  // integration
205  //------------------------------------------------
206  unsigned const n_integration_points =
207  _integration_method.getNumberOfPoints();
208 
209  SpatialPosition x_position;
210  x_position.setElementID(_element.getID());
211 
212  for (unsigned ip = 0; ip < n_integration_points; ip++)
213  {
214  x_position.setIntegrationPoint(ip);
215 
216  auto& ip_data = _ip_data[ip];
217  auto const& w = _ip_data[ip].integration_weight;
218 
219  auto const& N = ip_data.N;
220  auto const& dNdx = ip_data.dNdx;
221 
222  // levelset functions
223  Eigen::Vector3d const ip_physical_coords(
224  computePhysicalCoordinates(_element, N).getCoords());
225  std::vector<double> const levelsets(
227  _fracID_to_local, ip_physical_coords));
228 
229  // u = u^hat + sum_i(enrich^br_i(x) * [u]_i) + sum_i(enrich^junc_i(x) *
230  // [u]_i)
231  NodalDisplacementVectorType nodal_total_u = nodal_u;
232  for (unsigned i = 0; i < n_enrich_var; i++)
233  {
234  nodal_total_u += levelsets[i] * vec_nodal_g[i];
235  }
236 
237  auto const x_coord =
238  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(_element,
239  N);
240  auto const B =
241  LinearBMatrix::computeBMatrix<DisplacementDim,
242  ShapeFunction::NPOINTS,
243  typename BMatricesType::BMatrixType>(
244  dNdx, N, x_coord, _is_axially_symmetric);
245 
246  // strain, stress
247  auto const& eps_prev = ip_data._eps_prev;
248  auto const& sigma_prev = ip_data._sigma_prev;
249 
250  auto& eps = ip_data._eps;
251  auto& sigma = ip_data._sigma;
252  auto& state = ip_data._material_state_variables;
253 
254  eps.noalias() = B * nodal_total_u;
255 
256  auto&& solution = _ip_data[ip]._solid_material.integrateStress(
257  t, x_position, _process_data.dt, eps_prev, eps, sigma_prev, *state,
258  _process_data._reference_temperature);
259 
260  if (!solution)
261  {
262  OGS_FATAL("Computation of local constitutive relation failed.");
263  }
264 
266  std::tie(sigma, state, C) = std::move(*solution);
267 
268  // r_u = B^T * Sigma = B^T * C * B * (u+phi*[u])
269  // r_[u] = (phi*B)^T * Sigma = (phi*B)^T * C * B * (u+phi*[u])
270  local_b_u.noalias() -= B.transpose() * sigma * w;
271  for (unsigned i = 0; i < n_enrich_var; i++)
272  {
273  vec_local_b_g[i].noalias() -=
274  levelsets[i] * B.transpose() * sigma * w;
275  }
276 
277  // J_uu += B^T * C * B
278  local_J_uu.noalias() += B.transpose() * C * B * w;
279 
280  for (unsigned i = 0; i < n_enrich_var; i++)
281  {
282  // J_u[u] += B^T * C * (levelset * B)
283  vec_local_J_ug[i].noalias() +=
284  B.transpose() * C * (levelsets[i] * B) * w;
285 
286  // J_[u]u += (levelset * B)^T * C * B
287  vec_local_J_gu[i].noalias() +=
288  (levelsets[i] * B.transpose()) * C * B * w;
289 
290  for (unsigned j = 0; j < n_enrich_var; j++)
291  {
292  // J_[u][u] += (levelset * B)^T * C * (levelset * B)
293  vec_local_J_gg[i][j].noalias() +=
294  (levelsets[i] * B.transpose()) * C * (levelsets[j] * B) * w;
295  }
296  }
297  }
298 }
299 
300 template <typename ShapeFunction,
301  typename IntegrationMethod,
302  int DisplacementDim>
304  IntegrationMethod,
305  DisplacementDim>::
307  double const /*t*/, Eigen::VectorXd const& /*local_x*/)
308 {
309  // Compute average value per element
310  const int n = DisplacementDim == 2 ? 4 : 6;
311  Eigen::VectorXd ele_stress = Eigen::VectorXd::Zero(n);
312  Eigen::VectorXd ele_strain = Eigen::VectorXd::Zero(n);
313 
314  unsigned const n_integration_points =
315  _integration_method.getNumberOfPoints();
316  for (unsigned ip = 0; ip < n_integration_points; ip++)
317  {
318  auto& ip_data = _ip_data[ip];
319 
320  ele_stress += ip_data._sigma;
321  ele_strain += ip_data._eps;
322  }
323  ele_stress /= n_integration_points;
324  ele_strain /= n_integration_points;
325 
326  (*_process_data._mesh_prop_stress_xx)[_element.getID()] = ele_stress[0];
327  (*_process_data._mesh_prop_stress_yy)[_element.getID()] = ele_stress[1];
328  (*_process_data._mesh_prop_stress_zz)[_element.getID()] = ele_stress[2];
329  (*_process_data._mesh_prop_stress_xy)[_element.getID()] = ele_stress[3];
330  if (DisplacementDim == 3)
331  {
332  (*_process_data._mesh_prop_stress_yz)[_element.getID()] = ele_stress[4];
333  (*_process_data._mesh_prop_stress_xz)[_element.getID()] = ele_stress[5];
334  }
335 
336  (*_process_data._mesh_prop_strain_xx)[_element.getID()] = ele_strain[0];
337  (*_process_data._mesh_prop_strain_yy)[_element.getID()] = ele_strain[1];
338  (*_process_data._mesh_prop_strain_zz)[_element.getID()] = ele_strain[2];
339  (*_process_data._mesh_prop_strain_xy)[_element.getID()] = ele_strain[3];
340  if (DisplacementDim == 3)
341  {
342  (*_process_data._mesh_prop_strain_yz)[_element.getID()] = ele_strain[4];
343  (*_process_data._mesh_prop_strain_xz)[_element.getID()] = ele_strain[5];
344  }
345 }
346 
347 } // namespace SmallDeformation
348 } // namespace LIE
349 } // namespace ProcessLib
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
Definition: SecondaryData.h:29
MathLib::Point3d computePhysicalCoordinates(MeshLib::Element const &e, Eigen::MatrixBase< Derived > const &shape)
Definition: Utils.h:46
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool is_axially_symmetric, IntegrationMethod const &integration_method)
Definition of the Node class.
MechanicsBase< DisplacementDim > & selectSolidConstitutiveRelation(std::map< int, std::unique_ptr< MechanicsBase< DisplacementDim >>> const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
Definition of the Point3d class.
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, KelvinVectorDimensions< DisplacementDim >::value, Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:59
void setElementID(std::size_t element_id)
virtual unsigned getDimension() const =0
Get dimension of the mesh element.
std::vector< double > uGlobalEnrichments(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)
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition: LinearBMatrix.h:41
void assembleWithJacobian(double const t, Eigen::VectorXd const &local_u, Eigen::VectorXd &local_b, Eigen::MatrixXd &local_J) override
std::vector< IntegrationPointDataMatrix< ShapeMatricesType, BMatricesType, DisplacementDim >, Eigen::aligned_allocator< IntegrationPointDataMatrix< ShapeMatricesType, BMatricesType, DisplacementDim > > > _ip_data
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:90
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:54