OGS
MaterialForces.h
Go to the documentation of this file.
1
10#pragma once
11
12#include <cassert>
13#include <vector>
14
21
22namespace ProcessLib
23{
24namespace SmallDeformation
25{
27{
28 virtual std::vector<double> const& getMaterialForces(
29 std::vector<double> const& local_x,
30 std::vector<double>& nodal_values) = 0;
31
32 virtual ~MaterialForcesInterface() = default;
33};
34
35template <int DisplacementDim, typename ShapeFunction,
36 typename ShapeMatricesType, typename NodalForceVectorType,
37 typename NodalDisplacementVectorType, typename GradientVectorType,
38 typename GradientMatrixType, typename IPData, typename StressData,
39 typename OutputData, typename IntegrationMethod>
40std::vector<double> const& getMaterialForces(
41 std::vector<double> const& local_x, std::vector<double>& nodal_values,
42 IntegrationMethod const& _integration_method, IPData const& _ip_data,
43 StressData const& stress_data, OutputData const& output_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();
51 nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
52
53 for (unsigned ip = 0; ip < n_integration_points; ip++)
54 {
55 auto const& sigma = stress_data[ip].stress_data.sigma;
56 auto const& N = _ip_data[ip].N_u;
57 auto const& dNdx = _ip_data[ip].dNdx_u;
58
59 auto const& psi =
60 output_data[ip].free_energy_density_data.free_energy_density;
61
62 auto const x_coord =
64 element, _ip_data[ip].N_u);
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);
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
158template <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.
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 {
186 auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
187 std::vector<double> local_data;
188 auto const local_x = x.get(indices);
189
190 local_assembler.getMaterialForces(local_x, local_data);
191
192 assert(local_data.size() == indices.size());
193 node_values.add(indices, local_data);
194 },
195 local_assemblers, local_to_global_index_map, x, *material_forces);
196 MathLib::LinAlg::finalizeAssembly(*material_forces);
197}
198
199} // namespace SmallDeformation
200} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Global vector based on Eigen vector.
Definition EigenVector.h:25
void add(IndexType rowId, double v)
add entry
Definition EigenVector.h:76
double get(IndexType rowId) const
get entry
Definition EigenVector.h:58
void finalizeAssembly(PETScMatrix &A)
Definition LinAlg.cpp:198
void setLocalAccessibleVector(PETScVector const &x)
Definition LinAlg.cpp:27
void set(PETScVector &x, PetscScalar const a)
Definition LinAlg.cpp:32
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
void computeGMatrix(DNDX_Type const &dNdx, GMatrixType &g_matrix, const bool is_axially_symmetric, N_Type const &N, const double radius)
Fills a G-matrix based on given shape function dN/dx values.
Definition GMatrix.h:25
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, StressData const &stress_data, OutputData const &output_data, MeshLib::Element const &element, bool const is_axially_symmetric)
static void executeDereferenced(F const &f, C const &c, Args_ &&... args)
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > sigma
Definition StressData.h:20
virtual std::vector< double > const & getMaterialForces(std::vector< double > const &local_x, std::vector< double > &nodal_values)=0