OGS 6.2.0-97-g4a610c866
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  double free_energy_density = 0;
52 
54  std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
55  DisplacementDim>::MaterialStateVariables>
57 
59  typename ShapeMatricesType::NodalRowVectorType N;
60  typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
61 
63  {
64  eps_prev = eps;
65  sigma_prev = sigma;
66  material_state_variables->pushBackState();
67  }
68 
70 };
71 
74 template <typename ShapeMatrixType>
76 {
77  std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
78 };
79 
80 template <typename ShapeFunction, typename IntegrationMethod,
81  int DisplacementDim>
83  : public SmallDeformationLocalAssemblerInterface<DisplacementDim>
84 {
85 public:
86  using ShapeMatricesType =
92 
98 
102 
104  delete;
106 
108  MeshLib::Element const& e,
109  std::size_t const /*local_matrix_size*/,
110  bool const is_axially_symmetric,
111  unsigned const integration_order,
113  : _process_data(process_data),
114  _integration_method(integration_order),
115  _element(e),
116  _is_axially_symmetric(is_axially_symmetric)
117  {
118  unsigned const n_integration_points =
119  _integration_method.getNumberOfPoints();
120 
121  _ip_data.reserve(n_integration_points);
122  _secondary_data.N.resize(n_integration_points);
123 
124  auto const shape_matrices =
125  initShapeMatrices<ShapeFunction, ShapeMatricesType,
126  IntegrationMethod, DisplacementDim>(
127  e, is_axially_symmetric, _integration_method);
128 
129  auto& solid_material =
131  _process_data.solid_materials,
132  _process_data.material_ids,
133  e.getID());
134 
135  for (unsigned ip = 0; ip < n_integration_points; ip++)
136  {
137  _ip_data.emplace_back(solid_material);
138  auto& ip_data = _ip_data[ip];
139  auto const& sm = shape_matrices[ip];
140  _ip_data[ip].integration_weight =
141  _integration_method.getWeightedPoint(ip).getWeight() *
142  sm.integralMeasure * sm.detJ;
143 
144  ip_data.N = sm.N;
145  ip_data.dNdx = sm.dNdx;
146 
147  static const int kelvin_vector_size =
149  DisplacementDim>::value;
150  // Initialize current time step values
151  ip_data.sigma.setZero(kelvin_vector_size);
152  ip_data.eps.setZero(kelvin_vector_size);
153 
154  // Previous time step values are not initialized and are set later.
155  ip_data.sigma_prev.resize(kelvin_vector_size);
156  ip_data.eps_prev.resize(kelvin_vector_size);
157 
158  _secondary_data.N[ip] = shape_matrices[ip].N;
159  }
160  }
161 
163  std::size_t setIPDataInitialConditions(std::string const& name,
164  double const* values,
165  int const integration_order) override
166  {
167  if (integration_order !=
168  static_cast<int>(_integration_method.getIntegrationOrder()))
169  {
170  OGS_FATAL(
171  "Setting integration point initial conditions; The integration "
172  "order of the local assembler for element %d is different from "
173  "the integration order in the initial condition.",
174  _element.getID());
175  }
176 
177  if (name == "sigma_ip")
178  {
179  return setSigma(values);
180  }
181 
182  return 0;
183  }
184 
185  void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
186  std::vector<double>& /*local_M_data*/,
187  std::vector<double>& /*local_K_data*/,
188  std::vector<double>& /*local_b_data*/) override
189  {
190  OGS_FATAL(
191  "SmallDeformationLocalAssembler: assembly without jacobian is not "
192  "implemented.");
193  }
194 
195  void assembleWithJacobian(double const t,
196  std::vector<double> const& local_x,
197  std::vector<double> const& /*local_xdot*/,
198  const double /*dxdot_dx*/, const double /*dx_dx*/,
199  std::vector<double>& /*local_M_data*/,
200  std::vector<double>& /*local_K_data*/,
201  std::vector<double>& local_b_data,
202  std::vector<double>& local_Jac_data) override
203  {
204  auto const local_matrix_size = local_x.size();
205 
206  auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
207  local_Jac_data, local_matrix_size, local_matrix_size);
208 
209  auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
210  local_b_data, local_matrix_size);
211 
212  unsigned const n_integration_points =
213  _integration_method.getNumberOfPoints();
214 
216  x_position.setElementID(_element.getID());
217 
218  for (unsigned ip = 0; ip < n_integration_points; ip++)
219  {
220  x_position.setIntegrationPoint(ip);
221  auto const& w = _ip_data[ip].integration_weight;
222  auto const& N = _ip_data[ip].N;
223  auto const& dNdx = _ip_data[ip].dNdx;
224 
225  typename ShapeMatricesType::template MatrixType<DisplacementDim,
226  displacement_size>
227  N_u_op = ShapeMatricesType::template MatrixType<
228  DisplacementDim,
229  displacement_size>::Zero(DisplacementDim,
230  displacement_size);
231  for (int i = 0; i < DisplacementDim; ++i)
232  {
233  N_u_op
234  .template block<1, displacement_size / DisplacementDim>(
235  i, i * displacement_size / DisplacementDim)
236  .noalias() = N;
237  }
238 
239  auto const x_coord =
240  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
241  _element, N);
242  auto const B = LinearBMatrix::computeBMatrix<
243  DisplacementDim, ShapeFunction::NPOINTS,
244  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
245  _is_axially_symmetric);
246 
247  auto const& eps_prev = _ip_data[ip].eps_prev;
248  auto const& sigma_prev = _ip_data[ip].sigma_prev;
249 
250  auto& eps = _ip_data[ip].eps;
251  auto& sigma = _ip_data[ip].sigma;
252  auto& state = _ip_data[ip].material_state_variables;
253 
254  eps.noalias() =
255  B *
256  Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
257  local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
258 
259  auto&& solution = _ip_data[ip].solid_material.integrateStress(
260  t, x_position, _process_data.dt, eps_prev, eps, sigma_prev,
261  *state, _process_data.reference_temperature);
262 
263  if (!solution)
264  {
265  OGS_FATAL("Computation of local constitutive relation failed.");
266  }
267 
269  std::tie(sigma, state, C) = std::move(*solution);
270 
271  auto const rho = _process_data.solid_density(t, x_position)[0];
272  auto const& b = _process_data.specific_body_force;
273  local_b.noalias() -=
274  (B.transpose() * sigma - N_u_op.transpose() * rho * b) * w;
275  local_Jac.noalias() += B.transpose() * C * B * w;
276  }
277  }
278 
279  void preTimestepConcrete(std::vector<double> const& /*local_x*/,
280  double const /*t*/,
281  double const /*delta_t*/) override
282  {
283  unsigned const n_integration_points =
284  _integration_method.getNumberOfPoints();
285 
286  for (unsigned ip = 0; ip < n_integration_points; ip++)
287  {
288  _ip_data[ip].pushBackState();
289  }
290  }
291 
292  void postTimestepConcrete(std::vector<double> const& /*local_x*/) override
293  {
294  unsigned const n_integration_points =
295  _integration_method.getNumberOfPoints();
296 
298  x_position.setElementID(_element.getID());
299 
300  for (unsigned ip = 0; ip < n_integration_points; ip++)
301  {
302  x_position.setIntegrationPoint(ip);
303 
304  auto& ip_data = _ip_data[ip];
305 
306  // Update free energy density needed for material forces.
307  ip_data.free_energy_density =
308  ip_data.solid_material.computeFreeEnergyDensity(
309  _process_data.t, x_position, _process_data.dt, ip_data.eps,
310  ip_data.sigma, *ip_data.material_state_variables);
311  }
312  }
313 
314  std::vector<double> const& getMaterialForces(
315  std::vector<double> const& local_x,
316  std::vector<double>& nodal_values) override
317  {
319  DisplacementDim, ShapeFunction, ShapeMatricesType,
322  GradientMatrixType>(local_x, nodal_values, _integration_method,
323  _ip_data, _element, _is_axially_symmetric);
324  }
325 
326  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
327  const unsigned integration_point) const override
328  {
329  auto const& N = _secondary_data.N[integration_point];
330 
331  // assumes N is stored contiguously in memory
332  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
333  }
334 
335  std::vector<double> const& getIntPtFreeEnergyDensity(
336  const double /*t*/,
337  GlobalVector const& /*current_solution*/,
338  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
339  std::vector<double>& cache) const override
340  {
341  cache.clear();
342  cache.reserve(_ip_data.size());
343 
344  for (auto const& ip_data : _ip_data)
345  {
346  cache.push_back(ip_data.free_energy_density);
347  }
348 
349  return cache;
350  }
351 
352  std::size_t setSigma(double const* values)
353  {
354  auto const kelvin_vector_size =
356  DisplacementDim>::value;
357  auto const n_integration_points = _ip_data.size();
358 
359  std::vector<double> ip_sigma_values;
360  auto sigma_values =
361  Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
362  Eigen::ColMajor> const>(
363  values, kelvin_vector_size, n_integration_points);
364 
365  for (unsigned ip = 0; ip < n_integration_points; ++ip)
366  {
367  _ip_data[ip].sigma =
369  sigma_values.col(ip));
370  }
371 
372  return n_integration_points;
373  }
374 
375  // TODO (naumov) This method is same as getIntPtSigma but for arguments and
376  // the ordering of the cache_mat.
377  // There should be only one.
378  std::vector<double> getSigma() const override
379  {
380  auto const kelvin_vector_size =
382  DisplacementDim>::value;
383  auto const n_integration_points = _ip_data.size();
384 
385  std::vector<double> ip_sigma_values;
386  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
387  double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
388  ip_sigma_values, n_integration_points, kelvin_vector_size);
389 
390  for (unsigned ip = 0; ip < n_integration_points; ++ip)
391  {
392  auto const& sigma = _ip_data[ip].sigma;
393  cache_mat.row(ip) =
395  }
396 
397  return ip_sigma_values;
398  }
399 
400  std::vector<double> const& getIntPtSigma(
401  const double /*t*/,
402  GlobalVector const& /*current_solution*/,
403  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
404  std::vector<double>& cache) const override
405  {
406  static const int kelvin_vector_size =
408  DisplacementDim>::value;
409  auto const num_intpts = _ip_data.size();
410 
411  cache.clear();
412  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
413  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
414  cache, kelvin_vector_size, num_intpts);
415 
416  for (unsigned ip = 0; ip < num_intpts; ++ip)
417  {
418  auto const& sigma = _ip_data[ip].sigma;
419  cache_mat.col(ip) =
421  }
422 
423  return cache;
424  }
425 
426  std::vector<double> const& getIntPtEpsilon(
427  const double /*t*/,
428  GlobalVector const& /*current_solution*/,
429  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
430  std::vector<double>& cache) const override
431  {
432  auto const kelvin_vector_size =
434  DisplacementDim>::value;
435  auto const num_intpts = _ip_data.size();
436 
437  cache.clear();
438  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
439  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
440  cache, kelvin_vector_size, num_intpts);
441 
442  for (unsigned ip = 0; ip < num_intpts; ++ip)
443  {
444  auto const& eps = _ip_data[ip].eps;
445  cache_mat.col(ip) =
447  }
448 
449  return cache;
450  }
451 
452  unsigned getNumberOfIntegrationPoints() const override
453  {
454  return _integration_method.getNumberOfPoints();
455  }
456 
458  DisplacementDim>::MaterialStateVariables const&
459  getMaterialStateVariablesAt(unsigned integration_point) const override
460  {
461  return *_ip_data[integration_point].material_state_variables;
462  }
463 
464 private:
466 
467  std::vector<
469  Eigen::aligned_allocator<IntegrationPointData<
470  BMatricesType, ShapeMatricesType, DisplacementDim>>>
472 
473  IntegrationMethod _integration_method;
477 
478  static const int displacement_size =
479  ShapeFunction::NPOINTS * DisplacementDim;
480 };
481 
482 } // namespace SmallDeformation
483 } // 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)
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
void preTimestepConcrete(std::vector< double > const &, double const, double const) override
#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