1 /// @file 2 /// Test curl with a 2D Simplex non-tensor H(curl) basis 3 /// \test Test curl with a 2D Simplex non-tensor H(curl) basis 4 #include <ceed.h> 5 #include <math.h> 6 7 #include "t340-basis.h" 8 9 int main(int argc, char **argv) { 10 Ceed ceed; 11 CeedVector u, v; 12 const CeedInt p = 8, q = 4, dim = 2; 13 CeedBasis basis; 14 CeedScalar q_ref[dim * q], q_weight[q]; 15 CeedScalar interp[dim * p * q], curl[p * q]; 16 CeedScalar row_sum[dim * q], column_sum[p]; 17 18 CeedInit(argv[1], &ceed); 19 20 BuildHcurl2DSimplex(q_ref, q_weight, interp, curl); 21 CeedBasisCreateHcurl(ceed, CEED_TOPOLOGY_TRIANGLE, 1, p, q, interp, curl, q_ref, q_weight, &basis); 22 23 // Test curl for H(curl) 24 { 25 const CeedScalar *curl_in_basis; 26 27 CeedBasisGetCurl(basis, &curl_in_basis); 28 for (CeedInt i = 0; i < p * q; i++) { 29 if (fabs(curl[i] - curl_in_basis[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", curl[i], curl_in_basis[i]); 30 } 31 } 32 33 for (int i = 0; i < q; i++) { 34 row_sum[i] = 0; 35 for (int j = 0; j < p; j++) { 36 row_sum[i] += curl[j + i * p]; 37 } 38 } 39 for (int i = 0; i < p; i++) { 40 column_sum[i] = 0; 41 for (int j = 0; j < q; j++) { 42 column_sum[i] += curl[i + j * p]; 43 } 44 } 45 46 CeedVectorCreate(ceed, p, &u); 47 CeedVectorSetValue(u, 1.0); 48 CeedVectorCreate(ceed, q, &v); 49 CeedVectorSetValue(v, 0.0); 50 51 CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_CURL, u, v); 52 53 { 54 const CeedScalar *v_array; 55 56 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); 57 for (CeedInt i = 0; i < q; i++) { 58 if (fabs(row_sum[i] - v_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", row_sum[i], v_array[i]); 59 } 60 CeedVectorRestoreArrayRead(v, &v_array); 61 } 62 63 CeedVectorSetValue(v, 1.0); 64 CeedVectorSetValue(u, 0.0); 65 66 CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_CURL, v, u); 67 68 { 69 const CeedScalar *u_array; 70 71 CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); 72 for (CeedInt i = 0; i < p; i++) { 73 if (fabs(column_sum[i] - u_array[i]) > 100. * CEED_EPSILON) printf("%f != %f\n", column_sum[i], u_array[i]); 74 } 75 CeedVectorRestoreArrayRead(u, &u_array); 76 } 77 78 CeedBasisDestroy(&basis); 79 CeedVectorDestroy(&u); 80 CeedVectorDestroy(&v); 81 CeedDestroy(&ceed); 82 return 0; 83 } 84