1 /// @file 2 /// Test polynomial interpolation to arbitrary points in multiple dimensions 3 /// \test Test polynomial interpolation to arbitrary points 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_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, num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); 73 74 for (CeedInt i = 0; i < num_points; i++) { 75 CeedScalar fx = 0.0; 76 CeedScalar coord[dim]; 77 const CeedScalar *x_array, *u_array, *v_array, *u_point_array; 78 79 CeedVectorGetArrayRead(x_points, CEED_MEM_HOST, &x_array); 80 CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); 81 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); 82 for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d + i * dim]; 83 CeedVectorSetArray(x_point, CEED_MEM_HOST, CEED_COPY_VALUES, coord); 84 CeedBasisApplyAtPoints(basis_u, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, x_point, v_point, u_point); 85 CeedVectorGetArrayRead(u_point, CEED_MEM_HOST, &u_point_array); 86 for (CeedInt j = 0; j < p_dim; j++) fx += u_array[j] * u_point_array[j]; 87 if (fabs(v_array[i] - fx) > 100. * CEED_EPSILON) { 88 // LCOV_EXCL_START 89 printf("[%" CeedInt_FMT "] %f != %f = f(%f", dim, v_array[i], fx, coord[0]); 90 for (CeedInt d = 1; d < dim; d++) printf(", %f", coord[d]); 91 puts(")"); 92 // LCOV_EXCL_STOP 93 } 94 CeedVectorRestoreArrayRead(u_point, &u_point_array); 95 CeedVectorRestoreArrayRead(x_points, &x_array); 96 CeedVectorRestoreArrayRead(u, &u_array); 97 CeedVectorRestoreArrayRead(v, &v_array); 98 } 99 100 CeedVectorDestroy(&x); 101 CeedVectorDestroy(&x_nodes); 102 CeedVectorDestroy(&x_points); 103 CeedVectorDestroy(&x_point); 104 CeedVectorDestroy(&u); 105 CeedVectorDestroy(&v); 106 CeedVectorDestroy(&u_point); 107 CeedVectorDestroy(&v_point); 108 CeedBasisDestroy(&basis_x); 109 CeedBasisDestroy(&basis_u); 110 } 111 112 CeedDestroy(&ceed); 113 return 0; 114 } 115