1 /// @file 2 /// Test grad in multiple dimensions 3 /// \test Test grad in multiple dimensions 4 #include <ceed.h> 5 #include <math.h> 6 7 static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { 8 CeedScalar result = tanh(x[0] + 0.1); 9 if (dim > 1) result += atan(x[1] + 0.2); 10 if (dim > 2) result += exp(-(x[2] + 0.3)*(x[2] + 0.3)); 11 return result; 12 } 13 14 int main(int argc, char **argv) { 15 Ceed ceed; 16 17 CeedInit(argv[1], &ceed); 18 19 for (CeedInt dim=1; dim<=3; dim++) { 20 CeedVector X, X_q, U, U_q, ones, grad_T_ones; 21 CeedBasis basis_x_lobatto, basis_u_gauss; 22 CeedInt P = 8, Q = 10, P_dim = CeedIntPow(P, dim), Q_dim = CeedIntPow(Q, dim), 23 X_dim = CeedIntPow(2, dim); 24 CeedScalar x[X_dim*dim], u[P_dim]; 25 const CeedScalar *x_q, *u_q, *grad_t_ones_array; 26 CeedScalar sum_1 = 0, sum_2 = 0; 27 28 for (CeedInt d=0; d<dim; d++) 29 for (CeedInt i=0; i<X_dim; i++) 30 x[d*X_dim + i] = (i % CeedIntPow(2, dim-d)) / 31 CeedIntPow(2, dim-d-1) ? 1 : -1; 32 33 CeedVectorCreate(ceed, X_dim*dim, &X); 34 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 35 CeedVectorCreate(ceed, P_dim*dim, &X_q); 36 CeedVectorSetValue(X_q, 0); 37 CeedVectorCreate(ceed, P_dim, &U); 38 CeedVectorCreate(ceed, Q_dim*dim, &U_q); 39 CeedVectorSetValue(U_q, 0); 40 CeedVectorCreate(ceed, Q_dim*dim, &ones); 41 CeedVectorSetValue(ones, 1); 42 CeedVectorCreate(ceed, P_dim, &grad_T_ones); 43 CeedVectorSetValue(grad_T_ones, 0); 44 45 // Get function values at quadrature points 46 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, P, 47 CEED_GAUSS_LOBATTO, &basis_x_lobatto); 48 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, X_q); 49 50 CeedVectorGetArrayRead(X_q, CEED_MEM_HOST, &x_q); 51 for (CeedInt i=0; i<P_dim; i++) { 52 CeedScalar xx[dim]; 53 for (CeedInt d=0; d<dim; d++) 54 xx[d] = x_q[d*P_dim + i]; 55 u[i] = Eval(dim, xx); 56 } 57 CeedVectorRestoreArrayRead(X_q, &x_q); 58 CeedVectorSetArray(U, CEED_MEM_HOST, CEED_USE_POINTER, u); 59 60 // Calculate G u at quadrature points, G' * 1 at dofs 61 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, CEED_GAUSS, &basis_u_gauss); 62 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, U, U_q); 63 CeedBasisApply(basis_u_gauss, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, ones, 64 grad_T_ones); 65 66 // Check if 1' * G * u = u' * (G' * 1) 67 CeedVectorGetArrayRead(grad_T_ones, CEED_MEM_HOST, &grad_t_ones_array); 68 CeedVectorGetArrayRead(U_q, CEED_MEM_HOST, &u_q); 69 for (CeedInt i=0; i<P_dim; i++) 70 sum_1 += grad_t_ones_array[i]*u[i]; 71 for (CeedInt i=0; i<dim*Q_dim; i++) 72 sum_2 += u_q[i]; 73 CeedVectorRestoreArrayRead(grad_T_ones, &grad_t_ones_array); 74 CeedVectorRestoreArrayRead(U_q, &u_q); 75 if (fabs(sum_1 - sum_2) > 1E-10) 76 // LCOV_EXCL_START 77 printf("[%d] %f != %f\n", dim, sum_1, sum_2); 78 // LCOV_EXCL_STOP 79 80 CeedVectorDestroy(&X); 81 CeedVectorDestroy(&X_q); 82 CeedVectorDestroy(&U); 83 CeedVectorDestroy(&U_q); 84 CeedVectorDestroy(&ones); 85 CeedVectorDestroy(&grad_T_ones); 86 CeedBasisDestroy(&basis_x_lobatto); 87 CeedBasisDestroy(&basis_u_gauss); 88 } 89 CeedDestroy(&ceed); 90 return 0; 91 } 92