xref: /libCEED/tests/t415-qfunction.c (revision a2e5d304d0c7d96eecfcbbd32f1ea5194beb84ca)
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, dx, w, du, dv;
12   CeedQFunction qf_setup, qf_diff;
13   CeedInt       q        = 8;
14   const CeedInt num_comp = 3;
15 
16   CeedInit(argv[1], &ceed);
17 
18   for (CeedInt dim = 1; dim <= 3; dim++) {
19     CeedInt num_qpts = CeedIntPow(q, dim);
20 
21     CeedVectorCreate(ceed, num_qpts * dim * dim, &dx);
22     CeedVectorCreate(ceed, num_qpts, &w);
23     CeedVectorCreate(ceed, num_qpts * dim * num_comp, &du);
24     {
25       CeedScalar dx_array[num_qpts * dim * dim], w_array[num_qpts], du_array[num_qpts * dim * num_comp];
26 
27       for (CeedInt i = 0; i < num_qpts; i++) {
28         w_array[i] = 1.0 / num_qpts;
29       }
30       for (CeedInt d = 0; d < dim; d++) {
31         for (CeedInt g = 0; g < dim; g++) {
32           for (CeedInt i = 0; i < num_qpts; i++) {
33             dx_array[i + (g * dim + d) * num_qpts] = d == g;
34           }
35         }
36       }
37       for (CeedInt c = 0; c < num_comp; c++) {
38         for (CeedInt g = 0; g < dim; g++) {
39           for (CeedInt i = 0; i < num_qpts; i++) {
40             du_array[i + (g * num_comp + c) * num_qpts] = c + 1;
41           }
42         }
43       }
44 
45       CeedVectorSetArray(dx, CEED_MEM_HOST, CEED_COPY_VALUES, dx_array);
46       CeedVectorSetArray(w, CEED_MEM_HOST, CEED_COPY_VALUES, w_array);
47       CeedVectorSetArray(du, CEED_MEM_HOST, CEED_COPY_VALUES, du_array);
48     }
49     CeedVectorCreate(ceed, num_qpts * dim * (dim + 1) / 2, &q_data);
50     CeedVectorSetValue(q_data, 0.0);
51     CeedVectorCreate(ceed, num_qpts * dim * num_comp, &dv);
52     CeedVectorSetValue(dv, 0.0);
53 
54     char name_setup[26] = "", name_apply[26] = "";
55     snprintf(name_setup, sizeof name_setup, "Poisson%" CeedInt_FMT "DBuild", dim);
56     CeedQFunctionCreateInteriorByName(ceed, name_setup, &qf_setup);
57     {
58       in[0]  = dx;
59       in[1]  = w;
60       out[0] = q_data;
61       CeedQFunctionApply(qf_setup, num_qpts, in, out);
62     }
63 
64     snprintf(name_apply, sizeof name_apply, "Vector3Poisson%" CeedInt_FMT "DApply", dim);
65     CeedQFunctionCreateInteriorByName(ceed, name_apply, &qf_diff);
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     // Verify results
74     {
75       const CeedScalar *v_array;
76 
77       CeedVectorGetArrayRead(dv, CEED_MEM_HOST, &v_array);
78       for (CeedInt c = 0; c < num_comp; c++) {
79         CeedScalar sum = 0;
80         for (CeedInt i = 0; i < num_qpts; i++) {
81           for (CeedInt g = 0; g < dim; g++) sum += v_array[i + (g * num_comp + c) * num_qpts];
82         }
83         if (fabs(sum - dim * (c + 1)) > 10 * CEED_EPSILON) {
84           // LCOV_EXCL_START
85           printf("%" CeedInt_FMT "D volume error in component %" CeedInt_FMT ": %f != %f\n", dim, c, sum, dim * (c + 1.0));
86           // LCOV_EXCL_STOP
87         }
88       }
89       CeedVectorRestoreArrayRead(dv, &v_array);
90     }
91 
92     CeedVectorDestroy(&dx);
93     CeedVectorDestroy(&w);
94     CeedVectorDestroy(&du);
95     CeedVectorDestroy(&dv);
96     CeedVectorDestroy(&q_data);
97     CeedQFunctionDestroy(&qf_setup);
98     CeedQFunctionDestroy(&qf_diff);
99   }
100 
101   CeedDestroy(&ceed);
102   return 0;
103 }
104