1*5aed82e4SJeremy L Thompson // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 2a0154adeSJed Brown // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3a0154adeSJed Brown // 4a0154adeSJed Brown // SPDX-License-Identifier: BSD-2-Clause 5a0154adeSJed Brown // 6a0154adeSJed Brown // This file is part of CEED: http://github.com/ceed 7a0154adeSJed Brown 8a0154adeSJed Brown /** 9ea61e9acSJeremy L Thompson @brief Ceed QFunction for building the geometric data for the 3D Poisson operator 10a0154adeSJed Brown **/ 11a0154adeSJed Brown 1294b7b29bSJeremy L Thompson #ifndef CEED_POISSON3DBUILD_H 1394b7b29bSJeremy L Thompson #define CEED_POISSON3DBUILD_H 14a0154adeSJed Brown 15c9c2c079SJeremy L Thompson #include <ceed.h> 16c9c2c079SJeremy L Thompson 172b730f8bSJeremy L Thompson CEED_QFUNCTION(Poisson3DBuild)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 18ea61e9acSJeremy L Thompson // At every quadrature point, compute w/det(J).adj(J).adj(J)^T and store the symmetric part of the result. 19a0154adeSJed Brown // in[0] is Jacobians with shape [3, nc=3, Q] 20a0154adeSJed Brown // in[1] is quadrature weights, size (Q) 212b730f8bSJeremy L Thompson const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[0], *w = in[1]; 22a0154adeSJed Brown // out[0] is qdata, size (6*Q) 23a0154adeSJed Brown CeedScalar(*q_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 24a0154adeSJed Brown 25a0154adeSJed Brown const CeedInt dim = 3; 26a0154adeSJed Brown 27a0154adeSJed Brown // Quadrature point loop 282b730f8bSJeremy L Thompson CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 29a0154adeSJed Brown // Compute the adjoint 30a0154adeSJed Brown CeedScalar A[3][3]; 31a0154adeSJed Brown for (CeedInt j = 0; j < dim; j++) 32a0154adeSJed Brown for (CeedInt k = 0; k < dim; k++) 33a0154adeSJed Brown // Equivalent code with no mod operations: 34a0154adeSJed Brown // A[k][j] = J[k+1][j+1]*J[k+2][j+2] - J[k+2][j+1]*J[k+1][j+2] 35a0154adeSJed Brown A[k][j] = J[(k + 1) % dim][(j + 1) % dim][i] * J[(k + 2) % dim][(j + 2) % dim][i] - 36a0154adeSJed Brown J[(k + 2) % dim][(j + 1) % dim][i] * J[(k + 1) % dim][(j + 2) % dim][i]; 37a0154adeSJed Brown 38a0154adeSJed Brown // Compute quadrature weight / det(J) 392b730f8bSJeremy L Thompson const CeedScalar qw = w[i] / (J[0][0][i] * A[0][0] + J[0][1][i] * A[0][1] + J[0][2][i] * A[0][2]); 40a0154adeSJed Brown 41a0154adeSJed Brown // Compute geometric factors 42a0154adeSJed Brown // Stored in Voigt convention 43a0154adeSJed Brown // 0 5 4 44a0154adeSJed Brown // 5 1 3 45a0154adeSJed Brown // 4 3 2 46a0154adeSJed Brown q_data[0][i] = qw * (A[0][0] * A[0][0] + A[0][1] * A[0][1] + A[0][2] * A[0][2]); 47a0154adeSJed Brown q_data[1][i] = qw * (A[1][0] * A[1][0] + A[1][1] * A[1][1] + A[1][2] * A[1][2]); 48a0154adeSJed Brown q_data[2][i] = qw * (A[2][0] * A[2][0] + A[2][1] * A[2][1] + A[2][2] * A[2][2]); 49a0154adeSJed Brown q_data[3][i] = qw * (A[1][0] * A[2][0] + A[1][1] * A[2][1] + A[1][2] * A[2][2]); 50a0154adeSJed Brown q_data[4][i] = qw * (A[0][0] * A[2][0] + A[0][1] * A[2][1] + A[0][2] * A[2][2]); 51a0154adeSJed Brown q_data[5][i] = qw * (A[0][0] * A[1][0] + A[0][1] * A[1][1] + A[0][2] * A[1][2]); 52a0154adeSJed Brown } // End of Quadrature Point Loop 53a0154adeSJed Brown 54a0154adeSJed Brown return CEED_ERROR_SUCCESS; 55a0154adeSJed Brown } 56a0154adeSJed Brown 5794b7b29bSJeremy L Thompson #endif // CEED_POISSON3DBUILD_H 58