OGS 6.2.0-97-g4a610c866
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 preTimestepConcrete(std::vector<double> const& /*local_x*/,
518  double const /*t*/,
519  double const /*delta_t*/) override
520  {
521  unsigned const n_integration_points =
522  _integration_method.getNumberOfPoints();
523 
524  for (unsigned ip = 0; ip < n_integration_points; ip++)
525  {
526  _ip_data[ip].pushBackState();
527  }
528  }
529 
530  void computeCrackIntegral(std::size_t mesh_item_id,
531  NumLib::LocalToGlobalIndexMap const& dof_table,
532  GlobalVector const& x,
533  double& crack_volume) override
534  {
535  auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
536  auto local_x = x.get(indices);
537 
538  auto u = Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
539  local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
540 
541  int const n_integration_points =
542  _integration_method.getNumberOfPoints();
543 
544  for (int ip = 0; ip < n_integration_points; ip++)
545  {
546  auto const& dNdx = _ip_data[ip].dNdx;
547  auto const& d = _ip_data[ip].damage;
548  auto const& w = _ip_data[ip].integration_weight;
549 
550  double const div_u =
551  Deformation::divergence<DisplacementDim,
552  ShapeFunction::NPOINTS>(u, dNdx);
553  crack_volume += div_u * d * w;
554  }
555  }
556 
557  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
558  const unsigned integration_point) const override
559  {
560  auto const& N = _secondary_data.N[integration_point];
561 
562  // assumes N is stored contiguously in memory
563  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
564  }
565 
566  std::vector<double> const& getNodalValues(
567  std::vector<double>& nodal_values) const override
568  {
569  nodal_values.clear();
570  auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
571  nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
572 
573  unsigned const n_integration_points =
574  _integration_method.getNumberOfPoints();
575 
576  for (unsigned ip = 0; ip < n_integration_points; ip++)
577  {
578  auto const& w = _ip_data[ip].integration_weight;
579 
580  auto const& N = _ip_data[ip].N;
581  auto const& dNdx = _ip_data[ip].dNdx;
582 
583  auto const x_coord =
584  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
585  _element, N);
586  auto const B = LinearBMatrix::computeBMatrix<
587  DisplacementDim, ShapeFunction::NPOINTS,
588  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
589  _is_axially_symmetric);
590  auto& sigma = _ip_data[ip].sigma;
591 
592  local_b.noalias() += B.transpose() * sigma * w;
593  }
594 
595  return nodal_values;
596  }
597 
598  std::vector<double> const& getIntPtFreeEnergyDensity(
599  const double /*t*/,
600  GlobalVector const& /*current_solution*/,
601  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
602  std::vector<double>& cache) const override
603  {
604  cache.clear();
605  cache.reserve(_ip_data.size());
606 
607  for (auto const& ip_data : _ip_data)
608  {
609  cache.push_back(ip_data.free_energy_density);
610  }
611 
612  return cache;
613  }
614 
615  std::vector<double> const& getIntPtEpsPV(
616  const double /*t*/,
617  GlobalVector const& /*current_solution*/,
618  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
619  std::vector<double>& cache) const override
620  {
621  cache.clear();
622  cache.reserve(_ip_data.size());
623 
624  for (auto const& ip_data : _ip_data)
625  {
626  cache.push_back(*ip_data.eps_p_V);
627  }
628 
629  return cache;
630  }
631 
632  std::vector<double> const& getIntPtEpsPDXX(
633  const double /*t*/,
634  GlobalVector const& /*current_solution*/,
635  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
636  std::vector<double>& cache) const override
637  {
638  cache.clear();
639  cache.reserve(_ip_data.size());
640 
641  for (auto const& ip_data : _ip_data)
642  {
643  cache.push_back(*ip_data.eps_p_D_xx);
644  }
645 
646  return cache;
647  }
648 
649  std::vector<double> const& getIntPtSigma(
650  const double /*t*/,
651  GlobalVector const& /*current_solution*/,
652  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
653  std::vector<double>& cache) const override
654  {
655  static const int kelvin_vector_size =
657  DisplacementDim>::value;
658  auto const num_intpts = _ip_data.size();
659 
660  cache.clear();
661  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
662  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
663  cache, kelvin_vector_size, num_intpts);
664 
665  for (unsigned ip = 0; ip < num_intpts; ++ip)
666  {
667  auto const& sigma = _ip_data[ip].sigma;
668  cache_mat.col(ip) =
670  }
671 
672  return cache;
673  }
674 
675  std::vector<double> const& getIntPtEpsilon(
676  const double /*t*/,
677  GlobalVector const& /*current_solution*/,
678  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
679  std::vector<double>& cache) const override
680  {
681  auto const kelvin_vector_size =
683  DisplacementDim>::value;
684  auto const num_intpts = _ip_data.size();
685 
686  cache.clear();
687  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
688  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
689  cache, kelvin_vector_size, num_intpts);
690 
691  for (unsigned ip = 0; ip < num_intpts; ++ip)
692  {
693  auto const& eps = _ip_data[ip].eps;
694  cache_mat.col(ip) =
696  }
697 
698  return cache;
699  }
700 
701  std::size_t setSigma(double const* values)
702  {
703  auto const kelvin_vector_size =
705  DisplacementDim>::value;
706  auto const n_integration_points = _ip_data.size();
707 
708  std::vector<double> ip_sigma_values;
709  auto sigma_values =
710  Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
711  Eigen::ColMajor> const>(
712  values, kelvin_vector_size, n_integration_points);
713 
714  for (unsigned ip = 0; ip < n_integration_points; ++ip)
715  {
716  _ip_data[ip].sigma =
718  sigma_values.col(ip));
719  }
720 
721  return n_integration_points;
722  }
723 
724  // TODO (naumov) This method is same as getIntPtSigma but for arguments and
725  // the ordering of the cache_mat.
726  // There should be only one.
727  std::vector<double> getSigma() const override
728  {
729  auto const kelvin_vector_size =
731  DisplacementDim>::value;
732  auto const n_integration_points = _ip_data.size();
733 
734  std::vector<double> ip_sigma_values;
735  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
736  double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
737  ip_sigma_values, n_integration_points, kelvin_vector_size);
738 
739  for (unsigned ip = 0; ip < n_integration_points; ++ip)
740  {
741  auto const& sigma = _ip_data[ip].sigma;
742  cache_mat.row(ip) =
744  }
745 
746  return ip_sigma_values;
747  }
748 
749  std::size_t setKappaD(double const* values)
750  {
751  unsigned const n_integration_points =
752  _integration_method.getNumberOfPoints();
753 
754  for (unsigned ip = 0; ip < n_integration_points; ++ip)
755  {
756  _ip_data[ip].kappa_d = values[ip];
757  }
758  return n_integration_points;
759  }
760 
761  void setKappaD(double value)
762  {
763  for (auto& ip_data : _ip_data)
764  {
765  ip_data.kappa_d = value;
766  }
767  }
768  std::vector<double> getKappaD() const override
769  {
770  unsigned const n_integration_points =
771  _integration_method.getNumberOfPoints();
772 
773  std::vector<double> result_values;
774  result_values.resize(n_integration_points);
775 
776  for (unsigned ip = 0; ip < n_integration_points; ++ip)
777  {
778  result_values[ip] = _ip_data[ip].kappa_d;
779  }
780 
781  return result_values;
782  }
783 
784  std::vector<double> const& getIntPtDamage(
785  const double /*t*/,
786  GlobalVector const& /*current_solution*/,
787  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
788  std::vector<double>& cache) const override
789  {
790  cache.clear();
791  cache.reserve(_ip_data.size());
792 
793  for (auto const& ip_data : _ip_data)
794  {
795  cache.push_back(ip_data.damage);
796  }
797 
798  return cache;
799  }
800 
801  unsigned getNumberOfIntegrationPoints() const override
802  {
803  return _integration_method.getNumberOfPoints();
804  }
805 
807  DisplacementDim>::MaterialStateVariables const&
808  getMaterialStateVariablesAt(int const integration_point) const override
809  {
810  return *_ip_data[integration_point].material_state_variables;
811  }
812 
813 private:
814  std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
815  std::size_t const component) const
816  {
817  cache.clear();
818  cache.reserve(_ip_data.size());
819 
820  for (auto const& ip_data : _ip_data)
821  {
822  if (component < 3)
823  { // xx, yy, zz components
824  cache.push_back(ip_data.sigma[component]);
825  }
826  else
827  { // mixed xy, yz, xz components
828  cache.push_back(ip_data.sigma[component] / std::sqrt(2));
829  }
830  }
831 
832  return cache;
833  }
834 
835  std::vector<double> const& getIntPtEpsilon(
836  std::vector<double>& cache, std::size_t const component) const
837  {
838  cache.clear();
839  cache.reserve(_ip_data.size());
840 
841  for (auto const& ip_data : _ip_data)
842  {
843  if (component < 3) // xx, yy, zz components
844  cache.push_back(ip_data.eps[component]);
845  else // mixed xy, yz, xz components
846  cache.push_back(ip_data.eps[component] / std::sqrt(2));
847  }
848 
849  return cache;
850  }
851 
853  getIPDataPtr(int const ip) override
854  {
855  return &_ip_data[ip];
856  }
857 
858 private:
860 
861  std::vector<
863  Eigen::aligned_allocator<IntegrationPointData<
864  BMatricesType, ShapeMatricesType, DisplacementDim>>>
866 
867  IntegrationMethod const _integration_method;
871 
872  static const int displacement_size =
873  ShapeFunction::NPOINTS * DisplacementDim;
874 };
875 
876 } // namespace SmallDeformationNonlocal
877 } // 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:248
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:342
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
void preTimestepConcrete(std::vector< double > const &, double const, double const) 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