1 /// @file
2 /// Test high order interpolation in multiple dimensions
3 /// \test Test high order interpolation 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 for (CeedInt dim = 1; dim <= 3; dim++) {
22 CeedVector x, x_q, u, u_q;
23 CeedBasis basis_x_lobatto, basis_u_lobatto, basis_x_gauss, basis_u_gauss;
24 CeedInt q = 11, q_dim = CeedIntPow(q, dim), x_dim = CeedIntPow(2, dim);
25
26 CeedVectorCreate(ceed, x_dim * dim, &x);
27 {
28 CeedScalar x_array[x_dim * dim];
29
30 for (CeedInt d = 0; d < dim; d++) {
31 for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, d + 1)) / CeedIntPow(2, d) ? 1 : -1;
32 }
33 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
34 }
35 CeedVectorCreate(ceed, q_dim * dim, &x_q);
36 CeedVectorSetValue(x_q, 0);
37 CeedVectorCreate(ceed, q_dim, &u);
38 CeedVectorSetValue(u, 0);
39 CeedVectorCreate(ceed, q_dim, &u_q);
40
41 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, q, CEED_GAUSS_LOBATTO, &basis_x_lobatto);
42 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, q, q, CEED_GAUSS_LOBATTO, &basis_u_lobatto);
43
44 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q);
45
46 {
47 const CeedScalar *x_q_array;
48 CeedScalar u_q_array[q_dim];
49
50 CeedVectorGetArrayRead(x_q, CEED_MEM_HOST, &x_q_array);
51 for (CeedInt i = 0; i < q_dim; i++) {
52 CeedScalar coord[dim];
53 for (CeedInt d = 0; d < dim; d++) coord[d] = x_q_array[d * q_dim + i];
54 u_q_array[i] = Eval(dim, coord);
55 }
56 CeedVectorRestoreArrayRead(x_q, &x_q_array);
57 CeedVectorSetArray(u_q, CEED_MEM_HOST, CEED_COPY_VALUES, u_q_array);
58 }
59
60 // This operation is the identity because the quadrature is collocated
61 CeedBasisApply(basis_u_lobatto, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, u_q, u);
62
63 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, q, CEED_GAUSS, &basis_x_gauss);
64 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, q, q, CEED_GAUSS, &basis_u_gauss);
65
66 CeedBasisApply(basis_x_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q);
67 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, u_q);
68
69 {
70 const CeedScalar *x_q_array, *u_array;
71
72 CeedVectorGetArrayRead(x_q, CEED_MEM_HOST, &x_q_array);
73 CeedVectorGetArrayRead(u_q, CEED_MEM_HOST, &u_array);
74 for (CeedInt i = 0; i < q_dim; i++) {
75 CeedScalar coord[dim];
76 for (CeedInt d = 0; d < dim; d++) coord[d] = x_q_array[d * q_dim + i];
77 CeedScalar fx = Eval(dim, coord);
78 if (fabs(u_array[i] - fx) > 1E-4) {
79 // LCOV_EXCL_START
80 printf("[%" CeedInt_FMT "] %f != %f = f(%f", dim, u_array[i], fx, coord[0]);
81 for (CeedInt d = 1; d < dim; d++) printf(",%f", coord[d]);
82 puts(")");
83 // LCOV_EXCL_STOP
84 }
85 }
86 CeedVectorRestoreArrayRead(x_q, &x_q_array);
87 CeedVectorRestoreArrayRead(u_q, &u_array);
88 }
89
90 CeedVectorDestroy(&x);
91 CeedVectorDestroy(&x_q);
92 CeedVectorDestroy(&u);
93 CeedVectorDestroy(&u_q);
94 CeedBasisDestroy(&basis_x_lobatto);
95 CeedBasisDestroy(&basis_u_lobatto);
96 CeedBasisDestroy(&basis_x_gauss);
97 CeedBasisDestroy(&basis_u_gauss);
98 }
99 CeedDestroy(&ceed);
100 return 0;
101 }
102