1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Geometric factors for solid mechanics example using PETSc 19 20 #ifndef COMMON_H 21 #define COMMON_H 22 23 // ----------------------------------------------------------------------------- 24 // This QFunction sets up the geometric factors required for integration and 25 // coordinate transformations 26 // 27 // Reference (parent) coordinates: X 28 // Physical (current) coordinates: x 29 // Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 30 // Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 31 // 32 // All quadrature data is stored in 10 field vector of quadrature data. 33 // 34 // We require the transpose of the inverse of the Jacobian to properly compute 35 // integrals of the form: int( gradv u ) 36 // 37 // Inverse of Jacobian: 38 // dXdx_i,j = Aij / detJ 39 // 40 // Stored: Aij / detJ 41 // in q_data[1:9] as 42 // [A11 A12 A13] 43 // (detJ^-1) * [A21 A22 A23] 44 // [A31 A32 A33] 45 // 46 // ----------------------------------------------------------------------------- 47 CEED_QFUNCTION(SetupGeo)(void *ctx, CeedInt Q, const CeedScalar *const *in, 48 CeedScalar *const *out) { 49 // *INDENT-OFF* 50 // Inputs 51 const CeedScalar (*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], 52 (*w) = in[1]; 53 54 // Outputs 55 CeedScalar (*q_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 56 // *INDENT-ON* 57 58 CeedPragmaSIMD 59 // Quadrature Point Loop 60 for (CeedInt i=0; i<Q; i++) { 61 // Setup 62 const CeedScalar J11 = J[0][0][i]; 63 const CeedScalar J21 = J[0][1][i]; 64 const CeedScalar J31 = J[0][2][i]; 65 const CeedScalar J12 = J[1][0][i]; 66 const CeedScalar J22 = J[1][1][i]; 67 const CeedScalar J32 = J[1][2][i]; 68 const CeedScalar J13 = J[2][0][i]; 69 const CeedScalar J23 = J[2][1][i]; 70 const CeedScalar J33 = J[2][2][i]; 71 const CeedScalar A11 = J22*J33 - J23*J32; 72 const CeedScalar A12 = J13*J32 - J12*J33; 73 const CeedScalar A13 = J12*J23 - J13*J22; 74 const CeedScalar A21 = J23*J31 - J21*J33; 75 const CeedScalar A22 = J11*J33 - J13*J31; 76 const CeedScalar A23 = J13*J21 - J11*J23; 77 const CeedScalar A31 = J21*J32 - J22*J31; 78 const CeedScalar A32 = J12*J31 - J11*J32; 79 const CeedScalar A33 = J11*J22 - J12*J21; 80 const CeedScalar detJ = J11*A11 + J21*A12 + J31*A13; 81 82 // Qdata 83 // -- Interp-to-Interp q_data 84 q_data[0][i] = w[i] * detJ; 85 86 // -- Interp-to-Grad q_data 87 // Inverse of change of coordinate matrix: X_i,j 88 q_data[1][i] = A11 / detJ; 89 q_data[2][i] = A12 / detJ; 90 q_data[3][i] = A13 / detJ; 91 q_data[4][i] = A21 / detJ; 92 q_data[5][i] = A22 / detJ; 93 q_data[6][i] = A23 / detJ; 94 q_data[7][i] = A31 / detJ; 95 q_data[8][i] = A32 / detJ; 96 q_data[9][i] = A33 / detJ; 97 98 } // End of Quadrature Point Loop 99 100 return 0; 101 } 102 // ----------------------------------------------------------------------------- 103 104 // ----------------------------------------------------------------------------- 105 // This QFunction computes the surface integral of the user traction vector on 106 // the constrained faces. 107 // 108 // Reference (parent) 2D coordinates: X 109 // Physical (current) 3D coordinates: x 110 // Change of coordinate matrix: 111 // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 112 // 113 // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j} 114 // 115 // detJb is the magnitude of (J1,J2,J3) 116 // 117 // Computed: 118 // t * (w detJb) 119 // 120 // ----------------------------------------------------------------------------- 121 CEED_QFUNCTION(SetupTractionBCs)(void *ctx, CeedInt Q, 122 const CeedScalar *const *in, CeedScalar *const *out) { 123 // *INDENT-OFF* 124 // Inputs 125 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], 126 (*w) = in[1]; 127 // Outputs 128 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 129 // *INDENT-ON* 130 131 // User stress tensor 132 const CeedScalar (*traction) = (const CeedScalar(*))ctx; 133 134 CeedPragmaSIMD 135 // Quadrature Point Loop 136 for (CeedInt i = 0; i < Q; i++) { 137 // Setup 138 // *INDENT-OFF* 139 const CeedScalar dxdX[3][2] = {{J[0][0][i], 140 J[1][0][i]}, 141 {J[0][1][i], 142 J[1][1][i]}, 143 {J[0][2][i], 144 J[1][2][i]}}; 145 // *INDENT-ON* 146 // J1, J2, and J3 are given by the cross product of the columns of dxdX 147 const CeedScalar J1 = dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1]; 148 const CeedScalar J2 = dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1]; 149 const CeedScalar J3 = dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]; 150 151 // Qdata 152 // -- Interp-to-Interp q_data 153 CeedScalar wdetJb = w[i] * sqrt(J1 * J1 + J2 * J2 + J3 * J3); 154 155 // Traction surface integral 156 for (CeedInt j = 0; j < 3; j++) 157 v[j][i] = traction[j] * wdetJb; 158 159 } // End of Quadrature Point Loop 160 161 // Return 162 return 0; 163 } 164 // ----------------------------------------------------------------------------- 165 166 #endif // End of COMMON_H 167