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 manufactured solution forcing term for solid mechanics example using PETSc 10 11 #ifndef MANUFACTURED_H 12 #define MANUFACTURED_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 // Forcing term for linear elasticity manufactured solution 28 // ----------------------------------------------------------------------------- 29 CEED_QFUNCTION(SetupMMSForce)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 30 // Inputs 31 const CeedScalar *coords = in[0], *q_data = in[1]; 32 33 // Outputs 34 CeedScalar *force = out[0]; 35 36 // Context 37 const Physics context = (Physics)ctx; 38 const CeedScalar E = context->E; 39 const CeedScalar nu = context->nu; 40 41 // Quadrature Point Loop 42 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 43 // Setup 44 CeedScalar x = coords[i + 0 * Q], y = coords[i + 1 * Q], z = coords[i + 2 * Q]; 45 CeedScalar wdetJ = q_data[i]; 46 47 // Forcing function 48 // -- Component 1 49 force[i + 0 * Q] = (-(E * (cos(x * 2.0) * cos(y * 3.0) * exp(z * 4.0) * 4.0 - cos(z * 4.0) * sin(y * 3.0) * exp(x * 2.0) * 8.0) * (nu - 0.5)) / 50 ((nu * 2.0 - 1.0) * (nu + 1.0)) + 51 (E * (cos(z * 4.0) * sin(y * 3.0) * exp(x * 2.0) * (4.5) + sin(x * 2.0) * sin(z * 4.0) * exp(y * 3.0) * 3.0) * (nu - 0.5)) / 52 ((nu * 2.0 - 1.0) * (nu + 1.0)) + 53 (E * nu * cos(x * 2.0) * cos(y * 3.0) * exp(z * 4.0) * 8.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) - 54 (E * nu * sin(x * 2.0) * sin(z * 4.0) * exp(y * 3.0) * 6.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) - 55 (E * cos(z * 4.0) * sin(y * 3.0) * exp(x * 2.0) * (nu - 1.0) * 4.0) / ((nu * 2.0 - 1.0) * (nu + 1.0))) * 56 wdetJ / 1e8; 57 58 // -- Component 2 59 force[i + 1 * Q] = (-(E * (cos(y * 3.0) * cos(z * 4.0) * exp(x * 2.0) * 3.0 - cos(x * 2.0) * sin(z * 4.0) * exp(y * 3.0) * 2.0) * (nu - 0.5)) / 60 ((nu * 2.0 - 1.0) * (nu + 1.0)) + 61 (E * (cos(x * 2.0) * sin(z * 4.0) * exp(y * 3.0) * 8.0 + sin(x * 2.0) * sin(y * 3.0) * exp(z * 4.0) * 6.0) * (nu - 0.5)) / 62 ((nu * 2.0 - 1.0) * (nu + 1.0)) + 63 (E * nu * cos(y * 3.0) * cos(z * 4.0) * exp(x * 2.0) * 6.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) - 64 (E * nu * sin(x * 2.0) * sin(y * 3.0) * exp(z * 4.0) * 12.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) - 65 (E * cos(x * 2.0) * sin(z * 4.0) * exp(y * 3.0) * (nu - 1.0) * 9.0) / ((nu * 2.0 - 1.0) * (nu + 1.0))) * 66 wdetJ / 1e8; 67 68 // -- Component 3 69 force[i + 2 * Q] = (-(E * (cos(x * 2.0) * cos(z * 4.0) * exp(y * 3.0) * 6.0 - cos(y * 3.0) * sin(x * 2.0) * exp(z * 4.0) * (4.5)) * (nu - 0.5)) / 70 ((nu * 2.0 - 1.0) * (nu + 1.0)) + 71 (E * (cos(y * 3.0) * sin(x * 2.0) * exp(z * 4.0) * 2.0 + sin(y * 3.0) * sin(z * 4.0) * exp(x * 2.0) * 4.0) * (nu - 0.5)) / 72 ((nu * 2.0 - 1.0) * (nu + 1.0)) + 73 (E * nu * cos(x * 2.0) * cos(z * 4.0) * exp(y * 3.0) * 12.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) - 74 (E * nu * sin(y * 3.0) * sin(z * 4.0) * exp(x * 2.0) * 8.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) - 75 (E * cos(y * 3.0) * sin(x * 2.0) * exp(z * 4.0) * (nu - 1.0) * 16.0) / ((nu * 2.0 - 1.0) * (nu + 1.0))) * 76 wdetJ / 1e8; 77 78 } // End of Quadrature Point Loop 79 80 return 0; 81 } 82 // ----------------------------------------------------------------------------- 83 84 #endif // End MANUFACTURED_H 85