14411cf47Sjeremylt /// @file 252bfb9bbSJeremy L Thompson /// Test Collocated Grad calculated matches basis with Lobatto points 352bfb9bbSJeremy L Thompson /// \test Test Collocated Grad calculated matches basis with Lobatto points 457c64913Sjeremylt #include <ceed.h> 5ec3da8bcSJed Brown #include <ceed/backend.h> 657c64913Sjeremylt #include <math.h> 757c64913Sjeremylt 857c64913Sjeremylt int main(int argc, char **argv) { 957c64913Sjeremylt Ceed ceed; 10*9ac7b42eSJeremy L Thompson CeedInt P = 4; 11*9ac7b42eSJeremy L Thompson CeedScalar collo_grad_1d[(P+2)*(P+2)], x_2[P+2]; 12*9ac7b42eSJeremy L Thompson const CeedScalar *grad_1d, *q_ref; 13*9ac7b42eSJeremy L Thompson CeedScalar sum = 0.0; 1452bfb9bbSJeremy L Thompson CeedBasis b; 1557c64913Sjeremylt 1657c64913Sjeremylt CeedInit(argv[1], &ceed); 17aedaa0e5Sjeremylt 18d1d35e2fSjeremylt // Already collocated, GetCollocatedGrad will return grad_1d 19*9ac7b42eSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, P, CEED_GAUSS_LOBATTO, &b); 20d1d35e2fSjeremylt CeedBasisGetCollocatedGrad(b, collo_grad_1d); 21*9ac7b42eSJeremy L Thompson CeedBasisGetGrad(b, &grad_1d); 22aedaa0e5Sjeremylt 23*9ac7b42eSJeremy L Thompson for (CeedInt i=0; i<P; i++) 24*9ac7b42eSJeremy L Thompson for (CeedInt j=0; j<P; j++) 25*9ac7b42eSJeremy L Thompson if (fabs(collo_grad_1d[j+P*i] - grad_1d[j+P*i]) > 100*CEED_EPSILON) 26*9ac7b42eSJeremy L Thompson // LCOV_EXCL_START 27*9ac7b42eSJeremy L Thompson printf("Error in collocated gradient %f != %f\n", collo_grad_1d[j+P*i], 28*9ac7b42eSJeremy L Thompson grad_1d[j+P*i]); 29*9ac7b42eSJeremy L Thompson // LCOV_EXCL_START 3052bfb9bbSJeremy L Thompson CeedBasisDestroy(&b); 3157c64913Sjeremylt 3252bfb9bbSJeremy L Thompson // Q = P, not already collocated 33*9ac7b42eSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, P, CEED_GAUSS, &b); 34d1d35e2fSjeremylt CeedBasisGetCollocatedGrad(b, collo_grad_1d); 3552bfb9bbSJeremy L Thompson 36*9ac7b42eSJeremy L Thompson CeedBasisGetQRef(b, &q_ref); 37*9ac7b42eSJeremy L Thompson for (CeedInt i=0; i<P; i++) 38*9ac7b42eSJeremy L Thompson x_2[i] = q_ref[i]*q_ref[i]; 39*9ac7b42eSJeremy L Thompson 40*9ac7b42eSJeremy L Thompson // Verify collo_grad * x^2 = 2x for quadrature points 41*9ac7b42eSJeremy L Thompson for (CeedInt i=0; i<P; i++) { 42*9ac7b42eSJeremy L Thompson sum = 0.0; 43*9ac7b42eSJeremy L Thompson for (CeedInt j=0; j<P; j++) 44*9ac7b42eSJeremy L Thompson sum += collo_grad_1d[j+P*i]*x_2[j]; 45*9ac7b42eSJeremy L Thompson if (fabs(sum - 2*q_ref[i]) > 100*CEED_EPSILON) 46*9ac7b42eSJeremy L Thompson // LCOV_EXCL_START 47*9ac7b42eSJeremy L Thompson printf("Error in collocated gradient %f != %f\n", sum, 2*q_ref[i]); 48*9ac7b42eSJeremy L Thompson // LCOV_EXCL_START 4952bfb9bbSJeremy L Thompson } 5052bfb9bbSJeremy L Thompson CeedBasisDestroy(&b); 5152bfb9bbSJeremy L Thompson 5252bfb9bbSJeremy L Thompson // Q = P + 2, not already collocated 5352bfb9bbSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, P+2, CEED_GAUSS, &b); 54*9ac7b42eSJeremy L Thompson CeedBasisGetCollocatedGrad(b, collo_grad_1d); 5552bfb9bbSJeremy L Thompson 56*9ac7b42eSJeremy L Thompson CeedBasisGetQRef(b, &q_ref); 57*9ac7b42eSJeremy L Thompson for (CeedInt i=0; i<P+2; i++) 58*9ac7b42eSJeremy L Thompson x_2[i] = q_ref[i]*q_ref[i]; 59*9ac7b42eSJeremy L Thompson 60*9ac7b42eSJeremy L Thompson // Verify collo_grad * x^2 = 2x for quadrature points 61*9ac7b42eSJeremy L Thompson for (CeedInt i=0; i<P+2; i++) { 62*9ac7b42eSJeremy L Thompson sum = 0.0; 63*9ac7b42eSJeremy L Thompson for (CeedInt j=0; j<P+2; j++) 64*9ac7b42eSJeremy L Thompson sum += collo_grad_1d[j+(P+2)*i]*x_2[j]; 65*9ac7b42eSJeremy L Thompson if (fabs(sum - 2*q_ref[i]) > 100*CEED_EPSILON) 66*9ac7b42eSJeremy L Thompson // LCOV_EXCL_START 67*9ac7b42eSJeremy L Thompson printf("Error in collocated gradient %f != %f\n", sum, 2*q_ref[i]); 68*9ac7b42eSJeremy L Thompson // LCOV_EXCL_START 6952bfb9bbSJeremy L Thompson } 7052bfb9bbSJeremy L Thompson CeedBasisDestroy(&b); 7152bfb9bbSJeremy L Thompson 7257c64913Sjeremylt CeedDestroy(&ceed); 7357c64913Sjeremylt return 0; 7457c64913Sjeremylt } 75