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