1 /// @file 2 /// Test polynomial interpolation to arbirtary points with multiple components in multiple dimensions 3 /// \test Test polynomial interpolation to arbitrary points with multiple components in multiple dimensions 4 #include <ceed.h> 5 #include <math.h> 6 #include <stdio.h> 7 8 static CeedScalar Eval(CeedInt dim, CeedScalar scale, 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 scale * 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, u, v; 24 CeedBasis basis_x, basis_u; 25 const CeedInt p = 9, q = 9, num_comp = 3, 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, num_comp * p_dim, &u); 31 CeedVectorCreate(ceed, num_comp * num_points, &v); 32 33 // Get nodal coordinates 34 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p, CEED_GAUSS_LOBATTO, &basis_x); 35 { 36 CeedScalar x_array[x_dim * dim]; 37 38 for (CeedInt d = 0; d < dim; d++) { 39 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; 40 } 41 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 42 } 43 CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_nodes); 44 45 // Set values of u at nodes 46 { 47 const CeedScalar *x_array; 48 CeedScalar u_array[num_comp * p_dim]; 49 50 CeedVectorGetArrayRead(x_nodes, CEED_MEM_HOST, &x_array); 51 for (CeedInt i = 0; i < p_dim; i++) { 52 CeedScalar coord[dim]; 53 54 for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * p_dim + i]; 55 for (CeedInt c = 0; c < num_comp; c++) u_array[i + c * p_dim] = Eval(dim, c, coord); 56 } 57 CeedVectorRestoreArrayRead(x_nodes, &x_array); 58 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array); 59 } 60 61 // Interpolate to arbitrary points 62 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, p, q, CEED_GAUSS, &basis_u); 63 { 64 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}; 65 66 CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 67 } 68 CeedBasisApplyAtPoints(basis_u, num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); 69 70 { 71 const CeedScalar *x_array, *v_array; 72 73 CeedVectorGetArrayRead(x_points, CEED_MEM_HOST, &x_array); 74 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); 75 for (CeedInt i = 0; i < num_points; i++) { 76 CeedScalar coord[dim]; 77 78 for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d + i * dim]; 79 for (CeedInt c = 0; c < num_comp; c++) { 80 CeedScalar fx = Eval(dim, c, coord); 81 if (fabs(v_array[c + i * num_comp] - fx) > 1E-4) { 82 // LCOV_EXCL_START 83 printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] %f != %f = f(%f", dim, c, v_array[c + i * num_comp], fx, coord[0]); 84 for (CeedInt d = 1; d < dim; d++) printf(", %f", coord[d]); 85 puts(")"); 86 // LCOV_EXCL_STOP 87 } 88 } 89 } 90 CeedVectorRestoreArrayRead(x_points, &x_array); 91 CeedVectorRestoreArrayRead(v, &v_array); 92 } 93 94 CeedVectorDestroy(&x); 95 CeedVectorDestroy(&x_nodes); 96 CeedVectorDestroy(&x_points); 97 CeedVectorDestroy(&u); 98 CeedVectorDestroy(&v); 99 CeedBasisDestroy(&basis_x); 100 CeedBasisDestroy(&basis_u); 101 } 102 103 CeedDestroy(&ceed); 104 return 0; 105 } 106