1*3d8e8822SJeremy L Thompson // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2*3d8e8822SJeremy L Thompson // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 350c301a5SRezgar Shakeri // 4*3d8e8822SJeremy L Thompson // SPDX-License-Identifier: BSD-2-Clause 550c301a5SRezgar Shakeri // 6*3d8e8822SJeremy L Thompson // This file is part of CEED: http://github.com/ceed 750c301a5SRezgar Shakeri 850c301a5SRezgar Shakeri // Hdiv basis for quadrilateral linear BDMelement in 2D 950c301a5SRezgar Shakeri // Local numbering is as follow (each edge has 2 vector dof) 1050c301a5SRezgar Shakeri // b4 b5 1150c301a5SRezgar Shakeri // 2---------3 1250c301a5SRezgar Shakeri // b7| |b3 1350c301a5SRezgar Shakeri // | | 1450c301a5SRezgar Shakeri // b6| |b2 1550c301a5SRezgar Shakeri // 0---------1 1650c301a5SRezgar Shakeri // b0 b1 1750c301a5SRezgar Shakeri // Bx[0-->7] = b0_x-->b7_x, By[0-->7] = b0_y-->b7_y 1850c301a5SRezgar Shakeri // To see how the nodal basis is constructed visit: 1950c301a5SRezgar Shakeri // https://github.com/rezgarshakeri/H-div-Tests 2050c301a5SRezgar Shakeri int NodalHdivBasisQuad(CeedScalar *X, CeedScalar *Bx, CeedScalar *By) { 2150c301a5SRezgar Shakeri CeedScalar x_hat = X[0]; 2250c301a5SRezgar Shakeri CeedScalar y_hat = X[1]; 2350c301a5SRezgar Shakeri Bx[0] = -0.125 + 0.125*x_hat*x_hat; 2450c301a5SRezgar Shakeri By[0] = -0.25 + 0.25*x_hat + 0.25*y_hat + -0.25*x_hat*y_hat; 2550c301a5SRezgar Shakeri Bx[1] = 0.125 + -0.125*x_hat*x_hat; 2650c301a5SRezgar Shakeri By[1] = -0.25 + -0.25*x_hat + 0.25*y_hat + 0.25*x_hat*y_hat; 2750c301a5SRezgar Shakeri Bx[2] = 0.25 + 0.25*x_hat + -0.25*y_hat + -0.25*x_hat*y_hat; 2850c301a5SRezgar Shakeri By[2] = -0.125 + 0.125*y_hat*y_hat; 2950c301a5SRezgar Shakeri Bx[3] = 0.25 + 0.25*x_hat + 0.25*y_hat + 0.25*x_hat*y_hat; 3050c301a5SRezgar Shakeri By[3] = 0.125 + -0.125*y_hat*y_hat; 3150c301a5SRezgar Shakeri Bx[4] = -0.125 + 0.125*x_hat*x_hat; 3250c301a5SRezgar Shakeri By[4] = 0.25 + -0.25*x_hat + 0.25*y_hat + -0.25*x_hat*y_hat; 3350c301a5SRezgar Shakeri Bx[5] = 0.125 + -0.125*x_hat*x_hat; 3450c301a5SRezgar Shakeri By[5] = 0.25 + 0.25*x_hat + 0.25*y_hat + 0.25*x_hat*y_hat; 3550c301a5SRezgar Shakeri Bx[6] = -0.25 + 0.25*x_hat + 0.25*y_hat + -0.25*x_hat*y_hat; 3650c301a5SRezgar Shakeri By[6] = -0.125 + 0.125*y_hat*y_hat; 3750c301a5SRezgar Shakeri Bx[7] = -0.25 + 0.25*x_hat + -0.25*y_hat + 0.25*x_hat*y_hat; 3850c301a5SRezgar Shakeri By[7] = 0.125 + -0.125*y_hat*y_hat; 3950c301a5SRezgar Shakeri return 0; 4050c301a5SRezgar Shakeri } 4150c301a5SRezgar Shakeri static void HdivBasisQuad(CeedInt Q, CeedScalar *q_ref, CeedScalar *q_weights, 4250c301a5SRezgar Shakeri CeedScalar *interp, CeedScalar *div, CeedQuadMode quad_mode) { 4350c301a5SRezgar Shakeri 4450c301a5SRezgar Shakeri // Get 1D quadrature on [-1,1] 4550c301a5SRezgar Shakeri CeedScalar q_ref_1d[Q], q_weight_1d[Q]; 4650c301a5SRezgar Shakeri switch (quad_mode) { 4750c301a5SRezgar Shakeri case CEED_GAUSS: 4850c301a5SRezgar Shakeri CeedGaussQuadrature(Q, q_ref_1d, q_weight_1d); 4950c301a5SRezgar Shakeri break; 5059058f14SJeremy L Thompson // LCOV_EXCL_START 5150c301a5SRezgar Shakeri case CEED_GAUSS_LOBATTO: 5250c301a5SRezgar Shakeri CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d); 5350c301a5SRezgar Shakeri break; 5450c301a5SRezgar Shakeri } 5559058f14SJeremy L Thompson // LCOV_EXCL_STOP 5650c301a5SRezgar Shakeri 5750c301a5SRezgar Shakeri // Divergence operator; Divergence of nodal basis for ref element 5850c301a5SRezgar Shakeri CeedScalar D[8] = {0.25,0.25,0.25,0.25,0.25,0.25,0.25,0.25}; 5950c301a5SRezgar Shakeri // Loop over quadrature points 6050c301a5SRezgar Shakeri CeedScalar Bx[8], By[8]; 6150c301a5SRezgar Shakeri CeedScalar X[2]; 6250c301a5SRezgar Shakeri 6350c301a5SRezgar Shakeri for (CeedInt i=0; i<Q; i++) { 6450c301a5SRezgar Shakeri for (CeedInt j=0; j<Q; j++) { 6550c301a5SRezgar Shakeri CeedInt k1 = Q*i+j; 6650c301a5SRezgar Shakeri q_ref[k1] = q_ref_1d[j]; 6750c301a5SRezgar Shakeri q_ref[k1 + Q*Q] = q_ref_1d[i]; 6850c301a5SRezgar Shakeri q_weights[k1] = q_weight_1d[j]*q_weight_1d[i]; 6950c301a5SRezgar Shakeri X[0] = q_ref_1d[j]; 7050c301a5SRezgar Shakeri X[1] = q_ref_1d[i]; 7150c301a5SRezgar Shakeri NodalHdivBasisQuad(X, Bx, By); 7250c301a5SRezgar Shakeri for (CeedInt k=0; k<8; k++) { 7350c301a5SRezgar Shakeri interp[k1*8+k] = Bx[k]; 7450c301a5SRezgar Shakeri interp[k1*8+k+8*Q*Q] = By[k]; 7550c301a5SRezgar Shakeri div[k1*8+k] = D[k]; 7650c301a5SRezgar Shakeri } 7750c301a5SRezgar Shakeri } 7850c301a5SRezgar Shakeri } 7950c301a5SRezgar Shakeri }