1 /// @file 2 /// Test interpolation ApplyAdd in multiple dimensions 3 /// \test Test interpolation ApplyAdd in multiple dimensions 4 #include <ceed.h> 5 #include <math.h> 6 #include <stdio.h> 7 8 int main(int argc, char **argv) { 9 Ceed ceed; 10 11 CeedInit(argv[1], &ceed); 12 13 for (CeedInt dim = 1; dim <= 3; dim++) { 14 CeedVector u, u_q, v, v_q, w_q; 15 CeedBasis basis; 16 CeedInt p = 4, q = 5, p_dim = CeedIntPow(p, dim), q_dim = CeedIntPow(q, dim); 17 18 CeedVectorCreate(ceed, p_dim, &u); 19 CeedVectorCreate(ceed, p_dim, &v); 20 CeedVectorSetValue(u, 1.0); 21 CeedVectorSetValue(v, 0.0); 22 CeedVectorCreate(ceed, q_dim, &u_q); 23 CeedVectorCreate(ceed, q_dim, &v_q); 24 CeedVectorCreate(ceed, q_dim, &w_q); 25 26 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis); 27 28 // Compute area 29 CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, u_q); 30 CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, w_q); 31 CeedVectorPointwiseMult(v_q, u_q, w_q); 32 CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, v_q, v); 33 // Double area computed 34 CeedBasisApplyAdd(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, v_q, v); 35 36 // Check area 37 { 38 const CeedScalar *v_array; 39 CeedScalar area = 0.0; 40 41 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); 42 for (CeedInt i = 0; i < p_dim; i++) area += v_array[i]; 43 if (fabs(area - 2.0 * CeedIntPow(2, dim)) > 5E-6) printf("Incorrect area computed %f != %f\n", area, 2.0 * CeedIntPow(2, dim)); 44 CeedVectorRestoreArrayRead(v, &v_array); 45 } 46 47 CeedVectorDestroy(&u); 48 CeedVectorDestroy(&v); 49 CeedVectorDestroy(&u_q); 50 CeedVectorDestroy(&v_q); 51 CeedVectorDestroy(&w_q); 52 CeedBasisDestroy(&basis); 53 } 54 CeedDestroy(&ceed); 55 return 0; 56 } 57