OGS 6.1.0-1721-g6382411ad
SmallDeformationFEM.h
Go to the documentation of this file.
1 
10 #pragma once
11 
12 #include <memory>
13 #include <vector>
14 
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 
215  SpatialPosition x_position;
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  N_u_op
233  .template block<1, displacement_size / DisplacementDim>(
234  i, i * displacement_size / DisplacementDim)
235  .noalias() = N;
236 
237  auto const x_coord =
238  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
239  _element, N);
240  auto const B = LinearBMatrix::computeBMatrix<
241  DisplacementDim, ShapeFunction::NPOINTS,
242  typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
243  _is_axially_symmetric);
244 
245  auto const& eps_prev = _ip_data[ip].eps_prev;
246  auto const& sigma_prev = _ip_data[ip].sigma_prev;
247 
248  auto& eps = _ip_data[ip].eps;
249  auto& sigma = _ip_data[ip].sigma;
250  auto& state = _ip_data[ip].material_state_variables;
251 
252  eps.noalias() =
253  B *
254  Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
255  local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
256 
257  auto&& solution = _ip_data[ip].solid_material.integrateStress(
258  t, x_position, _process_data.dt, eps_prev, eps, sigma_prev,
259  *state, _process_data.reference_temperature);
260 
261  if (!solution)
262  OGS_FATAL("Computation of local constitutive relation failed.");
263 
265  std::tie(sigma, state, C) = std::move(*solution);
266 
267  auto const rho = _process_data.solid_density(t, x_position)[0];
268  auto const& b = _process_data.specific_body_force;
269  local_b.noalias() -=
270  (B.transpose() * sigma - N_u_op.transpose() * rho * b) * w;
271  local_Jac.noalias() += B.transpose() * C * B * w;
272  }
273  }
274 
275  void preTimestepConcrete(std::vector<double> const& /*local_x*/,
276  double const /*t*/,
277  double const /*delta_t*/) override
278  {
279  unsigned const n_integration_points =
280  _integration_method.getNumberOfPoints();
281 
282  for (unsigned ip = 0; ip < n_integration_points; ip++)
283  {
284  _ip_data[ip].pushBackState();
285  }
286  }
287 
288  void postTimestepConcrete(std::vector<double> const& /*local_x*/) override
289  {
290  unsigned const n_integration_points =
291  _integration_method.getNumberOfPoints();
292 
293  SpatialPosition x_position;
294  x_position.setElementID(_element.getID());
295 
296  for (unsigned ip = 0; ip < n_integration_points; ip++)
297  {
298  x_position.setIntegrationPoint(ip);
299 
300  auto& ip_data = _ip_data[ip];
301 
302  // Update free energy density needed for material forces.
303  ip_data.free_energy_density =
304  ip_data.solid_material.computeFreeEnergyDensity(
305  _process_data.t, x_position, _process_data.dt, ip_data.eps,
306  ip_data.sigma, *ip_data.material_state_variables);
307  }
308  }
309 
310  std::vector<double> const& getMaterialForces(
311  std::vector<double> const& local_x,
312  std::vector<double>& nodal_values) override
313  {
315  DisplacementDim, ShapeFunction, ShapeMatricesType,
318  GradientMatrixType>(local_x, nodal_values, _integration_method,
319  _ip_data, _element, _is_axially_symmetric);
320  }
321 
322  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
323  const unsigned integration_point) const override
324  {
325  auto const& N = _secondary_data.N[integration_point];
326 
327  // assumes N is stored contiguously in memory
328  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
329  }
330 
331  std::vector<double> const& getIntPtFreeEnergyDensity(
332  const double /*t*/,
333  GlobalVector const& /*current_solution*/,
334  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
335  std::vector<double>& cache) const override
336  {
337  cache.clear();
338  cache.reserve(_ip_data.size());
339 
340  for (auto const& ip_data : _ip_data)
341  {
342  cache.push_back(ip_data.free_energy_density);
343  }
344 
345  return cache;
346  }
347 
348  std::size_t setSigma(double const* values)
349  {
350  auto const kelvin_vector_size =
352  DisplacementDim>::value;
353  auto const n_integration_points = _ip_data.size();
354 
355  std::vector<double> ip_sigma_values;
356  auto sigma_values =
357  Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
358  Eigen::ColMajor> const>(
359  values, kelvin_vector_size, n_integration_points);
360 
361  for (unsigned ip = 0; ip < n_integration_points; ++ip)
362  {
363  _ip_data[ip].sigma =
365  sigma_values.col(ip));
366  }
367 
368  return n_integration_points;
369  }
370 
371  // TODO (naumov) This method is same as getIntPtSigma but for arguments and
372  // the ordering of the cache_mat.
373  // There should be only one.
374  std::vector<double> getSigma() const override
375  {
376  auto const kelvin_vector_size =
378  DisplacementDim>::value;
379  auto const n_integration_points = _ip_data.size();
380 
381  std::vector<double> ip_sigma_values;
382  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
383  double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
384  ip_sigma_values, n_integration_points, kelvin_vector_size);
385 
386  for (unsigned ip = 0; ip < n_integration_points; ++ip)
387  {
388  auto const& sigma = _ip_data[ip].sigma;
389  cache_mat.row(ip) =
391  }
392 
393  return ip_sigma_values;
394  }
395 
396  std::vector<double> const& getIntPtSigma(
397  const double /*t*/,
398  GlobalVector const& /*current_solution*/,
399  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
400  std::vector<double>& cache) const override
401  {
402  static const int kelvin_vector_size =
404  DisplacementDim>::value;
405  auto const num_intpts = _ip_data.size();
406 
407  cache.clear();
408  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
409  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
410  cache, kelvin_vector_size, num_intpts);
411 
412  for (unsigned ip = 0; ip < num_intpts; ++ip)
413  {
414  auto const& sigma = _ip_data[ip].sigma;
415  cache_mat.col(ip) =
417  }
418 
419  return cache;
420  }
421 
422  virtual std::vector<double> const& getIntPtEpsilon(
423  const double /*t*/,
424  GlobalVector const& /*current_solution*/,
425  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
426  std::vector<double>& cache) const override
427  {
428  auto const kelvin_vector_size =
430  DisplacementDim>::value;
431  auto const num_intpts = _ip_data.size();
432 
433  cache.clear();
434  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
435  double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
436  cache, kelvin_vector_size, num_intpts);
437 
438  for (unsigned ip = 0; ip < num_intpts; ++ip)
439  {
440  auto const& eps = _ip_data[ip].eps;
441  cache_mat.col(ip) =
443  }
444 
445  return cache;
446  }
447 
448  unsigned getNumberOfIntegrationPoints() const override
449  {
450  return _integration_method.getNumberOfPoints();
451  }
452 
454  DisplacementDim>::MaterialStateVariables const&
455  getMaterialStateVariablesAt(unsigned integration_point) const override
456  {
457  return *_ip_data[integration_point].material_state_variables;
458  }
459 
460 private:
462 
463  std::vector<
465  Eigen::aligned_allocator<IntegrationPointData<
466  BMatricesType, ShapeMatricesType, DisplacementDim>>>
468 
469  IntegrationMethod _integration_method;
473 
474  static const int displacement_size =
475  ShapeFunction::NPOINTS * DisplacementDim;
476 };
477 
478 } // namespace SmallDeformation
479 } // 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.
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 setElementID(std::size_t element_id)
void postTimestepConcrete(std::vector< double > const &) override
typename GMatricesType::GradientVectorType GradientVectorType
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
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:71
virtual std::vector< double > const & getIntPtEpsilon(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const override
void setIntegrationPoint(unsigned integration_point)
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