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