1 /// @file
2 /// Test polynomial gradient to arbitrary points in multiple dimensions
3 /// \test Test polynomial gradient to arbitrary points in multiple dimensions
4 #include <ceed.h>
5 #include <math.h>
6 #include <stdio.h>
7
Eval(CeedInt dim,const CeedScalar x[])8 static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) {
9 CeedScalar result = 1, center = 0.1;
10 for (CeedInt d = 0; d < dim; d++) {
11 result *= tanh(x[d] - center);
12 center += 0.1;
13 }
14 return result;
15 }
16
EvalGrad(CeedInt dim,const CeedScalar x[],CeedInt direction)17 static CeedScalar EvalGrad(CeedInt dim, const CeedScalar x[], CeedInt direction) {
18 CeedScalar result = 1, center = 0.1;
19 for (CeedInt d = 0; d < dim; d++) {
20 if (d == direction) result *= 1.0 / cosh(x[d] - center) / cosh(x[d] - center);
21 else result *= tanh(x[d] - center);
22 center += 0.1;
23 }
24 return result;
25 }
26
main(int argc,char ** argv)27 int main(int argc, char **argv) {
28 Ceed ceed;
29
30 CeedInit(argv[1], &ceed);
31
32 for (CeedInt dim = 1; dim <= 3; dim++) {
33 CeedVector x, x_nodes, x_points, u, v;
34 CeedBasis basis_x, basis_u;
35 const CeedInt p = 9, q = 9, num_points = 4, x_dim = CeedIntPow(2, dim), p_dim = CeedIntPow(p, dim);
36
37 CeedVectorCreate(ceed, x_dim * dim, &x);
38 CeedVectorCreate(ceed, p_dim * dim, &x_nodes);
39 CeedVectorCreate(ceed, num_points * dim, &x_points);
40 CeedVectorCreate(ceed, p_dim, &u);
41 CeedVectorCreate(ceed, num_points * dim, &v);
42
43 // Get nodal coordinates
44 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p, CEED_GAUSS_LOBATTO, &basis_x);
45 {
46 CeedScalar x_array[x_dim * dim];
47
48 for (CeedInt d = 0; d < dim; d++) {
49 for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, d + 1)) / CeedIntPow(2, d) ? 1 : -1;
50 }
51 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
52 }
53 CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_nodes);
54
55 // Set values of u at nodes
56 {
57 const CeedScalar *x_array;
58 CeedScalar u_array[p_dim];
59
60 CeedVectorGetArrayRead(x_nodes, CEED_MEM_HOST, &x_array);
61 for (CeedInt i = 0; i < p_dim; i++) {
62 CeedScalar coord[dim];
63
64 for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * p_dim + i];
65 u_array[i] = Eval(dim, coord);
66 }
67 CeedVectorRestoreArrayRead(x_nodes, &x_array);
68 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array);
69 }
70
71 // Interpolate to arbitrary points
72 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u);
73 {
74 CeedScalar x_array[12] = {-0.33, -0.65, 0.16, 0.99, -0.65, 0.16, 0.99, -0.33, 0.16, 0.99, -0.33, -0.65};
75
76 CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
77 }
78 CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, x_points, u, v);
79
80 {
81 const CeedScalar *x_array, *v_array;
82
83 CeedVectorGetArrayRead(x_points, CEED_MEM_HOST, &x_array);
84 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
85 for (CeedInt i = 0; i < num_points; i++) {
86 CeedScalar coord[dim];
87
88 for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * num_points + i];
89 for (CeedInt d = 0; d < dim; d++) {
90 const CeedScalar dfx = EvalGrad(dim, coord, d);
91
92 if (fabs(v_array[d * num_points + i] - dfx) > 1E-3) {
93 // LCOV_EXCL_START
94 printf("[%" CeedInt_FMT "] %f != %f = df(%f", dim, v_array[d * num_points + i], dfx, coord[0]);
95 for (CeedInt d = 1; d < dim; d++) printf(", %f", coord[d]);
96 printf(")\n");
97 // LCOV_EXCL_STOP
98 }
99 }
100 }
101 CeedVectorRestoreArrayRead(x_points, &x_array);
102 CeedVectorRestoreArrayRead(v, &v_array);
103 }
104
105 CeedVectorDestroy(&x);
106 CeedVectorDestroy(&x_nodes);
107 CeedVectorDestroy(&x_points);
108 CeedVectorDestroy(&u);
109 CeedVectorDestroy(&v);
110 CeedBasisDestroy(&basis_x);
111 CeedBasisDestroy(&basis_u);
112 }
113
114 CeedDestroy(&ceed);
115 return 0;
116 }
117