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