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