1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Linear elasticity for solid mechanics example using PETSc 19 20 #ifndef ELAS_LINEAR_H 21 #define ELAS_LINEAR_H 22 23 #include <math.h> 24 25 #ifndef PHYSICS_STRUCT 26 #define PHYSICS_STRUCT 27 typedef struct Physics_private *Physics; 28 struct Physics_private { 29 CeedScalar nu; // Poisson's ratio 30 CeedScalar E; // Young's Modulus 31 }; 32 #endif 33 34 // ----------------------------------------------------------------------------- 35 // Residual evaluation for linear elasticity 36 // ----------------------------------------------------------------------------- 37 CEED_QFUNCTION(ElasLinearF)(void *ctx, CeedInt Q, const CeedScalar *const *in, 38 CeedScalar *const *out) { 39 // *INDENT-OFF* 40 // Inputs 41 const CeedScalar (*ug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], 42 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 43 44 // Outputs 45 CeedScalar (*dvdX)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[0]; 46 // grad_u not used for linear elasticity 47 // (*grad_u)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[1]; 48 // *INDENT-ON* 49 50 // Context 51 const Physics context = (Physics)ctx; 52 const CeedScalar E = context->E; 53 const CeedScalar nu = context->nu; 54 55 // Quadrature Point Loop 56 CeedPragmaSIMD 57 for (CeedInt i=0; i<Q; i++) { 58 // Read spatial derivatives of u 59 // *INDENT-OFF* 60 const CeedScalar du[3][3] = {{ug[0][0][i], 61 ug[1][0][i], 62 ug[2][0][i]}, 63 {ug[0][1][i], 64 ug[1][1][i], 65 ug[2][1][i]}, 66 {ug[0][2][i], 67 ug[1][2][i], 68 ug[2][2][i]} 69 }; 70 // -- Qdata 71 const CeedScalar wdetJ = q_data[0][i]; 72 const CeedScalar dXdx[3][3] = {{q_data[1][i], 73 q_data[2][i], 74 q_data[3][i]}, 75 {q_data[4][i], 76 q_data[5][i], 77 q_data[6][i]}, 78 {q_data[7][i], 79 q_data[8][i], 80 q_data[9][i]} 81 }; 82 // *INDENT-ON* 83 84 // Compute grad_u 85 // dXdx = (dx/dX)^(-1) 86 // Apply dXdx to du = grad_u 87 CeedScalar grad_u[3][3]; 88 for (CeedInt j = 0; j < 3; j++) // Component 89 for (CeedInt k = 0; k < 3; k++) { // Derivative 90 grad_u[j][k] = 0; 91 for (CeedInt m = 0; m < 3; m++) 92 grad_u[j][k] += dXdx[m][k] * du[j][m]; 93 } 94 95 // Compute Strain : e (epsilon) 96 // e = 1/2 (grad u + (grad u)^T) 97 98 // *INDENT-OFF* 99 const CeedScalar e[3][3] = {{(grad_u[0][0] + grad_u[0][0])/2., 100 (grad_u[0][1] + grad_u[1][0])/2., 101 (grad_u[0][2] + grad_u[2][0])/2.}, 102 {(grad_u[1][0] + grad_u[0][1])/2., 103 (grad_u[1][1] + grad_u[1][1])/2., 104 (grad_u[1][2] + grad_u[2][1])/2.}, 105 {(grad_u[2][0] + grad_u[0][2])/2., 106 (grad_u[2][1] + grad_u[1][2])/2., 107 (grad_u[2][2] + grad_u[2][2])/2.} 108 }; 109 // *INDENT-ON* 110 111 // 112 // Formulation uses Voigt notation: 113 // stress (sigma) strain (epsilon) 114 // 115 // [sigma00] [e00] 116 // [sigma11] [e11] 117 // [sigma22] = S * [e22] 118 // [sigma12] [e12] 119 // [sigma02] [e02] 120 // [sigma01] [e01] 121 // 122 // where 123 // [1-nu nu nu ] 124 // [ nu 1-nu nu ] 125 // [ nu nu 1-nu ] 126 // S = E/((1+nu)*(1-2*nu)) [ (1-2*nu)/2 ] 127 // [ (1-2*nu)/2 ] 128 // [ (1-2*nu)/2 ] 129 130 // Above Voigt Notation is placed in a 3x3 matrix: 131 const CeedScalar ss = E / ((1 + nu)*(1 - 2*nu)); 132 const CeedScalar sigma00 = ss*((1 - nu)*e[0][0] + nu*e[1][1] + nu*e[2][2]), 133 sigma11 = ss*(nu*e[0][0] + (1 - nu)*e[1][1] + nu*e[2][2]), 134 sigma22 = ss*(nu*e[0][0] + nu*e[1][1] + (1 - nu)*e[2][2]), 135 sigma12 = ss*(1 - 2*nu)*e[1][2]*0.5, 136 sigma02 = ss*(1 - 2*nu)*e[0][2]*0.5, 137 sigma01 = ss*(1 - 2*nu)*e[0][1]*0.5; 138 // *INDENT-OFF* 139 const CeedScalar sigma[3][3] = {{sigma00, sigma01, sigma02}, 140 {sigma01, sigma11, sigma12}, 141 {sigma02, sigma12, sigma22} 142 }; 143 // *INDENT-ON* 144 145 // Apply dXdx^T and weight to sigma 146 for (CeedInt j = 0; j < 3; j++) // Component 147 for (CeedInt k = 0; k < 3; k++) { // Derivative 148 dvdX[k][j][i] = 0; 149 for (CeedInt m = 0; m < 3; m++) 150 dvdX[k][j][i] += dXdx[k][m] * sigma[j][m] * wdetJ; 151 } 152 153 } // End of Quadrature Point Loop 154 155 return 0; 156 } 157 158 // ----------------------------------------------------------------------------- 159 // Jacobian evaluation for linear elasticity 160 // ----------------------------------------------------------------------------- 161 CEED_QFUNCTION(ElasLineardF)(void *ctx, CeedInt Q, const CeedScalar *const *in, 162 CeedScalar *const *out) { 163 // *INDENT-OFF* 164 // Inputs 165 const CeedScalar (*deltaug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], 166 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 167 // grad_u not used for linear elasticity 168 // (*grad_u)[3][Q] = (CeedScalar(*)[3][Q])in[2]; 169 170 // Outputs 171 CeedScalar (*deltadvdX)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[0]; 172 // *INDENT-ON* 173 174 // Context 175 const Physics context = (Physics)ctx; 176 const CeedScalar E = context->E; 177 const CeedScalar nu = context->nu; 178 179 // Quadrature Point Loop 180 CeedPragmaSIMD 181 for (CeedInt i=0; i<Q; i++) { 182 // Read spatial derivatives of u 183 // *INDENT-OFF* 184 const CeedScalar deltadu[3][3] = {{deltaug[0][0][i], 185 deltaug[1][0][i], 186 deltaug[2][0][i]}, 187 {deltaug[0][1][i], 188 deltaug[1][1][i], 189 deltaug[2][1][i]}, 190 {deltaug[0][2][i], 191 deltaug[1][2][i], 192 deltaug[2][2][i]} 193 }; 194 // -- Qdata 195 const CeedScalar wdetJ = q_data[0][i]; 196 const CeedScalar dXdx[3][3] = {{q_data[1][i], 197 q_data[2][i], 198 q_data[3][i]}, 199 {q_data[4][i], 200 q_data[5][i], 201 q_data[6][i]}, 202 {q_data[7][i], 203 q_data[8][i], 204 q_data[9][i]} 205 }; 206 // *INDENT-ON* 207 208 // Compute graddeltau 209 // dXdx = (dx/dX)^(-1) 210 // Apply dXdx to deltadu = graddeltau 211 CeedScalar graddeltau[3][3]; 212 for (CeedInt j = 0; j < 3; j++) // Component 213 for (CeedInt k = 0; k < 3; k++) { // Derivative 214 graddeltau[j][k] = 0; 215 for (CeedInt m = 0; m < 3; m++) 216 graddeltau[j][k] += dXdx[m][k] * deltadu[j][m]; 217 } 218 219 // Compute Strain : e (epsilon) 220 // e = 1/2 (grad u + (grad u)^T) 221 // *INDENT-OFF* 222 const CeedScalar de[3][3] = {{(graddeltau[0][0] + graddeltau[0][0])/2., 223 (graddeltau[0][1] + graddeltau[1][0])/2., 224 (graddeltau[0][2] + graddeltau[2][0])/2.}, 225 {(graddeltau[1][0] + graddeltau[0][1])/2., 226 (graddeltau[1][1] + graddeltau[1][1])/2., 227 (graddeltau[1][2] + graddeltau[2][1])/2.}, 228 {(graddeltau[2][0] + graddeltau[0][2])/2., 229 (graddeltau[2][1] + graddeltau[1][2])/2., 230 (graddeltau[2][2] + graddeltau[2][2])/2.} 231 }; 232 233 // *INDENT-ON* 234 // Formulation uses Voigt notation: 235 // stress (sigma) strain (epsilon) 236 // 237 // [dsigma00] [de00] 238 // [dsigma11] [de11] 239 // [dsigma22] = S * [de22] 240 // [dsigma12] [de12] 241 // [dsigma02] [de02] 242 // [dsigma01] [de01] 243 // 244 // where 245 // [1-nu nu nu ] 246 // [ nu 1-nu nu ] 247 // [ nu nu 1-nu ] 248 // S = E/((1+nu)*(1-2*nu)) [ (1-2*nu)/2 ] 249 // [ (1-2*nu)/2 ] 250 // [ (1-2*nu)/2 ] 251 252 // Above Voigt Notation is placed in a 3x3 matrix: 253 const CeedScalar ss = E / ((1 + nu)*(1 - 2*nu)); 254 const CeedScalar dsigma00 = ss*((1 - nu)*de[0][0]+nu*de[1][1]+nu*de[2][2]), 255 dsigma11 = ss*(nu*de[0][0]+(1 - nu)*de[1][1]+nu*de[2][2]), 256 dsigma22 = ss*(nu*de[0][0]+nu*de[1][1]+(1 - nu)*de[2][2]), 257 dsigma12 = ss*(1 - 2*nu)*de[1][2] / 2, 258 dsigma02 = ss*(1 - 2*nu)*de[0][2] / 2, 259 dsigma01 = ss*(1 - 2*nu)*de[0][1] / 2; 260 // *INDENT-OFF* 261 const CeedScalar dsigma[3][3] = {{dsigma00, dsigma01, dsigma02}, 262 {dsigma01, dsigma11, dsigma12}, 263 {dsigma02, dsigma12, dsigma22} 264 }; 265 // *INDENT-ON* 266 267 // Apply dXdx^T and weight 268 for (CeedInt j = 0; j < 3; j++) // Component 269 for (CeedInt k = 0; k < 3 ; k++) { // Derivative 270 deltadvdX[k][j][i] = 0; 271 for (CeedInt m = 0; m < 3; m++) 272 deltadvdX[k][j][i] += dXdx[k][m] * dsigma[j][m] * wdetJ; 273 } 274 275 } // End of Quadrature Point Loop 276 277 return 0; 278 } 279 280 // ----------------------------------------------------------------------------- 281 // Strain energy computation for linear elasticity 282 // ----------------------------------------------------------------------------- 283 CEED_QFUNCTION(ElasLinearEnergy)(void *ctx, CeedInt Q, 284 const CeedScalar *const *in, 285 CeedScalar *const *out) { 286 // *INDENT-OFF* 287 // Inputs 288 const CeedScalar (*ug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], 289 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 290 291 // Outputs 292 CeedScalar (*energy) = (CeedScalar(*))out[0]; 293 // *INDENT-ON* 294 295 // Context 296 const Physics context = (Physics)ctx; 297 const CeedScalar E = context->E; 298 const CeedScalar nu = context->nu; 299 300 // Constants 301 const CeedScalar TwoMu = E / (1 + nu); 302 const CeedScalar mu = TwoMu / 2; 303 const CeedScalar Kbulk = E / (3*(1 - 2*nu)); // Bulk Modulus 304 const CeedScalar lambda = (3*Kbulk - TwoMu) / 3; 305 306 // Quadrature Point Loop 307 CeedPragmaSIMD 308 for (CeedInt i=0; i<Q; i++) { 309 // Read spatial derivatives of u 310 // *INDENT-OFF* 311 const CeedScalar du[3][3] = {{ug[0][0][i], 312 ug[1][0][i], 313 ug[2][0][i]}, 314 {ug[0][1][i], 315 ug[1][1][i], 316 ug[2][1][i]}, 317 {ug[0][2][i], 318 ug[1][2][i], 319 ug[2][2][i]} 320 }; 321 // -- Qdata 322 const CeedScalar wdetJ = q_data[0][i]; 323 const CeedScalar dXdx[3][3] = {{q_data[1][i], 324 q_data[2][i], 325 q_data[3][i]}, 326 {q_data[4][i], 327 q_data[5][i], 328 q_data[6][i]}, 329 {q_data[7][i], 330 q_data[8][i], 331 q_data[9][i]} 332 }; 333 // *INDENT-ON* 334 335 // Compute grad_u 336 // dXdx = (dx/dX)^(-1) 337 // Apply dXdx to du = grad_u 338 CeedScalar grad_u[3][3]; 339 for (CeedInt j = 0; j < 3; j++) // Component 340 for (CeedInt k = 0; k < 3; k++) { // Derivative 341 grad_u[j][k] = 0; 342 for (CeedInt m = 0; m < 3; m++) 343 grad_u[j][k] += dXdx[m][k] * du[j][m]; 344 } 345 346 // Compute Strain : e (epsilon) 347 // e = 1/2 (grad u + (grad u)^T) 348 349 // *INDENT-OFF* 350 const CeedScalar e[3][3] = {{(grad_u[0][0] + grad_u[0][0])/2., 351 (grad_u[0][1] + grad_u[1][0])/2., 352 (grad_u[0][2] + grad_u[2][0])/2.}, 353 {(grad_u[1][0] + grad_u[0][1])/2., 354 (grad_u[1][1] + grad_u[1][1])/2., 355 (grad_u[1][2] + grad_u[2][1])/2.}, 356 {(grad_u[2][0] + grad_u[0][2])/2., 357 (grad_u[2][1] + grad_u[1][2])/2., 358 (grad_u[2][2] + grad_u[2][2])/2.} 359 }; 360 // *INDENT-ON* 361 362 // Strain energy 363 const CeedScalar strain_vol = e[0][0] + e[1][1] + e[2][2]; 364 energy[i] = (lambda*strain_vol*strain_vol/2. + strain_vol*mu + 365 (e[0][1]*e[0][1]+e[0][2]*e[0][2]+e[1][2]*e[1][2])*2*mu)*wdetJ; 366 367 } // End of Quadrature Point Loop 368 369 return 0; 370 } 371 372 // ----------------------------------------------------------------------------- 373 // Nodal diagnostic quantities for linear elasticity 374 // ----------------------------------------------------------------------------- 375 CEED_QFUNCTION(ElasLinearDiagnostic)(void *ctx, CeedInt Q, 376 const CeedScalar *const *in, 377 CeedScalar *const *out) { 378 // *INDENT-OFF* 379 // Inputs 380 const CeedScalar (*u)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 381 (*ug)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1], 382 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 383 384 // Outputs 385 CeedScalar (*diagnostic)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 386 // *INDENT-ON* 387 388 // Context 389 const Physics context = (Physics)ctx; 390 const CeedScalar E = context->E; 391 const CeedScalar nu = context->nu; 392 393 // Constants 394 const CeedScalar TwoMu = E / (1 + nu); 395 const CeedScalar mu = TwoMu / 2; 396 const CeedScalar Kbulk = E / (3*(1 - 2*nu)); // Bulk Modulus 397 const CeedScalar lambda = (3*Kbulk - TwoMu) / 3; 398 399 // Quadrature Point Loop 400 CeedPragmaSIMD 401 for (CeedInt i=0; i<Q; i++) { 402 // Read spatial derivatives of u 403 // *INDENT-OFF* 404 const CeedScalar du[3][3] = {{ug[0][0][i], 405 ug[1][0][i], 406 ug[2][0][i]}, 407 {ug[0][1][i], 408 ug[1][1][i], 409 ug[2][1][i]}, 410 {ug[0][2][i], 411 ug[1][2][i], 412 ug[2][2][i]} 413 }; 414 // -- Qdata 415 const CeedScalar dXdx[3][3] = {{q_data[1][i], 416 q_data[2][i], 417 q_data[3][i]}, 418 {q_data[4][i], 419 q_data[5][i], 420 q_data[6][i]}, 421 {q_data[7][i], 422 q_data[8][i], 423 q_data[9][i]} 424 }; 425 // *INDENT-ON* 426 427 // Compute grad_u 428 // dXdx = (dx/dX)^(-1) 429 // Apply dXdx to du = grad_u 430 CeedScalar grad_u[3][3]; 431 for (CeedInt j = 0; j < 3; j++) // Component 432 for (CeedInt k = 0; k < 3; k++) { // Derivative 433 grad_u[j][k] = 0; 434 for (CeedInt m = 0; m < 3; m++) 435 grad_u[j][k] += dXdx[m][k] * du[j][m]; 436 } 437 438 // Compute Strain : e (epsilon) 439 // e = 1/2 (grad u + (grad u)^T) 440 441 // *INDENT-OFF* 442 const CeedScalar e[3][3] = {{(grad_u[0][0] + grad_u[0][0])/2., 443 (grad_u[0][1] + grad_u[1][0])/2., 444 (grad_u[0][2] + grad_u[2][0])/2.}, 445 {(grad_u[1][0] + grad_u[0][1])/2., 446 (grad_u[1][1] + grad_u[1][1])/2., 447 (grad_u[1][2] + grad_u[2][1])/2.}, 448 {(grad_u[2][0] + grad_u[0][2])/2., 449 (grad_u[2][1] + grad_u[1][2])/2., 450 (grad_u[2][2] + grad_u[2][2])/2.} 451 }; 452 // *INDENT-ON* 453 454 // Displacement 455 diagnostic[0][i] = u[0][i]; 456 diagnostic[1][i] = u[1][i]; 457 diagnostic[2][i] = u[2][i]; 458 459 // Pressure 460 const CeedScalar strain_vol = e[0][0] + e[1][1] + e[2][2]; 461 diagnostic[3][i] = -lambda*strain_vol; 462 463 // Stress tensor invariants 464 diagnostic[4][i] = strain_vol; 465 diagnostic[5][i] = 0.; 466 for (CeedInt j = 0; j < 3; j++) 467 for (CeedInt m = 0; m < 3; m++) 468 diagnostic[5][i] += e[j][m] * e[m][j]; 469 diagnostic[6][i] = 1 + strain_vol; 470 471 // Strain energy 472 diagnostic[7][i] = (lambda*strain_vol*strain_vol/2. + strain_vol*mu + 473 (e[0][1]*e[0][1]+e[0][2]*e[0][2]+e[1][2]*e[1][2])*2*mu); 474 475 } // End of Quadrature Point Loop 476 477 return 0; 478 } 479 // ----------------------------------------------------------------------------- 480 481 #endif //End of ELAS_LINEAR_H 482