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