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