1 // Copyright (c) 2017-2025, 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 CEED_QFUNCTION(setup_diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 11 // in[0] is Jacobians with shape [2, nc=2, Q] 12 // in[1] is quadrature weights, size (Q) 13 const CeedScalar *J = in[0], *w = in[1]; 14 15 // out[0] is qdata, size (Q) 16 CeedScalar *q_data = out[0]; 17 18 // Quadrature point loop 19 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 20 // Qdata stored in Voigt convention 21 // J: 0 2 q_data: 0 2 adj(J): J22 -J12 22 // 1 3 2 1 -J21 J11 23 const CeedScalar J11 = J[i + Q * 0]; 24 const CeedScalar J21 = J[i + Q * 1]; 25 const CeedScalar J12 = J[i + Q * 2]; 26 const CeedScalar J22 = J[i + Q * 3]; 27 const CeedScalar qw = w[i] / (J11 * J22 - J21 * J12); 28 q_data[i + Q * 0] = qw * (J12 * J12 + J22 * J22); 29 q_data[i + Q * 1] = qw * (J11 * J11 + J21 * J21); 30 q_data[i + Q * 2] = -qw * (J11 * J12 + J21 * J22); 31 } // End of Quadrature Point Loop 32 return 0; 33 } 34 35 CEED_QFUNCTION(apply)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 36 // in[0] is gradient u, shape [2, nc=1, Q] 37 // in[1] is quadrature data, size (3*Q) 38 const CeedScalar *ug = in[0], *q_data = in[1]; 39 40 // out[0] is output to multiply against gradient v, shape [2, nc=1, Q] 41 CeedScalar *vg = out[0]; 42 43 // Quadrature point loop 44 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 45 // Read spatial derivatives of u 46 const CeedScalar du[2] = {ug[i + Q * 0], ug[i + Q * 1]}; 47 48 // Read qdata (dXdxdXdxT symmetric matrix) 49 // Stored in Voigt convention 50 // 0 2 51 // 2 1 52 const CeedScalar dXdxdXdxT[2][2] = { 53 {q_data[i + 0 * Q], q_data[i + 2 * Q]}, 54 {q_data[i + 2 * Q], q_data[i + 1 * Q]} 55 }; 56 57 // Apply Poisson operator 58 // j = direction of vg 59 for (int j = 0; j < 2; j++) vg[i + j * Q] = (du[0] * dXdxdXdxT[0][j] + du[1] * dXdxdXdxT[1][j]); 60 } // End of Quadrature Point Loop 61 return 0; 62 } 63