1 /// @file 2 /// Test polynomial gradient to arbitrary points in multiple dimensions 3 /// \test Test polynomial gradient 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 static CeedScalar EvalGrad(CeedInt dim, const CeedScalar x[], CeedInt direction) { 18 CeedScalar result = 1, center = 0.1; 19 for (CeedInt d = 0; d < dim; d++) { 20 if (d == direction) result *= 1.0 / cosh(x[d] - center) / cosh(x[d] - center); 21 else result *= tanh(x[d] - center); 22 center += 0.1; 23 } 24 return result; 25 } 26 27 int main(int argc, char **argv) { 28 Ceed ceed; 29 30 CeedInit(argv[1], &ceed); 31 32 for (CeedInt dim = 1; dim <= 3; dim++) { 33 CeedVector x, x_nodes, x_points, u, v; 34 CeedBasis basis_x, basis_u; 35 const CeedInt p = 9, q = 9, num_points = 4, x_dim = CeedIntPow(2, dim), p_dim = CeedIntPow(p, dim); 36 37 CeedVectorCreate(ceed, x_dim * dim, &x); 38 CeedVectorCreate(ceed, p_dim * dim, &x_nodes); 39 CeedVectorCreate(ceed, num_points * dim, &x_points); 40 CeedVectorCreate(ceed, p_dim, &u); 41 CeedVectorCreate(ceed, num_points * dim, &v); 42 43 // Get nodal coordinates 44 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p, CEED_GAUSS_LOBATTO, &basis_x); 45 { 46 CeedScalar x_array[x_dim * dim]; 47 48 for (CeedInt d = 0; d < dim; d++) { 49 for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, d + 1)) / CeedIntPow(2, d) ? 1 : -1; 50 } 51 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 52 } 53 CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_nodes); 54 55 // Set values of u at nodes 56 { 57 const CeedScalar *x_array; 58 CeedScalar u_array[p_dim]; 59 60 CeedVectorGetArrayRead(x_nodes, CEED_MEM_HOST, &x_array); 61 for (CeedInt i = 0; i < p_dim; i++) { 62 CeedScalar coord[dim]; 63 64 for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * p_dim + i]; 65 u_array[i] = Eval(dim, coord); 66 } 67 CeedVectorRestoreArrayRead(x_nodes, &x_array); 68 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array); 69 } 70 71 // Interpolate to arbitrary points 72 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u); 73 { 74 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}; 75 76 CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 77 } 78 CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, x_points, u, v); 79 80 { 81 const CeedScalar *x_array, *v_array; 82 83 CeedVectorGetArrayRead(x_points, CEED_MEM_HOST, &x_array); 84 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); 85 for (CeedInt i = 0; i < num_points; i++) { 86 CeedScalar coord[dim]; 87 88 for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * num_points + i]; 89 for (CeedInt d = 0; d < dim; d++) { 90 const CeedScalar dfx = EvalGrad(dim, coord, d); 91 92 if (fabs(v_array[d * num_points + i] - dfx) > 1E-3) { 93 // LCOV_EXCL_START 94 printf("[%" CeedInt_FMT "] %f != %f = df(%f", dim, v_array[d * num_points + i], dfx, coord[0]); 95 for (CeedInt d = 1; d < dim; d++) printf(", %f", coord[d]); 96 printf(")\n"); 97 // LCOV_EXCL_STOP 98 } 99 } 100 } 101 CeedVectorRestoreArrayRead(x_points, &x_array); 102 CeedVectorRestoreArrayRead(v, &v_array); 103 } 104 105 CeedVectorDestroy(&x); 106 CeedVectorDestroy(&x_nodes); 107 CeedVectorDestroy(&x_points); 108 CeedVectorDestroy(&u); 109 CeedVectorDestroy(&v); 110 CeedBasisDestroy(&basis_x); 111 CeedBasisDestroy(&basis_u); 112 } 113 114 CeedDestroy(&ceed); 115 return 0; 116 } 117