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