1*3d8e8822SJeremy L Thompson // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2*3d8e8822SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 315ed3f7dSjeremylt // 4*3d8e8822SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause 515ed3f7dSjeremylt // 6*3d8e8822SJeremy L Thompson // This file is part of CEED: http://github.com/ceed 715ed3f7dSjeremylt 815ed3f7dSjeremylt CEED_QFUNCTION(setup_diff)(void *ctx, const CeedInt Q, 915ed3f7dSjeremylt const CeedScalar *const *in, 1015ed3f7dSjeremylt CeedScalar *const *out) { 1115ed3f7dSjeremylt // in[0] is Jacobians with shape [2, nc=2, Q] 1215ed3f7dSjeremylt // in[1] is quadrature weights, size (Q) 1315ed3f7dSjeremylt const CeedScalar *J = in[0], *w = in[1]; 1415ed3f7dSjeremylt 1515ed3f7dSjeremylt // out[0] is qdata, size (Q) 1615ed3f7dSjeremylt CeedScalar *q_data = out[0]; 1715ed3f7dSjeremylt 1815ed3f7dSjeremylt // Quadrature point loop 1915ed3f7dSjeremylt CeedPragmaSIMD 2015ed3f7dSjeremylt for (CeedInt i=0; i<Q; i++) { 2115ed3f7dSjeremylt // Qdata stored in Voigt convention 2215ed3f7dSjeremylt // J: 0 2 q_data: 0 2 adj(J): J22 -J12 2315ed3f7dSjeremylt // 1 3 2 1 -J21 J11 2415ed3f7dSjeremylt const CeedScalar J11 = J[i+Q*0]; 2515ed3f7dSjeremylt const CeedScalar J21 = J[i+Q*1]; 2615ed3f7dSjeremylt const CeedScalar J12 = J[i+Q*2]; 2715ed3f7dSjeremylt const CeedScalar J22 = J[i+Q*3]; 2815ed3f7dSjeremylt const CeedScalar qw = w[i] / (J11*J22 - J21*J12); 2915ed3f7dSjeremylt q_data[i+Q*0] = qw * (J12*J12 + J22*J22); 3015ed3f7dSjeremylt q_data[i+Q*1] = qw * (J11*J11 + J21*J21); 3115ed3f7dSjeremylt q_data[i+Q*2] = - qw * (J11*J12 + J21*J22); 3215ed3f7dSjeremylt } // End of Quadrature Point Loop 3315ed3f7dSjeremylt return 0; 3415ed3f7dSjeremylt } 3515ed3f7dSjeremylt 3615ed3f7dSjeremylt CEED_QFUNCTION(apply)(void *ctx, const CeedInt Q, const CeedScalar *const *in, 3715ed3f7dSjeremylt CeedScalar *const *out) { 3815ed3f7dSjeremylt // in[0] is gradient u, shape [2, nc=1, Q] 3915ed3f7dSjeremylt // in[1] is quadrature data, size (3*Q) 4015ed3f7dSjeremylt const CeedScalar *ug = in[0], *q_data = in[1]; 4115ed3f7dSjeremylt 4215ed3f7dSjeremylt // out[0] is output to multiply against gradient v, shape [2, nc=1, Q] 4315ed3f7dSjeremylt CeedScalar *vg = out[0]; 4415ed3f7dSjeremylt 4515ed3f7dSjeremylt // Quadrature point loop 4615ed3f7dSjeremylt CeedPragmaSIMD 4715ed3f7dSjeremylt for (CeedInt i=0; i<Q; i++) { 4815ed3f7dSjeremylt // Read spatial derivatives of u 4915ed3f7dSjeremylt const CeedScalar du[2] = {ug[i+Q*0], 5015ed3f7dSjeremylt ug[i+Q*1] 5115ed3f7dSjeremylt }; 5215ed3f7dSjeremylt 5315ed3f7dSjeremylt // Read qdata (dXdxdXdxT symmetric matrix) 5415ed3f7dSjeremylt // Stored in Voigt convention 5515ed3f7dSjeremylt // 0 2 5615ed3f7dSjeremylt // 2 1 5715ed3f7dSjeremylt // *INDENT-OFF* 5815ed3f7dSjeremylt const CeedScalar dXdxdXdxT[2][2] = {{q_data[i+0*Q], 5915ed3f7dSjeremylt q_data[i+2*Q]}, 6015ed3f7dSjeremylt {q_data[i+2*Q], 6115ed3f7dSjeremylt q_data[i+1*Q]} 6215ed3f7dSjeremylt }; 6315ed3f7dSjeremylt // *INDENT-ON* 6415ed3f7dSjeremylt 6515ed3f7dSjeremylt // Apply Poisson operator 6615ed3f7dSjeremylt // j = direction of vg 6715ed3f7dSjeremylt for (int j=0; j<2; j++) 6815ed3f7dSjeremylt vg[i+j*Q] = (du[0] * dXdxdXdxT[0][j] + 6915ed3f7dSjeremylt du[1] * dXdxdXdxT[1][j]); 7015ed3f7dSjeremylt } // End of Quadrature Point Loop 7115ed3f7dSjeremylt return 0; 7215ed3f7dSjeremylt } 73