1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 /// @file
9 /// Linear elasticity for solid mechanics example using PETSc
10
11 #include <ceed/types.h>
12 #ifndef CEED_RUNNING_JIT_PASS
13 #include <math.h>
14 #endif
15
16 #ifndef PHYSICS_STRUCT
17 #define PHYSICS_STRUCT
18 typedef struct Physics_private *Physics;
19 struct Physics_private {
20 CeedScalar nu; // Poisson's ratio
21 CeedScalar E; // Young's Modulus
22 };
23 #endif
24
25 // -----------------------------------------------------------------------------
26 // Residual evaluation for linear elasticity
27 // -----------------------------------------------------------------------------
ElasResidual_Linear(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)28 CEED_QFUNCTION(ElasResidual_Linear)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
29 // Inputs
30 const CeedScalar(*ug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];
31
32 // Outputs
33 CeedScalar(*dvdX)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[0];
34 // grad_u not used for linear elasticity
35 // (*grad_u)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[1];
36
37 // Context
38 const Physics context = (Physics)ctx;
39 const CeedScalar E = context->E;
40 const CeedScalar nu = context->nu;
41
42 // Quadrature Point Loop
43 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
44 // Read spatial derivatives of u
45 const CeedScalar du[3][3] = {
46 {ug[0][0][i], ug[1][0][i], ug[2][0][i]},
47 {ug[0][1][i], ug[1][1][i], ug[2][1][i]},
48 {ug[0][2][i], ug[1][2][i], ug[2][2][i]}
49 };
50 // -- Qdata
51 const CeedScalar wdetJ = q_data[0][i];
52 const CeedScalar dXdx[3][3] = {
53 {q_data[1][i], q_data[2][i], q_data[3][i]},
54 {q_data[4][i], q_data[5][i], q_data[6][i]},
55 {q_data[7][i], q_data[8][i], q_data[9][i]}
56 };
57
58 // Compute grad_u
59 // dXdx = (dx/dX)^(-1)
60 // Apply dXdx to du = grad_u
61 CeedScalar grad_u[3][3];
62 for (CeedInt j = 0; j < 3; j++) { // Component
63 for (CeedInt k = 0; k < 3; k++) { // Derivative
64 grad_u[j][k] = 0;
65 for (CeedInt m = 0; m < 3; m++) grad_u[j][k] += dXdx[m][k] * du[j][m];
66 }
67 }
68
69 // Compute Strain : e (epsilon)
70 // e = 1/2 (grad u + (grad u)^T)
71
72 const CeedScalar e[3][3] = {
73 {(grad_u[0][0] + grad_u[0][0]) / 2., (grad_u[0][1] + grad_u[1][0]) / 2., (grad_u[0][2] + grad_u[2][0]) / 2.},
74 {(grad_u[1][0] + grad_u[0][1]) / 2., (grad_u[1][1] + grad_u[1][1]) / 2., (grad_u[1][2] + grad_u[2][1]) / 2.},
75 {(grad_u[2][0] + grad_u[0][2]) / 2., (grad_u[2][1] + grad_u[1][2]) / 2., (grad_u[2][2] + grad_u[2][2]) / 2.}
76 };
77
78 //
79 // Formulation uses Voigt notation:
80 // stress (sigma) strain (epsilon)
81 //
82 // [sigma00] [e00]
83 // [sigma11] [e11]
84 // [sigma22] = S * [e22]
85 // [sigma12] [e12]
86 // [sigma02] [e02]
87 // [sigma01] [e01]
88 //
89 // where
90 // [1-nu nu nu ]
91 // [ nu 1-nu nu ]
92 // [ nu nu 1-nu ]
93 // S = E/((1+nu)*(1-2*nu)) [ (1-2*nu)/2 ]
94 // [ (1-2*nu)/2 ]
95 // [ (1-2*nu)/2 ]
96
97 // Above Voigt Notation is placed in a 3x3 matrix:
98 const CeedScalar ss = E / ((1 + nu) * (1 - 2 * nu));
99 const CeedScalar sigma00 = ss * ((1 - nu) * e[0][0] + nu * e[1][1] + nu * e[2][2]),
100 sigma11 = ss * (nu * e[0][0] + (1 - nu) * e[1][1] + nu * e[2][2]),
101 sigma22 = ss * (nu * e[0][0] + nu * e[1][1] + (1 - nu) * e[2][2]), sigma12 = ss * (1 - 2 * nu) * e[1][2] * 0.5,
102 sigma02 = ss * (1 - 2 * nu) * e[0][2] * 0.5, sigma01 = ss * (1 - 2 * nu) * e[0][1] * 0.5;
103 const CeedScalar sigma[3][3] = {
104 {sigma00, sigma01, sigma02},
105 {sigma01, sigma11, sigma12},
106 {sigma02, sigma12, sigma22}
107 };
108
109 // Apply dXdx^T and weight to sigma
110 for (CeedInt j = 0; j < 3; j++) { // Component
111 for (CeedInt k = 0; k < 3; k++) { // Derivative
112 dvdX[k][j][i] = 0;
113 for (CeedInt m = 0; m < 3; m++) dvdX[k][j][i] += dXdx[k][m] * sigma[j][m] * wdetJ;
114 }
115 }
116 } // End of Quadrature Point Loop
117
118 return CEED_ERROR_SUCCESS;
119 }
120
121 // -----------------------------------------------------------------------------
122 // Jacobian evaluation for linear elasticity
123 // -----------------------------------------------------------------------------
ElasJacobian_Linear(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)124 CEED_QFUNCTION(ElasJacobian_Linear)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
125 // Inputs
126 const CeedScalar(*deltaug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0],
127 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];
128 // grad_u not used for linear elasticity
129 // (*grad_u)[3][Q] = (CeedScalar(*)[3][Q])in[2];
130
131 // Outputs
132 CeedScalar(*deltadvdX)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[0];
133
134 // Context
135 const Physics context = (Physics)ctx;
136 const CeedScalar E = context->E;
137 const CeedScalar nu = context->nu;
138
139 // Quadrature Point Loop
140 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
141 // Read spatial derivatives of u
142 const CeedScalar deltadu[3][3] = {
143 {deltaug[0][0][i], deltaug[1][0][i], deltaug[2][0][i]},
144 {deltaug[0][1][i], deltaug[1][1][i], deltaug[2][1][i]},
145 {deltaug[0][2][i], deltaug[1][2][i], deltaug[2][2][i]}
146 };
147 // -- Qdata
148 const CeedScalar wdetJ = q_data[0][i];
149 const CeedScalar dXdx[3][3] = {
150 {q_data[1][i], q_data[2][i], q_data[3][i]},
151 {q_data[4][i], q_data[5][i], q_data[6][i]},
152 {q_data[7][i], q_data[8][i], q_data[9][i]}
153 };
154
155 // Compute graddeltau
156 // dXdx = (dx/dX)^(-1)
157 // Apply dXdx to deltadu = graddeltau
158 CeedScalar graddeltau[3][3];
159 for (CeedInt j = 0; j < 3; j++) { // Component
160 for (CeedInt k = 0; k < 3; k++) { // Derivative
161 graddeltau[j][k] = 0;
162 for (CeedInt m = 0; m < 3; m++) graddeltau[j][k] += dXdx[m][k] * deltadu[j][m];
163 }
164 }
165
166 // Compute Strain : e (epsilon)
167 // e = 1/2 (grad u + (grad u)^T)
168 const CeedScalar de[3][3] = {
169 {(graddeltau[0][0] + graddeltau[0][0]) / 2., (graddeltau[0][1] + graddeltau[1][0]) / 2., (graddeltau[0][2] + graddeltau[2][0]) / 2.},
170 {(graddeltau[1][0] + graddeltau[0][1]) / 2., (graddeltau[1][1] + graddeltau[1][1]) / 2., (graddeltau[1][2] + graddeltau[2][1]) / 2.},
171 {(graddeltau[2][0] + graddeltau[0][2]) / 2., (graddeltau[2][1] + graddeltau[1][2]) / 2., (graddeltau[2][2] + graddeltau[2][2]) / 2.}
172 };
173
174 // Formulation uses Voigt notation:
175 // stress (sigma) strain (epsilon)
176 //
177 // [dsigma00] [de00]
178 // [dsigma11] [de11]
179 // [dsigma22] = S * [de22]
180 // [dsigma12] [de12]
181 // [dsigma02] [de02]
182 // [dsigma01] [de01]
183 //
184 // where
185 // [1-nu nu nu ]
186 // [ nu 1-nu nu ]
187 // [ nu nu 1-nu ]
188 // S = E/((1+nu)*(1-2*nu)) [ (1-2*nu)/2 ]
189 // [ (1-2*nu)/2 ]
190 // [ (1-2*nu)/2 ]
191
192 // Above Voigt Notation is placed in a 3x3 matrix:
193 const CeedScalar ss = E / ((1 + nu) * (1 - 2 * nu));
194 const CeedScalar dsigma00 = ss * ((1 - nu) * de[0][0] + nu * de[1][1] + nu * de[2][2]),
195 dsigma11 = ss * (nu * de[0][0] + (1 - nu) * de[1][1] + nu * de[2][2]),
196 dsigma22 = ss * (nu * de[0][0] + nu * de[1][1] + (1 - nu) * de[2][2]), dsigma12 = ss * (1 - 2 * nu) * de[1][2] / 2,
197 dsigma02 = ss * (1 - 2 * nu) * de[0][2] / 2, dsigma01 = ss * (1 - 2 * nu) * de[0][1] / 2;
198 const CeedScalar dsigma[3][3] = {
199 {dsigma00, dsigma01, dsigma02},
200 {dsigma01, dsigma11, dsigma12},
201 {dsigma02, dsigma12, dsigma22}
202 };
203
204 // Apply dXdx^T and weight
205 for (CeedInt j = 0; j < 3; j++) { // Component
206 for (CeedInt k = 0; k < 3; k++) { // Derivative
207 deltadvdX[k][j][i] = 0;
208 for (CeedInt m = 0; m < 3; m++) deltadvdX[k][j][i] += dXdx[k][m] * dsigma[j][m] * wdetJ;
209 }
210 }
211 } // End of Quadrature Point Loop
212
213 return CEED_ERROR_SUCCESS;
214 }
215
216 // -----------------------------------------------------------------------------
217 // Strain energy computation for linear elasticity
218 // -----------------------------------------------------------------------------
ElasEnergy_Linear(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)219 CEED_QFUNCTION(ElasEnergy_Linear)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
220 // Inputs
221 const CeedScalar(*ug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];
222
223 // Outputs
224 CeedScalar(*energy) = (CeedScalar(*))out[0];
225
226 // Context
227 const Physics context = (Physics)ctx;
228 const CeedScalar E = context->E;
229 const CeedScalar nu = context->nu;
230
231 // Constants
232 const CeedScalar TwoMu = E / (1 + nu);
233 const CeedScalar mu = TwoMu / 2;
234 const CeedScalar Kbulk = E / (3 * (1 - 2 * nu)); // Bulk Modulus
235 const CeedScalar lambda = (3 * Kbulk - TwoMu) / 3;
236
237 // Quadrature Point Loop
238 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
239 // Read spatial derivatives of u
240 const CeedScalar du[3][3] = {
241 {ug[0][0][i], ug[1][0][i], ug[2][0][i]},
242 {ug[0][1][i], ug[1][1][i], ug[2][1][i]},
243 {ug[0][2][i], ug[1][2][i], ug[2][2][i]}
244 };
245 // -- Qdata
246 const CeedScalar wdetJ = q_data[0][i];
247 const CeedScalar dXdx[3][3] = {
248 {q_data[1][i], q_data[2][i], q_data[3][i]},
249 {q_data[4][i], q_data[5][i], q_data[6][i]},
250 {q_data[7][i], q_data[8][i], q_data[9][i]}
251 };
252
253 // Compute grad_u
254 // dXdx = (dx/dX)^(-1)
255 // Apply dXdx to du = grad_u
256 CeedScalar grad_u[3][3];
257 for (CeedInt j = 0; j < 3; j++) { // Component
258 for (CeedInt k = 0; k < 3; k++) { // Derivative
259 grad_u[j][k] = 0;
260 for (CeedInt m = 0; m < 3; m++) grad_u[j][k] += dXdx[m][k] * du[j][m];
261 }
262 }
263
264 // Compute Strain : e (epsilon)
265 // e = 1/2 (grad u + (grad u)^T)
266
267 const CeedScalar e[3][3] = {
268 {(grad_u[0][0] + grad_u[0][0]) / 2., (grad_u[0][1] + grad_u[1][0]) / 2., (grad_u[0][2] + grad_u[2][0]) / 2.},
269 {(grad_u[1][0] + grad_u[0][1]) / 2., (grad_u[1][1] + grad_u[1][1]) / 2., (grad_u[1][2] + grad_u[2][1]) / 2.},
270 {(grad_u[2][0] + grad_u[0][2]) / 2., (grad_u[2][1] + grad_u[1][2]) / 2., (grad_u[2][2] + grad_u[2][2]) / 2.}
271 };
272
273 // Strain energy
274 const CeedScalar strain_vol = e[0][0] + e[1][1] + e[2][2];
275 energy[i] =
276 (lambda * strain_vol * strain_vol / 2. + strain_vol * mu + (e[0][1] * e[0][1] + e[0][2] * e[0][2] + e[1][2] * e[1][2]) * 2 * mu) * wdetJ;
277
278 } // End of Quadrature Point Loop
279
280 return CEED_ERROR_SUCCESS;
281 }
282
283 // -----------------------------------------------------------------------------
284 // Nodal diagnostic quantities for linear elasticity
285 // -----------------------------------------------------------------------------
ElasDiagnostic_Linear(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)286 CEED_QFUNCTION(ElasDiagnostic_Linear)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
287 // Inputs
288 const CeedScalar(*u)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], (*ug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1],
289 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2];
290
291 // Outputs
292 CeedScalar(*diagnostic)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
293
294 // Context
295 const Physics context = (Physics)ctx;
296 const CeedScalar E = context->E;
297 const CeedScalar nu = context->nu;
298
299 // Constants
300 const CeedScalar TwoMu = E / (1 + nu);
301 const CeedScalar mu = TwoMu / 2;
302 const CeedScalar Kbulk = E / (3 * (1 - 2 * nu)); // Bulk Modulus
303 const CeedScalar lambda = (3 * Kbulk - TwoMu) / 3;
304
305 // Quadrature Point Loop
306 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
307 // Read spatial derivatives of u
308 const CeedScalar du[3][3] = {
309 {ug[0][0][i], ug[1][0][i], ug[2][0][i]},
310 {ug[0][1][i], ug[1][1][i], ug[2][1][i]},
311 {ug[0][2][i], ug[1][2][i], ug[2][2][i]}
312 };
313 // -- Qdata
314 const CeedScalar dXdx[3][3] = {
315 {q_data[1][i], q_data[2][i], q_data[3][i]},
316 {q_data[4][i], q_data[5][i], q_data[6][i]},
317 {q_data[7][i], q_data[8][i], q_data[9][i]}
318 };
319
320 // Compute grad_u
321 // dXdx = (dx/dX)^(-1)
322 // Apply dXdx to du = grad_u
323 CeedScalar grad_u[3][3];
324 for (CeedInt j = 0; j < 3; j++) { // Component
325 for (CeedInt k = 0; k < 3; k++) { // Derivative
326 grad_u[j][k] = 0;
327 for (CeedInt m = 0; m < 3; m++) grad_u[j][k] += dXdx[m][k] * du[j][m];
328 }
329 }
330
331 // Compute Strain : e (epsilon)
332 // e = 1/2 (grad u + (grad u)^T)
333
334 const CeedScalar e[3][3] = {
335 {(grad_u[0][0] + grad_u[0][0]) / 2., (grad_u[0][1] + grad_u[1][0]) / 2., (grad_u[0][2] + grad_u[2][0]) / 2.},
336 {(grad_u[1][0] + grad_u[0][1]) / 2., (grad_u[1][1] + grad_u[1][1]) / 2., (grad_u[1][2] + grad_u[2][1]) / 2.},
337 {(grad_u[2][0] + grad_u[0][2]) / 2., (grad_u[2][1] + grad_u[1][2]) / 2., (grad_u[2][2] + grad_u[2][2]) / 2.}
338 };
339
340 // Displacement
341 diagnostic[0][i] = u[0][i];
342 diagnostic[1][i] = u[1][i];
343 diagnostic[2][i] = u[2][i];
344
345 // Pressure
346 const CeedScalar strain_vol = e[0][0] + e[1][1] + e[2][2];
347 diagnostic[3][i] = -lambda * strain_vol;
348
349 // Stress tensor invariants
350 diagnostic[4][i] = strain_vol;
351 diagnostic[5][i] = 0.;
352 for (CeedInt j = 0; j < 3; j++) {
353 for (CeedInt m = 0; m < 3; m++) diagnostic[5][i] += e[j][m] * e[m][j];
354 }
355 diagnostic[6][i] = 1 + strain_vol;
356
357 // Strain energy
358 diagnostic[7][i] =
359 (lambda * strain_vol * strain_vol / 2. + strain_vol * mu + (e[0][1] * e[0][1] + e[0][2] * e[0][2] + e[1][2] * e[1][2]) * 2 * mu);
360 } // End of Quadrature Point Loop
361
362 return CEED_ERROR_SUCCESS;
363 }
364 // -----------------------------------------------------------------------------
365