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 #include <ceed.h> 9 10 CEED_QFUNCTION(setup)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 11 // *INDENT-OFF* 12 const CeedScalar *w = in[0], (*J)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[1]; 13 CeedScalar(*q_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 14 // *INDENT-ON* 15 16 // Quadrature point loop 17 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 18 // Qdata stored in Voigt convention 19 // J: 0 2 q_data: 0 2 adj(J): J11 -J01 20 // 1 3 2 1 -J10 J00 21 const CeedScalar J00 = J[0][0][i]; 22 const CeedScalar J10 = J[0][1][i]; 23 const CeedScalar J01 = J[1][0][i]; 24 const CeedScalar J11 = J[1][1][i]; 25 const CeedScalar qw = w[i] / (J00 * J11 - J10 * J01); 26 q_data[0][i] = qw * (J01 * J01 + J11 * J11); 27 q_data[1][i] = qw * (J00 * J00 + J10 * J10); 28 q_data[2][i] = -qw * (J00 * J01 + J10 * J11); 29 } // End of Quadrature Point Loop 30 31 return 0; 32 } 33 34 CEED_QFUNCTION(diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 35 // *INDENT-OFF* 36 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], (*ug)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[1]; 37 CeedScalar(*vg)[2][CEED_Q_VLA] = (CeedScalar(*)[2][CEED_Q_VLA])out[0]; 38 // *INDENT-ON* 39 40 const CeedInt dim = 2; 41 const CeedScalar num_comp = 2; 42 const CeedScalar scale[2][2] = { 43 {1.0, 2.0}, 44 {3.0, 4.0}, 45 }; 46 47 // Quadrature point loop 48 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 49 // Read qdata (dXdxdXdxT symmetric matrix) 50 // Stored in Voigt convention 51 // 0 2 52 // 2 1 53 // *INDENT-OFF* 54 const CeedScalar dXdxdXdxT[2][2] = { 55 {q_data[0][i], q_data[2][i]}, 56 {q_data[2][i], q_data[1][i]} 57 }; 58 // *INDENT-ON* 59 60 // Apply Poisson operator 61 // j = direction of vg 62 for (CeedInt j = 0; j < dim; j++) { 63 for (CeedInt k = 0; k < num_comp; k++) { 64 vg[j][k][i] = (ug[0][k][i] * dXdxdXdxT[0][j] * scale[0][j] + ug[1][k][i] * dXdxdXdxT[1][j] * scale[1][j]); 65 } 66 } 67 } // End of Quadrature Point Loop 68 69 return 0; 70 } 71