1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC. 2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. 3 // All Rights reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 CEED_QFUNCTION(setup_diff)(void *ctx, const CeedInt Q, 18 const CeedScalar *const *in, 19 CeedScalar *const *out) { 20 // in[0] is Jacobians with shape [2, nc=2, Q] 21 // in[1] is quadrature weights, size (Q) 22 const CeedScalar *J = in[0], *w = in[1]; 23 24 // out[0] is qdata, size (Q) 25 CeedScalar *q_data = out[0]; 26 27 // Quadrature point loop 28 CeedPragmaSIMD 29 for (CeedInt i=0; i<Q; i++) { 30 // Qdata stored in Voigt convention 31 // J: 0 2 q_data: 0 2 adj(J): J22 -J12 32 // 1 3 2 1 -J21 J11 33 const CeedScalar J11 = J[i+Q*0]; 34 const CeedScalar J21 = J[i+Q*1]; 35 const CeedScalar J12 = J[i+Q*2]; 36 const CeedScalar J22 = J[i+Q*3]; 37 const CeedScalar qw = w[i] / (J11*J22 - J21*J12); 38 q_data[i+Q*0] = qw * (J12*J12 + J22*J22); 39 q_data[i+Q*1] = qw * (J11*J11 + J21*J21); 40 q_data[i+Q*2] = - qw * (J11*J12 + J21*J22); 41 } // End of Quadrature Point Loop 42 return 0; 43 } 44 45 CEED_QFUNCTION(apply)(void *ctx, const CeedInt Q, const CeedScalar *const *in, 46 CeedScalar *const *out) { 47 // in[0] is gradient u, shape [2, nc=1, Q] 48 // in[1] is quadrature data, size (3*Q) 49 const CeedScalar *ug = in[0], *q_data = in[1]; 50 51 // out[0] is output to multiply against gradient v, shape [2, nc=1, Q] 52 CeedScalar *vg = out[0]; 53 54 // Quadrature point loop 55 CeedPragmaSIMD 56 for (CeedInt i=0; i<Q; i++) { 57 // Read spatial derivatives of u 58 const CeedScalar du[2] = {ug[i+Q*0], 59 ug[i+Q*1] 60 }; 61 62 // Read qdata (dXdxdXdxT symmetric matrix) 63 // Stored in Voigt convention 64 // 0 2 65 // 2 1 66 // *INDENT-OFF* 67 const CeedScalar dXdxdXdxT[2][2] = {{q_data[i+0*Q], 68 q_data[i+2*Q]}, 69 {q_data[i+2*Q], 70 q_data[i+1*Q]} 71 }; 72 // *INDENT-ON* 73 74 // Apply Poisson operator 75 // j = direction of vg 76 for (int j=0; j<2; j++) 77 vg[i+j*Q] = (du[0] * dXdxdXdxT[0][j] + 78 du[1] * dXdxdXdxT[1][j]); 79 } // End of Quadrature Point Loop 80 return 0; 81 } 82