OGS 6.2.1-97-g73d1aeda3
SmallDeformationNonlocalFEM.h
Go to the documentation of this file.
1 
12 #pragma once
13 
14 #include <algorithm>
15 #include <limits>
16 #include <memory>
17 #include <vector>
18 
26 #include "ParameterLib/Parameter.h"
32 
33 #include "Damage.h"
34 #include "IntegrationPointData.h"
37 
38 namespace ProcessLib
39 {
40 namespace SmallDeformationNonlocal
41 {
44 template <typename ShapeMatrixType>
46 {
47  std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
48 };
49 
50 template <typename ShapeFunction, typename IntegrationMethod,
51  int DisplacementDim>
53  : public SmallDeformationNonlocalLocalAssemblerInterface<DisplacementDim>
54 {
55 public:
56  using ShapeMatricesType =
63 
69 
74 
76  MeshLib::Element const& e,
77  std::size_t const /*local_matrix_size*/,
78  bool const is_axially_symmetric,
79  unsigned const integration_order,
81  : _process_data(process_data),
82  _integration_method(integration_order),
83  _element(e),
84  _is_axially_symmetric(is_axially_symmetric)
85  {
86  unsigned const n_integration_points =
87  _integration_method.getNumberOfPoints();
88 
89  _ip_data.reserve(n_integration_points);
90  _secondary_data.N.resize(n_integration_points);
91 
92  auto const shape_matrices =
93  initShapeMatrices<ShapeFunction, ShapeMatricesType,
94  IntegrationMethod, DisplacementDim>(
95  e, is_axially_symmetric, _integration_method);
96 
97  auto& solid_material =
99  _process_data.solid_materials,
100  _process_data.material_ids,
101  e.getID());
102  auto* ehlers_solid_material = dynamic_cast<
104  &solid_material);
105  if (ehlers_solid_material == nullptr)
106  {
107  OGS_FATAL(
108  "The SmallDeformationNonlocalLocal process supports only "
109  "Ehlers material at the moment. For other materials the "
110  "interface must be extended first.");
111  }
112 
113  for (unsigned ip = 0; ip < n_integration_points; ip++)
114  {
115  _ip_data.emplace_back(*ehlers_solid_material);
116  auto& ip_data = _ip_data[ip];
117  auto const& sm = shape_matrices[ip];
118  _ip_data[ip].integration_weight =
119  _integration_method.getWeightedPoint(ip).getWeight() *
120  sm.integralMeasure * sm.detJ;
121 
122  ip_data.N = sm.N;
123  ip_data.dNdx = sm.dNdx;
124 
125  // Initialize current time step values
126  ip_data.sigma.setZero(MathLib::KelvinVector::KelvinVectorDimensions<
127  DisplacementDim>::value);
129  DisplacementDim>::value);
130 
131  // Previous time step values are not initialized and are set later.
132  ip_data.sigma_prev.resize(
134  DisplacementDim>::value);
135  ip_data.eps_prev.resize(
137  DisplacementDim>::value);
138 
139  _secondary_data.N[ip] = shape_matrices[ip].N;
140 
141  ip_data.coordinates = getSingleIntegrationPointCoordinates(ip);
142  }
143  }
144 
145  std::size_t setIPDataInitialConditions(std::string const& name,
146  double const* values,
147  int const integration_order) override
148  {
149  if (integration_order !=
150  static_cast<int>(_integration_method.getIntegrationOrder()))
151  {
152  OGS_FATAL(
153  "Setting integration point initial conditions; The integration "
154  "order of the local assembler for element %d is different from "
155  "the integration order in the initial condition.",
156  _element.getID());
157  }
158 
159  if (name == "sigma_ip")
160  {
161  return setSigma(values);
162  }
163 
164  if (name == "kappa_d_ip")
165  {
166  return setKappaD(values);
167  }
168 
169  return 0;
170  }
171 
173  std::string const& name, std::vector<double> const& value) override
174  {
175  if (name == "kappa_d_ip")
176  {
177  if (value.size() != 1)
178  {
179  OGS_FATAL(
180  "CellData for kappa_d initial conditions has wrong number "
181  "of components. 1 expected, got %d.",
182  value.size());
183  }
184  setKappaD(value[0]);
185  }
186  }
187 
188  double alpha_0(double const distance2) const
189  {
190  double const internal_length2 = _process_data.internal_length_squared;
191  return (distance2 > internal_length2)
192  ? 0
193  : (1 - distance2 / (internal_length2)) *
194  (1 - distance2 / (internal_length2));
195  }
196 
197  void nonlocal(
198  std::size_t const /*mesh_item_id*/,
199  std::vector<
201  DisplacementDim>>> const& local_assemblers) override
202  {
203  auto const search_element_ids = MeshLib::findElementsWithinRadius(
204  _element, _process_data.internal_length_squared);
205 
206  unsigned const n_integration_points =
207  _integration_method.getNumberOfPoints();
208 
209  std::vector<double> distances; // Cache for ip-ip distances.
210  //
211  // For every integration point in this element collect the neighbouring
212  // integration points falling in given radius (internal length) and
213  // compute the alpha_kl weight.
214  //
215  for (unsigned k = 0; k < n_integration_points; k++)
216  {
217  //
218  // Collect the integration points.
219  //
220 
221  auto const& xyz = _ip_data[k].coordinates;
222 
223  // For all neighbors of element
224  for (auto const search_element_id : search_element_ids)
225  {
226  auto const& la = local_assemblers[search_element_id];
227  la->getIntegrationPointCoordinates(xyz, distances);
228  for (int ip = 0; ip < static_cast<int>(distances.size()); ++ip)
229  {
230  if (distances[ip] >= _process_data.internal_length_squared)
231  {
232  continue;
233  }
234  // save into current ip_k
235  _ip_data[k].non_local_assemblers.push_back(
236  {la->getIPDataPtr(ip),
237  std::numeric_limits<double>::quiet_NaN(),
238  distances[ip]});
239  }
240  }
241  if (_ip_data[k].non_local_assemblers.size() == 0)
242  {
243  OGS_FATAL("no neighbours found!");
244  }
245 
246  double a_k_sum_m = 0;
247  for (auto const& tuple : _ip_data[k].non_local_assemblers)
248  {
249  double const distance2_m = tuple.distance2;
250 
251  auto const& w_m = tuple.ip_l_pointer->integration_weight;
252 
253  a_k_sum_m += w_m * alpha_0(distance2_m);
254  }
255 
256  //
257  // Calculate alpha_kl =
258  // alpha_0(|x_k - x_l|) / int_{m \in ip} alpha_0(|x_k - x_m|)
259  //
260  for (auto& tuple : _ip_data[k].non_local_assemblers)
261  {
262  double const distance2_l = tuple.distance2;
263  double const a_kl = alpha_0(distance2_l) / a_k_sum_m;
264 
265  // Store the a_kl already multiplied with the integration
266  // weight of that l integration point.
267  auto const w_l = tuple.ip_l_pointer->integration_weight;
268  tuple.alpha_kl_times_w_l = a_kl * w_l;
269  }
270  }
271  }
272 
274  int integration_point) const
275  {
276  auto const& N = _secondary_data.N[integration_point];
277 
278  Eigen::Vector3d xyz = Eigen::Vector3d::Zero(); // Resulting coordinates
279  auto* nodes = _element.getNodes();
280  for (int i = 0; i < N.size(); ++i)
281  {
282  Eigen::Map<Eigen::Vector3d const> node_coordinates{
283  nodes[i]->getCoords(), 3};
284  xyz += node_coordinates * N[i];
285  }
286  return xyz;
287  }
288 
293  Eigen::Vector3d const& coords,
294  std::vector<double>& distances) const override
295  {
296  unsigned const n_integration_points =
297  _integration_method.getNumberOfPoints();
298 
299  distances.resize(n_integration_points);
300 
301  for (unsigned ip = 0; ip < n_integration_points; ip++)
302  {
303  auto const& xyz = _ip_data[ip].coordinates;
304  distances[ip] = (xyz - coords).squaredNorm();
305  }
306  }
307 
308  void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
309  std::vector<double>& /*local_M_data*/,
310  std::vector<double>& /*local_K_data*/,
311  std::vector<double>& /*local_b_data*/) override
312  {
313  OGS_FATAL(
314  "SmallDeformationNonlocalLocalAssembler: assembly without jacobian "
315  "is not "
316  "implemented.");
317  }
318 
319  void preAssemble(double const t,
320  std::vector<double> const& local_x) override
321  {
322  auto const n_integration_points =
323  _integration_method.getNumberOfPoints();
324 
326  x_position.setElementID(_element.getID());
327 
328  for (unsigned ip = 0; ip < n_integration_points; ip++)
329  {
330  x_position.setIntegrationPoint(ip);
331 
332  auto const& N = _ip_data[ip].N;
333  auto const& dNdx = _ip_data[ip].dNdx;
334 
335  auto const x_coord =
336  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
337  _element, N);
338  auto const B = LinearBMatrix::computeBMatrix<
339  DisplacementDim, ShapeFunction::NPOINTS,
340  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
341  _is_axially_symmetric);
342  auto const& eps_prev = _ip_data[ip].eps_prev;
343  auto const& sigma_prev = _ip_data[ip].sigma_prev;
344 
345  auto& eps = _ip_data[ip].eps;
346  auto& sigma = _ip_data[ip].sigma;
347  auto& C = _ip_data[ip].C;
348  auto& state = _ip_data[ip].material_state_variables;
349  double const& damage_prev = _ip_data[ip].damage_prev;
350 
351  eps.noalias() =
352  B *
353  Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
354  local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
355 
356  // sigma is for plastic part only.
357  std::unique_ptr<
359  new_C;
360  std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
361  DisplacementDim>::MaterialStateVariables>
362  new_state;
363 
364  // Compute sigma_eff from damage total stress sigma
366  KelvinVectorType const sigma_eff_prev =
367  sigma_prev /
368  (1. - damage_prev); // damage_prev is in [0,1) range. See
369  // calculateDamage() function.
370 
371  auto&& solution = _ip_data[ip].solid_material.integrateStress(
372  t, x_position, _process_data.dt, eps_prev, eps, sigma_eff_prev,
373  *state, _process_data.reference_temperature);
374 
375  if (!solution)
376  {
377  OGS_FATAL("Computation of local constitutive relation failed.");
378  }
379 
380  std::tie(sigma, state, C) = std::move(*solution);
381 
383  {
384  auto const& ehlers_material =
386  DisplacementDim> const&>(_ip_data[ip].solid_material);
387  auto const damage_properties =
388  ehlers_material.evaluatedDamageProperties(t, x_position);
389  auto const material_properties =
390  ehlers_material.evaluatedMaterialProperties(t, x_position);
391 
392  // Ehlers material state variables
393  auto& state_vars =
395  DisplacementDim>&>(
396  *_ip_data[ip].material_state_variables);
397 
398  double const eps_p_eff_diff =
399  state_vars.eps_p.eff - state_vars.eps_p_prev.eff;
400 
401  _ip_data[ip].kappa_d = calculateDamageKappaD<DisplacementDim>(
402  eps_p_eff_diff, sigma, _ip_data[ip].kappa_d_prev,
403  damage_properties.h_d, material_properties);
404 
405  if (!_ip_data[ip].active_self)
406  {
407  _ip_data[ip].active_self |= _ip_data[ip].kappa_d > 0;
408  if (_ip_data[ip].active_self)
409  {
410  for (auto const& tuple :
411  _ip_data[ip].non_local_assemblers)
412  {
413  // Activate the integration point.
414  tuple.ip_l_pointer->activated = true;
415  }
416  }
417  }
418  }
419  }
420  }
421 
422  void assembleWithJacobian(double const t,
423  std::vector<double> const& local_x,
424  std::vector<double> const& /*local_xdot*/,
425  const double /*dxdot_dx*/, const double /*dx_dx*/,
426  std::vector<double>& /*local_M_data*/,
427  std::vector<double>& /*local_K_data*/,
428  std::vector<double>& local_b_data,
429  std::vector<double>& local_Jac_data) override
430  {
431  auto const local_matrix_size = local_x.size();
432 
433  auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
434  local_Jac_data, local_matrix_size, local_matrix_size);
435 
436  auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
437  local_b_data, local_matrix_size);
438 
439  unsigned const n_integration_points =
440  _integration_method.getNumberOfPoints();
441 
443  x_position.setElementID(_element.getID());
444 
445  // Non-local integration.
446  for (unsigned ip = 0; ip < n_integration_points; ip++)
447  {
448  x_position.setIntegrationPoint(ip);
449  auto const& w = _ip_data[ip].integration_weight;
450 
451  auto const& N = _ip_data[ip].N;
452  auto const& dNdx = _ip_data[ip].dNdx;
453 
454  auto const x_coord =
455  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
456  _element, N);
457  auto const B = LinearBMatrix::computeBMatrix<
458  DisplacementDim, ShapeFunction::NPOINTS,
459  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
460  _is_axially_symmetric);
461 
462  auto& sigma = _ip_data[ip].sigma;
463  auto& C = _ip_data[ip].C;
464  double& damage = _ip_data[ip].damage;
465 
466  {
467  double nonlocal_kappa_d = 0;
468 
469  if (_ip_data[ip].active_self || _ip_data[ip].activated)
470  {
471  for (auto const& tuple : _ip_data[ip].non_local_assemblers)
472  {
473  // Get local variable for the integration point l.
474  double const kappa_d_l = tuple.ip_l_pointer->kappa_d;
475  double const a_kl_times_w_l = tuple.alpha_kl_times_w_l;
476  nonlocal_kappa_d += a_kl_times_w_l * kappa_d_l;
477  }
478  }
479 
480  auto const& ehlers_material =
482  DisplacementDim> const&>(_ip_data[ip].solid_material);
483 
484  //
485  // Overnonlocal formulation
486  //
487  // See (Di Luzio & Bazant 2005, IJSS) for details.
488  // The implementation would go here and would be for a given
489  // gamma_nonlocal:
490  //
491  // Update nonlocal damage with local damage (scaled with 1 -
492  // \gamma_{nonlocal}) for the current integration point and the
493  // nonlocal integral part.
494  // nonlocal_kappa_d = (1. - gamma_nonlocal) * kappa_d +
495  // gamma_nonlocal * nonlocal_kappa_d;
496 
497  nonlocal_kappa_d = std::max(0., nonlocal_kappa_d);
498 
499  // Update damage based on nonlocal kappa_d
500  {
501  auto const damage_properties =
502  ehlers_material.evaluatedDamageProperties(t,
503  x_position);
504  damage = calculateDamage(nonlocal_kappa_d,
505  damage_properties.alpha_d,
506  damage_properties.beta_d);
507  damage = std::max(0., damage);
508  }
509  sigma = sigma * (1. - damage);
510  }
511 
512  local_b.noalias() -= B.transpose() * sigma * w;
513  local_Jac.noalias() += B.transpose() * C * (1. - damage) * B * w;
514  }
515  }
516 
517  void initializeConcrete() override
518  {
519  unsigned const n_integration_points =
520  _integration_method.getNumberOfPoints();
521 
522  for (unsigned ip = 0; ip < n_integration_points; ip++)
523  {
524  _ip_data[ip].pushBackState();
525  }
526  }
527 
528  void postTimestepConcrete(std::vector<double> const& /*local_x*/) override
529  {
530  unsigned const n_integration_points =
531  _integration_method.getNumberOfPoints();
532 
533  for (unsigned ip = 0; ip < n_integration_points; ip++)
534  {
535  _ip_data[ip].pushBackState();
536  }
537  }
538 
539  void computeCrackIntegral(std::size_t mesh_item_id,
540  NumLib::LocalToGlobalIndexMap const& dof_table,
541  GlobalVector const& x,
542  double& crack_volume) override
543  {
544  auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
545  auto local_x = x.get(indices);
546 
547  auto u = Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
548  local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
549 
550  int const n_integration_points =
551  _integration_method.getNumberOfPoints();
552 
553  for (int ip = 0; ip < n_integration_points; ip++)
554  {
555  auto const& dNdx = _ip_data[ip].dNdx;
556  auto const& d = _ip_data[ip].damage;
557  auto const& w = _ip_data[ip].integration_weight;
558 
559  double const div_u =
560  Deformation::divergence<DisplacementDim,
561  ShapeFunction::NPOINTS>(u, dNdx);
562  crack_volume += div_u * d * w;
563  }
564  }
565 
566  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
567  const unsigned integration_point) const override
568  {
569  auto const& N = _secondary_data.N[integration_point];
570 
571  // assumes N is stored contiguously in memory
572  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
573  }
574 
575  std::vector<double> const& getNodalValues(
576  std::vector<double>& nodal_values) const override
577  {
578  nodal_values.clear();
579  auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
580  nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
581 
582  unsigned const n_integration_points =
583  _integration_method.getNumberOfPoints();
584 
585  for (unsigned ip = 0; ip < n_integration_points; ip++)
586  {
587  auto const& w = _ip_data[ip].integration_weight;
588 
589  auto const& N = _ip_data[ip].N;
590  auto const& dNdx = _ip_data[ip].dNdx;
591 
592  auto const x_coord =
593  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
594  _element, N);
595  auto const B = LinearBMatrix::computeBMatrix<
596  DisplacementDim, ShapeFunction::NPOINTS,
597  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
598  _is_axially_symmetric);
599  auto& sigma = _ip_data[ip].sigma;
600 
601  local_b.noalias() += B.transpose() * sigma * w;
602  }
603 
604  return nodal_values;
605  }
606 
607  std::vector<double> const& getIntPtFreeEnergyDensity(
608  const double /*t*/,
609  GlobalVector const& /*current_solution*/,
610  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
611  std::vector<double>& cache) const override
612  {
613  cache.clear();
614  cache.reserve(_ip_data.size());
615 
616  for (auto const& ip_data : _ip_data)
617  {
618  cache.push_back(ip_data.free_energy_density);
619  }
620 
621  return cache;
622  }
623 
624  std::vector<double> const& getIntPtEpsPV(
625  const double /*t*/,
626  GlobalVector const& /*current_solution*/,
627  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
628  std::vector<double>& cache) const override
629  {
630  cache.clear();
631  cache.reserve(_ip_data.size());
632 
633  for (auto const& ip_data : _ip_data)
634  {
635  cache.push_back(*ip_data.eps_p_V);
636  }
637 
638  return cache;
639  }
640 
641  std::vector<double> const& getIntPtEpsPDXX(
642  const double /*t*/,
643  GlobalVector const& /*current_solution*/,
644  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
645  std::vector<double>& cache) const override
646  {
647  cache.clear();
648  cache.reserve(_ip_data.size());
649 
650  for (auto const& ip_data : _ip_data)
651  {
652  cache.push_back(*ip_data.eps_p_D_xx);
653  }
654 
655  return cache;
656  }
657 
658  std::vector<double> const& getIntPtSigma(
659  const double /*t*/,
660  GlobalVector const& /*current_solution*/,
661  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
662  std::vector<double>& cache) const override
663  {
664  static const int kelvin_vector_size =
666  DisplacementDim>::value;
667  auto const num_intpts = _ip_data.size();
668 
669  cache.clear();
670  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
671  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
672  cache, kelvin_vector_size, num_intpts);
673 
674  for (unsigned ip = 0; ip < num_intpts; ++ip)
675  {
676  auto const& sigma = _ip_data[ip].sigma;
677  cache_mat.col(ip) =
679  }
680 
681  return cache;
682  }
683 
684  std::vector<double> const& getIntPtEpsilon(
685  const double /*t*/,
686  GlobalVector const& /*current_solution*/,
687  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
688  std::vector<double>& cache) const override
689  {
690  auto const kelvin_vector_size =
692  DisplacementDim>::value;
693  auto const num_intpts = _ip_data.size();
694 
695  cache.clear();
696  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
697  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
698  cache, kelvin_vector_size, num_intpts);
699 
700  for (unsigned ip = 0; ip < num_intpts; ++ip)
701  {
702  auto const& eps = _ip_data[ip].eps;
703  cache_mat.col(ip) =
705  }
706 
707  return cache;
708  }
709 
710  std::size_t setSigma(double const* values)
711  {
712  auto const kelvin_vector_size =
714  DisplacementDim>::value;
715  auto const n_integration_points = _ip_data.size();
716 
717  auto sigma_values =
718  Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
719  Eigen::ColMajor> const>(
720  values, kelvin_vector_size, n_integration_points);
721 
722  for (unsigned ip = 0; ip < n_integration_points; ++ip)
723  {
724  _ip_data[ip].sigma =
726  sigma_values.col(ip));
727  }
728 
729  return n_integration_points;
730  }
731 
732  // TODO (naumov) This method is same as getIntPtSigma but for arguments and
733  // the ordering of the cache_mat.
734  // There should be only one.
735  std::vector<double> getSigma() const override
736  {
737  auto const kelvin_vector_size =
739  DisplacementDim>::value;
740  auto const n_integration_points = _ip_data.size();
741 
742  std::vector<double> ip_sigma_values;
743  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
744  double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
745  ip_sigma_values, n_integration_points, kelvin_vector_size);
746 
747  for (unsigned ip = 0; ip < n_integration_points; ++ip)
748  {
749  auto const& sigma = _ip_data[ip].sigma;
750  cache_mat.row(ip) =
752  }
753 
754  return ip_sigma_values;
755  }
756 
757  std::size_t setKappaD(double const* values)
758  {
759  unsigned const n_integration_points =
760  _integration_method.getNumberOfPoints();
761 
762  for (unsigned ip = 0; ip < n_integration_points; ++ip)
763  {
764  _ip_data[ip].kappa_d = values[ip];
765  }
766  return n_integration_points;
767  }
768 
769  void setKappaD(double value)
770  {
771  for (auto& ip_data : _ip_data)
772  {
773  ip_data.kappa_d = value;
774  }
775  }
776  std::vector<double> getKappaD() const override
777  {
778  unsigned const n_integration_points =
779  _integration_method.getNumberOfPoints();
780 
781  std::vector<double> result_values;
782  result_values.resize(n_integration_points);
783 
784  for (unsigned ip = 0; ip < n_integration_points; ++ip)
785  {
786  result_values[ip] = _ip_data[ip].kappa_d;
787  }
788 
789  return result_values;
790  }
791 
792  std::vector<double> const& getIntPtDamage(
793  const double /*t*/,
794  GlobalVector const& /*current_solution*/,
795  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
796  std::vector<double>& cache) const override
797  {
798  cache.clear();
799  cache.reserve(_ip_data.size());
800 
801  for (auto const& ip_data : _ip_data)
802  {
803  cache.push_back(ip_data.damage);
804  }
805 
806  return cache;
807  }
808 
809  unsigned getNumberOfIntegrationPoints() const override
810  {
811  return _integration_method.getNumberOfPoints();
812  }
813 
815  DisplacementDim>::MaterialStateVariables const&
816  getMaterialStateVariablesAt(int const integration_point) const override
817  {
818  return *_ip_data[integration_point].material_state_variables;
819  }
820 
821 private:
822  std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
823  std::size_t const component) const
824  {
825  cache.clear();
826  cache.reserve(_ip_data.size());
827 
828  for (auto const& ip_data : _ip_data)
829  {
830  if (component < 3)
831  { // xx, yy, zz components
832  cache.push_back(ip_data.sigma[component]);
833  }
834  else
835  { // mixed xy, yz, xz components
836  cache.push_back(ip_data.sigma[component] / std::sqrt(2));
837  }
838  }
839 
840  return cache;
841  }
842 
843  std::vector<double> const& getIntPtEpsilon(
844  std::vector<double>& cache, std::size_t const component) const
845  {
846  cache.clear();
847  cache.reserve(_ip_data.size());
848 
849  for (auto const& ip_data : _ip_data)
850  {
851  if (component < 3) // xx, yy, zz components
852  cache.push_back(ip_data.eps[component]);
853  else // mixed xy, yz, xz components
854  cache.push_back(ip_data.eps[component] / std::sqrt(2));
855  }
856 
857  return cache;
858  }
859 
861  getIPDataPtr(int const ip) override
862  {
863  return &_ip_data[ip];
864  }
865 
866 private:
868 
869  std::vector<
871  Eigen::aligned_allocator<IntegrationPointData<
872  BMatricesType, ShapeMatricesType, DisplacementDim>>>
874 
875  IntegrationMethod const _integration_method;
879 
880  static const int displacement_size =
881  ShapeFunction::NPOINTS * DisplacementDim;
882 };
883 
884 } // namespace SmallDeformationNonlocal
885 } // namespace ProcessLib
IntegrationPointDataNonlocalInterface * getIPDataPtr(int const ip) override
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
Definition: KelvinVector.h:173
std::vector< double > const & getNodalValues(std::vector< double > &nodal_values) const override
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
double divergence(const Eigen::Ref< Eigen::Matrix< double, NPOINTS *DisplacementDim, 1 > const > &u, DNDX_Type const &dNdx)
Divergence of displacement, the volumetric strain.
Definition: Divergence.h:20
void setElementID(std::size_t element_id)
MatrixType< _number_of_dof, _number_of_dof > StiffnessMatrixType
Definition: BMatrixPolicy.h:40
VectorType< _number_of_dof > NodalForceVectorType
Rhs residual.
Definition: BMatrixPolicy.h:43
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool is_axially_symmetric, IntegrationMethod const &integration_method)
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)
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, KelvinVectorDimensions< DisplacementDim >::value, Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:59
void setIPDataInitialConditionsFromCellData(std::string const &name, std::vector< double > const &value) override
PlasticStrain< KelvinVector > eps_p
plastic part of the state.
Definition: Ehlers.h:239
void preAssemble(double const t, std::vector< double > const &local_x) override
void getIntegrationPointCoordinates(Eigen::Vector3d const &coords, std::vector< double > &distances) const override
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
void setIntegrationPoint(unsigned integration_point)
std::vector< double > const & getIntPtEpsPV(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
void assembleWithJacobian(double const t, std::vector< double > const &local_x, std::vector< double > const &, const double, const double, std::vector< double > &, std::vector< double > &, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:31
std::vector< double > const & getIntPtEpsPDXX(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
std::vector< double > const & getIntPtFreeEnergyDensity(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
DamageProperties evaluatedDamageProperties(double const t, ParameterLib::SpatialPosition const &x) const
Definition: Ehlers.h:339
std::vector< double > const & getIntPtDamage(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt(int const integration_point) const override
double calculateDamage(double const kappa_d, double const alpha_d, double const beta_d)
Definition: Damage.h:25
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
Coordinates mapping matrices at particular location.
Definition: ShapeMatrices.h:45
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
VectorType< ShapeFunction::NPOINTS > NodalVectorType
void nonlocal(std::size_t const, std::vector< std::unique_ptr< SmallDeformationNonlocalLocalAssemblerInterface< DisplacementDim >>> const &local_assemblers) override
VectorType< GlobalDim > GlobalDimVectorType
std::vector< std::size_t > findElementsWithinRadius(Element const &start_element, double const radius_squared)
SmallDeformationNonlocalLocalAssembler(MeshLib::Element const &e, std::size_t const, bool const is_axially_symmetric, unsigned const integration_order, SmallDeformationNonlocalProcessData< DisplacementDim > &process_data)
std::vector< double > const & getIntPtSigma(std::vector< double > &cache, std::size_t const component) const
void assemble(double const, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
std::vector< double > const & getIntPtEpsilon(std::vector< double > &cache, std::size_t const component) const
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
void computeCrackIntegral(std::size_t mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double &crack_volume) override
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
std::size_t setIPDataInitialConditions(std::string const &name, double const *values, int const integration_order) override
std::vector< double > const & getIntPtSigma(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:54
std::vector< double > const & getIntPtEpsilon(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
std::vector< IntegrationPointData< BMatricesType, ShapeMatricesType, DisplacementDim >, Eigen::aligned_allocator< IntegrationPointData< BMatricesType, ShapeMatricesType, DisplacementDim > > > _ip_data
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49