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