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