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; 109ac7b42eSJeremy L Thompson CeedInt P = 4; 119ac7b42eSJeremy L Thompson CeedScalar collo_grad_1d[(P + 2) * (P + 2)], x_2[P + 2]; 129ac7b42eSJeremy L Thompson const CeedScalar *grad_1d, *q_ref; 139ac7b42eSJeremy 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 199ac7b42eSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, P, CEED_GAUSS_LOBATTO, &b); 20d1d35e2fSjeremylt CeedBasisGetCollocatedGrad(b, collo_grad_1d); 219ac7b42eSJeremy L Thompson CeedBasisGetGrad(b, &grad_1d); 22aedaa0e5Sjeremylt 23*2b730f8bSJeremy L Thompson for (CeedInt i = 0; i < P; i++) { 24*2b730f8bSJeremy L Thompson for (CeedInt j = 0; j < P; j++) { 25*2b730f8bSJeremy L Thompson if (fabs(collo_grad_1d[j + P * i] - grad_1d[j + P * i]) > 100 * CEED_EPSILON) { 269ac7b42eSJeremy L Thompson // LCOV_EXCL_START 27*2b730f8bSJeremy L Thompson printf("Error in collocated gradient %f != %f\n", collo_grad_1d[j + P * i], grad_1d[j + P * i]); 289ac7b42eSJeremy L Thompson // LCOV_EXCL_START 29*2b730f8bSJeremy L Thompson } 30*2b730f8bSJeremy L Thompson } 31*2b730f8bSJeremy L Thompson } 3252bfb9bbSJeremy L Thompson CeedBasisDestroy(&b); 3357c64913Sjeremylt 3452bfb9bbSJeremy L Thompson // Q = P, not already collocated 359ac7b42eSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, P, CEED_GAUSS, &b); 36d1d35e2fSjeremylt CeedBasisGetCollocatedGrad(b, collo_grad_1d); 3752bfb9bbSJeremy L Thompson 389ac7b42eSJeremy L Thompson CeedBasisGetQRef(b, &q_ref); 39*2b730f8bSJeremy L Thompson for (CeedInt i = 0; i < P; i++) x_2[i] = q_ref[i] * q_ref[i]; 409ac7b42eSJeremy L Thompson 419ac7b42eSJeremy L Thompson // Verify collo_grad * x^2 = 2x for quadrature points 429ac7b42eSJeremy L Thompson for (CeedInt i = 0; i < P; i++) { 439ac7b42eSJeremy L Thompson sum = 0.0; 44*2b730f8bSJeremy L Thompson for (CeedInt j = 0; j < P; j++) sum += collo_grad_1d[j + P * i] * x_2[j]; 45*2b730f8bSJeremy L Thompson if (fabs(sum - 2 * q_ref[i]) > 100 * CEED_EPSILON) printf("Error in collocated gradient %f != %f\n", sum, 2 * q_ref[i]); 4652bfb9bbSJeremy L Thompson } 4752bfb9bbSJeremy L Thompson CeedBasisDestroy(&b); 4852bfb9bbSJeremy L Thompson 4952bfb9bbSJeremy L Thompson // Q = P + 2, not already collocated 5052bfb9bbSJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, P + 2, CEED_GAUSS, &b); 519ac7b42eSJeremy L Thompson CeedBasisGetCollocatedGrad(b, collo_grad_1d); 5252bfb9bbSJeremy L Thompson 539ac7b42eSJeremy L Thompson CeedBasisGetQRef(b, &q_ref); 54*2b730f8bSJeremy L Thompson for (CeedInt i = 0; i < P + 2; i++) x_2[i] = q_ref[i] * q_ref[i]; 559ac7b42eSJeremy L Thompson 569ac7b42eSJeremy L Thompson // Verify collo_grad * x^2 = 2x for quadrature points 579ac7b42eSJeremy L Thompson for (CeedInt i = 0; i < P + 2; i++) { 589ac7b42eSJeremy L Thompson sum = 0.0; 59*2b730f8bSJeremy L Thompson for (CeedInt j = 0; j < P + 2; j++) sum += collo_grad_1d[j + (P + 2) * i] * x_2[j]; 60*2b730f8bSJeremy L Thompson if (fabs(sum - 2 * q_ref[i]) > 100 * CEED_EPSILON) printf("Error in collocated gradient %f != %f\n", sum, 2 * q_ref[i]); 6152bfb9bbSJeremy L Thompson } 6252bfb9bbSJeremy L Thompson CeedBasisDestroy(&b); 6352bfb9bbSJeremy L Thompson 6457c64913Sjeremylt CeedDestroy(&ceed); 6557c64913Sjeremylt return 0; 6657c64913Sjeremylt } 67