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 /// Geometric factors for solid mechanics example using PETSc
10
11 #include <ceed/types.h>
12
13 // -----------------------------------------------------------------------------
14 // This QFunction computes the surface integral of the user traction vector on the constrained faces.
15 //
16 // Reference (parent) 2D coordinates: X
17 // Physical (current) 3D coordinates: x
18 // Change of coordinate matrix:
19 // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2]
20 //
21 // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j}
22 //
23 // detJb is the magnitude of (J1,J2,J3)
24 //
25 // Computed:
26 // t * (w detJb)
27 // -----------------------------------------------------------------------------
SetupTractionBCs(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)28 CEED_QFUNCTION(SetupTractionBCs)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
29 // Inputs
30 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], (*w) = in[1];
31 // Outputs
32 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
33
34 // User stress tensor
35 const CeedScalar(*traction) = (const CeedScalar(*))ctx;
36
37 CeedPragmaSIMD
38 // Quadrature Point Loop
39 for (CeedInt i = 0; i < Q; i++) {
40 // Setup
41 const CeedScalar dxdX[3][2] = {
42 {J[0][0][i], J[1][0][i]},
43 {J[0][1][i], J[1][1][i]},
44 {J[0][2][i], J[1][2][i]}
45 };
46 // J1, J2, and J3 are given by the cross product of the columns of dxdX
47 const CeedScalar J1 = dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1];
48 const CeedScalar J2 = dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1];
49 const CeedScalar J3 = dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1];
50
51 // Qdata
52 // -- Interp-to-Interp q_data
53 CeedScalar wdetJb = w[i] * sqrt(J1 * J1 + J2 * J2 + J3 * J3);
54
55 // Traction surface integral
56 for (CeedInt j = 0; j < 3; j++) v[j][i] = traction[j] * wdetJb;
57
58 } // End of Quadrature Point Loop
59
60 // Return
61 return 0;
62 }
63 // -----------------------------------------------------------------------------
64