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 /// Hyperelasticity, finite strain 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 // Series approximation of log1p()
27 // log1p() is not vectorized in libc
28 //
29 // The series expansion is accurate to 1e-7 in the range sqrt(2)/2 < J < sqrt(2), with machine precision accuracy near J=1.
30 // The initialization extends this range to 0.35 ~= sqrt(2)/4 < J < sqrt(2)*2 ~= 2.83, which should be sufficient for applications of the Neo-Hookean
31 // model.
32 // -----------------------------------------------------------------------------
33 #ifndef LOG1P_SERIES_SHIFTED
34 #define LOG1P_SERIES_SHIFTED
log1p_series_shifted(CeedScalar x)35 CEED_QFUNCTION_HELPER CeedScalar log1p_series_shifted(CeedScalar x) {
36 const CeedScalar left = sqrt(2.) / 2 - 1, right = sqrt(2.) - 1;
37 CeedScalar sum = 0;
38 if (1) { // Disable if the smaller range sqrt(2)/2 < J < sqrt(2) is sufficient
39 if (x < left) { // Replace if with while for arbitrary range (may hurt vectorization)
40 sum -= log(2.) / 2;
41 x = 1 + 2 * x;
42 } else if (right < x) {
43 sum += log(2.) / 2;
44 x = (x - 1) / 2;
45 }
46 }
47 CeedScalar y = x / (2. + x);
48 const CeedScalar y2 = y * y;
49 sum += y;
50 y *= y2;
51 sum += y / 3;
52 y *= y2;
53 sum += y / 5;
54 y *= y2;
55 sum += y / 7;
56 return 2 * sum;
57 }
58 #endif
59
60 // -----------------------------------------------------------------------------
61 // Compute det F - 1
62 // -----------------------------------------------------------------------------
63 #ifndef DETJM1
64 #define DETJM1
computeJM1(const CeedScalar grad_u[3][3])65 CEED_QFUNCTION_HELPER CeedScalar computeJM1(const CeedScalar grad_u[3][3]) {
66 return grad_u[0][0] * (grad_u[1][1] * grad_u[2][2] - grad_u[1][2] * grad_u[2][1]) +
67 grad_u[0][1] * (grad_u[1][2] * grad_u[2][0] - grad_u[1][0] * grad_u[2][2]) +
68 grad_u[0][2] * (grad_u[1][0] * grad_u[2][1] - grad_u[2][0] * grad_u[1][1]) + grad_u[0][0] + grad_u[1][1] + grad_u[2][2] +
69 grad_u[0][0] * grad_u[1][1] + grad_u[0][0] * grad_u[2][2] + grad_u[1][1] * grad_u[2][2] - grad_u[0][1] * grad_u[1][0] -
70 grad_u[0][2] * grad_u[2][0] - grad_u[1][2] * grad_u[2][1];
71 }
72 #endif
73
74 // -----------------------------------------------------------------------------
75 // Compute matrix^(-1), where matrix is symetric, returns array of 6
76 // -----------------------------------------------------------------------------
77 #ifndef MatinvSym
78 #define MatinvSym
computeMatinvSym(const CeedScalar A[3][3],const CeedScalar detA,CeedScalar Ainv[6])79 CEED_QFUNCTION_HELPER int computeMatinvSym(const CeedScalar A[3][3], const CeedScalar detA, CeedScalar Ainv[6]) {
80 // Compute A^(-1) : A-Inverse
81 CeedScalar B[6] = {
82 A[1][1] * A[2][2] - A[1][2] * A[2][1], /* *NOPAD* */
83 A[0][0] * A[2][2] - A[0][2] * A[2][0], /* *NOPAD* */
84 A[0][0] * A[1][1] - A[0][1] * A[1][0], /* *NOPAD* */
85 A[0][2] * A[1][0] - A[0][0] * A[1][2], /* *NOPAD* */
86 A[0][1] * A[1][2] - A[0][2] * A[1][1], /* *NOPAD* */
87 A[0][2] * A[2][1] - A[0][1] * A[2][2] /* *NOPAD* */
88 };
89 for (CeedInt m = 0; m < 6; m++) Ainv[m] = B[m] / (detA);
90
91 return CEED_ERROR_SUCCESS;
92 }
93 #endif
94
95 // -----------------------------------------------------------------------------
96 // Common computations between FS and dFS
97 // -----------------------------------------------------------------------------
commonFS(const CeedScalar lambda,const CeedScalar mu,const CeedScalar grad_u[3][3],CeedScalar Swork[6],CeedScalar Cinvwork[6],CeedScalar * logJ)98 CEED_QFUNCTION_HELPER int commonFS(const CeedScalar lambda, const CeedScalar mu, const CeedScalar grad_u[3][3], CeedScalar Swork[6],
99 CeedScalar Cinvwork[6], CeedScalar *logJ) {
100 // E - Green-Lagrange strain tensor
101 // E = 1/2 (grad_u + grad_u^T + grad_u^T*grad_u)
102 const CeedInt indj[6] = {0, 1, 2, 1, 0, 0}, indk[6] = {0, 1, 2, 2, 2, 1};
103 CeedScalar E2work[6];
104 for (CeedInt m = 0; m < 6; m++) {
105 E2work[m] = grad_u[indj[m]][indk[m]] + grad_u[indk[m]][indj[m]];
106 for (CeedInt n = 0; n < 3; n++) E2work[m] += grad_u[n][indj[m]] * grad_u[n][indk[m]];
107 }
108 CeedScalar E2[3][3] = {
109 {E2work[0], E2work[5], E2work[4]},
110 {E2work[5], E2work[1], E2work[3]},
111 {E2work[4], E2work[3], E2work[2]}
112 };
113 // J-1
114 const CeedScalar Jm1 = computeJM1(grad_u);
115
116 // C : right Cauchy-Green tensor
117 // C = I + 2E
118 const CeedScalar C[3][3] = {
119 {1 + E2[0][0], E2[0][1], E2[0][2] },
120 {E2[0][1], 1 + E2[1][1], E2[1][2] },
121 {E2[0][2], E2[1][2], 1 + E2[2][2]}
122 };
123
124 // Compute C^(-1) : C-Inverse
125 const CeedScalar detC = (Jm1 + 1.) * (Jm1 + 1.);
126 computeMatinvSym(C, detC, Cinvwork);
127
128 const CeedScalar C_inv[3][3] = {
129 {Cinvwork[0], Cinvwork[5], Cinvwork[4]},
130 {Cinvwork[5], Cinvwork[1], Cinvwork[3]},
131 {Cinvwork[4], Cinvwork[3], Cinvwork[2]}
132 };
133
134 // Compute the Second Piola-Kirchhoff (S)
135 *logJ = log1p_series_shifted(Jm1);
136 for (CeedInt m = 0; m < 6; m++) {
137 Swork[m] = (lambda * (*logJ)) * Cinvwork[m];
138 for (CeedInt n = 0; n < 3; n++) Swork[m] += mu * C_inv[indj[m]][n] * E2[n][indk[m]];
139 }
140
141 return CEED_ERROR_SUCCESS;
142 }
143
144 // -----------------------------------------------------------------------------
145 // Residual evaluation for hyperelasticity, finite strain
146 // -----------------------------------------------------------------------------
ElasFSResidual_NH(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)147 CEED_QFUNCTION(ElasFSResidual_NH)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
148 // Inputs
149 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];
150
151 // Outputs
152 CeedScalar(*dvdX)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[0];
153 // Store grad_u for HyperFSdF (Jacobian of HyperFSF)
154 CeedScalar(*grad_u)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[1];
155
156 // Context
157 const Physics context = (Physics)ctx;
158 const CeedScalar E = context->E;
159 const CeedScalar nu = context->nu;
160 const CeedScalar TwoMu = E / (1 + nu);
161 const CeedScalar mu = TwoMu / 2;
162 const CeedScalar Kbulk = E / (3 * (1 - 2 * nu)); // Bulk Modulus
163 const CeedScalar lambda = (3 * Kbulk - TwoMu) / 3;
164
165 // Formulation Terminology:
166 // I3 : 3x3 Identity matrix
167 // C : right Cauchy-Green tensor
168 // C_inv : inverse of C
169 // F : deformation gradient
170 // S : 2nd Piola-Kirchhoff (in current config)
171 // P : 1st Piola-Kirchhoff (in referential config)
172 // Formulation:
173 // F = I3 + grad_ue
174 // J = det(F)
175 // C = F(^T)*F
176 // S = mu*I3 + (lambda*log(J)-mu)*C_inv;
177 // P = F*S
178
179 // Quadrature Point Loop
180 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
181 // Read spatial derivatives of u
182 const CeedScalar du[3][3] = {
183 {ug[0][0][i], ug[1][0][i], ug[2][0][i]},
184 {ug[0][1][i], ug[1][1][i], ug[2][1][i]},
185 {ug[0][2][i], ug[1][2][i], ug[2][2][i]}
186 };
187 // -- Qdata
188 const CeedScalar wdetJ = q_data[0][i];
189 const CeedScalar dXdx[3][3] = {
190 {q_data[1][i], q_data[2][i], q_data[3][i]},
191 {q_data[4][i], q_data[5][i], q_data[6][i]},
192 {q_data[7][i], q_data[8][i], q_data[9][i]}
193 };
194
195 // Compute grad_u
196 // dXdx = (dx/dX)^(-1)
197 // Apply dXdx to du = grad_u
198 for (CeedInt j = 0; j < 3; j++) { // Component
199 for (CeedInt k = 0; k < 3; k++) { // Derivative
200 grad_u[j][k][i] = 0;
201 for (CeedInt m = 0; m < 3; m++) grad_u[j][k][i] += dXdx[m][k] * du[j][m];
202 }
203 }
204
205 // I3 : 3x3 Identity matrix
206 // Compute The Deformation Gradient : F = I3 + grad_u
207 const CeedScalar F[3][3] = {
208 {grad_u[0][0][i] + 1, grad_u[0][1][i], grad_u[0][2][i] },
209 {grad_u[1][0][i], grad_u[1][1][i] + 1, grad_u[1][2][i] },
210 {grad_u[2][0][i], grad_u[2][1][i], grad_u[2][2][i] + 1}
211 };
212
213 // Common components of finite strain calculations
214 CeedScalar Swork[6], Cinvwork[6], logJ;
215 const CeedScalar tempgradu[3][3] = {
216 {grad_u[0][0][i], grad_u[0][1][i], grad_u[0][2][i]},
217 {grad_u[1][0][i], grad_u[1][1][i], grad_u[1][2][i]},
218 {grad_u[2][0][i], grad_u[2][1][i], grad_u[2][2][i]}
219 };
220 commonFS(lambda, mu, tempgradu, Swork, Cinvwork, &logJ);
221
222 // Second Piola-Kirchhoff (S)
223 const CeedScalar S[3][3] = {
224 {Swork[0], Swork[5], Swork[4]},
225 {Swork[5], Swork[1], Swork[3]},
226 {Swork[4], Swork[3], Swork[2]}
227 };
228
229 // Compute the First Piola-Kirchhoff : P = F*S
230 CeedScalar P[3][3];
231 for (CeedInt j = 0; j < 3; j++) {
232 for (CeedInt k = 0; k < 3; k++) {
233 P[j][k] = 0;
234 for (CeedInt m = 0; m < 3; m++) P[j][k] += F[j][m] * S[m][k];
235 }
236 }
237
238 // Apply dXdx^T and weight to P (First Piola-Kirchhoff)
239 for (CeedInt j = 0; j < 3; j++) { // Component
240 for (CeedInt k = 0; k < 3; k++) { // Derivative
241 dvdX[k][j][i] = 0;
242 for (CeedInt m = 0; m < 3; m++) dvdX[k][j][i] += dXdx[k][m] * P[j][m] * wdetJ;
243 }
244 }
245 } // End of Quadrature Point Loop
246
247 return CEED_ERROR_SUCCESS;
248 }
249
250 // -----------------------------------------------------------------------------
251 // Jacobian evaluation for hyperelasticity, finite strain
252 // -----------------------------------------------------------------------------
ElasFSJacobian_NH(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)253 CEED_QFUNCTION(ElasFSJacobian_NH)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
254 // Inputs
255 const CeedScalar(*deltaug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0],
256 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];
257 // grad_u is used for hyperelasticity (non-linear)
258 const CeedScalar(*grad_u)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[2];
259
260 // Outputs
261 CeedScalar(*deltadvdX)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[0];
262
263 // Context
264 const Physics context = (Physics)ctx;
265 const CeedScalar E = context->E;
266 const CeedScalar nu = context->nu;
267
268 // Constants
269 const CeedScalar TwoMu = E / (1 + nu);
270 const CeedScalar mu = TwoMu / 2;
271 const CeedScalar Kbulk = E / (3 * (1 - 2 * nu)); // Bulk Modulus
272 const CeedScalar lambda = (3 * Kbulk - TwoMu) / 3;
273
274 // Quadrature Point Loop
275 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
276 // Read spatial derivatives of delta_u
277 const CeedScalar deltadu[3][3] = {
278 {deltaug[0][0][i], deltaug[1][0][i], deltaug[2][0][i]},
279 {deltaug[0][1][i], deltaug[1][1][i], deltaug[2][1][i]},
280 {deltaug[0][2][i], deltaug[1][2][i], deltaug[2][2][i]}
281 };
282 // -- Qdata
283 const CeedScalar wdetJ = q_data[0][i];
284 const CeedScalar dXdx[3][3] = {
285 {q_data[1][i], q_data[2][i], q_data[3][i]},
286 {q_data[4][i], q_data[5][i], q_data[6][i]},
287 {q_data[7][i], q_data[8][i], q_data[9][i]}
288 };
289
290 // Compute graddeltau
291 // dXdx = (dx/dX)^(-1)
292 // Apply dXdx to deltadu = graddelta
293 CeedScalar graddeltau[3][3];
294 for (CeedInt j = 0; j < 3; j++) { // Component
295 for (CeedInt k = 0; k < 3; k++) { // Derivative
296 graddeltau[j][k] = 0;
297 for (CeedInt m = 0; m < 3; m++) graddeltau[j][k] += dXdx[m][k] * deltadu[j][m];
298 }
299 }
300
301 // I3 : 3x3 Identity matrix
302 // Deformation Gradient : F = I3 + grad_u
303 const CeedScalar F[3][3] = {
304 {grad_u[0][0][i] + 1, grad_u[0][1][i], grad_u[0][2][i] },
305 {grad_u[1][0][i], grad_u[1][1][i] + 1, grad_u[1][2][i] },
306 {grad_u[2][0][i], grad_u[2][1][i], grad_u[2][2][i] + 1}
307 };
308
309 // Common components of finite strain calculations
310 CeedScalar Swork[6], Cinvwork[6], logJ;
311 const CeedScalar tempgradu[3][3] = {
312 {grad_u[0][0][i], grad_u[0][1][i], grad_u[0][2][i]},
313 {grad_u[1][0][i], grad_u[1][1][i], grad_u[1][2][i]},
314 {grad_u[2][0][i], grad_u[2][1][i], grad_u[2][2][i]}
315 };
316 commonFS(lambda, mu, tempgradu, Swork, Cinvwork, &logJ);
317
318 // deltaE - Green-Lagrange strain tensor
319 const CeedInt indj[6] = {0, 1, 2, 1, 0, 0}, indk[6] = {0, 1, 2, 2, 2, 1};
320 CeedScalar deltaEwork[6];
321 for (CeedInt m = 0; m < 6; m++) {
322 deltaEwork[m] = 0;
323 for (CeedInt n = 0; n < 3; n++) deltaEwork[m] += (graddeltau[n][indj[m]] * F[n][indk[m]] + F[n][indj[m]] * graddeltau[n][indk[m]]) / 2.;
324 }
325 CeedScalar deltaE[3][3] = {
326 {deltaEwork[0], deltaEwork[5], deltaEwork[4]},
327 {deltaEwork[5], deltaEwork[1], deltaEwork[3]},
328 {deltaEwork[4], deltaEwork[3], deltaEwork[2]}
329 };
330
331 // C : right Cauchy-Green tensor
332 // C^(-1) : C-Inverse
333 const CeedScalar C_inv[3][3] = {
334 {Cinvwork[0], Cinvwork[5], Cinvwork[4]},
335 {Cinvwork[5], Cinvwork[1], Cinvwork[3]},
336 {Cinvwork[4], Cinvwork[3], Cinvwork[2]}
337 };
338
339 // Second Piola-Kirchhoff (S)
340 const CeedScalar S[3][3] = {
341 {Swork[0], Swork[5], Swork[4]},
342 {Swork[5], Swork[1], Swork[3]},
343 {Swork[4], Swork[3], Swork[2]}
344 };
345
346 // deltaS = dSdE:deltaE
347 // = lambda(C_inv:deltaE)C_inv + 2(mu-lambda*log(J))C_inv*deltaE*C_inv
348 // -- C_inv:deltaE
349 CeedScalar Cinv_contract_E = 0;
350 for (CeedInt j = 0; j < 3; j++) {
351 for (CeedInt k = 0; k < 3; k++) Cinv_contract_E += C_inv[j][k] * deltaE[j][k];
352 }
353 // -- deltaE*C_inv
354 CeedScalar deltaECinv[3][3];
355 for (CeedInt j = 0; j < 3; j++) {
356 for (CeedInt k = 0; k < 3; k++) {
357 deltaECinv[j][k] = 0;
358 for (CeedInt m = 0; m < 3; m++) deltaECinv[j][k] += deltaE[j][m] * C_inv[m][k];
359 }
360 }
361 // -- intermediate deltaS = C_inv*deltaE*C_inv
362 CeedScalar deltaS[3][3];
363 for (CeedInt j = 0; j < 3; j++) {
364 for (CeedInt k = 0; k < 3; k++) {
365 deltaS[j][k] = 0;
366 for (CeedInt m = 0; m < 3; m++) deltaS[j][k] += C_inv[j][m] * deltaECinv[m][k];
367 }
368 }
369 // -- deltaS = lambda(C_inv:deltaE)C_inv - 2(lambda*log(J)-mu)*(intermediate)
370 for (CeedInt j = 0; j < 3; j++) {
371 for (CeedInt k = 0; k < 3; k++) deltaS[j][k] = lambda * Cinv_contract_E * C_inv[j][k] - 2. * (lambda * logJ - mu) * deltaS[j][k];
372 }
373
374 // deltaP = dPdF:deltaF = deltaF*S + F*deltaS
375 CeedScalar deltaP[3][3];
376 for (CeedInt j = 0; j < 3; j++) {
377 for (CeedInt k = 0; k < 3; k++) {
378 deltaP[j][k] = 0;
379 for (CeedInt m = 0; m < 3; m++) deltaP[j][k] += graddeltau[j][m] * S[m][k] + F[j][m] * deltaS[m][k];
380 }
381 }
382
383 // Apply dXdx^T and weight
384 for (CeedInt j = 0; j < 3; j++) { // Component
385 for (CeedInt k = 0; k < 3; k++) { // Derivative
386 deltadvdX[k][j][i] = 0;
387 for (CeedInt m = 0; m < 3; m++) deltadvdX[k][j][i] += dXdx[k][m] * deltaP[j][m] * wdetJ;
388 }
389 }
390 } // End of Quadrature Point Loop
391
392 return CEED_ERROR_SUCCESS;
393 }
394
395 // -----------------------------------------------------------------------------
396 // Strain energy computation for hyperelasticity, finite strain
397 // -----------------------------------------------------------------------------
ElasFSEnergy_NH(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)398 CEED_QFUNCTION(ElasFSEnergy_NH)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
399 // Inputs
400 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];
401
402 // Outputs
403 CeedScalar(*energy) = (CeedScalar(*))out[0];
404
405 // Context
406 const Physics context = (Physics)ctx;
407 const CeedScalar E = context->E;
408 const CeedScalar nu = context->nu;
409 const CeedScalar TwoMu = E / (1 + nu);
410 const CeedScalar mu = TwoMu / 2;
411 const CeedScalar Kbulk = E / (3 * (1 - 2 * nu)); // Bulk Modulus
412 const CeedScalar lambda = (3 * Kbulk - TwoMu) / 3;
413
414 // Quadrature Point Loop
415 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
416 // Read spatial derivatives of u
417 const CeedScalar du[3][3] = {
418 {ug[0][0][i], ug[1][0][i], ug[2][0][i]},
419 {ug[0][1][i], ug[1][1][i], ug[2][1][i]},
420 {ug[0][2][i], ug[1][2][i], ug[2][2][i]}
421 };
422 // -- Qdata
423 const CeedScalar wdetJ = q_data[0][i];
424 const CeedScalar dXdx[3][3] = {
425 {q_data[1][i], q_data[2][i], q_data[3][i]},
426 {q_data[4][i], q_data[5][i], q_data[6][i]},
427 {q_data[7][i], q_data[8][i], q_data[9][i]}
428 };
429
430 // Compute grad_u
431 // dXdx = (dx/dX)^(-1)
432 // Apply dXdx to du = grad_u
433 CeedScalar grad_u[3][3];
434 for (int j = 0; j < 3; j++) { // Component
435 for (int k = 0; k < 3; k++) { // Derivative
436 grad_u[j][k] = 0;
437 for (int m = 0; m < 3; m++) grad_u[j][k] += dXdx[m][k] * du[j][m];
438 }
439 }
440
441 // E - Green-Lagrange strain tensor
442 // E = 1/2 (grad_u + grad_u^T + grad_u^T*grad_u)
443 const CeedInt indj[6] = {0, 1, 2, 1, 0, 0}, indk[6] = {0, 1, 2, 2, 2, 1};
444 CeedScalar E2work[6];
445 for (CeedInt m = 0; m < 6; m++) {
446 E2work[m] = grad_u[indj[m]][indk[m]] + grad_u[indk[m]][indj[m]];
447 for (CeedInt n = 0; n < 3; n++) E2work[m] += grad_u[n][indj[m]] * grad_u[n][indk[m]];
448 }
449 CeedScalar E2[3][3] = {
450 {E2work[0], E2work[5], E2work[4]},
451 {E2work[5], E2work[1], E2work[3]},
452 {E2work[4], E2work[3], E2work[2]}
453 };
454 const CeedScalar Jm1 = computeJM1(grad_u);
455 const CeedScalar logJ = log1p_series_shifted(Jm1);
456
457 // Strain energy Phi(E) for compressible Neo-Hookean
458 energy[i] = (lambda * logJ * logJ / 2. - mu * logJ + mu * (E2[0][0] + E2[1][1] + E2[2][2]) / 2.) * wdetJ;
459
460 } // End of Quadrature Point Loop
461
462 return CEED_ERROR_SUCCESS;
463 }
464
465 // -----------------------------------------------------------------------------
466 // Nodal diagnostic quantities for hyperelasticity, finite strain
467 // -----------------------------------------------------------------------------
ElasFSDiagnostic_NH(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)468 CEED_QFUNCTION(ElasFSDiagnostic_NH)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
469 // Inputs
470 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],
471 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2];
472
473 // Outputs
474 CeedScalar(*diagnostic)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
475
476 // Context
477 const Physics context = (Physics)ctx;
478 const CeedScalar E = context->E;
479 const CeedScalar nu = context->nu;
480 const CeedScalar TwoMu = E / (1 + nu);
481 const CeedScalar mu = TwoMu / 2;
482 const CeedScalar Kbulk = E / (3 * (1 - 2 * nu)); // Bulk Modulus
483 const CeedScalar lambda = (3 * Kbulk - TwoMu) / 3;
484
485 // Quadrature Point Loop
486 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
487 // Read spatial derivatives of u
488 const CeedScalar du[3][3] = {
489 {ug[0][0][i], ug[1][0][i], ug[2][0][i]},
490 {ug[0][1][i], ug[1][1][i], ug[2][1][i]},
491 {ug[0][2][i], ug[1][2][i], ug[2][2][i]}
492 };
493 // -- Qdata
494 const CeedScalar dXdx[3][3] = {
495 {q_data[1][i], q_data[2][i], q_data[3][i]},
496 {q_data[4][i], q_data[5][i], q_data[6][i]},
497 {q_data[7][i], q_data[8][i], q_data[9][i]}
498 };
499
500 // Compute grad_u
501 // dXdx = (dx/dX)^(-1)
502 // Apply dXdx to du = grad_u
503 CeedScalar grad_u[3][3];
504 for (int j = 0; j < 3; j++) { // Component
505 for (int k = 0; k < 3; k++) { // Derivative
506 grad_u[j][k] = 0;
507 for (int m = 0; m < 3; m++) grad_u[j][k] += dXdx[m][k] * du[j][m];
508 }
509 }
510
511 // E - Green-Lagrange strain tensor
512 // E = 1/2 (grad_u + grad_u^T + grad_u^T*grad_u)
513 const CeedInt indj[6] = {0, 1, 2, 1, 0, 0}, indk[6] = {0, 1, 2, 2, 2, 1};
514 CeedScalar E2work[6];
515 for (CeedInt m = 0; m < 6; m++) {
516 E2work[m] = grad_u[indj[m]][indk[m]] + grad_u[indk[m]][indj[m]];
517 for (CeedInt n = 0; n < 3; n++) E2work[m] += grad_u[n][indj[m]] * grad_u[n][indk[m]];
518 }
519 CeedScalar E2[3][3] = {
520 {E2work[0], E2work[5], E2work[4]},
521 {E2work[5], E2work[1], E2work[3]},
522 {E2work[4], E2work[3], E2work[2]}
523 };
524
525 // Displacement
526 diagnostic[0][i] = u[0][i];
527 diagnostic[1][i] = u[1][i];
528 diagnostic[2][i] = u[2][i];
529
530 // Pressure
531 const CeedScalar Jm1 = computeJM1(grad_u);
532 const CeedScalar logJ = log1p_series_shifted(Jm1);
533 diagnostic[3][i] = -lambda * logJ;
534
535 // Stress tensor invariants
536 diagnostic[4][i] = (E2[0][0] + E2[1][1] + E2[2][2]) / 2.;
537 diagnostic[5][i] = 0.;
538 for (CeedInt j = 0; j < 3; j++) {
539 for (CeedInt m = 0; m < 3; m++) diagnostic[5][i] += E2[j][m] * E2[m][j] / 4.;
540 }
541 diagnostic[6][i] = Jm1 + 1.;
542
543 // Strain energy
544 diagnostic[7][i] = (lambda * logJ * logJ / 2. - mu * logJ + mu * (E2[0][0] + E2[1][1] + E2[2][2]) / 2.);
545 } // End of Quadrature Point Loop
546
547 return CEED_ERROR_SUCCESS;
548 }
549 // -----------------------------------------------------------------------------
550