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