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 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 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, dim - d)) / CeedIntPow(2, dim - d - 1) ? 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