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 #include <ceed.h> 9 10 // Compute det(A) 11 CEED_QFUNCTION_HELPER CeedScalar MatDet2x2(const CeedScalar A[2][2]) { return A[0][0] * A[1][1] - A[1][0] * A[0][1]; } 12 13 // Compute alpha * A^T * B = C 14 CEED_QFUNCTION_HELPER int AlphaMatTransposeMatMult2x2(const CeedScalar alpha, const CeedScalar A[2][2], const CeedScalar B[2][2], 15 CeedScalar C[2][2]) { 16 for (CeedInt j = 0; j < 2; j++) { 17 for (CeedInt k = 0; k < 2; k++) { 18 C[j][k] = 0; 19 for (CeedInt m = 0; m < 2; m++) { 20 C[j][k] += alpha * A[m][j] * B[m][k]; 21 } 22 } 23 } 24 25 return 0; 26 } 27 28 CEED_QFUNCTION(mass)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 29 // Inputs 30 const CeedScalar(*w) = in[0], (*dxdX)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[1], 31 (*u)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 32 33 // Outputs 34 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 35 36 // Quadrature Point Loop 37 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 38 // Setup, J = dx/dX 39 const CeedScalar J[2][2] = { 40 {dxdX[0][0][i], dxdX[1][0][i]}, 41 {dxdX[0][1][i], dxdX[1][1][i]} 42 }; 43 const CeedScalar det_J = MatDet2x2(J); 44 45 // Piola map: J^T*J*u*w/detJ 46 // 1) Compute J^T * J 47 CeedScalar JT_J[2][2]; 48 AlphaMatTransposeMatMult2x2(1., J, J, JT_J); 49 50 // 2) Compute J^T*J*u * w /detJ 51 for (CeedInt k = 0; k < 2; k++) { 52 v[k][i] = 0; 53 for (CeedInt m = 0; m < 2; m++) v[k][i] += JT_J[k][m] * u[m][i] * w[i] / det_J; 54 } 55 } // End of Quadrature Point Loop 56 57 return 0; 58 } 59