OGS 6.2.1-97-g73d1aeda3
MaterialForces.h
Go to the documentation of this file.
1 
12 #pragma once
13 
14 #include <cassert>
15 #include <vector>
16 
18 #include "MathLib/LinAlg/LinAlg.h"
23 
24 namespace ProcessLib
25 {
26 namespace SmallDeformation
27 {
29 {
30  virtual std::vector<double> const& getMaterialForces(
31  std::vector<double> const& local_x,
32  std::vector<double>& nodal_values) = 0;
33 
34  virtual ~MaterialForcesInterface() = default;
35 };
36 
37 template <int DisplacementDim, typename ShapeFunction,
38  typename ShapeMatricesType, typename NodalForceVectorType,
39  typename NodalDisplacementVectorType, typename GradientVectorType,
40  typename GradientMatrixType, typename IPData,
41  typename IntegrationMethod>
42 std::vector<double> const& getMaterialForces(
43  std::vector<double> const& local_x, std::vector<double>& nodal_values,
44  IntegrationMethod const& _integration_method, IPData const& _ip_data,
45  MeshLib::Element const& element, bool const is_axially_symmetric)
46 {
47  unsigned const n_integration_points =
48  _integration_method.getNumberOfPoints();
49 
50  nodal_values.clear();
51  auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
52  nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
53 
54  for (unsigned ip = 0; ip < n_integration_points; ip++)
55  {
56  auto const& sigma = _ip_data[ip].sigma;
57  auto const& N = _ip_data[ip].N;
58  auto const& dNdx = _ip_data[ip].dNdx;
59 
60  auto const& psi = _ip_data[ip].free_energy_density;
61 
62  auto const x_coord =
63  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
64  element, _ip_data[ip].N);
65 
66  // For the 2D case the 33-component is needed (and the four entries
67  // of the non-symmetric matrix); In 3d there are nine entries.
68  GradientMatrixType G(
69  DisplacementDim * DisplacementDim + (DisplacementDim == 2 ? 1 : 0),
70  ShapeFunction::NPOINTS * DisplacementDim);
71  Deformation::computeGMatrix<DisplacementDim, ShapeFunction::NPOINTS>(
72  dNdx, G, is_axially_symmetric, N, x_coord);
73 
74  GradientVectorType const grad_u =
75  G * Eigen::Map<NodalForceVectorType const>(
76  local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
77 
78  GradientVectorType eshelby_stress;
79  eshelby_stress.setZero(DisplacementDim * DisplacementDim +
80  (DisplacementDim == 2 ? 1 : 0));
81 
82  if (DisplacementDim == 3)
83  {
84  eshelby_stress[0] = eshelby_stress[DisplacementDim + 1] =
85  eshelby_stress[8] = psi;
86 
87  eshelby_stress[0] -= sigma[0] * grad_u[0] +
88  sigma[3] / std::sqrt(2) * grad_u[3] +
89  sigma[5] / std::sqrt(2) * grad_u[6];
90 
91  eshelby_stress[1] -= sigma[3] / std::sqrt(2) * grad_u[0] +
92  sigma[1] * grad_u[3] +
93  sigma[4] / std::sqrt(2) * grad_u[6];
94 
95  eshelby_stress[2] -= sigma[5] / std::sqrt(2) * grad_u[0] +
96  sigma[4] / std::sqrt(2) * grad_u[3] +
97  sigma[2] * grad_u[6];
98 
99  eshelby_stress[3] -= sigma[0] * grad_u[1] +
100  sigma[3] / std::sqrt(2) * grad_u[4] +
101  sigma[5] / std::sqrt(2) * grad_u[7];
102 
103  eshelby_stress[4] -= sigma[3] / std::sqrt(2) * grad_u[1] +
104  sigma[1] * grad_u[4] +
105  sigma[4] / std::sqrt(2) * grad_u[7];
106 
107  eshelby_stress[5] -= sigma[5] / std::sqrt(2) * grad_u[1] +
108  sigma[4] / std::sqrt(2) * grad_u[4] +
109  sigma[2] * grad_u[7];
110 
111  eshelby_stress[6] -= sigma[0] * grad_u[2] +
112  sigma[3] / std::sqrt(2) * grad_u[5] +
113  sigma[5] / std::sqrt(2) * grad_u[8];
114 
115  eshelby_stress[7] -= sigma[3] / std::sqrt(2) * grad_u[2] +
116  sigma[1] * grad_u[5] +
117  sigma[4] / std::sqrt(2) * grad_u[8];
118 
119  eshelby_stress[8] -= sigma[5] / std::sqrt(2) * grad_u[2] +
120  sigma[4] / std::sqrt(2) * grad_u[5] +
121  sigma[2] * grad_u[8];
122  }
123  else if (DisplacementDim == 2)
124  {
125  eshelby_stress[0] = eshelby_stress[DisplacementDim + 1] =
126  eshelby_stress[4] = psi;
127 
128  eshelby_stress[0] -=
129  sigma[0] * grad_u[0] + sigma[3] / std::sqrt(2) * grad_u[2];
130 
131  eshelby_stress[1] -=
132  sigma[3] / std::sqrt(2) * grad_u[0] + sigma[1] * grad_u[2];
133 
134  eshelby_stress[2] -=
135  sigma[0] * grad_u[1] + sigma[3] / std::sqrt(2) * grad_u[3];
136 
137  eshelby_stress[3] -=
138  sigma[3] / std::sqrt(2) * grad_u[1] + sigma[1] * grad_u[3];
139 
140  // for axial-symmetric case the following is not zero in general
141  eshelby_stress[4] -= sigma[2] * grad_u[4];
142  }
143  else
144  {
145  OGS_FATAL(
146  "Material forces not implemented for displacement "
147  "dimension "
148  "other than 2 and 3.");
149  }
150 
151  auto const& w = _ip_data[ip].integration_weight;
152  local_b += G.transpose() * eshelby_stress * w;
153  }
154 
155  return nodal_values;
156 }
157 
158 template <typename LocalAssemblerInterface>
160  std::unique_ptr<GlobalVector>& material_forces,
161  std::vector<std::unique_ptr<LocalAssemblerInterface>> const&
162  local_assemblers,
163  NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
164  GlobalVector const& x)
165 {
166  DBUG("Compute material forces for small deformation process.");
167 
168  // Prepare local storage and fetch values.
169  MathLib::LinAlg::setLocalAccessibleVector(x);
170 
171  if (!material_forces)
172  {
173  material_forces =
175  }
176 
177  MathLib::LinAlg::set(*material_forces, 0);
178 
180  [](const std::size_t mesh_item_id,
181  LocalAssemblerInterface& local_assembler,
182  const NumLib::LocalToGlobalIndexMap& dof_table,
183  GlobalVector const& x,
184  GlobalVector& node_values) {
185  auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
186  std::vector<double> local_data;
187  auto const local_x = x.get(indices);
188 
189  local_assembler.getMaterialForces(local_x, local_data);
190 
191  assert(local_data.size() == indices.size());
192  node_values.add(indices, local_data);
193  },
194  local_assemblers, local_to_global_index_map, x, *material_forces);
195  MathLib::LinAlg::finalizeAssembly(*material_forces);
196 }
197 
198 } // namespace SmallDeformation
199 } // namespace ProcessLib
virtual std::vector< double > const & getMaterialForces(std::vector< double > const &local_x, std::vector< double > &nodal_values)=0
void finalizeAssembly(Matrix &)
static void executeDereferenced(F const &f, C const &c, Args_ &&... args)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
void writeMaterialForces(std::unique_ptr< GlobalVector > &material_forces, std::vector< std::unique_ptr< LocalAssemblerInterface >> const &local_assemblers, NumLib::LocalToGlobalIndexMap const &local_to_global_index_map, GlobalVector const &x)
#define OGS_FATAL(fmt,...)
Definition: Error.h:63