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 sets up the geometric factors required for integration and coordinate transformations
15 //
16 // Reference (parent) coordinates: X
17 // Physical (current) coordinates: x
18 // Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation)
19 // Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j}
20 //
21 // All quadrature data is stored in 10 field vector of quadrature data.
22 //
23 // We require the transpose of the inverse of the Jacobian to properly compute integrals of the form: int( gradv u )
24 //
25 // Inverse of Jacobian:
26 // dXdx_i,j = Aij / detJ
27 //
28 // Stored: Aij / detJ
29 // in q_data[1:9] as
30 // [A11 A12 A13]
31 // (detJ^-1) * [A21 A22 A23]
32 // [A31 A32 A33]
33 // -----------------------------------------------------------------------------
SetupGeo(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)34 CEED_QFUNCTION(SetupGeo)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
35 // Inputs
36 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], (*w) = in[1];
37
38 // Outputs
39 CeedScalar(*q_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
40
41 CeedPragmaSIMD
42 // Quadrature Point Loop
43 for (CeedInt i = 0; i < Q; i++) {
44 // Setup
45 const CeedScalar J11 = J[0][0][i];
46 const CeedScalar J21 = J[0][1][i];
47 const CeedScalar J31 = J[0][2][i];
48 const CeedScalar J12 = J[1][0][i];
49 const CeedScalar J22 = J[1][1][i];
50 const CeedScalar J32 = J[1][2][i];
51 const CeedScalar J13 = J[2][0][i];
52 const CeedScalar J23 = J[2][1][i];
53 const CeedScalar J33 = J[2][2][i];
54 const CeedScalar A11 = J22 * J33 - J23 * J32;
55 const CeedScalar A12 = J13 * J32 - J12 * J33;
56 const CeedScalar A13 = J12 * J23 - J13 * J22;
57 const CeedScalar A21 = J23 * J31 - J21 * J33;
58 const CeedScalar A22 = J11 * J33 - J13 * J31;
59 const CeedScalar A23 = J13 * J21 - J11 * J23;
60 const CeedScalar A31 = J21 * J32 - J22 * J31;
61 const CeedScalar A32 = J12 * J31 - J11 * J32;
62 const CeedScalar A33 = J11 * J22 - J12 * J21;
63 const CeedScalar detJ = J11 * A11 + J21 * A12 + J31 * A13;
64
65 // Qdata
66 // -- Interp-to-Interp q_data
67 q_data[0][i] = w[i] * detJ;
68
69 // -- Interp-to-Grad q_data
70 // Inverse of change of coordinate matrix: X_i,j
71 q_data[1][i] = A11 / detJ;
72 q_data[2][i] = A12 / detJ;
73 q_data[3][i] = A13 / detJ;
74 q_data[4][i] = A21 / detJ;
75 q_data[5][i] = A22 / detJ;
76 q_data[6][i] = A23 / detJ;
77 q_data[7][i] = A31 / detJ;
78 q_data[8][i] = A32 / detJ;
79 q_data[9][i] = A33 / detJ;
80
81 } // End of Quadrature Point Loop
82
83 return 0;
84 }
85 // -----------------------------------------------------------------------------
86