13d8e8822SJeremy L Thompson // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 23d8e8822SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 35754ecacSJeremy L Thompson // 43d8e8822SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause 55754ecacSJeremy L Thompson // 63d8e8822SJeremy L Thompson // This file is part of CEED: http://github.com/ceed 75754ecacSJeremy L Thompson 85754ecacSJeremy L Thompson /// @file 95754ecacSJeremy L Thompson /// Linear elasticity manufactured solution forcing term for solid mechanics example using PETSc 105754ecacSJeremy L Thompson 115754ecacSJeremy L Thompson #ifndef MANUFACTURED_H 125754ecacSJeremy L Thompson #define MANUFACTURED_H 135754ecacSJeremy L Thompson 14*c9c2c079SJeremy L Thompson #include <ceed.h> 155754ecacSJeremy L Thompson #include <math.h> 165754ecacSJeremy L Thompson 175754ecacSJeremy L Thompson #ifndef PHYSICS_STRUCT 185754ecacSJeremy L Thompson #define PHYSICS_STRUCT 195754ecacSJeremy L Thompson typedef struct Physics_private *Physics; 205754ecacSJeremy L Thompson struct Physics_private { 215754ecacSJeremy L Thompson CeedScalar nu; // Poisson's ratio 225754ecacSJeremy L Thompson CeedScalar E; // Young's Modulus 235754ecacSJeremy L Thompson }; 245754ecacSJeremy L Thompson #endif 255754ecacSJeremy L Thompson 265754ecacSJeremy L Thompson // ----------------------------------------------------------------------------- 275754ecacSJeremy L Thompson // Forcing term for linear elasticity manufactured solution 285754ecacSJeremy L Thompson // ----------------------------------------------------------------------------- 295754ecacSJeremy L Thompson CEED_QFUNCTION(SetupMMSForce)(void *ctx, const CeedInt Q, 305754ecacSJeremy L Thompson const CeedScalar *const *in, 315754ecacSJeremy L Thompson CeedScalar *const *out) { 325754ecacSJeremy L Thompson // Inputs 335754ecacSJeremy L Thompson const CeedScalar *coords = in[0], *q_data = in[1]; 345754ecacSJeremy L Thompson 355754ecacSJeremy L Thompson // Outputs 365754ecacSJeremy L Thompson CeedScalar *force = out[0]; 375754ecacSJeremy L Thompson 385754ecacSJeremy L Thompson // Context 395754ecacSJeremy L Thompson const Physics context = (Physics)ctx; 405754ecacSJeremy L Thompson const CeedScalar E = context->E; 415754ecacSJeremy L Thompson const CeedScalar nu = context->nu; 425754ecacSJeremy L Thompson 435754ecacSJeremy L Thompson // Quadrature Point Loop 445754ecacSJeremy L Thompson CeedPragmaSIMD 455754ecacSJeremy L Thompson for (CeedInt i=0; i<Q; i++) { 465754ecacSJeremy L Thompson // Setup 475754ecacSJeremy L Thompson CeedScalar x = coords[i+0*Q], y = coords[i+1*Q], z = coords[i+2*Q]; 485754ecacSJeremy L Thompson CeedScalar wdetJ = q_data[i]; 495754ecacSJeremy L Thompson 505754ecacSJeremy L Thompson // Forcing function 515754ecacSJeremy L Thompson // -- Component 1 525754ecacSJeremy L Thompson force[i+0*Q] = (-(E*(cos(x*2.0)*cos(y*3.0)*exp(z*4.0)*4.0 - 535754ecacSJeremy L Thompson cos(z*4.0)*sin( y*3.0)*exp(x*2.0)*8.0)*(nu-0.5))/ 545754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0)) + 555754ecacSJeremy L Thompson (E*(cos(z*4.0)*sin(y*3.0)*exp(x*2.0)*(4.5) + 565754ecacSJeremy L Thompson sin(x*2.0)*sin(z*4.0)*exp( y*3.0)*3.0)*(nu-0.5))/ 575754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0)) + 585754ecacSJeremy L Thompson (E*nu*cos(x*2.0)*cos(y*3.0)*exp(z*4.0)*8.0)/ 595754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0)) - 605754ecacSJeremy L Thompson (E*nu*sin(x*2.0)*sin(z*4.0)*exp(y*3.0)*6.0)/ 615754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0)) - 625754ecacSJeremy L Thompson (E*cos(z*4.0)*sin(y*3.0)*exp(x*2.0)*(nu-1.0)*4.0)/ 635754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0))) * wdetJ / 1e8; 645754ecacSJeremy L Thompson 655754ecacSJeremy L Thompson // -- Component 2 665754ecacSJeremy L Thompson force[i+1*Q] = (-(E*(cos(y*3.0)*cos(z*4.0)*exp(x*2.0)*3.0 - 675754ecacSJeremy L Thompson cos(x*2.0)*sin( z*4.0)*exp(y*3.0)*2.0)*(nu-0.5))/ 685754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0)) + 695754ecacSJeremy L Thompson (E*(cos(x*2.0)*sin(z*4.0)*exp(y*3.0)*8.0 + 705754ecacSJeremy L Thompson sin(x*2.0)*sin(y*3.0)*exp(z*4.0)*6.0)*(nu-0.5))/ 715754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0)) + 725754ecacSJeremy L Thompson (E*nu*cos(y*3.0)*cos(z*4.0)*exp(x*2.0)*6.0)/ 735754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0)) - 745754ecacSJeremy L Thompson (E*nu*sin( x*2.0)*sin(y*3.0)*exp(z*4.0)*12.0)/ 755754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0)) - 765754ecacSJeremy L Thompson (E*cos(x*2.0)*sin(z*4.0)*exp(y*3.0)*(nu-1.0)*9.0)/ 775754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0))) * wdetJ / 1e8; 785754ecacSJeremy L Thompson 795754ecacSJeremy L Thompson // -- Component 3 805754ecacSJeremy L Thompson force[i+2*Q] = (-(E*(cos(x*2.0)*cos(z*4.0)*exp(y*3.0)*6.0 - 815754ecacSJeremy L Thompson cos(y*3.0)*sin( x*2.0)*exp(z*4.0)*(4.5))*(nu-0.5))/ 825754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0)) + 835754ecacSJeremy L Thompson (E*(cos(y*3.0)*sin(x*2.0)*exp(z*4.0)*2.0 + 845754ecacSJeremy L Thompson sin(y*3.0)*sin(z*4.0)*exp(x*2.0)*4.0)*(nu-0.5))/ 855754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0)) + 865754ecacSJeremy L Thompson (E*nu*cos(x*2.0)*cos(z*4.0)*exp(y*3.0)*12.0)/ 875754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0)) - 885754ecacSJeremy L Thompson (E*nu*sin( y*3.0)*sin(z*4.0)*exp(x*2.0)*8.0)/ 895754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0)) - 905754ecacSJeremy L Thompson (E*cos(y*3.0)*sin(x*2.0)*exp(z*4.0)*(nu-1.0)*16.0)/ 915754ecacSJeremy L Thompson ((nu*2.0-1.0)*(nu+1.0))) * wdetJ / 1e8; 925754ecacSJeremy L Thompson 935754ecacSJeremy L Thompson } // End of Quadrature Point Loop 945754ecacSJeremy L Thompson 955754ecacSJeremy L Thompson return 0; 965754ecacSJeremy L Thompson } 975754ecacSJeremy L Thompson // ----------------------------------------------------------------------------- 985754ecacSJeremy L Thompson 995754ecacSJeremy L Thompson #endif // End MANUFACTURED_H 100