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