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 #ifndef __CUDACC__ 24 # include <math.h> 25 #endif 26 27 #ifndef PHYSICS_STRUCT 28 #define PHYSICS_STRUCT 29 typedef struct Physics_private *Physics; 30 struct Physics_private { 31 CeedScalar nu; // Poisson's ratio 32 CeedScalar E; // Young's Modulus 33 }; 34 #endif 35 36 // ----------------------------------------------------------------------------- 37 // Forcing term for linear elasticity manufactured solution 38 // ----------------------------------------------------------------------------- 39 CEED_QFUNCTION(SetupMMSForce)(void *ctx, const CeedInt Q, 40 const CeedScalar *const *in, 41 CeedScalar *const *out) { 42 // Inputs 43 const CeedScalar *coords = in[0], *q_data = in[1]; 44 45 // Outputs 46 CeedScalar *force = out[0]; 47 48 // Context 49 const Physics context = (Physics)ctx; 50 const CeedScalar E = context->E; 51 const CeedScalar nu = context->nu; 52 53 // Quadrature Point Loop 54 CeedPragmaSIMD 55 for (CeedInt i=0; i<Q; i++) { 56 // Setup 57 CeedScalar x = coords[i+0*Q], y = coords[i+1*Q], z = coords[i+2*Q]; 58 CeedScalar wdetJ = q_data[i]; 59 60 // Forcing function 61 // -- Component 1 62 force[i+0*Q] = (-(E*(cos(x*2.0)*cos(y*3.0)*exp(z*4.0)*4.0 - 63 cos(z*4.0)*sin( y*3.0)*exp(x*2.0)*8.0)*(nu-0.5))/ 64 ((nu*2.0-1.0)*(nu+1.0)) + 65 (E*(cos(z*4.0)*sin(y*3.0)*exp(x*2.0)*(4.5) + 66 sin(x*2.0)*sin(z*4.0)*exp( y*3.0)*3.0)*(nu-0.5))/ 67 ((nu*2.0-1.0)*(nu+1.0)) + 68 (E*nu*cos(x*2.0)*cos(y*3.0)*exp(z*4.0)*8.0)/ 69 ((nu*2.0-1.0)*(nu+1.0)) - 70 (E*nu*sin(x*2.0)*sin(z*4.0)*exp(y*3.0)*6.0)/ 71 ((nu*2.0-1.0)*(nu+1.0)) - 72 (E*cos(z*4.0)*sin(y*3.0)*exp(x*2.0)*(nu-1.0)*4.0)/ 73 ((nu*2.0-1.0)*(nu+1.0))) * wdetJ / 1e8; 74 75 // -- Component 2 76 force[i+1*Q] = (-(E*(cos(y*3.0)*cos(z*4.0)*exp(x*2.0)*3.0 - 77 cos(x*2.0)*sin( z*4.0)*exp(y*3.0)*2.0)*(nu-0.5))/ 78 ((nu*2.0-1.0)*(nu+1.0)) + 79 (E*(cos(x*2.0)*sin(z*4.0)*exp(y*3.0)*8.0 + 80 sin(x*2.0)*sin(y*3.0)*exp(z*4.0)*6.0)*(nu-0.5))/ 81 ((nu*2.0-1.0)*(nu+1.0)) + 82 (E*nu*cos(y*3.0)*cos(z*4.0)*exp(x*2.0)*6.0)/ 83 ((nu*2.0-1.0)*(nu+1.0)) - 84 (E*nu*sin( x*2.0)*sin(y*3.0)*exp(z*4.0)*12.0)/ 85 ((nu*2.0-1.0)*(nu+1.0)) - 86 (E*cos(x*2.0)*sin(z*4.0)*exp(y*3.0)*(nu-1.0)*9.0)/ 87 ((nu*2.0-1.0)*(nu+1.0))) * wdetJ / 1e8; 88 89 // -- Component 3 90 force[i+2*Q] = (-(E*(cos(x*2.0)*cos(z*4.0)*exp(y*3.0)*6.0 - 91 cos(y*3.0)*sin( x*2.0)*exp(z*4.0)*(4.5))*(nu-0.5))/ 92 ((nu*2.0-1.0)*(nu+1.0)) + 93 (E*(cos(y*3.0)*sin(x*2.0)*exp(z*4.0)*2.0 + 94 sin(y*3.0)*sin(z*4.0)*exp(x*2.0)*4.0)*(nu-0.5))/ 95 ((nu*2.0-1.0)*(nu+1.0)) + 96 (E*nu*cos(x*2.0)*cos(z*4.0)*exp(y*3.0)*12.0)/ 97 ((nu*2.0-1.0)*(nu+1.0)) - 98 (E*nu*sin( y*3.0)*sin(z*4.0)*exp(x*2.0)*8.0)/ 99 ((nu*2.0-1.0)*(nu+1.0)) - 100 (E*cos(y*3.0)*sin(x*2.0)*exp(z*4.0)*(nu-1.0)*16.0)/ 101 ((nu*2.0-1.0)*(nu+1.0))) * wdetJ / 1e8; 102 103 } // End of Quadrature Point Loop 104 105 return 0; 106 } 107 // ----------------------------------------------------------------------------- 108 109 #endif // End MANUFACTURED_H 110