1 // Copyright (c) 2017-2026, 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 #include <ceed/types.h>
12 #ifndef CEED_RUNNING_JIT_PASS
13 #include <math.h>
14 #endif
15
16 #ifndef PHYSICS_STRUCT
17 #define PHYSICS_STRUCT
18 typedef struct Physics_private *Physics;
19 struct Physics_private {
20 CeedScalar nu; // Poisson's ratio
21 CeedScalar E; // Young's Modulus
22 };
23 #endif
24
25 // -----------------------------------------------------------------------------
26 // Forcing term for linear elasticity manufactured solution
27 // -----------------------------------------------------------------------------
SetupMMSForce(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)28 CEED_QFUNCTION(SetupMMSForce)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
29 // Inputs
30 const CeedScalar *coords = in[0], *q_data = in[1];
31
32 // Outputs
33 CeedScalar *force = out[0];
34
35 // Context
36 const Physics context = (Physics)ctx;
37 const CeedScalar E = context->E;
38 const CeedScalar nu = context->nu;
39
40 // Quadrature Point Loop
41 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
42 // Setup
43 CeedScalar x = coords[i + 0 * Q], y = coords[i + 1 * Q], z = coords[i + 2 * Q];
44 CeedScalar wdetJ = q_data[i];
45
46 // Forcing function
47 // -- Component 1
48 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)) /
49 ((nu * 2.0 - 1.0) * (nu + 1.0)) +
50 (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)) /
51 ((nu * 2.0 - 1.0) * (nu + 1.0)) +
52 (E * nu * cos(x * 2.0) * cos(y * 3.0) * exp(z * 4.0) * 8.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) -
53 (E * nu * sin(x * 2.0) * sin(z * 4.0) * exp(y * 3.0) * 6.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) -
54 (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))) *
55 wdetJ / 1e8;
56
57 // -- Component 2
58 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)) /
59 ((nu * 2.0 - 1.0) * (nu + 1.0)) +
60 (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)) /
61 ((nu * 2.0 - 1.0) * (nu + 1.0)) +
62 (E * nu * cos(y * 3.0) * cos(z * 4.0) * exp(x * 2.0) * 6.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) -
63 (E * nu * sin(x * 2.0) * sin(y * 3.0) * exp(z * 4.0) * 12.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) -
64 (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))) *
65 wdetJ / 1e8;
66
67 // -- Component 3
68 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)) /
69 ((nu * 2.0 - 1.0) * (nu + 1.0)) +
70 (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)) /
71 ((nu * 2.0 - 1.0) * (nu + 1.0)) +
72 (E * nu * cos(x * 2.0) * cos(z * 4.0) * exp(y * 3.0) * 12.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) -
73 (E * nu * sin(y * 3.0) * sin(z * 4.0) * exp(x * 2.0) * 8.0) / ((nu * 2.0 - 1.0) * (nu + 1.0)) -
74 (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))) *
75 wdetJ / 1e8;
76
77 } // End of Quadrature Point Loop
78
79 return 0;
80 }
81 // -----------------------------------------------------------------------------
82