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