OGS 6.1.0-1721-g6382411ad
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 
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  continue;
232  // save into current ip_k
233  _ip_data[k].non_local_assemblers.push_back(
234  {la->getIPDataPtr(ip),
235  std::numeric_limits<double>::quiet_NaN(),
236  distances[ip]});
237  }
238  }
239  if (_ip_data[k].non_local_assemblers.size() == 0)
240  {
241  OGS_FATAL("no neighbours found!");
242  }
243 
244  double a_k_sum_m = 0;
245  for (auto const& tuple : _ip_data[k].non_local_assemblers)
246  {
247  double const distance2_m = tuple.distance2;
248 
249  auto const& w_m = tuple.ip_l_pointer->integration_weight;
250 
251  a_k_sum_m += w_m * alpha_0(distance2_m);
252  }
253 
254  //
255  // Calculate alpha_kl =
256  // alpha_0(|x_k - x_l|) / int_{m \in ip} alpha_0(|x_k - x_m|)
257  //
258  for (auto& tuple : _ip_data[k].non_local_assemblers)
259  {
260  double const distance2_l = tuple.distance2;
261  double const a_kl = alpha_0(distance2_l) / a_k_sum_m;
262 
263  // Store the a_kl already multiplied with the integration
264  // weight of that l integration point.
265  auto const w_l = tuple.ip_l_pointer->integration_weight;
266  tuple.alpha_kl_times_w_l = a_kl * w_l;
267  }
268  }
269  }
270 
272  int integration_point) const
273  {
274  auto const& N = _secondary_data.N[integration_point];
275 
276  Eigen::Vector3d xyz = Eigen::Vector3d::Zero(); // Resulting coordinates
277  auto* nodes = _element.getNodes();
278  for (int i = 0; i < N.size(); ++i)
279  {
280  Eigen::Map<Eigen::Vector3d const> node_coordinates{
281  nodes[i]->getCoords(), 3};
282  xyz += node_coordinates * N[i];
283  }
284  return xyz;
285  }
286 
291  Eigen::Vector3d const& coords,
292  std::vector<double>& distances) const override
293  {
294  unsigned const n_integration_points =
295  _integration_method.getNumberOfPoints();
296 
297  distances.resize(n_integration_points);
298 
299  for (unsigned ip = 0; ip < n_integration_points; ip++)
300  {
301  auto const& xyz = _ip_data[ip].coordinates;
302  distances[ip] = (xyz - coords).squaredNorm();
303  }
304  }
305 
306  void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
307  std::vector<double>& /*local_M_data*/,
308  std::vector<double>& /*local_K_data*/,
309  std::vector<double>& /*local_b_data*/) override
310  {
311  OGS_FATAL(
312  "SmallDeformationNonlocalLocalAssembler: assembly without jacobian "
313  "is not "
314  "implemented.");
315  }
316 
317  void preAssemble(double const t,
318  std::vector<double> const& local_x) override
319  {
320  auto const n_integration_points =
321  _integration_method.getNumberOfPoints();
322 
323  SpatialPosition x_position;
324  x_position.setElementID(_element.getID());
325 
326  for (unsigned ip = 0; ip < n_integration_points; ip++)
327  {
328  x_position.setIntegrationPoint(ip);
329 
330  auto const& N = _ip_data[ip].N;
331  auto const& dNdx = _ip_data[ip].dNdx;
332 
333  auto const x_coord =
334  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
335  _element, N);
336  auto const B = LinearBMatrix::computeBMatrix<
337  DisplacementDim, ShapeFunction::NPOINTS,
338  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
339  _is_axially_symmetric);
340  auto const& eps_prev = _ip_data[ip].eps_prev;
341  auto const& sigma_prev = _ip_data[ip].sigma_prev;
342 
343  auto& eps = _ip_data[ip].eps;
344  auto& sigma = _ip_data[ip].sigma;
345  auto& C = _ip_data[ip].C;
346  auto& state = _ip_data[ip].material_state_variables;
347  double const& damage_prev = _ip_data[ip].damage_prev;
348 
349  eps.noalias() =
350  B *
351  Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
352  local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
353 
354  // sigma is for plastic part only.
355  std::unique_ptr<
357  new_C;
358  std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
359  DisplacementDim>::MaterialStateVariables>
360  new_state;
361 
362  // Compute sigma_eff from damage total stress sigma
364  KelvinVectorType const sigma_eff_prev =
365  sigma_prev /
366  (1. - damage_prev); // damage_prev is in [0,1) range. See
367  // calculateDamage() function.
368 
369  auto&& solution = _ip_data[ip].solid_material.integrateStress(
370  t, x_position, _process_data.dt, eps_prev, eps, sigma_eff_prev,
371  *state, _process_data.reference_temperature);
372 
373  if (!solution)
374  OGS_FATAL("Computation of local constitutive relation failed.");
375 
376  std::tie(sigma, state, C) = std::move(*solution);
377 
379  {
380  auto const& ehlers_material =
382  DisplacementDim> const&>(_ip_data[ip].solid_material);
383  auto const damage_properties =
384  ehlers_material.evaluatedDamageProperties(t, x_position);
385  auto const material_properties =
386  ehlers_material.evaluatedMaterialProperties(t, x_position);
387 
388  // Ehlers material state variables
389  auto& state_vars =
391  DisplacementDim>&>(
392  *_ip_data[ip].material_state_variables);
393 
394  double const eps_p_eff_diff =
395  state_vars.eps_p.eff - state_vars.eps_p_prev.eff;
396 
397  _ip_data[ip].kappa_d = calculateDamageKappaD<DisplacementDim>(
398  eps_p_eff_diff, sigma, _ip_data[ip].kappa_d_prev,
399  damage_properties.h_d, material_properties);
400 
401  if (!_ip_data[ip].active_self)
402  {
403  _ip_data[ip].active_self |= _ip_data[ip].kappa_d > 0;
404  if (_ip_data[ip].active_self)
405  {
406  for (auto const& tuple :
407  _ip_data[ip].non_local_assemblers)
408  {
409  // Activate the integration point.
410  tuple.ip_l_pointer->activated = true;
411  }
412  }
413  }
414  }
415  }
416  }
417 
418  void assembleWithJacobian(double const t,
419  std::vector<double> const& local_x,
420  std::vector<double> const& /*local_xdot*/,
421  const double /*dxdot_dx*/, const double /*dx_dx*/,
422  std::vector<double>& /*local_M_data*/,
423  std::vector<double>& /*local_K_data*/,
424  std::vector<double>& local_b_data,
425  std::vector<double>& local_Jac_data) override
426  {
427  auto const local_matrix_size = local_x.size();
428 
429  auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
430  local_Jac_data, local_matrix_size, local_matrix_size);
431 
432  auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
433  local_b_data, local_matrix_size);
434 
435  unsigned const n_integration_points =
436  _integration_method.getNumberOfPoints();
437 
438  SpatialPosition x_position;
439  x_position.setElementID(_element.getID());
440 
441  // Non-local integration.
442  for (unsigned ip = 0; ip < n_integration_points; ip++)
443  {
444  x_position.setIntegrationPoint(ip);
445  auto const& w = _ip_data[ip].integration_weight;
446 
447  auto const& N = _ip_data[ip].N;
448  auto const& dNdx = _ip_data[ip].dNdx;
449 
450  auto const x_coord =
451  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
452  _element, N);
453  auto const B = LinearBMatrix::computeBMatrix<
454  DisplacementDim, ShapeFunction::NPOINTS,
455  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
456  _is_axially_symmetric);
457 
458  auto& sigma = _ip_data[ip].sigma;
459  auto& C = _ip_data[ip].C;
460  double& damage = _ip_data[ip].damage;
461 
462  {
463  double nonlocal_kappa_d = 0;
464 
465  if (_ip_data[ip].active_self || _ip_data[ip].activated)
466  {
467  for (auto const& tuple : _ip_data[ip].non_local_assemblers)
468  {
469  // Get local variable for the integration point l.
470  double const kappa_d_l = tuple.ip_l_pointer->kappa_d;
471  double const a_kl_times_w_l = tuple.alpha_kl_times_w_l;
472  nonlocal_kappa_d += a_kl_times_w_l * kappa_d_l;
473  }
474  }
475 
476  auto const& ehlers_material =
478  DisplacementDim> const&>(_ip_data[ip].solid_material);
479 
480  //
481  // Overnonlocal formulation
482  //
483  // See (Di Luzio & Bazant 2005, IJSS) for details.
484  // The implementation would go here and would be for a given
485  // gamma_nonlocal:
486  //
487  // Update nonlocal damage with local damage (scaled with 1 -
488  // \gamma_{nonlocal}) for the current integration point and the
489  // nonlocal integral part.
490  // nonlocal_kappa_d = (1. - gamma_nonlocal) * kappa_d +
491  // gamma_nonlocal * nonlocal_kappa_d;
492 
493  nonlocal_kappa_d = std::max(0., nonlocal_kappa_d);
494 
495  // Update damage based on nonlocal kappa_d
496  {
497  auto const damage_properties =
498  ehlers_material.evaluatedDamageProperties(t,
499  x_position);
500  damage = calculateDamage(nonlocal_kappa_d,
501  damage_properties.alpha_d,
502  damage_properties.beta_d);
503  damage = std::max(0., damage);
504  }
505  sigma = sigma * (1. - damage);
506  }
507 
508  local_b.noalias() -= B.transpose() * sigma * w;
509  local_Jac.noalias() += B.transpose() * C * (1. - damage) * B * w;
510  }
511  }
512 
513  void preTimestepConcrete(std::vector<double> const& /*local_x*/,
514  double const /*t*/,
515  double const /*delta_t*/) override
516  {
517  unsigned const n_integration_points =
518  _integration_method.getNumberOfPoints();
519 
520  for (unsigned ip = 0; ip < n_integration_points; ip++)
521  {
522  _ip_data[ip].pushBackState();
523  }
524  }
525 
526  void computeCrackIntegral(std::size_t mesh_item_id,
527  NumLib::LocalToGlobalIndexMap const& dof_table,
528  GlobalVector const& x,
529  double& crack_volume) override
530  {
531  auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
532  auto local_x = x.get(indices);
533 
534  auto u = Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
535  local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
536 
537  int const n_integration_points =
538  _integration_method.getNumberOfPoints();
539 
540  for (int ip = 0; ip < n_integration_points; ip++)
541  {
542  auto const& dNdx = _ip_data[ip].dNdx;
543  auto const& d = _ip_data[ip].damage;
544  auto const& w = _ip_data[ip].integration_weight;
545 
546  double const div_u =
547  Deformation::divergence<DisplacementDim,
548  ShapeFunction::NPOINTS>(u, dNdx);
549  crack_volume += div_u * d * w;
550  }
551  }
552 
553  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
554  const unsigned integration_point) const override
555  {
556  auto const& N = _secondary_data.N[integration_point];
557 
558  // assumes N is stored contiguously in memory
559  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
560  }
561 
562  std::vector<double> const& getNodalValues(
563  std::vector<double>& nodal_values) const override
564  {
565  nodal_values.clear();
566  auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
567  nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
568 
569  unsigned const n_integration_points =
570  _integration_method.getNumberOfPoints();
571 
572  for (unsigned ip = 0; ip < n_integration_points; ip++)
573  {
574  auto const& w = _ip_data[ip].integration_weight;
575 
576  auto const& N = _ip_data[ip].N;
577  auto const& dNdx = _ip_data[ip].dNdx;
578 
579  auto const x_coord =
580  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
581  _element, N);
582  auto const B = LinearBMatrix::computeBMatrix<
583  DisplacementDim, ShapeFunction::NPOINTS,
584  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
585  _is_axially_symmetric);
586  auto& sigma = _ip_data[ip].sigma;
587 
588  local_b.noalias() += B.transpose() * sigma * w;
589  }
590 
591  return nodal_values;
592  }
593 
594  std::vector<double> const& getIntPtFreeEnergyDensity(
595  const double /*t*/,
596  GlobalVector const& /*current_solution*/,
597  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
598  std::vector<double>& cache) const override
599  {
600  cache.clear();
601  cache.reserve(_ip_data.size());
602 
603  for (auto const& ip_data : _ip_data)
604  {
605  cache.push_back(ip_data.free_energy_density);
606  }
607 
608  return cache;
609  }
610 
611  std::vector<double> const& getIntPtEpsPV(
612  const double /*t*/,
613  GlobalVector const& /*current_solution*/,
614  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
615  std::vector<double>& cache) const override
616  {
617  cache.clear();
618  cache.reserve(_ip_data.size());
619 
620  for (auto const& ip_data : _ip_data)
621  {
622  cache.push_back(*ip_data.eps_p_V);
623  }
624 
625  return cache;
626  }
627 
628  std::vector<double> const& getIntPtEpsPDXX(
629  const double /*t*/,
630  GlobalVector const& /*current_solution*/,
631  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
632  std::vector<double>& cache) const override
633  {
634  cache.clear();
635  cache.reserve(_ip_data.size());
636 
637  for (auto const& ip_data : _ip_data)
638  {
639  cache.push_back(*ip_data.eps_p_D_xx);
640  }
641 
642  return cache;
643  }
644 
645  std::vector<double> const& getIntPtSigma(
646  const double /*t*/,
647  GlobalVector const& /*current_solution*/,
648  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
649  std::vector<double>& cache) const override
650  {
651  static const int kelvin_vector_size =
653  DisplacementDim>::value;
654  auto const num_intpts = _ip_data.size();
655 
656  cache.clear();
657  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
658  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
659  cache, kelvin_vector_size, num_intpts);
660 
661  for (unsigned ip = 0; ip < num_intpts; ++ip)
662  {
663  auto const& sigma = _ip_data[ip].sigma;
664  cache_mat.col(ip) =
666  }
667 
668  return cache;
669  }
670 
671  virtual std::vector<double> const& getIntPtEpsilon(
672  const double /*t*/,
673  GlobalVector const& /*current_solution*/,
674  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
675  std::vector<double>& cache) const override
676  {
677  auto const kelvin_vector_size =
679  DisplacementDim>::value;
680  auto const num_intpts = _ip_data.size();
681 
682  cache.clear();
683  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
684  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
685  cache, kelvin_vector_size, num_intpts);
686 
687  for (unsigned ip = 0; ip < num_intpts; ++ip)
688  {
689  auto const& eps = _ip_data[ip].eps;
690  cache_mat.col(ip) =
692  }
693 
694  return cache;
695  }
696 
697  std::size_t setSigma(double const* values)
698  {
699  auto const kelvin_vector_size =
701  DisplacementDim>::value;
702  auto const n_integration_points = _ip_data.size();
703 
704  std::vector<double> ip_sigma_values;
705  auto sigma_values =
706  Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
707  Eigen::ColMajor> const>(
708  values, kelvin_vector_size, n_integration_points);
709 
710  for (unsigned ip = 0; ip < n_integration_points; ++ip)
711  {
712  _ip_data[ip].sigma =
714  sigma_values.col(ip));
715  }
716 
717  return n_integration_points;
718  }
719 
720  // TODO (naumov) This method is same as getIntPtSigma but for arguments and
721  // the ordering of the cache_mat.
722  // There should be only one.
723  std::vector<double> getSigma() const override
724  {
725  auto const kelvin_vector_size =
727  DisplacementDim>::value;
728  auto const n_integration_points = _ip_data.size();
729 
730  std::vector<double> ip_sigma_values;
731  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
732  double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
733  ip_sigma_values, n_integration_points, kelvin_vector_size);
734 
735  for (unsigned ip = 0; ip < n_integration_points; ++ip)
736  {
737  auto const& sigma = _ip_data[ip].sigma;
738  cache_mat.row(ip) =
740  }
741 
742  return ip_sigma_values;
743  }
744 
745  std::size_t setKappaD(double const* values)
746  {
747  unsigned const n_integration_points =
748  _integration_method.getNumberOfPoints();
749 
750  for (unsigned ip = 0; ip < n_integration_points; ++ip)
751  {
752  _ip_data[ip].kappa_d = values[ip];
753  }
754  return n_integration_points;
755  }
756 
757  void setKappaD(double value)
758  {
759  for (auto& ip_data : _ip_data)
760  {
761  ip_data.kappa_d = value;
762  }
763  }
764  std::vector<double> getKappaD() const override
765  {
766  unsigned const n_integration_points =
767  _integration_method.getNumberOfPoints();
768 
769  std::vector<double> result_values;
770  result_values.resize(n_integration_points);
771 
772  for (unsigned ip = 0; ip < n_integration_points; ++ip)
773  {
774  result_values[ip] = _ip_data[ip].kappa_d;
775  }
776 
777  return result_values;
778  }
779 
780  std::vector<double> const& getIntPtDamage(
781  const double /*t*/,
782  GlobalVector const& /*current_solution*/,
783  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
784  std::vector<double>& cache) const override
785  {
786  cache.clear();
787  cache.reserve(_ip_data.size());
788 
789  for (auto const& ip_data : _ip_data)
790  {
791  cache.push_back(ip_data.damage);
792  }
793 
794  return cache;
795  }
796 
797  unsigned getNumberOfIntegrationPoints() const override
798  {
799  return _integration_method.getNumberOfPoints();
800  }
801 
803  DisplacementDim>::MaterialStateVariables const&
804  getMaterialStateVariablesAt(int const integration_point) const override
805  {
806  return *_ip_data[integration_point].material_state_variables;
807  }
808 
809 private:
810  std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
811  std::size_t const component) const
812  {
813  cache.clear();
814  cache.reserve(_ip_data.size());
815 
816  for (auto const& ip_data : _ip_data)
817  {
818  if (component < 3) // xx, yy, zz components
819  cache.push_back(ip_data.sigma[component]);
820  else // mixed xy, yz, xz components
821  cache.push_back(ip_data.sigma[component] / std::sqrt(2));
822  }
823 
824  return cache;
825  }
826 
827  std::vector<double> const& getIntPtEpsilon(
828  std::vector<double>& cache, std::size_t const component) const
829  {
830  cache.clear();
831  cache.reserve(_ip_data.size());
832 
833  for (auto const& ip_data : _ip_data)
834  {
835  if (component < 3) // xx, yy, zz components
836  cache.push_back(ip_data.eps[component]);
837  else // mixed xy, yz, xz components
838  cache.push_back(ip_data.eps[component] / std::sqrt(2));
839  }
840 
841  return cache;
842  }
843 
845  getIPDataPtr(int const ip) override
846  {
847  return &_ip_data[ip];
848  }
849 
850 private:
852 
853  std::vector<
855  Eigen::aligned_allocator<IntegrationPointData<
856  BMatricesType, ShapeMatricesType, DisplacementDim>>>
858 
859  IntegrationMethod const _integration_method;
863 
864  static const int displacement_size =
865  ShapeFunction::NPOINTS * DisplacementDim;
866 };
867 
868 } // namespace SmallDeformationNonlocal
869 } // 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
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
void setElementID(std::size_t element_id)
PlasticStrain< KelvinVector > eps_p
plastic part of the state.
Definition: Ehlers.h:242
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.
std::vector< double > const & getIntPtEpsPV(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
DamageProperties evaluatedDamageProperties(double const t, ProcessLib::SpatialPosition const &x) const
Definition: Ehlers.h:336
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
virtual std::vector< double > const & getIntPtEpsilon(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
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
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:71
void computeCrackIntegral(std::size_t mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double &crack_volume) override
void setIntegrationPoint(unsigned integration_point)
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< 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