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