xref: /libCEED/tests/t415-qfunction.c (revision 9e201c85545dd39529c090846df629a32c15659b)
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%" CeedInt_FMT "DBuild", dim);
26     CeedQFunctionCreateInteriorByName(ceed, name_setup, &qf_setup);
27     snprintf(name_apply, sizeof name_apply, "Vector3Poisson%" CeedInt_FMT "DApply",
28              dim);
29     CeedQFunctionCreateInteriorByName(ceed, name_apply, &qf_diff);
30 
31     for (CeedInt i=0; i<num_qpts; i++) {
32       w[i] = 1.0 / num_qpts;
33     }
34     for (CeedInt d=0; d<dim; d++) {
35       for (CeedInt g=0; g<dim; g++) {
36         for (CeedInt i=0; i<num_qpts; i++) {
37           j[i + (g * dim + d) * num_qpts] = d == g;
38         }
39       }
40     }
41     for (CeedInt c=0; c<num_comp; c++) {
42       for (CeedInt g=0; g<dim; g++) {
43         for (CeedInt i=0; i<num_qpts; i++) {
44           du[i + (g * num_comp + c) * num_qpts] = c + 1;
45         }
46       }
47     }
48 
49     CeedVectorCreate(ceed, num_qpts * dim * dim, &J);
50     CeedVectorSetArray(J, CEED_MEM_HOST, CEED_USE_POINTER, j);
51     CeedVectorCreate(ceed, num_qpts, &W);
52     CeedVectorSetArray(W, CEED_MEM_HOST, CEED_USE_POINTER, w);
53     CeedVectorCreate(ceed, num_qpts * dim * (dim + 1) / 2, &Q_data);
54     CeedVectorSetValue(Q_data, 0.0);
55     CeedVectorCreate(ceed, num_qpts * dim * num_comp, &dU);
56     CeedVectorSetArray(dU, CEED_MEM_HOST, CEED_USE_POINTER, du);
57     CeedVectorCreate(ceed, num_qpts * dim * num_comp, &dV);
58     CeedVectorSetValue(dV, 0.0);
59 
60     {
61       in[0] = J;
62       in[1] = W;
63       out[0] = Q_data;
64       CeedQFunctionApply(qf_setup, num_qpts, in, out);
65     }
66     {
67       in[0] = dU;
68       in[1] = Q_data;
69       out[0] = dV;
70       CeedQFunctionApply(qf_diff, num_qpts, in, out);
71     }
72 
73     CeedVectorGetArrayRead(dV, CEED_MEM_HOST, &vv);
74     for (CeedInt c=0; c<num_comp; c++) {
75       CeedScalar sum = 0;
76       for (CeedInt i=0; i<num_qpts; i++)
77         for (CeedInt g=0; g<dim; g++)
78           sum += vv[i + (g * num_comp + c) * num_qpts];
79       if (fabs(sum - dim * (c + 1)) > 10*CEED_EPSILON)
80         // LCOV_EXCL_START
81         printf("%" CeedInt_FMT "D volume error in component %" CeedInt_FMT
82                ": %f != %f\n", dim, c, sum, dim * (c + 1.0));
83       // LCOV_EXCL_STOP
84     }
85     CeedVectorRestoreArrayRead(dV, &vv);
86 
87     CeedVectorDestroy(&J);
88     CeedVectorDestroy(&W);
89     CeedVectorDestroy(&dU);
90     CeedVectorDestroy(&dV);
91     CeedVectorDestroy(&Q_data);
92     CeedQFunctionDestroy(&qf_setup);
93     CeedQFunctionDestroy(&qf_diff);
94   }
95 
96   CeedDestroy(&ceed);
97   return 0;
98 }
99