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(curl) basis for order 2 Nédélec (first kind) triangle element in 2D
11 // See: https://defelement.com/elements/examples/triangle-N1curl-2.html
BuildHcurl2DSimplex(CeedScalar * q_ref,CeedScalar * q_weight,CeedScalar * interp,CeedScalar * curl)12 static void BuildHcurl2DSimplex(CeedScalar *q_ref, CeedScalar *q_weight, CeedScalar *interp, CeedScalar *curl) {
13 CeedInt P = 8, Q = 4;
14
15 q_ref[0] = 1. / 3.;
16 q_ref[4] = 1. / 3.;
17 q_ref[1] = 0.2;
18 q_ref[5] = 0.2;
19 q_ref[2] = 0.2;
20 q_ref[6] = 0.6;
21 q_ref[3] = 0.6;
22 q_ref[7] = 0.2;
23 q_weight[0] = -25. / 96.;
24 q_weight[1] = 25. / 96.;
25 q_weight[2] = 27. / 96.;
26 q_weight[3] = 25. / 96.;
27
28 // Loop over quadrature points
29 for (int i = 0; i < Q; i++) {
30 CeedScalar x1 = q_ref[0 * Q + i], x2 = q_ref[1 * Q + i];
31 // Interp
32 interp[(i + 0) * P + 0] = 2. * x2 * (1. - 4. * x1);
33 interp[(i + Q) * P + 0] = 4. * x1 * (2. * x1 - 1.);
34 interp[(i + 0) * P + 1] = 4. * x2 * (1. - 2. * x2);
35 interp[(i + Q) * P + 1] = 2. * x1 * (4. * x2 - 1.);
36 interp[(i + 0) * P + 2] = 2. * x2 * (-4. * x1 - 4. * x2 + 3.);
37 interp[(i + Q) * P + 2] = 8. * x1 * x1 + 8. * x1 * x2 - 12. * x1 - 6. * x2 + 4.;
38 interp[(i + 0) * P + 3] = 4. * x2 * (2. * x2 - 1.);
39 interp[(i + Q) * P + 3] = -8. * x1 * x2 + 2. * x1 + 6. * x2 - 2.;
40 interp[(i + 0) * P + 4] = 8. * x1 * x2 - 6. * x1 + 8. * x2 * x2 - 12. * x2 + 4.;
41 interp[(i + Q) * P + 4] = 2. * x1 * (-4. * x1 - 4. * x2 + 3.);
42 interp[(i + 0) * P + 5] = -8. * x1 * x2 + 6. * x1 + 2. * x2 - 2.;
43 interp[(i + Q) * P + 5] = 4. * x1 * (2. * x1 - 1.);
44 interp[(i + 0) * P + 6] = 8. * x2 * (-x1 - 2. * x2 + 2.);
45 interp[(i + Q) * P + 6] = 8. * x1 * (x1 + 2. * x2 - 1.);
46 interp[(i + 0) * P + 7] = 8. * x2 * (2. * x1 + x2 - 1.);
47 interp[(i + Q) * P + 7] = 8. * x1 * (-2. * x1 - x2 + 2.);
48 // Curl
49 curl[i * P + 0] = 24. * x1 - 6.;
50 curl[i * P + 1] = 24. * x2 - 6.;
51 curl[i * P + 2] = 24. * x1 + 24. * x2 - 18.;
52 curl[i * P + 3] = -24. * x2 + 6.;
53 curl[i * P + 4] = -24. * x1 - 24. * x2 + 18.;
54 curl[i * P + 5] = 24. * x1 - 6.;
55 curl[i * P + 6] = 24. * x1 + 48. * x2 - 24.;
56 curl[i * P + 7] = -48. * x1 - 24. * x1 + 24.;
57 }
58 }
59