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