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 // Hdiv basis for quadrilateral linear BDM element in 2D 11 // Local numbering is as follow (each edge has 2 vector DoF) 12 // b4 b5 13 // 2---------3 14 // b7| |b3 15 // | | 16 // b6| |b2 17 // 0---------1 18 // b0 b1 19 // Bx[0-->7] = b0_x-->b7_x, By[0-->7] = b0_y-->b7_y 20 // To see how the nodal basis is constructed visit: 21 // https://github.com/rezgarshakeri/H-div-Tests 22 int BuildNodalHdivQuadrilateral(CeedScalar *x, CeedScalar *Bx, CeedScalar *By) { 23 CeedScalar x_hat = x[0]; 24 CeedScalar y_hat = x[1]; 25 Bx[0] = -0.125 + 0.125 * x_hat * x_hat; 26 By[0] = -0.25 + 0.25 * x_hat + 0.25 * y_hat + -0.25 * x_hat * y_hat; 27 Bx[1] = 0.125 + -0.125 * x_hat * x_hat; 28 By[1] = -0.25 + -0.25 * x_hat + 0.25 * y_hat + 0.25 * x_hat * y_hat; 29 Bx[2] = 0.25 + 0.25 * x_hat + -0.25 * y_hat + -0.25 * x_hat * y_hat; 30 By[2] = -0.125 + 0.125 * y_hat * y_hat; 31 Bx[3] = 0.25 + 0.25 * x_hat + 0.25 * y_hat + 0.25 * x_hat * y_hat; 32 By[3] = 0.125 + -0.125 * y_hat * y_hat; 33 Bx[4] = -0.125 + 0.125 * x_hat * x_hat; 34 By[4] = 0.25 + -0.25 * x_hat + 0.25 * y_hat + -0.25 * x_hat * y_hat; 35 Bx[5] = 0.125 + -0.125 * x_hat * x_hat; 36 By[5] = 0.25 + 0.25 * x_hat + 0.25 * y_hat + 0.25 * x_hat * y_hat; 37 Bx[6] = -0.25 + 0.25 * x_hat + 0.25 * y_hat + -0.25 * x_hat * y_hat; 38 By[6] = -0.125 + 0.125 * y_hat * y_hat; 39 Bx[7] = -0.25 + 0.25 * x_hat + -0.25 * y_hat + 0.25 * x_hat * y_hat; 40 By[7] = 0.125 + -0.125 * y_hat * y_hat; 41 return 0; 42 } 43 static void BuildHdivQuadrilateral(CeedInt q, CeedScalar *q_ref, CeedScalar *q_weights, CeedScalar *interp, CeedScalar *div, CeedQuadMode quad_mode) { 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 BuildNodalHdivQuadrilateral(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 }