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