1 // Copyright (c) 2017-2023, 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 (3D) for Navier-Stokes example using PETSc 10 11 #ifndef setup_geo_h 12 #define setup_geo_h 13 14 #include <ceed.h> 15 #include <math.h> 16 17 #include "setupgeo_helpers.h" 18 19 // ***************************************************************************** 20 // This QFunction sets up the geometric factors required for integration and coordinate transformations 21 // 22 // Reference (parent) coordinates: X 23 // Physical (current) coordinates: x 24 // Change of coordinate matrix: dxdX_{i,j} = x_{i,j} (indicial notation) 25 // Inverse of change of coordinate matrix: dXdx_{i,j} = (detJ^-1) * X_{i,j} 26 // 27 // All quadrature data is stored in 10 field vector of quadrature data. 28 // 29 // We require the determinant of the Jacobian to properly compute integrals of the form: int( v u ) 30 // 31 // Determinant of Jacobian: 32 // detJ = J11*A11 + J21*A12 + J31*A13 33 // Jij = Jacobian entry ij 34 // Aij = Adjugate ij 35 // 36 // Stored: w detJ 37 // in q_data[0] 38 // 39 // We require the transpose of the inverse of the Jacobian to properly compute integrals of the form: int( gradv u ) 40 // 41 // Inverse of Jacobian: 42 // dXdx_i,j = Aij / detJ 43 // 44 // Stored: Aij / detJ 45 // in q_data[1:9] as 46 // (detJ^-1) * [A11 A12 A13] 47 // [A21 A22 A23] 48 // [A31 A32 A33] 49 // ***************************************************************************** 50 CEED_QFUNCTION(Setup)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 51 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0]; 52 const CeedScalar(*w) = in[1]; 53 CeedScalar(*q_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 54 55 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 56 CeedScalar detJ, dXdx[3][3]; 57 InvertMappingJacobian_3D(Q, i, J, dXdx, &detJ); 58 q_data[0][i] = w[i] * detJ; 59 q_data[1][i] = dXdx[0][0]; 60 q_data[2][i] = dXdx[0][1]; 61 q_data[3][i] = dXdx[0][2]; 62 q_data[4][i] = dXdx[1][0]; 63 q_data[5][i] = dXdx[1][1]; 64 q_data[6][i] = dXdx[1][2]; 65 q_data[7][i] = dXdx[2][0]; 66 q_data[8][i] = dXdx[2][1]; 67 q_data[9][i] = dXdx[2][2]; 68 } 69 return 0; 70 } 71 72 // ***************************************************************************** 73 // This QFunction sets up the geometric factor required for integration when reference coordinates are in 2D and the physical coordinates are in 3D 74 // 75 // Reference (parent) 2D coordinates: X 76 // Physical (current) 3D coordinates: x 77 // Change of coordinate matrix: 78 // dxdX_{i,j} = dx_i/dX_j (indicial notation) [3 * 2] 79 // Inverse change of coordinate matrix: 80 // dXdx_{i,j} = dX_i/dx_j (indicial notation) [2 * 3] 81 // 82 // (J1,J2,J3) is given by the cross product of the columns of dxdX_{i,j} 83 // 84 // detJb is the magnitude of (J1,J2,J3) 85 // 86 // dXdx is calculated via Moore–Penrose inverse: 87 // 88 // dX_i/dx_j = (dxdX^T dxdX)^(-1) dxdX 89 // = (dx_l/dX_i * dx_l/dX_k)^(-1) dx_j/dX_k 90 // 91 // All quadrature data is stored in 10 field vector of quadrature data. 92 // 93 // We require the determinant of the Jacobian to properly compute integrals of 94 // the form: int( u v ) 95 // 96 // Stored: w detJb 97 // in q_data_sur[0] 98 // 99 // Normal vector = (J1,J2,J3) / detJb 100 // 101 // - TODO Could possibly remove normal vector, as it could be calculated in the Qfunction from dXdx 102 // See https://github.com/CEED/libCEED/pull/868#discussion_r871979484 103 // Stored: (J1,J2,J3) / detJb 104 // in q_data_sur[1:3] as 105 // (detJb^-1) * [ J1 ] 106 // [ J2 ] 107 // [ J3 ] 108 // 109 // Stored: dXdx_{i,j} 110 // in q_data_sur[4:9] as 111 // [dXdx_11 dXdx_12 dXdx_13] 112 // [dXdx_21 dXdx_22 dXdx_23] 113 // ***************************************************************************** 114 CEED_QFUNCTION(SetupBoundary)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 115 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0]; 116 const CeedScalar(*w) = in[1]; 117 CeedScalar(*q_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 118 119 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 120 CeedScalar detJb, normal[3], dXdx[2][3]; 121 122 NormalVectorFromdxdX_3D(Q, i, J, normal, &detJb); 123 q_data_sur[0][i] = w[i] * detJb; 124 q_data_sur[1][i] = normal[0]; 125 q_data_sur[2][i] = normal[1]; 126 q_data_sur[3][i] = normal[2]; 127 128 InvertBoundaryMappingJacobian_3D(Q, i, J, dXdx); 129 q_data_sur[4][i] = dXdx[0][0]; 130 q_data_sur[5][i] = dXdx[0][1]; 131 q_data_sur[6][i] = dXdx[0][2]; 132 q_data_sur[7][i] = dXdx[1][0]; 133 q_data_sur[8][i] = dXdx[1][1]; 134 q_data_sur[9][i] = dXdx[1][2]; 135 } 136 return 0; 137 } 138 139 // ***************************************************************************** 140 141 #endif // setup_geo_h 142