ExaDG
Loading...
Searching...
No Matches
continuum_mechanics.h
1/* ______________________________________________________________________
2 *
3 * ExaDG - High-Order Discontinuous Galerkin for the Exa-Scale
4 *
5 * Copyright (C) 2021 by the ExaDG authors
6 *
7 * This program is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program. If not, see <https://www.gnu.org/licenses/>.
19 * ______________________________________________________________________
20 */
21
22#ifndef EXADG_STRUCTURE_SPATIAL_DISCRETIZATION_OPERATORS_CONTINUUM_MECHANICS_H_
23#define EXADG_STRUCTURE_SPATIAL_DISCRETIZATION_OPERATORS_CONTINUUM_MECHANICS_H_
24
25// deal.II
26#include <deal.II/base/tensor.h>
27#include <deal.II/physics/transformations.h>
28
29namespace ExaDG
30{
31namespace Structure
32{
33template<int dim, typename Number, typename TensorType>
34inline DEAL_II_ALWAYS_INLINE //
35 TensorType
36 add_identity(TensorType tensor)
37{
38 for(unsigned int i = 0; i < dim; i++)
39 {
40 tensor[i][i] = tensor[i][i] + 1.0;
41 }
42
43 return tensor;
44}
45
46template<int dim, typename Number, typename TensorType>
47inline DEAL_II_ALWAYS_INLINE //
48 TensorType
49 subtract_identity(TensorType tensor)
50{
51 for(unsigned int i = 0; i < dim; i++)
52 {
53 tensor[i][i] = tensor[i][i] - 1.0;
54 }
55
56 return tensor;
57}
58
59// Create a symmetric tensor from a tensor H plus H^T.
60// Note that we have (H + H^T)^T = H^T + (H^T)^T = H^T + H = H + H^T.
61template<int dim, typename Number>
62inline DEAL_II_ALWAYS_INLINE //
63 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>>
64 compute_H_plus_HT(dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> const & H)
65{
66 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>> result;
67
68 for(unsigned int i = 0; i < dim; ++i)
69 {
70 for(unsigned int j = i; j < dim; ++j)
71 {
72 result[i][j] = H[i][j] + H[j][i];
73 }
74 }
75
76 // Debug output.
77 if constexpr(false)
78 {
79 if constexpr(std::is_same_v<Number, double>)
80 {
81 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> check = H + transpose(H);
82 dealii::VectorizedArray<Number> sum_err = 0.0;
83 for(unsigned int i = 0; i < dim; ++i)
84 {
85 for(unsigned int j = 0; j < dim; ++j)
86 {
87 sum_err += std::abs(check[i][j] - result[i][j]);
88 }
89 }
90 std::cout << "compute_H_plus_HT : sum_err = " << sum_err << "\n";
91 AssertThrow(sum_err.sum() < 1e-18,
92 dealii::ExcMessage("Check in `compute_H_plus_HT()` failed."));
93 }
94 }
95
96 return result;
97}
98
99// Create a symmetric tensor from a tensor H^T times H.
100// Note that we have (H^T * H)^T = H^T * (H^T)^T = H^T * H.
101template<int dim, typename Number>
102inline DEAL_II_ALWAYS_INLINE //
103 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>>
104 compute_HT_times_H(dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> const & H)
105{
106 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>> result;
107
108 for(unsigned int i = 0; i < dim; ++i)
109 {
110 for(unsigned int j = i; j < dim; ++j)
111 {
112 for(unsigned int k = 0; k < dim; ++k)
113 {
114 result[i][j] += /* HT[i][k] */ H[k][i] * H[k][j];
115 }
116 }
117 }
118
119 // Debug output.
120 if constexpr(false)
121 {
122 if constexpr(std::is_same_v<Number, double>)
123 {
124 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> check = transpose(H) * H;
125 dealii::VectorizedArray<Number> sum_err = 0.0;
126 for(unsigned int i = 0; i < dim; ++i)
127 {
128 for(unsigned int j = 0; j < dim; ++j)
129 {
130 sum_err += std::abs(check[i][j] - result[i][j]);
131 }
132 }
133 std::cout << "compute_HT_times_H : sum_err = " << sum_err << "\n";
134 AssertThrow(sum_err.sum() < 1e-18,
135 dealii::ExcMessage("Check in `compute_HT_times_H()` failed."));
136 }
137 }
138
139 return result;
140}
141
142// Create a symmetric tensor from a tensor H times H^T.
143// Note that we have (H * H^T)^T = (H^T)^T * H^T = H * H^T.
144template<int dim, typename Number>
145inline DEAL_II_ALWAYS_INLINE //
146 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>>
147 compute_H_times_HT(dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> const & H)
148{
149 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>> result;
150
151 for(unsigned int i = 0; i < dim; ++i)
152 {
153 for(unsigned int j = i; j < dim; ++j)
154 {
155 for(unsigned int k = 0; k < dim; ++k)
156 {
157 result[i][j] += H[i][k] * H[j][k] /* HT[k][j] */;
158 }
159 }
160 }
161
162 // Debug output.
163 if constexpr(false)
164 {
165 if constexpr(std::is_same_v<Number, double>)
166 {
167 dealii::VectorizedArray<Number> sum_err = 0.0;
168
169 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> check = H * transpose(H);
170 for(unsigned int i = 0; i < dim; ++i)
171 {
172 for(unsigned int j = 0; j < dim; ++j)
173 {
174 sum_err += std::abs(check[i][j] - result[i][j]);
175 }
176 }
177 std::cout << "compute_H_times_HT : sum_err = " << sum_err << "\n";
178 AssertThrow(sum_err.sum() < 1e-18,
179 dealii::ExcMessage("Check in `compute_H_times_HT()` failed."));
180 }
181 }
182
183 return result;
184}
185
186// Compute the push-forward of a symmetric tensor `S`, i.e.,
187// tau = (1/J) * F * S * F^T,
188// with `S` being symmetric, and therefore symmetric `tau`.
189template<int dim, typename Number>
190inline DEAL_II_ALWAYS_INLINE //
191 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>>
192 compute_push_forward(dealii::VectorizedArray<Number> const & J,
193 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>> const & S,
194 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> const & F)
195{
196 // Compute the non-symmetric S * F^T, since we do not want to recompute
197 // the components for the triple matrix product.
198 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> S_times_FT;
199 for(unsigned int i = 0; i < dim; ++i)
200 {
201 for(unsigned int j = 0; j < dim; ++j)
202 {
203 for(unsigned int k = 0; k < dim; ++k)
204 {
205 S_times_FT[i][j] += S[i][k] * F[j][k] /* FT[k][j] */;
206 }
207 }
208 }
209
210 // Compute the remaining product, exploit symmetry of the result.
211 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>> result;
212 for(unsigned int i = 0; i < dim; ++i)
213 {
214 for(unsigned int j = i; j < dim; ++j)
215 {
216 for(unsigned int k = 0; k < dim; ++k)
217 {
218 result[i][j] += F[i][k] * S_times_FT[k][j];
219 }
220 }
221 }
222
223 // Perform scaling with 1/J.
224 result /= J;
225
226 // Debug output.
227 if constexpr(false)
228 {
229 if constexpr(std::is_same_v<Number, double>)
230 {
231 dealii::VectorizedArray<Number> sum_err = 0.0;
232
233 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> check =
234 (1.0 / J) * F * S * transpose(F);
235 for(unsigned int i = 0; i < dim; ++i)
236 {
237 for(unsigned int j = 0; j < dim; ++j)
238 {
239 sum_err += std::abs(check[i][j] - result[i][j]);
240 }
241 }
242 std::cout << "compute_push_forward : sum_err = " << sum_err << "\n";
243 AssertThrow(sum_err.sum() < 1e-18,
244 dealii::ExcMessage("Check in `compute_push_forward()` failed."));
245 }
246 }
247
248 return result;
249}
250
251template<int dim, typename Number>
252inline DEAL_II_ALWAYS_INLINE //
253 dealii::Tensor<2, dim, dealii::VectorizedArray<Number>>
254 get_F(const dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> & gradient_displacement)
255{
256 return add_identity<dim, Number>(gradient_displacement);
257}
258
259template<int dim, typename Number>
260inline DEAL_II_ALWAYS_INLINE //
261 dealii::SymmetricTensor<2, dim, dealii::VectorizedArray<Number>>
262 get_E(const dealii::Tensor<2, dim, dealii::VectorizedArray<Number>> & gradient_displacement)
263{
264 // E = 0.5 (F^T * F) = 0.5 * (H + H^T + H^T * H),
265 // where F is the deformation gradient and H is `gradient_displacement`.
266 // Note that this version has also improved numerical stability.
267 return (0.5 *
268 (compute_H_plus_HT(gradient_displacement) + compute_HT_times_H(gradient_displacement)));
269}
270
271} // namespace Structure
272} // namespace ExaDG
273
274#endif /* EXADG_STRUCTURE_SPATIAL_DISCRETIZATION_OPERATORS_CONTINUUM_MECHANICS_H_ */
ExaDG
Definition driver.cpp:33
Generated by doxygen 1.13.2