1 // Copyright (c) 2017-2026, 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 // H(div) basis for quadrilateral linear BDM element in 2D
11 // Local numbering is as follows (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
BuildNodalHdivQuadrilateral(CeedScalar * x,CeedScalar * Bx,CeedScalar * By)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
BuildHdivQuadrilateral(CeedInt q,CeedScalar * q_ref,CeedScalar * q_weights,CeedScalar * interp,CeedScalar * div,CeedQuadMode quad_mode)44 static void BuildHdivQuadrilateral(CeedInt q, CeedScalar *q_ref, CeedScalar *q_weights, CeedScalar *interp, CeedScalar *div, CeedQuadMode quad_mode) {
45 // Get 1D quadrature on [-1,1]
46 CeedScalar q_ref_1d[q], q_weight_1d[q];
47 switch (quad_mode) {
48 case CEED_GAUSS:
49 CeedGaussQuadrature(q, q_ref_1d, q_weight_1d);
50 break;
51 // LCOV_EXCL_START
52 case CEED_GAUSS_LOBATTO:
53 CeedLobattoQuadrature(q, q_ref_1d, q_weight_1d);
54 break;
55 }
56 // LCOV_EXCL_STOP
57
58 // Divergence operator; Divergence of nodal basis for ref element
59 CeedScalar D[8] = {0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25};
60 CeedScalar Bx[8], By[8];
61 CeedScalar X[2];
62
63 // Loop over quadrature points
64 for (CeedInt i = 0; i < q; i++) {
65 for (CeedInt j = 0; j < q; j++) {
66 CeedInt k1 = q * i + j;
67 q_ref[k1] = q_ref_1d[j];
68 q_ref[k1 + q * q] = q_ref_1d[i];
69 q_weights[k1] = q_weight_1d[j] * q_weight_1d[i];
70 X[0] = q_ref_1d[j];
71 X[1] = q_ref_1d[i];
72 BuildNodalHdivQuadrilateral(X, Bx, By);
73 for (CeedInt k = 0; k < 8; k++) {
74 interp[k1 * 8 + k] = Bx[k];
75 interp[k1 * 8 + k + 8 * q * q] = By[k];
76 div[k1 * 8 + k] = D[k];
77 }
78 }
79 }
80 }
81