1 /// @file 2 /// Test Collocated Grad calculated matches basis with Lobatto points 3 /// \test Test Collocated Grad calculated matches basis with Lobatto points 4 #include <ceed.h> 5 #include <ceed/backend.h> 6 #include <math.h> 7 #include <stdio.h> 8 9 int main(int argc, char **argv) { 10 Ceed ceed; 11 CeedInt p = 4; 12 CeedScalar collocated_gradient_1d[(p + 2) * (p + 2)], x_2[p + 2]; 13 const CeedScalar *gradient_1d, *q_ref; 14 CeedScalar sum = 0.0; 15 CeedBasis basis; 16 17 CeedInit(argv[1], &ceed); 18 19 // Already collocated, GetCollocatedGrad will return grad_1d 20 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, p, CEED_GAUSS_LOBATTO, &basis); 21 CeedBasisGetCollocatedGrad(basis, collocated_gradient_1d); 22 CeedBasisGetGrad(basis, &gradient_1d); 23 24 for (CeedInt i = 0; i < p; i++) { 25 for (CeedInt j = 0; j < p; j++) { 26 if (fabs(collocated_gradient_1d[j + p * i] - gradient_1d[j + p * i]) > 100 * CEED_EPSILON) { 27 // LCOV_EXCL_START 28 printf("Error in collocated gradient %f != %f\n", collocated_gradient_1d[j + p * i], gradient_1d[j + p * i]); 29 // LCOV_EXCL_START 30 } 31 } 32 } 33 CeedBasisDestroy(&basis); 34 35 // Q = P, not already collocated 36 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, p, CEED_GAUSS, &basis); 37 CeedBasisGetCollocatedGrad(basis, collocated_gradient_1d); 38 39 CeedBasisGetQRef(basis, &q_ref); 40 for (CeedInt i = 0; i < p; i++) x_2[i] = q_ref[i] * q_ref[i]; 41 42 // Verify collo_grad * x^2 = 2x for quadrature points 43 for (CeedInt i = 0; i < p; i++) { 44 sum = 0.0; 45 for (CeedInt j = 0; j < p; j++) sum += collocated_gradient_1d[j + p * i] * x_2[j]; 46 if (fabs(sum - 2 * q_ref[i]) > 100 * CEED_EPSILON) printf("Error in collocated gradient %f != %f\n", sum, 2 * q_ref[i]); 47 } 48 CeedBasisDestroy(&basis); 49 50 // Q = P + 2, not already collocated 51 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, p + 2, CEED_GAUSS, &basis); 52 CeedBasisGetCollocatedGrad(basis, collocated_gradient_1d); 53 54 CeedBasisGetQRef(basis, &q_ref); 55 for (CeedInt i = 0; i < p + 2; i++) x_2[i] = q_ref[i] * q_ref[i]; 56 57 // Verify collocated_gradient * x^2 = 2x for quadrature points 58 for (CeedInt i = 0; i < p + 2; i++) { 59 sum = 0.0; 60 for (CeedInt j = 0; j < p + 2; j++) sum += collocated_gradient_1d[j + (p + 2) * i] * x_2[j]; 61 if (fabs(sum - 2 * q_ref[i]) > 100 * CEED_EPSILON) printf("Error in collocated gradient %f != %f\n", sum, 2 * q_ref[i]); 62 } 63 CeedBasisDestroy(&basis); 64 65 CeedDestroy(&ceed); 66 return 0; 67 } 68