1*9ba83ac0SJeremy L Thompson // Copyright (c) 2017-2026, 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
11c0b5abf0SJeremy L Thompson #include <ceed/types.h>
12c0b5abf0SJeremy L Thompson #ifndef CEED_RUNNING_JIT_PASS
135754ecacSJeremy L Thompson #include <math.h>
14c0b5abf0SJeremy 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 // -----------------------------------------------------------------------------
SetupMMSForce(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)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