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