15aed82e4SJeremy L Thompson // Copyright (c) 2017-2024, 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 11*c0b5abf0SJeremy L Thompson #include <ceed/types.h> 12*c0b5abf0SJeremy L Thompson #ifndef CEED_RUNNING_JIT_PASS 135754ecacSJeremy L Thompson #include <math.h> 14*c0b5abf0SJeremy L Thompson #endif 155754ecacSJeremy L Thompson 165754ecacSJeremy L Thompson #ifndef PHYSICS_STRUCT 175754ecacSJeremy L Thompson #define PHYSICS_STRUCT 185754ecacSJeremy L Thompson typedef struct Physics_private *Physics; 195754ecacSJeremy L Thompson struct Physics_private { 205754ecacSJeremy L Thompson CeedScalar nu; // Poisson's ratio 215754ecacSJeremy L Thompson CeedScalar E; // Young's Modulus 225754ecacSJeremy L Thompson }; 235754ecacSJeremy L Thompson #endif 245754ecacSJeremy L Thompson 255754ecacSJeremy L Thompson // ----------------------------------------------------------------------------- 265754ecacSJeremy L Thompson // Forcing term for linear elasticity manufactured solution 275754ecacSJeremy L Thompson // ----------------------------------------------------------------------------- 282b730f8bSJeremy L Thompson CEED_QFUNCTION(SetupMMSForce)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 295754ecacSJeremy L Thompson // Inputs 305754ecacSJeremy L Thompson const CeedScalar *coords = in[0], *q_data = in[1]; 315754ecacSJeremy L Thompson 325754ecacSJeremy L Thompson // Outputs 335754ecacSJeremy L Thompson CeedScalar *force = out[0]; 345754ecacSJeremy L Thompson 355754ecacSJeremy L Thompson // Context 365754ecacSJeremy L Thompson const Physics context = (Physics)ctx; 375754ecacSJeremy L Thompson const CeedScalar E = context->E; 385754ecacSJeremy L Thompson const CeedScalar nu = context->nu; 395754ecacSJeremy L Thompson 405754ecacSJeremy L Thompson // Quadrature Point Loop 412b730f8bSJeremy L Thompson CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 425754ecacSJeremy L Thompson // Setup 435754ecacSJeremy L Thompson CeedScalar x = coords[i + 0 * Q], y = coords[i + 1 * Q], z = coords[i + 2 * Q]; 445754ecacSJeremy L Thompson CeedScalar wdetJ = q_data[i]; 455754ecacSJeremy L Thompson 465754ecacSJeremy L Thompson // Forcing function 475754ecacSJeremy L Thompson // -- Component 1 482b730f8bSJeremy L Thompson force[i + 0 * Q] = (-(E * (cos(x * 2.0) * cos(y * 3.0) * exp(z * 4.0) * 4.0 - cos(z * 4.0) * sin(y * 3.0) * exp(x * 2.0) * 8.0) * (nu - 0.5)) / 495754ecacSJeremy L Thompson ((nu * 2.0 - 1.0) * (nu + 1.0)) + 502b730f8bSJeremy L Thompson (E * (cos(z * 4.0) * sin(y * 3.0) * exp(x * 2.0) * (4.5) + sin(x * 2.0) * sin(z * 4.0) * exp(y * 3.0) * 3.0) * (nu - 0.5)) / 515754ecacSJeremy L Thompson ((nu * 2.0 - 1.0) * (nu + 1.0)) + 522b730f8bSJeremy L Thompson (E * nu * cos(x * 2.0) * cos(y * 3.0) * exp(z * 4.0) * 8.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) - 532b730f8bSJeremy L Thompson (E * nu * sin(x * 2.0) * sin(z * 4.0) * exp(y * 3.0) * 6.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) - 542b730f8bSJeremy L Thompson (E * cos(z * 4.0) * sin(y * 3.0) * exp(x * 2.0) * (nu - 1.0) * 4.0) / ((nu * 2.0 - 1.0) * (nu + 1.0))) * 552b730f8bSJeremy L Thompson wdetJ / 1e8; 565754ecacSJeremy L Thompson 575754ecacSJeremy L Thompson // -- Component 2 582b730f8bSJeremy L Thompson force[i + 1 * Q] = (-(E * (cos(y * 3.0) * cos(z * 4.0) * exp(x * 2.0) * 3.0 - cos(x * 2.0) * sin(z * 4.0) * exp(y * 3.0) * 2.0) * (nu - 0.5)) / 595754ecacSJeremy L Thompson ((nu * 2.0 - 1.0) * (nu + 1.0)) + 602b730f8bSJeremy L Thompson (E * (cos(x * 2.0) * sin(z * 4.0) * exp(y * 3.0) * 8.0 + sin(x * 2.0) * sin(y * 3.0) * exp(z * 4.0) * 6.0) * (nu - 0.5)) / 615754ecacSJeremy L Thompson ((nu * 2.0 - 1.0) * (nu + 1.0)) + 622b730f8bSJeremy L Thompson (E * nu * cos(y * 3.0) * cos(z * 4.0) * exp(x * 2.0) * 6.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) - 632b730f8bSJeremy L Thompson (E * nu * sin(x * 2.0) * sin(y * 3.0) * exp(z * 4.0) * 12.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) - 642b730f8bSJeremy L Thompson (E * cos(x * 2.0) * sin(z * 4.0) * exp(y * 3.0) * (nu - 1.0) * 9.0) / ((nu * 2.0 - 1.0) * (nu + 1.0))) * 652b730f8bSJeremy L Thompson wdetJ / 1e8; 665754ecacSJeremy L Thompson 675754ecacSJeremy L Thompson // -- Component 3 682b730f8bSJeremy L Thompson force[i + 2 * Q] = (-(E * (cos(x * 2.0) * cos(z * 4.0) * exp(y * 3.0) * 6.0 - cos(y * 3.0) * sin(x * 2.0) * exp(z * 4.0) * (4.5)) * (nu - 0.5)) / 695754ecacSJeremy L Thompson ((nu * 2.0 - 1.0) * (nu + 1.0)) + 702b730f8bSJeremy L Thompson (E * (cos(y * 3.0) * sin(x * 2.0) * exp(z * 4.0) * 2.0 + sin(y * 3.0) * sin(z * 4.0) * exp(x * 2.0) * 4.0) * (nu - 0.5)) / 715754ecacSJeremy L Thompson ((nu * 2.0 - 1.0) * (nu + 1.0)) + 722b730f8bSJeremy L Thompson (E * nu * cos(x * 2.0) * cos(z * 4.0) * exp(y * 3.0) * 12.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) - 732b730f8bSJeremy L Thompson (E * nu * sin(y * 3.0) * sin(z * 4.0) * exp(x * 2.0) * 8.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) - 742b730f8bSJeremy L Thompson (E * cos(y * 3.0) * sin(x * 2.0) * exp(z * 4.0) * (nu - 1.0) * 16.0) / ((nu * 2.0 - 1.0) * (nu + 1.0))) * 752b730f8bSJeremy L Thompson wdetJ / 1e8; 765754ecacSJeremy L Thompson 775754ecacSJeremy L Thompson } // End of Quadrature Point Loop 785754ecacSJeremy L Thompson 795754ecacSJeremy L Thompson return 0; 805754ecacSJeremy L Thompson } 815754ecacSJeremy L Thompson // ----------------------------------------------------------------------------- 82