xref: /libCEED/tests/t316-basis.c (revision 52bfb9bbf17f17edbcd45876cdc8689a879bc683) 
1*52bfb9bbSJeremy L Thompson /// @file
2*52bfb9bbSJeremy L Thompson /// Test underintegrated grad in multiple dimensions
3*52bfb9bbSJeremy L Thompson /// \test Test underintegraded grad in multiple dimensions
4*52bfb9bbSJeremy L Thompson #include <ceed.h>
5*52bfb9bbSJeremy L Thompson #include <math.h>
6*52bfb9bbSJeremy L Thompson 
7*52bfb9bbSJeremy L Thompson static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) {
8*52bfb9bbSJeremy L Thompson   CeedScalar result = tanh(x[0] + 0.1);
9*52bfb9bbSJeremy L Thompson   if (dim > 1) result += atan(x[1] + 0.2);
10*52bfb9bbSJeremy L Thompson   if (dim > 2) result += exp(-(x[2] + 0.3)*(x[2] + 0.3));
11*52bfb9bbSJeremy L Thompson   return result;
12*52bfb9bbSJeremy L Thompson }
13*52bfb9bbSJeremy L Thompson 
14*52bfb9bbSJeremy L Thompson int main(int argc, char **argv) {
15*52bfb9bbSJeremy L Thompson   Ceed ceed;
16*52bfb9bbSJeremy L Thompson 
17*52bfb9bbSJeremy L Thompson   CeedInit(argv[1], &ceed);
18*52bfb9bbSJeremy L Thompson 
19*52bfb9bbSJeremy L Thompson   for (CeedInt dim=1; dim<=3; dim++) {
20*52bfb9bbSJeremy L Thompson     CeedVector X, Xq, U, Uq, Ones, Gtposeones;
21*52bfb9bbSJeremy L Thompson     CeedBasis bxl, bug;
22*52bfb9bbSJeremy L Thompson     CeedInt P = 8, Q = 7, Pdim = CeedIntPow(P, dim), Qdim = CeedIntPow(Q, dim),
23*52bfb9bbSJeremy L Thompson             Xdim = CeedIntPow(2, dim);
24*52bfb9bbSJeremy L Thompson     CeedScalar x[Xdim*dim], u[Pdim];
25*52bfb9bbSJeremy L Thompson     const CeedScalar *xq, *uq, *gtposeones;
26*52bfb9bbSJeremy L Thompson     CeedScalar sum1 = 0, sum2 = 0;
27*52bfb9bbSJeremy L Thompson 
28*52bfb9bbSJeremy L Thompson     for (CeedInt d=0; d<dim; d++)
29*52bfb9bbSJeremy L Thompson       for (CeedInt i=0; i<Xdim; i++)
30*52bfb9bbSJeremy L Thompson         x[d*Xdim + i] = (i % CeedIntPow(2, dim-d)) /
31*52bfb9bbSJeremy L Thompson                         CeedIntPow(2, dim-d-1) ? 1 : -1;
32*52bfb9bbSJeremy L Thompson 
33*52bfb9bbSJeremy L Thompson     CeedVectorCreate(ceed, Xdim*dim, &X);
34*52bfb9bbSJeremy L Thompson     CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x);
35*52bfb9bbSJeremy L Thompson     CeedVectorCreate(ceed, Pdim*dim, &Xq);
36*52bfb9bbSJeremy L Thompson     CeedVectorSetValue(Xq, 0);
37*52bfb9bbSJeremy L Thompson     CeedVectorCreate(ceed, Pdim, &U);
38*52bfb9bbSJeremy L Thompson     CeedVectorCreate(ceed, Qdim*dim, &Uq);
39*52bfb9bbSJeremy L Thompson     CeedVectorSetValue(Uq, 0);
40*52bfb9bbSJeremy L Thompson     CeedVectorCreate(ceed, Qdim*dim, &Ones);
41*52bfb9bbSJeremy L Thompson     CeedVectorSetValue(Ones, 1);
42*52bfb9bbSJeremy L Thompson     CeedVectorCreate(ceed, Pdim, &Gtposeones);
43*52bfb9bbSJeremy L Thompson     CeedVectorSetValue(Gtposeones, 0);
44*52bfb9bbSJeremy L Thompson 
45*52bfb9bbSJeremy L Thompson     // Get function values at quadrature points
46*52bfb9bbSJeremy L Thompson     CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, P,
47*52bfb9bbSJeremy L Thompson                                     CEED_GAUSS_LOBATTO, &bxl);
48*52bfb9bbSJeremy L Thompson     CeedBasisApply(bxl, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, X, Xq);
49*52bfb9bbSJeremy L Thompson 
50*52bfb9bbSJeremy L Thompson     CeedVectorGetArrayRead(Xq, CEED_MEM_HOST, &xq);
51*52bfb9bbSJeremy L Thompson     for (CeedInt i=0; i<Pdim; i++) {
52*52bfb9bbSJeremy L Thompson       CeedScalar xx[dim];
53*52bfb9bbSJeremy L Thompson       for (CeedInt d=0; d<dim; d++)
54*52bfb9bbSJeremy L Thompson         xx[d] = xq[d*Pdim + i];
55*52bfb9bbSJeremy L Thompson       u[i] = Eval(dim, xx);
56*52bfb9bbSJeremy L Thompson     }
57*52bfb9bbSJeremy L Thompson     CeedVectorRestoreArrayRead(Xq, &xq);
58*52bfb9bbSJeremy L Thompson     CeedVectorSetArray(U, CEED_MEM_HOST, CEED_USE_POINTER, u);
59*52bfb9bbSJeremy L Thompson 
60*52bfb9bbSJeremy L Thompson     // Calculate G u at quadrature points, G' * 1 at dofs
61*52bfb9bbSJeremy L Thompson     CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, CEED_GAUSS, &bug);
62*52bfb9bbSJeremy L Thompson     CeedBasisApply(bug, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, U, Uq);
63*52bfb9bbSJeremy L Thompson     CeedBasisApply(bug, 1, CEED_TRANSPOSE, CEED_EVAL_GRAD, Ones, Gtposeones);
64*52bfb9bbSJeremy L Thompson 
65*52bfb9bbSJeremy L Thompson     // Check if 1' * G * u = u' * (G' * 1)
66*52bfb9bbSJeremy L Thompson     CeedVectorGetArrayRead(Gtposeones, CEED_MEM_HOST, &gtposeones);
67*52bfb9bbSJeremy L Thompson     CeedVectorGetArrayRead(Uq, CEED_MEM_HOST, &uq);
68*52bfb9bbSJeremy L Thompson     for (CeedInt i=0; i<Pdim; i++)
69*52bfb9bbSJeremy L Thompson       sum1 += gtposeones[i]*u[i];
70*52bfb9bbSJeremy L Thompson     for (CeedInt i=0; i<dim*Qdim; i++)
71*52bfb9bbSJeremy L Thompson       sum2 += uq[i];
72*52bfb9bbSJeremy L Thompson     CeedVectorRestoreArrayRead(Gtposeones, &gtposeones);
73*52bfb9bbSJeremy L Thompson     CeedVectorRestoreArrayRead(Uq, &uq);
74*52bfb9bbSJeremy L Thompson     if (fabs(sum1 - sum2) > 1e-10) {
75*52bfb9bbSJeremy L Thompson       printf("[%d] %f != %f\n", dim, sum1, sum2);
76*52bfb9bbSJeremy L Thompson     }
77*52bfb9bbSJeremy L Thompson 
78*52bfb9bbSJeremy L Thompson     CeedVectorDestroy(&X);
79*52bfb9bbSJeremy L Thompson     CeedVectorDestroy(&Xq);
80*52bfb9bbSJeremy L Thompson     CeedVectorDestroy(&U);
81*52bfb9bbSJeremy L Thompson     CeedVectorDestroy(&Uq);
82*52bfb9bbSJeremy L Thompson     CeedVectorDestroy(&Ones);
83*52bfb9bbSJeremy L Thompson     CeedVectorDestroy(&Gtposeones);
84*52bfb9bbSJeremy L Thompson     CeedBasisDestroy(&bxl);
85*52bfb9bbSJeremy L Thompson     CeedBasisDestroy(&bug);
86*52bfb9bbSJeremy L Thompson   }
87*52bfb9bbSJeremy L Thompson   CeedDestroy(&ceed);
88*52bfb9bbSJeremy L Thompson   return 0;
89*52bfb9bbSJeremy L Thompson }
90