xref: /libCEED/tests/t415-qfunction.c (revision 13964f0727a62e5421e6d3b433e838b96a9ce891)
1 /// @file
2 /// Test creation, evaluation, and destruction for vector Poisson QFunction by name
3 /// \test Test creation, evaluation, and destruction for vector Poisson QFunction by name
4 #include <ceed.h>
5 #include <math.h>
6 #include <string.h>
7 
8 int main(int argc, char **argv) {
9   Ceed ceed;
10   CeedVector in[16], out[16];
11   CeedVector Q_data, J, W, dU, dV;
12   CeedQFunction qf_setup, qf_diff;
13   CeedInt Q = 8;
14   const CeedInt num_comp = 3;
15   const CeedScalar *vv;
16 
17   CeedInit(argv[1], &ceed);
18 
19   for (CeedInt dim = 1; dim <= 3; dim++) {
20     CeedInt num_qpts = CeedIntPow(Q, dim);
21     CeedScalar j[num_qpts * dim * dim], w[num_qpts],
22                du[num_qpts * dim * num_comp];
23 
24     char name_setup[26] = "", name_apply[26] = "";
25     snprintf(name_setup, sizeof name_setup, "Poisson%dDBuild", dim);
26     CeedQFunctionCreateInteriorByName(ceed, name_setup, &qf_setup);
27     snprintf(name_apply, sizeof name_apply, "Vector3Poisson%dDApply", dim);
28     CeedQFunctionCreateInteriorByName(ceed, name_apply, &qf_diff);
29 
30     for (CeedInt i=0; i<num_qpts; i++) {
31       w[i] = 1.0 / num_qpts;
32     }
33     for (CeedInt d=0; d<dim; d++) {
34       for (CeedInt g=0; g<dim; g++) {
35         for (CeedInt i=0; i<num_qpts; i++) {
36           j[i + (g * dim + d) * num_qpts] = d == g;
37         }
38       }
39     }
40     for (CeedInt c=0; c<num_comp; c++) {
41       for (CeedInt g=0; g<dim; g++) {
42         for (CeedInt i=0; i<num_qpts; i++) {
43           du[i + (g * num_comp + c) * num_qpts] = c + 1;
44         }
45       }
46     }
47 
48     CeedVectorCreate(ceed, num_qpts * dim * dim, &J);
49     CeedVectorSetArray(J, CEED_MEM_HOST, CEED_USE_POINTER, j);
50     CeedVectorCreate(ceed, num_qpts, &W);
51     CeedVectorSetArray(W, CEED_MEM_HOST, CEED_USE_POINTER, w);
52     CeedVectorCreate(ceed, num_qpts * dim * (dim + 1) / 2, &Q_data);
53     CeedVectorSetValue(Q_data, 0.0);
54     CeedVectorCreate(ceed, num_qpts * dim * num_comp, &dU);
55     CeedVectorSetArray(dU, CEED_MEM_HOST, CEED_USE_POINTER, du);
56     CeedVectorCreate(ceed, num_qpts * dim * num_comp, &dV);
57     CeedVectorSetValue(dV, 0.0);
58 
59     {
60       in[0] = J;
61       in[1] = W;
62       out[0] = Q_data;
63       CeedQFunctionApply(qf_setup, num_qpts, in, out);
64     }
65     {
66       in[0] = dU;
67       in[1] = Q_data;
68       out[0] = dV;
69       CeedQFunctionApply(qf_diff, num_qpts, in, out);
70     }
71 
72     CeedVectorGetArrayRead(dV, CEED_MEM_HOST, &vv);
73     for (CeedInt c=0; c<num_comp; c++) {
74       CeedScalar sum = 0;
75       for (CeedInt i=0; i<num_qpts; i++)
76         for (CeedInt g=0; g<dim; g++)
77           sum += vv[i + (g * num_comp + c) * num_qpts];
78       if (fabs(sum - dim * (c + 1)) > 10*CEED_EPSILON)
79         // LCOV_EXCL_START
80         printf("%dD volume error in component %d: %f != %f\n",
81                dim, c, sum, dim * (c + 1.0));
82       // LCOV_EXCL_STOP
83     }
84     CeedVectorRestoreArrayRead(dV, &vv);
85 
86     CeedVectorDestroy(&J);
87     CeedVectorDestroy(&W);
88     CeedVectorDestroy(&dU);
89     CeedVectorDestroy(&dV);
90     CeedVectorDestroy(&Q_data);
91     CeedQFunctionDestroy(&qf_setup);
92     CeedQFunctionDestroy(&qf_diff);
93   }
94 
95   CeedDestroy(&ceed);
96   return 0;
97 }
98