1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Linear elasticity manufactured solution forcing term for solid mechanics example using PETSc 19 20 #ifndef MANUFACTURED_H 21 #define MANUFACTURED_H 22 23 #include <math.h> 24 25 #ifndef PHYSICS_STRUCT 26 #define PHYSICS_STRUCT 27 typedef struct Physics_private *Physics; 28 struct Physics_private { 29 CeedScalar nu; // Poisson's ratio 30 CeedScalar E; // Young's Modulus 31 }; 32 #endif 33 34 // ----------------------------------------------------------------------------- 35 // Forcing term for linear elasticity manufactured solution 36 // ----------------------------------------------------------------------------- 37 CEED_QFUNCTION(SetupMMSForce)(void *ctx, const CeedInt Q, 38 const CeedScalar *const *in, 39 CeedScalar *const *out) { 40 // Inputs 41 const CeedScalar *coords = in[0], *q_data = in[1]; 42 43 // Outputs 44 CeedScalar *force = out[0]; 45 46 // Context 47 const Physics context = (Physics)ctx; 48 const CeedScalar E = context->E; 49 const CeedScalar nu = context->nu; 50 51 // Quadrature Point Loop 52 CeedPragmaSIMD 53 for (CeedInt i=0; i<Q; i++) { 54 // Setup 55 CeedScalar x = coords[i+0*Q], y = coords[i+1*Q], z = coords[i+2*Q]; 56 CeedScalar wdetJ = q_data[i]; 57 58 // Forcing function 59 // -- Component 1 60 force[i+0*Q] = (-(E*(cos(x*2.0)*cos(y*3.0)*exp(z*4.0)*4.0 - 61 cos(z*4.0)*sin( y*3.0)*exp(x*2.0)*8.0)*(nu-0.5))/ 62 ((nu*2.0-1.0)*(nu+1.0)) + 63 (E*(cos(z*4.0)*sin(y*3.0)*exp(x*2.0)*(4.5) + 64 sin(x*2.0)*sin(z*4.0)*exp( y*3.0)*3.0)*(nu-0.5))/ 65 ((nu*2.0-1.0)*(nu+1.0)) + 66 (E*nu*cos(x*2.0)*cos(y*3.0)*exp(z*4.0)*8.0)/ 67 ((nu*2.0-1.0)*(nu+1.0)) - 68 (E*nu*sin(x*2.0)*sin(z*4.0)*exp(y*3.0)*6.0)/ 69 ((nu*2.0-1.0)*(nu+1.0)) - 70 (E*cos(z*4.0)*sin(y*3.0)*exp(x*2.0)*(nu-1.0)*4.0)/ 71 ((nu*2.0-1.0)*(nu+1.0))) * wdetJ / 1e8; 72 73 // -- Component 2 74 force[i+1*Q] = (-(E*(cos(y*3.0)*cos(z*4.0)*exp(x*2.0)*3.0 - 75 cos(x*2.0)*sin( z*4.0)*exp(y*3.0)*2.0)*(nu-0.5))/ 76 ((nu*2.0-1.0)*(nu+1.0)) + 77 (E*(cos(x*2.0)*sin(z*4.0)*exp(y*3.0)*8.0 + 78 sin(x*2.0)*sin(y*3.0)*exp(z*4.0)*6.0)*(nu-0.5))/ 79 ((nu*2.0-1.0)*(nu+1.0)) + 80 (E*nu*cos(y*3.0)*cos(z*4.0)*exp(x*2.0)*6.0)/ 81 ((nu*2.0-1.0)*(nu+1.0)) - 82 (E*nu*sin( x*2.0)*sin(y*3.0)*exp(z*4.0)*12.0)/ 83 ((nu*2.0-1.0)*(nu+1.0)) - 84 (E*cos(x*2.0)*sin(z*4.0)*exp(y*3.0)*(nu-1.0)*9.0)/ 85 ((nu*2.0-1.0)*(nu+1.0))) * wdetJ / 1e8; 86 87 // -- Component 3 88 force[i+2*Q] = (-(E*(cos(x*2.0)*cos(z*4.0)*exp(y*3.0)*6.0 - 89 cos(y*3.0)*sin( x*2.0)*exp(z*4.0)*(4.5))*(nu-0.5))/ 90 ((nu*2.0-1.0)*(nu+1.0)) + 91 (E*(cos(y*3.0)*sin(x*2.0)*exp(z*4.0)*2.0 + 92 sin(y*3.0)*sin(z*4.0)*exp(x*2.0)*4.0)*(nu-0.5))/ 93 ((nu*2.0-1.0)*(nu+1.0)) + 94 (E*nu*cos(x*2.0)*cos(z*4.0)*exp(y*3.0)*12.0)/ 95 ((nu*2.0-1.0)*(nu+1.0)) - 96 (E*nu*sin( y*3.0)*sin(z*4.0)*exp(x*2.0)*8.0)/ 97 ((nu*2.0-1.0)*(nu+1.0)) - 98 (E*cos(y*3.0)*sin(x*2.0)*exp(z*4.0)*(nu-1.0)*16.0)/ 99 ((nu*2.0-1.0)*(nu+1.0))) * wdetJ / 1e8; 100 101 } // End of Quadrature Point Loop 102 103 return 0; 104 } 105 // ----------------------------------------------------------------------------- 106 107 #endif // End MANUFACTURED_H 108