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