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