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 /** 9 @brief Ceed QFunction for building the geometric data for the 3D Poisson 10 operator 11 **/ 12 13 #ifndef poisson3dbuild_h 14 #define poisson3dbuild_h 15 16 #include <ceed.h> 17 18 CEED_QFUNCTION(Poisson3DBuild)(void *ctx, const CeedInt Q, 19 const CeedScalar *const *in, 20 CeedScalar *const *out) { 21 // At every quadrature point, compute w/det(J).adj(J).adj(J)^T and store 22 // the symmetric part of the result. 23 // *INDENT-OFF* 24 // in[0] is Jacobians with shape [3, nc=3, Q] 25 // in[1] is quadrature weights, size (Q) 26 const CeedScalar (*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], 27 *w = in[1]; 28 // out[0] is qdata, size (6*Q) 29 CeedScalar (*q_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 30 // *INDENT-ON* 31 32 const CeedInt dim = 3; 33 34 // Quadrature point loop 35 CeedPragmaSIMD 36 for (CeedInt i=0; i<Q; i++) { 37 // Compute the adjoint 38 CeedScalar A[3][3]; 39 for (CeedInt j=0; j<dim; j++) 40 for (CeedInt k=0; k<dim; k++) 41 // Equivalent code with no mod operations: 42 // A[k][j] = J[k+1][j+1]*J[k+2][j+2] - J[k+2][j+1]*J[k+1][j+2] 43 A[k][j] = J[(k+1)%dim][(j+1)%dim][i]*J[(k+2)%dim][(j+2)%dim][i] - 44 J[(k+2)%dim][(j+1)%dim][i]*J[(k+1)%dim][(j+2)%dim][i]; 45 46 // Compute quadrature weight / det(J) 47 const CeedScalar qw = w[i] / (J[0][0][i]*A[0][0] + J[0][1][i]*A[0][1] + 48 J[0][2][i]*A[0][2]); 49 50 // Compute geometric factors 51 // Stored in Voigt convention 52 // 0 5 4 53 // 5 1 3 54 // 4 3 2 55 q_data[0][i] = qw * (A[0][0]*A[0][0] + A[0][1]*A[0][1] + A[0][2]*A[0][2]); 56 q_data[1][i] = qw * (A[1][0]*A[1][0] + A[1][1]*A[1][1] + A[1][2]*A[1][2]); 57 q_data[2][i] = qw * (A[2][0]*A[2][0] + A[2][1]*A[2][1] + A[2][2]*A[2][2]); 58 q_data[3][i] = qw * (A[1][0]*A[2][0] + A[1][1]*A[2][1] + A[1][2]*A[2][2]); 59 q_data[4][i] = qw * (A[0][0]*A[2][0] + A[0][1]*A[2][1] + A[0][2]*A[2][2]); 60 q_data[5][i] = qw * (A[0][0]*A[1][0] + A[0][1]*A[1][1] + A[0][2]*A[1][2]); 61 } // End of Quadrature Point Loop 62 63 return CEED_ERROR_SUCCESS; 64 } 65 66 #endif // poisson3dbuild_h 67