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 #include <ceed/types.h>
9
10 // Compute det(A)
MatDet2x2(const CeedScalar A[2][2])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
AlphaMatTransposeMatMult2x2(const CeedScalar alpha,const CeedScalar A[2][2],const CeedScalar B[2][2],CeedScalar C[2][2])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
mass(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)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