1 /// @file
2 /// Test Collocated Grad calculated matches basis with Lobatto points
3 /// \test Test Collocated Grad calculated matches basis with Lobatto points
4 #include <ceed.h>
5 #include <ceed/backend.h>
6 #include <math.h>
7 #include <stdio.h>
8
main(int argc,char ** argv)9 int main(int argc, char **argv) {
10 Ceed ceed;
11 CeedInt p = 4;
12 CeedScalar collocated_gradient_1d[(p + 2) * (p + 2)], x_2[p + 2];
13 const CeedScalar *gradient_1d, *q_ref;
14 CeedScalar sum = 0.0;
15 CeedBasis basis;
16
17 CeedInit(argv[1], &ceed);
18
19 // Already collocated, GetCollocatedGrad will return grad_1d
20 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, p, CEED_GAUSS_LOBATTO, &basis);
21 CeedBasisGetCollocatedGrad(basis, collocated_gradient_1d);
22 CeedBasisGetGrad(basis, &gradient_1d);
23
24 for (CeedInt i = 0; i < p; i++) {
25 for (CeedInt j = 0; j < p; j++) {
26 if (fabs(collocated_gradient_1d[j + p * i] - gradient_1d[j + p * i]) > 100 * CEED_EPSILON) {
27 // LCOV_EXCL_START
28 printf("Error in collocated gradient %f != %f\n", collocated_gradient_1d[j + p * i], gradient_1d[j + p * i]);
29 // LCOV_EXCL_STOP
30 }
31 }
32 }
33 CeedBasisDestroy(&basis);
34
35 // Q = P, not already collocated
36 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, p, CEED_GAUSS, &basis);
37 CeedBasisGetCollocatedGrad(basis, collocated_gradient_1d);
38
39 CeedBasisGetQRef(basis, &q_ref);
40 for (CeedInt i = 0; i < p; i++) x_2[i] = q_ref[i] * q_ref[i];
41
42 // Verify collo_grad * x^2 = 2x for quadrature points
43 for (CeedInt i = 0; i < p; i++) {
44 sum = 0.0;
45 for (CeedInt j = 0; j < p; j++) sum += collocated_gradient_1d[j + p * i] * x_2[j];
46 if (fabs(sum - 2 * q_ref[i]) > 100 * CEED_EPSILON) printf("Error in collocated gradient %f != %f\n", sum, 2 * q_ref[i]);
47 }
48 CeedBasisDestroy(&basis);
49
50 // Q = P + 2, not already collocated
51 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, p + 2, CEED_GAUSS, &basis);
52 CeedBasisGetCollocatedGrad(basis, collocated_gradient_1d);
53
54 CeedBasisGetQRef(basis, &q_ref);
55 for (CeedInt i = 0; i < p + 2; i++) x_2[i] = q_ref[i] * q_ref[i];
56
57 // Verify collocated_gradient * x^2 = 2x for quadrature points
58 for (CeedInt i = 0; i < p + 2; i++) {
59 sum = 0.0;
60 for (CeedInt j = 0; j < p + 2; j++) sum += collocated_gradient_1d[j + (p + 2) * i] * x_2[j];
61 if (fabs(sum - 2 * q_ref[i]) > 100 * CEED_EPSILON) printf("Error in collocated gradient %f != %f\n", sum, 2 * q_ref[i]);
62 }
63 CeedBasisDestroy(&basis);
64
65 CeedDestroy(&ceed);
66 return 0;
67 }
68