13d8e8822SJeremy L Thompson // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 23d8e8822SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3b6faaefaSjeremylt // 43d8e8822SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause 5b6faaefaSjeremylt // 63d8e8822SJeremy L Thompson // This file is part of CEED: http://github.com/ceed 7b6faaefaSjeremylt 8*c9c2c079SJeremy L Thompson #include <ceed.h> 9*c9c2c079SJeremy L Thompson 10b6faaefaSjeremylt CEED_QFUNCTION(setup)(void *ctx, const CeedInt Q, 11b6faaefaSjeremylt const CeedScalar *const *in, 12b6faaefaSjeremylt CeedScalar *const *out) { 13b6faaefaSjeremylt const CeedScalar *qw = in[0], *J = in[1]; 14b6faaefaSjeremylt CeedScalar *qd = out[0]; 15b6faaefaSjeremylt 16b6faaefaSjeremylt for (CeedInt i=0; i<Q; i++) { 179a53e505Sjeremylt // Qdata stored in Voigt convention 18b6faaefaSjeremylt // J: 0 2 qd: 0 2 adj(J): J22 -J12 19b6faaefaSjeremylt // 1 3 2 1 -J21 J11 20b6faaefaSjeremylt const CeedScalar J11 = J[i+Q*0]; 21b6faaefaSjeremylt const CeedScalar J21 = J[i+Q*1]; 22b6faaefaSjeremylt const CeedScalar J12 = J[i+Q*2]; 23b6faaefaSjeremylt const CeedScalar J22 = J[i+Q*3]; 24b6faaefaSjeremylt const CeedScalar w = qw[i] / (J11*J22 - J21*J12); 25b6faaefaSjeremylt qd[i+Q*0] = w * (J12*J12 + J22*J22); 26b6faaefaSjeremylt qd[i+Q*1] = w * (J11*J11 + J21*J21); 27b6faaefaSjeremylt qd[i+Q*2] = - w * (J11*J12 + J21*J22); 28b6faaefaSjeremylt } // End of Quadrature Point Loop 29b6faaefaSjeremylt return 0; 30b6faaefaSjeremylt } 31b6faaefaSjeremylt 32b6faaefaSjeremylt CEED_QFUNCTION(diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, 33b6faaefaSjeremylt CeedScalar *const *out) { 34b6faaefaSjeremylt const CeedScalar *qd = in[0], *ug = in[1]; 35b6faaefaSjeremylt CeedScalar *vg = out[0]; 36b6faaefaSjeremylt 37b6faaefaSjeremylt for (CeedInt i=0; i<Q; i++) { 38b6faaefaSjeremylt // Read spatial derivatives of u 39b6faaefaSjeremylt const CeedScalar du[2] = {ug[i+Q*0], 40b6faaefaSjeremylt ug[i+Q*1] 41b6faaefaSjeremylt }; 42b6faaefaSjeremylt 43b6faaefaSjeremylt // Read qdata (dXdxdXdxT symmetric matrix) 44b6faaefaSjeremylt // Stored in Voigt convention 45b6faaefaSjeremylt // 0 2 46b6faaefaSjeremylt // 2 1 47cb0b5415Sjeremylt // *INDENT-OFF* 48b6faaefaSjeremylt const CeedScalar dXdxdXdxT[2][2] = {{qd[i+0*Q], 49b6faaefaSjeremylt qd[i+2*Q]}, 50b6faaefaSjeremylt {qd[i+2*Q], 51b6faaefaSjeremylt qd[i+1*Q]} 52b6faaefaSjeremylt }; 53cb0b5415Sjeremylt // *INDENT-ON* 54b6faaefaSjeremylt // j = direction of vg 55b6faaefaSjeremylt for (int j=0; j<2; j++) 56b6faaefaSjeremylt vg[i+j*Q] = (du[0] * dXdxdXdxT[0][j] + 57b6faaefaSjeremylt du[1] * dXdxdXdxT[1][j]); 58b6faaefaSjeremylt } 59b6faaefaSjeremylt return 0; 60b6faaefaSjeremylt } 61