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  const double /*dxdot_dx*/, const double /*dx_dx*/,
451  std::vector<double>& /*local_M_data*/,
452  std::vector<double>& /*local_K_data*/,
453  std::vector<double>& local_b_data,
454  std::vector<double>& local_Jac_data) override
455  {
456  auto const local_matrix_size = local_x.size();
457 
458  auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
459  local_Jac_data, local_matrix_size, local_matrix_size);
460 
461  auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
462  local_b_data, local_matrix_size);
463 
464  unsigned const n_integration_points =
465  _integration_method.getNumberOfPoints();
466 
468  x_position.setElementID(_element.getID());
469 
470  // Non-local integration.
471  for (unsigned ip = 0; ip < n_integration_points; ip++)
472  {
473  x_position.setIntegrationPoint(ip);
474  auto const& w = _ip_data[ip].integration_weight;
475 
476  auto const& N = _ip_data[ip].N;
477  auto const& dNdx = _ip_data[ip].dNdx;
478 
479  auto const x_coord =
480  NumLib::interpolateXCoordinate<ShapeFunction,
482  auto const B = LinearBMatrix::computeBMatrix<
483  DisplacementDim, ShapeFunction::NPOINTS,
484  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
486 
487  auto& sigma = _ip_data[ip].sigma;
488  auto& C = _ip_data[ip].C;
489  double& damage = _ip_data[ip].damage;
490 
491  {
492  double nonlocal_kappa_d = 0;
493 
494  if (_ip_data[ip].active_self || _ip_data[ip].activated)
495  {
496  for (auto const& tuple : _ip_data[ip].non_local_assemblers)
497  {
498  // Get local variable for the integration point l.
499  double const kappa_d_l = tuple.ip_l_pointer->kappa_d;
500  double const a_kl_times_w_l = tuple.alpha_kl_times_w_l;
501  nonlocal_kappa_d += a_kl_times_w_l * kappa_d_l;
502  }
503  }
504 
505  auto const& ehlers_material =
507  DisplacementDim> const&>(_ip_data[ip].solid_material);
508 
509  //
510  // Overnonlocal formulation
511  //
512  // See (Di Luzio & Bazant 2005, IJSS) for details.
513  // The implementation would go here and would be for a given
514  // gamma_nonlocal:
515  //
516  // Update nonlocal damage with local damage (scaled with 1 -
517  // \gamma_{nonlocal}) for the current integration point and the
518  // nonlocal integral part.
519  // nonlocal_kappa_d = (1. - gamma_nonlocal) * kappa_d +
520  // gamma_nonlocal * nonlocal_kappa_d;
521 
522  nonlocal_kappa_d = std::max(0., nonlocal_kappa_d);
523 
524  // Update damage based on nonlocal kappa_d
525  {
526  auto const damage_properties =
527  ehlers_material.evaluatedDamageProperties(t,
528  x_position);
529  damage = calculateDamage(nonlocal_kappa_d,
530  damage_properties.alpha_d,
531  damage_properties.beta_d);
532  damage = std::max(0., damage);
533  }
534  sigma = sigma * (1. - damage);
535  }
536 
537  local_b.noalias() -= B.transpose() * sigma * w;
538  local_Jac.noalias() += B.transpose() * C * (1. - damage) * B * w;
539  }
540  }
541 
542  void initializeConcrete() override
543  {
544  unsigned const n_integration_points =
545  _integration_method.getNumberOfPoints();
546 
547  for (unsigned ip = 0; ip < n_integration_points; ip++)
548  {
549  _ip_data[ip].pushBackState();
550  }
551  }
552 
553  void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
554  double const /*t*/, double const /*dt*/) override
555  {
556  unsigned const n_integration_points =
557  _integration_method.getNumberOfPoints();
558 
559  for (unsigned ip = 0; ip < n_integration_points; ip++)
560  {
561  _ip_data[ip].pushBackState();
562  }
563  }
564 
565  void computeCrackIntegral(std::size_t mesh_item_id,
566  NumLib::LocalToGlobalIndexMap const& dof_table,
567  GlobalVector const& x,
568  double& crack_volume) override
569  {
570  auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
571  auto local_x = x.get(indices);
572 
573  auto u = Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
574  local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
575 
576  int const n_integration_points =
577  _integration_method.getNumberOfPoints();
578 
579  for (int ip = 0; ip < n_integration_points; ip++)
580  {
581  auto const& dNdx = _ip_data[ip].dNdx;
582  auto const& d = _ip_data[ip].damage;
583  auto const& w = _ip_data[ip].integration_weight;
584 
585  double const div_u =
586  Deformation::divergence<DisplacementDim,
587  ShapeFunction::NPOINTS>(u, dNdx);
588  crack_volume += div_u * d * w;
589  }
590  }
591 
592  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
593  const unsigned integration_point) const override
594  {
595  auto const& N = _secondary_data.N[integration_point];
596 
597  // assumes N is stored contiguously in memory
598  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
599  }
600 
601  std::vector<double> const& getNodalValues(
602  std::vector<double>& nodal_values) const override
603  {
604  nodal_values.clear();
605  auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
606  nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
607 
608  unsigned const n_integration_points =
609  _integration_method.getNumberOfPoints();
610 
611  for (unsigned ip = 0; ip < n_integration_points; ip++)
612  {
613  auto const& w = _ip_data[ip].integration_weight;
614 
615  auto const& N = _ip_data[ip].N;
616  auto const& dNdx = _ip_data[ip].dNdx;
617 
618  auto const x_coord =
619  NumLib::interpolateXCoordinate<ShapeFunction,
621  auto const B = LinearBMatrix::computeBMatrix<
622  DisplacementDim, ShapeFunction::NPOINTS,
623  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
625  auto& sigma = _ip_data[ip].sigma;
626 
627  local_b.noalias() += B.transpose() * sigma * w;
628  }
629 
630  return nodal_values;
631  }
632 
633  std::vector<double> const& getIntPtFreeEnergyDensity(
634  const double /*t*/,
635  std::vector<GlobalVector*> const& /*x*/,
636  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
637  std::vector<double>& cache) const override
638  {
639  cache.clear();
640  cache.reserve(_ip_data.size());
641 
642  transform(
643  cbegin(_ip_data), cend(_ip_data), back_inserter(cache),
644  [](auto const& ip_data) { return ip_data.free_energy_density; });
645 
646  return cache;
647  }
648 
649  std::vector<double> const& getIntPtEpsPV(
650  const double /*t*/,
651  std::vector<GlobalVector*> const& /*x*/,
652  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
653  std::vector<double>& cache) const override
654  {
655  cache.clear();
656  cache.reserve(_ip_data.size());
657 
658  transform(cbegin(_ip_data), cend(_ip_data), back_inserter(cache),
659  [](auto const& ip_data) { return *ip_data.eps_p_V; });
660 
661  return cache;
662  }
663 
664  std::vector<double> const& getIntPtEpsPDXX(
665  const double /*t*/,
666  std::vector<GlobalVector*> const& /*x*/,
667  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
668  std::vector<double>& cache) const override
669  {
670  cache.clear();
671  cache.reserve(_ip_data.size());
672 
673  transform(cbegin(_ip_data), cend(_ip_data), back_inserter(cache),
674  [](auto const& ip_data) { return *ip_data.eps_p_D_xx; });
675 
676  return cache;
677  }
678 
679  std::vector<double> const& getIntPtSigma(
680  const double /*t*/,
681  std::vector<GlobalVector*> const& /*x*/,
682  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
683  std::vector<double>& cache) const override
684  {
685  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
686  _ip_data, &IpData::sigma, cache);
687  }
688 
689  std::vector<double> const& getIntPtEpsilon(
690  const double /*t*/,
691  std::vector<GlobalVector*> const& /*x*/,
692  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
693  std::vector<double>& cache) const override
694  {
695  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
696  _ip_data, &IpData::eps, cache);
697  }
698 
699  std::size_t setSigma(double const* values)
700  {
701  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
702  values, _ip_data, &IpData::sigma);
703  }
704 
705  // TODO (naumov) This method is same as getIntPtSigma but for arguments and
706  // the ordering of the cache_mat.
707  // There should be only one.
708  std::vector<double> getSigma() const override
709  {
710  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
712  }
713 
714  void setKappaD(double value)
715  {
716  for (auto& ip_data : _ip_data)
717  {
718  ip_data.kappa_d = value;
719  }
720  }
721  std::vector<double> getKappaD() const override
722  {
723  unsigned const n_integration_points =
724  _integration_method.getNumberOfPoints();
725 
726  std::vector<double> result_values;
727  result_values.resize(n_integration_points);
728 
729  for (unsigned ip = 0; ip < n_integration_points; ++ip)
730  {
731  result_values[ip] = _ip_data[ip].kappa_d;
732  }
733 
734  return result_values;
735  }
736 
737  std::vector<double> const& getIntPtDamage(
738  const double /*t*/,
739  std::vector<GlobalVector*> const& /*x*/,
740  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
741  std::vector<double>& cache) const override
742  {
744  _ip_data, &IpData::damage, cache);
745  }
746 
747  unsigned getNumberOfIntegrationPoints() const override
748  {
749  return _integration_method.getNumberOfPoints();
750  }
751 
753  DisplacementDim>::MaterialStateVariables const&
754  getMaterialStateVariablesAt(int const integration_point) const override
755  {
756  return *_ip_data[integration_point].material_state_variables;
757  }
758 
759 private:
760  std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
761  std::size_t const component) const
762  {
763  cache.clear();
764  cache.reserve(_ip_data.size());
765 
766  for (auto const& ip_data : _ip_data)
767  {
768  if (component < 3)
769  { // xx, yy, zz components
770  cache.push_back(ip_data.sigma[component]);
771  }
772  else
773  { // mixed xy, yz, xz components
774  cache.push_back(ip_data.sigma[component] / std::sqrt(2));
775  }
776  }
777 
778  return cache;
779  }
780 
781  std::vector<double> const& getIntPtEpsilon(
782  std::vector<double>& cache, std::size_t const component) const
783  {
784  cache.clear();
785  cache.reserve(_ip_data.size());
786 
787  for (auto const& ip_data : _ip_data)
788  {
789  if (component < 3) // xx, yy, zz components
790  cache.push_back(ip_data.eps[component]);
791  else // mixed xy, yz, xz components
792  cache.push_back(ip_data.eps[component] / std::sqrt(2));
793  }
794 
795  return cache;
796  }
797 
799  getIPDataPtr(int const ip) override
800  {
801  return &_ip_data[ip];
802  }
803 
804 private:
806 
807  std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
808 
809  IntegrationMethod const _integration_method;
813 
814  static const int displacement_size =
815  ShapeFunction::NPOINTS * DisplacementDim;
816 };
817 
818 } // namespace SmallDeformationNonlocal
819 } // 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:26
double get(IndexType rowId) const
get entry
Definition: EigenVector.h:62
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
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
void assembleWithJacobian(double const t, double const, 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
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:108
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
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