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