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 #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 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 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 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 // ----------------------------------------------------------------------------- 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 // ----------------------------------------------------------------------------- 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 // ----------------------------------------------------------------------------- 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 // ----------------------------------------------------------------------------- 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 // ----------------------------------------------------------------------------- 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