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