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