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 // Hdiv basis for quadrilateral linear BDMelement in 2D 18 // Local numbering is as follow (each edge has 2 vector dof) 19 // b4 b5 20 // 2---------3 21 // b7| |b3 22 // | | 23 // b6| |b2 24 // 0---------1 25 // b0 b1 26 // Bx[0-->7] = b0_x-->b7_x, By[0-->7] = b0_y-->b7_y 27 // To see how the nodal basis is constructed visit: 28 // https://github.com/rezgarshakeri/H-div-Tests 29 int NodalHdivBasisQuad(CeedScalar *X, CeedScalar *Bx, CeedScalar *By) { 30 CeedScalar x_hat = X[0]; 31 CeedScalar y_hat = X[1]; 32 Bx[0] = -0.125 + 0.125*x_hat*x_hat; 33 By[0] = -0.25 + 0.25*x_hat + 0.25*y_hat + -0.25*x_hat*y_hat; 34 Bx[1] = 0.125 + -0.125*x_hat*x_hat; 35 By[1] = -0.25 + -0.25*x_hat + 0.25*y_hat + 0.25*x_hat*y_hat; 36 Bx[2] = 0.25 + 0.25*x_hat + -0.25*y_hat + -0.25*x_hat*y_hat; 37 By[2] = -0.125 + 0.125*y_hat*y_hat; 38 Bx[3] = 0.25 + 0.25*x_hat + 0.25*y_hat + 0.25*x_hat*y_hat; 39 By[3] = 0.125 + -0.125*y_hat*y_hat; 40 Bx[4] = -0.125 + 0.125*x_hat*x_hat; 41 By[4] = 0.25 + -0.25*x_hat + 0.25*y_hat + -0.25*x_hat*y_hat; 42 Bx[5] = 0.125 + -0.125*x_hat*x_hat; 43 By[5] = 0.25 + 0.25*x_hat + 0.25*y_hat + 0.25*x_hat*y_hat; 44 Bx[6] = -0.25 + 0.25*x_hat + 0.25*y_hat + -0.25*x_hat*y_hat; 45 By[6] = -0.125 + 0.125*y_hat*y_hat; 46 Bx[7] = -0.25 + 0.25*x_hat + -0.25*y_hat + 0.25*x_hat*y_hat; 47 By[7] = 0.125 + -0.125*y_hat*y_hat; 48 return 0; 49 } 50 static void HdivBasisQuad(CeedInt Q, CeedScalar *q_ref, CeedScalar *q_weights, 51 CeedScalar *interp, CeedScalar *div, CeedQuadMode quad_mode) { 52 53 // Get 1D quadrature on [-1,1] 54 CeedScalar q_ref_1d[Q], q_weight_1d[Q]; 55 switch (quad_mode) { 56 case CEED_GAUSS: 57 CeedGaussQuadrature(Q, q_ref_1d, q_weight_1d); 58 break; 59 // LCOV_EXCL_START 60 case CEED_GAUSS_LOBATTO: 61 CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d); 62 break; 63 } 64 // LCOV_EXCL_STOP 65 66 // Divergence operator; Divergence of nodal basis for ref element 67 CeedScalar D[8] = {0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25}; 68 // Loop over quadrature points 69 CeedScalar Bx[8], By[8]; 70 CeedScalar X[2]; 71 72 for (CeedInt i=0; i<Q; i++) { 73 for (CeedInt j=0; j<Q; j++) { 74 CeedInt k1 = Q*i+j; 75 q_ref[k1] = q_ref_1d[j]; 76 q_ref[k1 + Q*Q] = q_ref_1d[i]; 77 q_weights[k1] = q_weight_1d[j]*q_weight_1d[i]; 78 X[0] = q_ref_1d[j]; 79 X[1] = q_ref_1d[i]; 80 NodalHdivBasisQuad(X, Bx, By); 81 for (CeedInt k=0; k<8; k++) { 82 interp[k1*8+k] = Bx[k]; 83 interp[k1*8+k+8*Q*Q] = By[k]; 84 div[k1*8+k] = D[k]; 85 } 86 } 87 } 88 }