1 /// @file 2 /// Test underintegrated grad in multiple dimensions 3 /// \test Test underintegraded 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 static CeedScalar GetTolerance(CeedScalarType scalar_type, int dim) { 15 CeedScalar tol; 16 if (scalar_type == CEED_SCALAR_FP32) { 17 if (dim == 3) tol = 0.005; 18 else tol = 1.e-4; 19 } else { 20 tol = 1.e-11; 21 } 22 return tol; 23 } 24 25 int main(int argc, char **argv) { 26 Ceed ceed; 27 28 CeedInit(argv[1], &ceed); 29 30 for (CeedInt dim = 1; dim <= 3; dim++) { 31 CeedVector x, x_q, u, u_q, ones, v; 32 CeedBasis basis_x_lobatto, basis_u_gauss; 33 CeedInt p = 8, q = 7, p_dim = CeedIntPow(p, dim), q_dim = CeedIntPow(q, dim), x_dim = CeedIntPow(2, dim); 34 CeedScalar sum_1 = 0, sum_2 = 0; 35 36 CeedVectorCreate(ceed, x_dim * dim, &x); 37 { 38 CeedScalar x_array[x_dim * dim]; 39 40 for (CeedInt d = 0; d < dim; d++) { 41 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; 42 } 43 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 44 } 45 CeedVectorCreate(ceed, p_dim * dim, &x_q); 46 CeedVectorSetValue(x_q, 0); 47 CeedVectorCreate(ceed, p_dim, &u); 48 CeedVectorCreate(ceed, q_dim * dim, &u_q); 49 CeedVectorSetValue(u_q, 0); 50 CeedVectorCreate(ceed, q_dim * dim, &ones); 51 CeedVectorSetValue(ones, 1); 52 CeedVectorCreate(ceed, p_dim, &v); 53 CeedVectorSetValue(v, 0); 54 55 // Get function values at quadrature points 56 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p, CEED_GAUSS_LOBATTO, &basis_x_lobatto); 57 CeedBasisApply(basis_x_lobatto, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_q); 58 59 { 60 const CeedScalar *x_q_array; 61 CeedScalar u_array[p_dim]; 62 63 CeedVectorGetArrayRead(x_q, CEED_MEM_HOST, &x_q_array); 64 for (CeedInt i = 0; i < p_dim; i++) { 65 CeedScalar coord[dim]; 66 for (CeedInt d = 0; d < dim; d++) coord[d] = x_q_array[d * p_dim + i]; 67 u_array[i] = Eval(dim, coord); 68 } 69 CeedVectorRestoreArrayRead(x_q, &x_q_array); 70 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, u_array); 71 } 72 73 // Calculate G u at quadrature points, G' * 1 at dofs 74 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u_gauss); 75 CeedBasisApply(basis_u_gauss, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, u, u_q); 76 CeedBasisApply(basis_u_gauss, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, ones, v); 77 78 // Check if 1' * G * u = u' * (G' * 1) 79 { 80 const CeedScalar *u_array, *v_array, *u_q_array; 81 82 CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); 83 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); 84 CeedVectorGetArrayRead(u_q, CEED_MEM_HOST, &u_q_array); 85 for (CeedInt i = 0; i < p_dim; i++) sum_1 += v_array[i] * u_array[i]; 86 for (CeedInt i = 0; i < dim * q_dim; i++) sum_2 += u_q_array[i]; 87 CeedVectorRestoreArrayRead(u, &u_array); 88 CeedVectorRestoreArrayRead(v, &v_array); 89 CeedVectorRestoreArrayRead(u_q, &u_q_array); 90 } 91 CeedScalar tol = GetTolerance(CEED_SCALAR_TYPE, dim); 92 if (fabs(sum_1 - sum_2) > tol) printf("[%" CeedInt_FMT "] %f != %f\n", dim, sum_1, sum_2); 93 94 CeedVectorDestroy(&x); 95 CeedVectorDestroy(&x_q); 96 CeedVectorDestroy(&u); 97 CeedVectorDestroy(&u_q); 98 CeedVectorDestroy(&ones); 99 CeedVectorDestroy(&v); 100 CeedBasisDestroy(&basis_x_lobatto); 101 CeedBasisDestroy(&basis_u_gauss); 102 } 103 CeedDestroy(&ceed); 104 return 0; 105 } 106