OGS 6.2.1-97-g73d1aeda3
SmallDeformationFEM.h
Go to the documentation of this file.
1 
10 #pragma once
11 
12 #include <memory>
13 #include <vector>
14 
21 #include "ParameterLib/Parameter.h"
28 
31 
32 namespace ProcessLib
33 {
34 namespace SmallDeformation
35 {
36 template <typename BMatricesType, typename ShapeMatricesType,
37  int DisplacementDim>
39 {
43  : solid_material(solid_material),
45  solid_material.createMaterialStateVariables())
46  {
47  }
48 
51 
54  BMatricesType::KelvinVectorType::Zero(
56  DisplacementDim>::value);
57 
58  double free_energy_density = 0;
59 
61  std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
62  DisplacementDim>::MaterialStateVariables>
64 
66  typename ShapeMatricesType::NodalRowVectorType N;
67  typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
68 
70  {
71  eps_prev = eps;
72  sigma_prev = sigma;
73  material_state_variables->pushBackState();
74  }
75 
77 };
78 
81 template <typename ShapeMatrixType>
83 {
84  std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
85 };
86 
87 template <typename ShapeFunction, typename IntegrationMethod,
88  int DisplacementDim>
90  : public SmallDeformationLocalAssemblerInterface<DisplacementDim>
91 {
92 public:
93  using ShapeMatricesType =
99 
105 
109 
111  delete;
113 
115  MeshLib::Element const& e,
116  std::size_t const /*local_matrix_size*/,
117  bool const is_axially_symmetric,
118  unsigned const integration_order,
120  : _process_data(process_data),
121  _integration_method(integration_order),
122  _element(e),
123  _is_axially_symmetric(is_axially_symmetric)
124  {
125  unsigned const n_integration_points =
126  _integration_method.getNumberOfPoints();
127 
128  _ip_data.reserve(n_integration_points);
129  _secondary_data.N.resize(n_integration_points);
130 
131  auto const shape_matrices =
132  initShapeMatrices<ShapeFunction, ShapeMatricesType,
133  IntegrationMethod, DisplacementDim>(
134  e, is_axially_symmetric, _integration_method);
135 
136  auto& solid_material =
138  _process_data.solid_materials,
139  _process_data.material_ids,
140  e.getID());
141 
143 
144  for (unsigned ip = 0; ip < n_integration_points; ip++)
145  {
146  _ip_data.emplace_back(solid_material);
147  auto& ip_data = _ip_data[ip];
148  auto const& sm = shape_matrices[ip];
149  _ip_data[ip].integration_weight =
150  _integration_method.getWeightedPoint(ip).getWeight() *
151  sm.integralMeasure * sm.detJ;
152 
153  ip_data.N = sm.N;
154  ip_data.dNdx = sm.dNdx;
155 
156  static const int kelvin_vector_size =
158  DisplacementDim>::value;
159 
160  // Initialize current time step values
161  if (_process_data.nonequilibrium_stress)
162  {
163  // Computation of non-equilibrium stress.
164  x_position.setCoordinates(MathLib::Point3d(
165  interpolateCoordinates<ShapeFunction, ShapeMatricesType>(
166  e, ip_data.N)));
167  std::vector<double> sigma_neq_data =
168  (*_process_data.nonequilibrium_stress)(
169  std::numeric_limits<
170  double>::quiet_NaN() /* time independent */,
171  x_position);
172  ip_data.sigma_neq =
173  Eigen::Map<typename BMatricesType::KelvinVectorType>(
174  sigma_neq_data.data(),
176  DisplacementDim>::value,
177  1);
178  }
179  // Initialization from non-equilibrium sigma, which is zero by
180  // default, or is set to some value.
181  ip_data.sigma = ip_data.sigma_neq;
182 
183  ip_data.eps.setZero(kelvin_vector_size);
184 
185  // Previous time step values are not initialized and are set later.
186  ip_data.sigma_prev.resize(kelvin_vector_size);
187  ip_data.eps_prev.resize(kelvin_vector_size);
188 
189  _secondary_data.N[ip] = shape_matrices[ip].N;
190  }
191  }
192 
194  std::size_t setIPDataInitialConditions(std::string const& name,
195  double const* values,
196  int const integration_order) override
197  {
198  if (integration_order !=
199  static_cast<int>(_integration_method.getIntegrationOrder()))
200  {
201  OGS_FATAL(
202  "Setting integration point initial conditions; The integration "
203  "order of the local assembler for element %d is different from "
204  "the integration order in the initial condition.",
205  _element.getID());
206  }
207 
208  if (name == "sigma_ip")
209  {
210  return setSigma(values);
211  }
212 
213  return 0;
214  }
215 
216  void initializeConcrete() override
217  {
218  unsigned const n_integration_points =
219  _integration_method.getNumberOfPoints();
220 
221  for (unsigned ip = 0; ip < n_integration_points; ip++)
222  {
223  _ip_data[ip].pushBackState();
224  }
225  }
226 
227  void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
228  std::vector<double>& /*local_M_data*/,
229  std::vector<double>& /*local_K_data*/,
230  std::vector<double>& /*local_b_data*/) override
231  {
232  OGS_FATAL(
233  "SmallDeformationLocalAssembler: assembly without jacobian is not "
234  "implemented.");
235  }
236 
237  void assembleWithJacobian(double const t,
238  std::vector<double> const& local_x,
239  std::vector<double> const& /*local_xdot*/,
240  const double /*dxdot_dx*/, const double /*dx_dx*/,
241  std::vector<double>& /*local_M_data*/,
242  std::vector<double>& /*local_K_data*/,
243  std::vector<double>& local_b_data,
244  std::vector<double>& local_Jac_data) override
245  {
246  auto const local_matrix_size = local_x.size();
247 
248  auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
249  local_Jac_data, local_matrix_size, local_matrix_size);
250 
251  auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
252  local_b_data, local_matrix_size);
253 
254  unsigned const n_integration_points =
255  _integration_method.getNumberOfPoints();
256 
258  x_position.setElementID(_element.getID());
259 
260  for (unsigned ip = 0; ip < n_integration_points; ip++)
261  {
262  x_position.setIntegrationPoint(ip);
263  auto const& w = _ip_data[ip].integration_weight;
264  auto const& N = _ip_data[ip].N;
265  auto const& dNdx = _ip_data[ip].dNdx;
266 
267  typename ShapeMatricesType::template MatrixType<DisplacementDim,
268  displacement_size>
269  N_u_op = ShapeMatricesType::template MatrixType<
270  DisplacementDim,
271  displacement_size>::Zero(DisplacementDim,
272  displacement_size);
273  for (int i = 0; i < DisplacementDim; ++i)
274  {
275  N_u_op
276  .template block<1, displacement_size / DisplacementDim>(
277  i, i * displacement_size / DisplacementDim)
278  .noalias() = N;
279  }
280 
281  auto const x_coord =
282  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
283  _element, N);
284  auto const B = LinearBMatrix::computeBMatrix<
285  DisplacementDim, ShapeFunction::NPOINTS,
286  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
287  _is_axially_symmetric);
288 
289  auto const& eps_prev = _ip_data[ip].eps_prev;
290  auto const& sigma_prev = _ip_data[ip].sigma_prev;
291  auto const& sigma_neq = _ip_data[ip].sigma_neq;
292 
293  auto& eps = _ip_data[ip].eps;
294  auto& sigma = _ip_data[ip].sigma;
295  auto& state = _ip_data[ip].material_state_variables;
296 
297  eps.noalias() =
298  B *
299  Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
300  local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
301 
302  auto&& solution = _ip_data[ip].solid_material.integrateStress(
303  t, x_position, _process_data.dt, eps_prev, eps, sigma_prev,
304  *state, _process_data.reference_temperature);
305 
306  if (!solution)
307  {
308  OGS_FATAL("Computation of local constitutive relation failed.");
309  }
310 
312  std::tie(sigma, state, C) = std::move(*solution);
313 
314  auto const rho = _process_data.solid_density(t, x_position)[0];
315  auto const& b = _process_data.specific_body_force;
316  local_b.noalias() -= (B.transpose() * (sigma - sigma_neq) -
317  N_u_op.transpose() * rho * b) *
318  w;
319  local_Jac.noalias() += B.transpose() * C * B * w;
320  }
321  }
322 
323  void postTimestepConcrete(std::vector<double> const& /*local_x*/) override
324  {
325  unsigned const n_integration_points =
326  _integration_method.getNumberOfPoints();
327 
328  for (unsigned ip = 0; ip < n_integration_points; ip++)
329  {
330  _ip_data[ip].pushBackState();
331  }
332 
334  x_position.setElementID(_element.getID());
335 
336  for (unsigned ip = 0; ip < n_integration_points; ip++)
337  {
338  x_position.setIntegrationPoint(ip);
339 
340  auto& ip_data = _ip_data[ip];
341 
342  // Update free energy density needed for material forces.
343  ip_data.free_energy_density =
344  ip_data.solid_material.computeFreeEnergyDensity(
345  _process_data.t, x_position, _process_data.dt, ip_data.eps,
346  ip_data.sigma, *ip_data.material_state_variables);
347  }
348  }
349 
350  std::vector<double> const& getMaterialForces(
351  std::vector<double> const& local_x,
352  std::vector<double>& nodal_values) override
353  {
355  DisplacementDim, ShapeFunction, ShapeMatricesType,
358  GradientMatrixType>(local_x, nodal_values, _integration_method,
359  _ip_data, _element, _is_axially_symmetric);
360  }
361 
362  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
363  const unsigned integration_point) const override
364  {
365  auto const& N = _secondary_data.N[integration_point];
366 
367  // assumes N is stored contiguously in memory
368  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
369  }
370 
371  std::vector<double> const& getIntPtFreeEnergyDensity(
372  const double /*t*/,
373  GlobalVector const& /*current_solution*/,
374  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
375  std::vector<double>& cache) const override
376  {
377  cache.clear();
378  cache.reserve(_ip_data.size());
379 
380  for (auto const& ip_data : _ip_data)
381  {
382  cache.push_back(ip_data.free_energy_density);
383  }
384 
385  return cache;
386  }
387 
388  std::size_t setSigma(double const* values)
389  {
390  auto const kelvin_vector_size =
392  DisplacementDim>::value;
393  auto const n_integration_points = _ip_data.size();
394 
395  auto sigma_values =
396  Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
397  Eigen::ColMajor> const>(
398  values, kelvin_vector_size, n_integration_points);
399 
400  for (unsigned ip = 0; ip < n_integration_points; ++ip)
401  {
402  _ip_data[ip].sigma =
404  sigma_values.col(ip));
405  }
406 
407  return n_integration_points;
408  }
409 
410  // TODO (naumov) This method is same as getIntPtSigma but for arguments and
411  // the ordering of the cache_mat.
412  // There should be only one.
413  std::vector<double> getSigma() const override
414  {
415  auto const kelvin_vector_size =
417  DisplacementDim>::value;
418  auto const n_integration_points = _ip_data.size();
419 
420  std::vector<double> ip_sigma_values;
421  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
422  double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
423  ip_sigma_values, n_integration_points, kelvin_vector_size);
424 
425  for (unsigned ip = 0; ip < n_integration_points; ++ip)
426  {
427  auto const& sigma = _ip_data[ip].sigma;
428  cache_mat.row(ip) =
430  }
431 
432  return ip_sigma_values;
433  }
434 
435  std::vector<double> const& getIntPtSigma(
436  const double /*t*/,
437  GlobalVector const& /*current_solution*/,
438  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
439  std::vector<double>& cache) const override
440  {
441  static const int kelvin_vector_size =
443  DisplacementDim>::value;
444  auto const num_intpts = _ip_data.size();
445 
446  cache.clear();
447  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
448  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
449  cache, kelvin_vector_size, num_intpts);
450 
451  for (unsigned ip = 0; ip < num_intpts; ++ip)
452  {
453  auto const& sigma = _ip_data[ip].sigma;
454  cache_mat.col(ip) =
456  }
457 
458  return cache;
459  }
460 
461  std::vector<double> const& getIntPtEpsilon(
462  const double /*t*/,
463  GlobalVector const& /*current_solution*/,
464  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
465  std::vector<double>& cache) const override
466  {
467  auto const kelvin_vector_size =
469  DisplacementDim>::value;
470  auto const num_intpts = _ip_data.size();
471 
472  cache.clear();
473  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
474  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
475  cache, kelvin_vector_size, num_intpts);
476 
477  for (unsigned ip = 0; ip < num_intpts; ++ip)
478  {
479  auto const& eps = _ip_data[ip].eps;
480  cache_mat.col(ip) =
482  }
483 
484  return cache;
485  }
486 
487  unsigned getNumberOfIntegrationPoints() const override
488  {
489  return _integration_method.getNumberOfPoints();
490  }
491 
493  DisplacementDim>::MaterialStateVariables const&
494  getMaterialStateVariablesAt(unsigned integration_point) const override
495  {
496  return *_ip_data[integration_point].material_state_variables;
497  }
498 
499 private:
501 
502  std::vector<
504  Eigen::aligned_allocator<IntegrationPointData<
505  BMatricesType, ShapeMatricesType, DisplacementDim>>>
507 
508  IntegrationMethod _integration_method;
512 
513  static const int displacement_size =
514  ShapeFunction::NPOINTS * DisplacementDim;
515 };
516 
517 } // namespace SmallDeformation
518 } // namespace ProcessLib
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
Definition: KelvinVector.h:173
typename GMatricesType::GradientMatrixType GradientMatrixType
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
void setElementID(std::size_t element_id)
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
void assemble(double const, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
MatrixType< _number_of_dof, _number_of_dof > StiffnessMatrixType
Definition: BMatrixPolicy.h:40
typename BMatricesType::NodalForceVectorType NodalForceVectorType
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)
BMatricesType::KelvinVectorType sigma_neq
Non-equilibrium initial stress.
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)
std::vector< double > const & getIntPtFreeEnergyDensity(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, KelvinVectorDimensions< DisplacementDim >::value, Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:59
void postTimestepConcrete(std::vector< double > const &) override
typename GMatricesType::GradientVectorType GradientVectorType
void setIntegrationPoint(unsigned integration_point)
SmallDeformationProcessData< DisplacementDim > & _process_data
typename BMatricesType::NodalForceVectorType NodalDisplacementVectorType
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:31
std::vector< double > const & getIntPtSigma(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
ShapeMatricesType::GlobalDimNodalMatrixType dNdx
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables
std::size_t setIPDataInitialConditions(std::string const &name, double const *values, int const integration_order) override
Returns number of read integration points.
SmallDeformationLocalAssembler(MeshLib::Element const &e, std::size_t const, bool const is_axially_symmetric, unsigned const integration_order, SmallDeformationProcessData< DisplacementDim > &process_data)
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)
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
VectorType< ShapeFunction::NPOINTS > NodalVectorType
std::vector< double > const & getIntPtEpsilon(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
std::vector< IntegrationPointData< BMatricesType, ShapeMatricesType, DisplacementDim >, Eigen::aligned_allocator< IntegrationPointData< BMatricesType, ShapeMatricesType, DisplacementDim > > > _ip_data
IntegrationPointData(MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
VectorType< DisplacementDim *DisplacementDim+(DisplacementDim==2 ? 1 :0)> GradientVectorType
Definition: GMatrixPolicy.h:45
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
SecondaryData< typename ShapeMatrices::ShapeType > _secondary_data
typename ShapeMatricesType::NodalVectorType NodalVectorType
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:54
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material
std::vector< double > const & getMaterialForces(std::vector< double > const &local_x, std::vector< double > &nodal_values) override
MatrixType< DisplacementDim *DisplacementDim+(DisplacementDim==2 ? 1 :0), _number_of_dof > GradientMatrixType
Definition: GMatrixPolicy.h:43
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt(unsigned integration_point) const override
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
typename BMatricesType::StiffnessMatrixType StiffnessMatrixType