OGS
SmallDeformationFEM.h
Go to the documentation of this file.
1 
11 #pragma once
12 
13 #include <limits>
14 #include <memory>
15 #include <vector>
16 
25 #include "ParameterLib/Parameter.h"
33 
34 namespace ProcessLib
35 {
36 namespace SmallDeformation
37 {
38 namespace MPL = MaterialPropertyLib;
39 
40 template <typename BMatricesType, typename ShapeMatricesType,
41  int DisplacementDim>
43 {
49  solid_material.createMaterialStateVariables())
50  {
51  }
52 
55  double free_energy_density = 0;
56 
58  std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
59  DisplacementDim>::MaterialStateVariables>
61 
63  typename ShapeMatricesType::NodalRowVectorType N;
64  typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
65 
67  {
68  eps_prev = eps;
69  sigma_prev = sigma;
70  material_state_variables->pushBackState();
71  }
72 
74 };
75 
78 template <typename ShapeMatrixType>
80 {
81  std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
82 };
83 
84 template <typename ShapeFunction, typename IntegrationMethod,
85  int DisplacementDim>
87  : public SmallDeformationLocalAssemblerInterface<DisplacementDim>
88 {
89 public:
96 
102 
106  using IpData =
108 
110  delete;
112 
114  MeshLib::Element const& e,
115  std::size_t const /*local_matrix_size*/,
116  bool const is_axially_symmetric,
117  unsigned const integration_order,
119  : _process_data(process_data),
120  _integration_method(integration_order),
121  _element(e),
122  _is_axially_symmetric(is_axially_symmetric)
123  {
124  unsigned const n_integration_points =
125  _integration_method.getNumberOfPoints();
126 
127  _ip_data.reserve(n_integration_points);
128  _secondary_data.N.resize(n_integration_points);
129 
130  auto const shape_matrices =
132  DisplacementDim>(e, is_axially_symmetric,
134 
135  auto& solid_material =
137  _process_data.solid_materials,
138  _process_data.material_ids,
139  e.getID());
140 
141  for (unsigned ip = 0; ip < n_integration_points; ip++)
142  {
143  _ip_data.emplace_back(solid_material);
144  auto& ip_data = _ip_data[ip];
145  auto const& sm = shape_matrices[ip];
146  _ip_data[ip].integration_weight =
147  _integration_method.getWeightedPoint(ip).getWeight() *
148  sm.integralMeasure * sm.detJ;
149 
150  ip_data.N = sm.N;
151  ip_data.dNdx = sm.dNdx;
152 
153  static const int kelvin_vector_size =
155  DisplacementDim);
156  // Initialize current time step values
157  ip_data.sigma.setZero(kelvin_vector_size);
158  ip_data.eps.setZero(kelvin_vector_size);
159 
160  // Previous time step values are not initialized and are set later.
161  ip_data.sigma_prev.resize(kelvin_vector_size);
162  ip_data.eps_prev.resize(kelvin_vector_size);
163 
164  _secondary_data.N[ip] = shape_matrices[ip].N;
165  }
166  }
167 
169  std::size_t setIPDataInitialConditions(std::string const& name,
170  double const* values,
171  int const integration_order) override
172  {
173  if (integration_order !=
174  static_cast<int>(_integration_method.getIntegrationOrder()))
175  {
176  OGS_FATAL(
177  "Setting integration point initial conditions; The integration "
178  "order of the local assembler for element {:d} is different "
179  "from the integration order in the initial condition.",
180  _element.getID());
181  }
182 
183  if (name == "sigma_ip")
184  {
185  if (_process_data.initial_stress != nullptr)
186  {
187  OGS_FATAL(
188  "Setting initial conditions for stress from integration "
189  "point data and from a parameter '{:s}' is not possible "
190  "simultaneously.",
191  _process_data.initial_stress->name);
192  }
193  return setSigma(values);
194  }
195 
196  return 0;
197  }
198 
199  void initializeConcrete() override
200  {
201  unsigned const n_integration_points =
202  _integration_method.getNumberOfPoints();
203  for (unsigned ip = 0; ip < n_integration_points; ip++)
204  {
205  auto& ip_data = _ip_data[ip];
206 
208  if (_process_data.initial_stress != nullptr)
209  {
210  ParameterLib::SpatialPosition const x_position{
211  std::nullopt, _element.getID(), ip,
213  NumLib::interpolateCoordinates<ShapeFunction,
215  _element, ip_data.N))};
216 
217  ip_data.sigma =
219  DisplacementDim>((*_process_data.initial_stress)(
220  std::numeric_limits<
221  double>::quiet_NaN() /* time independent */,
222  x_position));
223  }
224 
225  ip_data.pushBackState();
226  }
227  }
228 
229  void assemble(double const /*t*/, double const /*dt*/,
230  std::vector<double> const& /*local_x*/,
231  std::vector<double> const& /*local_xdot*/,
232  std::vector<double>& /*local_M_data*/,
233  std::vector<double>& /*local_K_data*/,
234  std::vector<double>& /*local_b_data*/) override
235  {
236  OGS_FATAL(
237  "SmallDeformationLocalAssembler: assembly without jacobian is not "
238  "implemented.");
239  }
240 
241  void assembleWithJacobian(double const t, double const dt,
242  std::vector<double> const& local_x,
243  std::vector<double> const& /*local_xdot*/,
244  const double /*dxdot_dx*/, const double /*dx_dx*/,
245  std::vector<double>& /*local_M_data*/,
246  std::vector<double>& /*local_K_data*/,
247  std::vector<double>& local_b_data,
248  std::vector<double>& local_Jac_data) override
249  {
250  auto const local_matrix_size = local_x.size();
251 
252  auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
253  local_Jac_data, local_matrix_size, local_matrix_size);
254 
255  auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
256  local_b_data, local_matrix_size);
257 
258  unsigned const n_integration_points =
259  _integration_method.getNumberOfPoints();
260 
261  MPL::VariableArray variables_prev;
262  MPL::VariableArray variables;
264  x_position.setElementID(_element.getID());
265 
266  auto const& b = _process_data.specific_body_force;
267 
268  for (unsigned ip = 0; ip < n_integration_points; ip++)
269  {
270  x_position.setIntegrationPoint(ip);
271  auto const& w = _ip_data[ip].integration_weight;
272  auto const& N = _ip_data[ip].N;
273  auto const& dNdx = _ip_data[ip].dNdx;
274 
275  typename ShapeMatricesType::template MatrixType<DisplacementDim,
277  N_u_op = ShapeMatricesType::template MatrixType<
278  DisplacementDim,
279  displacement_size>::Zero(DisplacementDim,
281  for (int i = 0; i < DisplacementDim; ++i)
282  {
283  N_u_op
284  .template block<1, displacement_size / DisplacementDim>(
285  i, i * displacement_size / DisplacementDim)
286  .noalias() = N;
287  }
288 
289  auto const x_coord =
290  NumLib::interpolateXCoordinate<ShapeFunction,
292  auto const B = LinearBMatrix::computeBMatrix<
293  DisplacementDim, ShapeFunction::NPOINTS,
294  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
296 
297  auto const& eps_prev = _ip_data[ip].eps_prev;
298  auto const& sigma_prev = _ip_data[ip].sigma_prev;
299 
300  auto& eps = _ip_data[ip].eps;
301  auto& sigma = _ip_data[ip].sigma;
302  auto& state = _ip_data[ip].material_state_variables;
303 
304  eps.noalias() =
305  B *
306  Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
307  local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
308 
309  variables_prev[static_cast<int>(MPL::Variable::stress)]
310  .emplace<
312  sigma_prev);
313  variables_prev[static_cast<int>(MPL::Variable::mechanical_strain)]
314  .emplace<
316  eps_prev);
317 
318  double const T_ref =
319  _process_data.reference_temperature
320  ? (*_process_data.reference_temperature)(t, x_position)[0]
321  : std::numeric_limits<double>::quiet_NaN();
322 
323  variables_prev[static_cast<int>(MPL::Variable::temperature)]
324  .emplace<double>(T_ref);
325  variables[static_cast<int>(MPL::Variable::mechanical_strain)]
326  .emplace<
328  eps);
329  variables[static_cast<int>(MPL::Variable::temperature)]
330  .emplace<double>(T_ref);
331 
332  auto&& solution = _ip_data[ip].solid_material.integrateStress(
333  variables_prev, variables, t, x_position, dt, *state);
334 
335  if (!solution)
336  {
337  OGS_FATAL("Computation of local constitutive relation failed.");
338  }
339 
341  std::tie(sigma, state, C) = std::move(*solution);
342 
343  auto const rho = _process_data.solid_density(t, x_position)[0];
344  local_b.noalias() -=
345  (B.transpose() * sigma - N_u_op.transpose() * rho * b) * w;
346  local_Jac.noalias() += B.transpose() * C * B * w;
347  }
348  }
349 
350  void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
351  double const t, double const dt) override
352  {
353  unsigned const n_integration_points =
354  _integration_method.getNumberOfPoints();
355 
356  for (unsigned ip = 0; ip < n_integration_points; ip++)
357  {
358  _ip_data[ip].pushBackState();
359  }
360 
362  x_position.setElementID(_element.getID());
363 
364  for (unsigned ip = 0; ip < n_integration_points; ip++)
365  {
366  x_position.setIntegrationPoint(ip);
367 
368  auto& ip_data = _ip_data[ip];
369 
370  // Update free energy density needed for material forces.
371  ip_data.free_energy_density =
372  ip_data.solid_material.computeFreeEnergyDensity(
373  t, x_position, dt, ip_data.eps, ip_data.sigma,
374  *ip_data.material_state_variables);
375  }
376  }
377 
378  std::vector<double> const& getMaterialForces(
379  std::vector<double> const& local_x,
380  std::vector<double>& nodal_values) override
381  {
383  DisplacementDim, ShapeFunction, ShapeMatricesType,
386  GradientMatrixType>(local_x, nodal_values, _integration_method,
388  }
389 
390  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
391  const unsigned integration_point) const override
392  {
393  auto const& N = _secondary_data.N[integration_point];
394 
395  // assumes N is stored contiguously in memory
396  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
397  }
398 
399  std::vector<double> const& getIntPtFreeEnergyDensity(
400  const double /*t*/,
401  std::vector<GlobalVector*> const& /*x*/,
402  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
403  std::vector<double>& cache) const override
404  {
405  cache.clear();
406  cache.reserve(_ip_data.size());
407 
408  transform(
409  cbegin(_ip_data), cend(_ip_data), back_inserter(cache),
410  [](auto const& ip_data) { return ip_data.free_energy_density; });
411 
412  return cache;
413  }
414 
415  std::size_t setSigma(double const* values)
416  {
417  return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
418  values, _ip_data, &IpData::sigma);
419  }
420 
421  // TODO (naumov) This method is same as getIntPtSigma but for arguments and
422  // the ordering of the cache_mat.
423  // There should be only one.
424  std::vector<double> getSigma() const override
425  {
426  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
428  }
429 
430  std::vector<double> const& getIntPtSigma(
431  const double /*t*/,
432  std::vector<GlobalVector*> const& /*x*/,
433  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
434  std::vector<double>& cache) const override
435  {
436  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
437  _ip_data, &IpData::sigma, cache);
438  }
439 
440  std::vector<double> const& getIntPtEpsilon(
441  const double /*t*/,
442  std::vector<GlobalVector*> const& /*x*/,
443  std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
444  std::vector<double>& cache) const override
445  {
446  return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
447  _ip_data, &IpData::eps, cache);
448  }
449 
450  unsigned getNumberOfIntegrationPoints() const override
451  {
452  return _integration_method.getNumberOfPoints();
453  }
454 
456  DisplacementDim>::MaterialStateVariables const&
457  getMaterialStateVariablesAt(unsigned integration_point) const override
458  {
459  return *_ip_data[integration_point].material_state_variables;
460  }
461 
463  double const /*t*/, double const /*dt*/, Eigen::VectorXd const& /*x*/,
464  Eigen::VectorXd const& /*x_dot*/) override
465  {
466  int const elem_id = _element.getID();
468  x_position.setElementID(elem_id);
469  unsigned const n_integration_points =
470  _integration_method.getNumberOfPoints();
471 
472  auto sigma_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
473  Eigen::Matrix<double, 3, 3>::Zero());
474 
475  for (unsigned ip = 0; ip < n_integration_points; ip++)
476  {
477  x_position.setIntegrationPoint(ip);
478  auto const& sigma = _ip_data[ip].sigma;
479  sigma_sum += sigma;
480  }
481 
482  Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
484  n_integration_points;
485 
486  Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(
487  sigma_avg);
488 
489  Eigen::Map<Eigen::Vector3d>(
490  &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
491  e_s.eigenvalues();
492 
493  auto eigen_vectors = e_s.eigenvectors();
494 
495  for (auto i = 0; i < 3; i++)
496  {
497  Eigen::Map<Eigen::Vector3d>(
498  &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
499  eigen_vectors.col(i);
500  }
501  }
502 
503 private:
505 
506  std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
507 
508  IntegrationMethod _integration_method;
512 
513  static const int displacement_size =
514  ShapeFunction::NPOINTS * DisplacementDim;
515 };
516 
517 } // namespace SmallDeformation
518 } // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
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
MatrixType< DisplacementDim *DisplacementDim+(DisplacementDim==2 ? 1 :0), _number_of_dof > GradientMatrixType
Definition: GMatrixPolicy.h:42
VectorType< DisplacementDim *DisplacementDim+(DisplacementDim==2 ? 1 :0)> GradientVectorType
Definition: GMatrixPolicy.h:44
std::vector< double > const & getIntPtEpsilon(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
SmallDeformationProcessData< DisplacementDim > & _process_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
typename GMatricesType::GradientMatrixType GradientMatrixType
SecondaryData< typename ShapeMatrices::ShapeType > _secondary_data
void assemble(double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
SmallDeformationLocalAssembler(SmallDeformationLocalAssembler const &)=delete
typename ShapeMatricesType::NodalVectorType NodalVectorType
typename BMatricesType::NodalForceVectorType NodalForceVectorType
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
typename GMatricesType::GradientVectorType GradientVectorType
std::size_t setIPDataInitialConditions(std::string const &name, double const *values, int const integration_order) override
Returns number of read integration points.
std::vector< double > const & getIntPtSigma(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
std::vector< double > const & getMaterialForces(std::vector< double > const &local_x, std::vector< double > &nodal_values) override
SmallDeformationLocalAssembler(MeshLib::Element const &e, std::size_t const, bool const is_axially_symmetric, unsigned const integration_order, SmallDeformationProcessData< DisplacementDim > &process_data)
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt(unsigned integration_point) const override
typename BMatricesType::StiffnessMatrixType StiffnessMatrixType
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
SmallDeformationLocalAssembler(SmallDeformationLocalAssembler &&)=delete
void computeSecondaryVariableConcrete(double const, double const, Eigen::VectorXd const &, Eigen::VectorXd const &) override
typename BMatricesType::NodalForceVectorType NodalDisplacementVectorType
void postTimestepConcrete(Eigen::VectorXd const &, double const t, double const dt) override
std::vector< double > const & getIntPtFreeEnergyDensity(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 dt, 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
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, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:56
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
Definition: KelvinVector.h:170
MathLib::TemplatePoint< double, 3 > Point3d
static const double t
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
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)
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
std::vector< double > const & getMaterialForces(std::vector< double > const &local_x, std::vector< double > &nodal_values, IntegrationMethod const &_integration_method, IPData const &_ip_data, MeshLib::Element const &element, bool const is_axially_symmetric)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
VectorType< ShapeFunction::NPOINTS > NodalVectorType
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
Definition: SecondaryData.h:28
IntegrationPointData(MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
ShapeMatricesType::GlobalDimNodalMatrixType dNdx
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N