1 // Copyright (c) 2017-2024, 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 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 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 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 // ----------------------------------------------------------------------------- 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 // ----------------------------------------------------------------------------- 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 // ----------------------------------------------------------------------------- 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 // ----------------------------------------------------------------------------- 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 // ----------------------------------------------------------------------------- 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