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 <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, 30 const CeedScalar *const *in, 31 CeedScalar *const *out) { 32 // Inputs 33 const CeedScalar *coords = in[0], *q_data = in[1]; 34 35 // Outputs 36 CeedScalar *force = out[0]; 37 38 // Context 39 const Physics context = (Physics)ctx; 40 const CeedScalar E = context->E; 41 const CeedScalar nu = context->nu; 42 43 // Quadrature Point Loop 44 CeedPragmaSIMD 45 for (CeedInt i=0; i<Q; i++) { 46 // Setup 47 CeedScalar x = coords[i+0*Q], y = coords[i+1*Q], z = coords[i+2*Q]; 48 CeedScalar wdetJ = q_data[i]; 49 50 // Forcing function 51 // -- Component 1 52 force[i+0*Q] = (-(E*(cos(x*2.0)*cos(y*3.0)*exp(z*4.0)*4.0 - 53 cos(z*4.0)*sin( y*3.0)*exp(x*2.0)*8.0)*(nu-0.5))/ 54 ((nu*2.0-1.0)*(nu+1.0)) + 55 (E*(cos(z*4.0)*sin(y*3.0)*exp(x*2.0)*(4.5) + 56 sin(x*2.0)*sin(z*4.0)*exp( y*3.0)*3.0)*(nu-0.5))/ 57 ((nu*2.0-1.0)*(nu+1.0)) + 58 (E*nu*cos(x*2.0)*cos(y*3.0)*exp(z*4.0)*8.0)/ 59 ((nu*2.0-1.0)*(nu+1.0)) - 60 (E*nu*sin(x*2.0)*sin(z*4.0)*exp(y*3.0)*6.0)/ 61 ((nu*2.0-1.0)*(nu+1.0)) - 62 (E*cos(z*4.0)*sin(y*3.0)*exp(x*2.0)*(nu-1.0)*4.0)/ 63 ((nu*2.0-1.0)*(nu+1.0))) * wdetJ / 1e8; 64 65 // -- Component 2 66 force[i+1*Q] = (-(E*(cos(y*3.0)*cos(z*4.0)*exp(x*2.0)*3.0 - 67 cos(x*2.0)*sin( z*4.0)*exp(y*3.0)*2.0)*(nu-0.5))/ 68 ((nu*2.0-1.0)*(nu+1.0)) + 69 (E*(cos(x*2.0)*sin(z*4.0)*exp(y*3.0)*8.0 + 70 sin(x*2.0)*sin(y*3.0)*exp(z*4.0)*6.0)*(nu-0.5))/ 71 ((nu*2.0-1.0)*(nu+1.0)) + 72 (E*nu*cos(y*3.0)*cos(z*4.0)*exp(x*2.0)*6.0)/ 73 ((nu*2.0-1.0)*(nu+1.0)) - 74 (E*nu*sin( x*2.0)*sin(y*3.0)*exp(z*4.0)*12.0)/ 75 ((nu*2.0-1.0)*(nu+1.0)) - 76 (E*cos(x*2.0)*sin(z*4.0)*exp(y*3.0)*(nu-1.0)*9.0)/ 77 ((nu*2.0-1.0)*(nu+1.0))) * wdetJ / 1e8; 78 79 // -- Component 3 80 force[i+2*Q] = (-(E*(cos(x*2.0)*cos(z*4.0)*exp(y*3.0)*6.0 - 81 cos(y*3.0)*sin( x*2.0)*exp(z*4.0)*(4.5))*(nu-0.5))/ 82 ((nu*2.0-1.0)*(nu+1.0)) + 83 (E*(cos(y*3.0)*sin(x*2.0)*exp(z*4.0)*2.0 + 84 sin(y*3.0)*sin(z*4.0)*exp(x*2.0)*4.0)*(nu-0.5))/ 85 ((nu*2.0-1.0)*(nu+1.0)) + 86 (E*nu*cos(x*2.0)*cos(z*4.0)*exp(y*3.0)*12.0)/ 87 ((nu*2.0-1.0)*(nu+1.0)) - 88 (E*nu*sin( y*3.0)*sin(z*4.0)*exp(x*2.0)*8.0)/ 89 ((nu*2.0-1.0)*(nu+1.0)) - 90 (E*cos(y*3.0)*sin(x*2.0)*exp(z*4.0)*(nu-1.0)*16.0)/ 91 ((nu*2.0-1.0)*(nu+1.0))) * wdetJ / 1e8; 92 93 } // End of Quadrature Point Loop 94 95 return 0; 96 } 97 // ----------------------------------------------------------------------------- 98 99 #endif // End MANUFACTURED_H 100