1 /// @file
2 /// Test creation, action, and destruction for diffusion matrix operator
3 /// \test Test creation, action, and destruction for diffusion matrix operator
4 #include "t522-operator.h"
5
6 #include <ceed.h>
7 #include <math.h>
8 #include <stdio.h>
9 #include <stdlib.h>
10
11 #include "t320-basis.h"
12
13 /* The mesh comprises of two rows of 3 quadrilaterals followed by one row
14 of 6 triangles:
15 _ _ _
16 |_|_|_|
17 |_|_|_|
18 |/|/|/|
19
20 */
21
main(int argc,char ** argv)22 int main(int argc, char **argv) {
23 Ceed ceed;
24 CeedElemRestriction elem_restriction_x_tet, elem_restriction_u_tet, elem_restriction_q_data_tet, elem_restriction_x_hex, elem_restriction_u_hex,
25 elem_restriction_q_data_hex;
26 CeedBasis basis_x_tet, basis_u_tet, basis_x_hex, basis_u_hex;
27 CeedQFunction qf_setup_tet, qf_diff_tet, qf_setup_hex, qf_diff_hex;
28 CeedOperator op_setup_tet, op_diff_tet, op_setup_hex, op_diff_hex, op_setup, op_diff;
29 CeedVector q_data_tet, q_data_hex, x, u, v;
30 CeedInt num_elem_tet = 6, p_tet = 6, q_tet = 4, num_elem_hex = 6, p_hex = 3, q_hex = 4, dim = 2;
31 CeedInt n_x = 3, n_y = 3, n_x_tet = 3, n_y_tet = 1, n_x_hex = 3;
32 CeedInt row, col, offset;
33 CeedInt num_dofs = (n_x * 2 + 1) * (n_y * 2 + 1), num_qpts_tet = num_elem_tet * q_tet, num_qpts_hex = num_elem_hex * q_hex * q_hex;
34 CeedInt ind_x_tet[num_elem_tet * p_tet], ind_x_hex[num_elem_hex * p_hex * p_hex];
35 CeedScalar q_ref[dim * q_tet], q_weight[q_tet];
36 CeedScalar interp[p_tet * q_tet], grad[dim * p_tet * q_tet];
37
38 CeedInit(argv[1], &ceed);
39
40 // Vectors
41 CeedVectorCreate(ceed, dim * num_dofs, &x);
42 {
43 CeedScalar x_array[dim * num_dofs];
44
45 for (CeedInt i = 0; i < n_y * 2 + 1; i++) {
46 for (CeedInt j = 0; j < n_x * 2 + 1; j++) {
47 x_array[i + j * (n_y * 2 + 1) + 0 * num_dofs] = (CeedScalar)i / (2 * n_y);
48 x_array[i + j * (n_y * 2 + 1) + 1 * num_dofs] = (CeedScalar)j / (2 * n_x);
49 }
50 }
51 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
52 }
53 CeedVectorCreate(ceed, num_dofs, &u);
54 CeedVectorCreate(ceed, num_dofs, &v);
55 CeedVectorCreate(ceed, num_qpts_tet * dim * (dim + 1) / 2, &q_data_tet);
56 CeedVectorCreate(ceed, num_qpts_hex * dim * (dim + 1) / 2, &q_data_hex);
57
58 // Tet Elements
59 // -- Restrictions
60 for (CeedInt i = 0; i < num_elem_tet / 2; i++) {
61 col = i % n_x_tet;
62 row = i / n_x_tet;
63 offset = col * 2 + row * (n_x_tet * 2 + 1) * 2;
64
65 ind_x_tet[i * 2 * p_tet + 0] = 2 + offset;
66 ind_x_tet[i * 2 * p_tet + 1] = 9 + offset;
67 ind_x_tet[i * 2 * p_tet + 2] = 16 + offset;
68 ind_x_tet[i * 2 * p_tet + 3] = 1 + offset;
69 ind_x_tet[i * 2 * p_tet + 4] = 8 + offset;
70 ind_x_tet[i * 2 * p_tet + 5] = 0 + offset;
71
72 ind_x_tet[i * 2 * p_tet + 6] = 14 + offset;
73 ind_x_tet[i * 2 * p_tet + 7] = 7 + offset;
74 ind_x_tet[i * 2 * p_tet + 8] = 0 + offset;
75 ind_x_tet[i * 2 * p_tet + 9] = 15 + offset;
76 ind_x_tet[i * 2 * p_tet + 10] = 8 + offset;
77 ind_x_tet[i * 2 * p_tet + 11] = 16 + offset;
78 }
79 CeedElemRestrictionCreate(ceed, num_elem_tet, p_tet, dim, num_dofs, dim * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x_tet,
80 &elem_restriction_x_tet);
81 CeedElemRestrictionCreate(ceed, num_elem_tet, p_tet, 1, 1, num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x_tet, &elem_restriction_u_tet);
82
83 CeedInt strides_q_data_tet[3] = {1, q_tet, q_tet * dim * (dim + 1) / 2};
84 CeedElemRestrictionCreateStrided(ceed, num_elem_tet, q_tet, dim * (dim + 1) / 2, dim * (dim + 1) / 2 * num_qpts_tet, strides_q_data_tet,
85 &elem_restriction_q_data_tet);
86
87 // -- Bases
88 Build2DSimplex(q_ref, q_weight, interp, grad);
89 CeedBasisCreateH1(ceed, CEED_TOPOLOGY_TRIANGLE, dim, p_tet, q_tet, interp, grad, q_ref, q_weight, &basis_x_tet);
90
91 Build2DSimplex(q_ref, q_weight, interp, grad);
92 CeedBasisCreateH1(ceed, CEED_TOPOLOGY_TRIANGLE, 1, p_tet, q_tet, interp, grad, q_ref, q_weight, &basis_u_tet);
93
94 // -- QFunctions
95 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup_tet);
96 CeedQFunctionAddInput(qf_setup_tet, "weight", 1, CEED_EVAL_WEIGHT);
97 CeedQFunctionAddInput(qf_setup_tet, "dx", dim * dim, CEED_EVAL_GRAD);
98 CeedQFunctionAddOutput(qf_setup_tet, "rho", dim * (dim + 1) / 2, CEED_EVAL_NONE);
99
100 CeedQFunctionCreateInterior(ceed, 1, diff, diff_loc, &qf_diff_tet);
101 CeedQFunctionAddInput(qf_diff_tet, "rho", dim * (dim + 1) / 2, CEED_EVAL_NONE);
102 CeedQFunctionAddInput(qf_diff_tet, "u", dim, CEED_EVAL_GRAD);
103 CeedQFunctionAddOutput(qf_diff_tet, "v", dim, CEED_EVAL_GRAD);
104
105 // -- Operators
106 // ---- Setup Tet
107 CeedOperatorCreate(ceed, qf_setup_tet, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_tet);
108 CeedOperatorSetField(op_setup_tet, "weight", CEED_ELEMRESTRICTION_NONE, basis_x_tet, CEED_VECTOR_NONE);
109 CeedOperatorSetField(op_setup_tet, "dx", elem_restriction_x_tet, basis_x_tet, CEED_VECTOR_ACTIVE);
110 CeedOperatorSetField(op_setup_tet, "rho", elem_restriction_q_data_tet, CEED_BASIS_NONE, q_data_tet);
111 // ---- Diff Tet
112 CeedOperatorCreate(ceed, qf_diff_tet, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_diff_tet);
113 CeedOperatorSetField(op_diff_tet, "rho", elem_restriction_q_data_tet, CEED_BASIS_NONE, q_data_tet);
114 CeedOperatorSetField(op_diff_tet, "u", elem_restriction_u_tet, basis_u_tet, CEED_VECTOR_ACTIVE);
115 CeedOperatorSetField(op_diff_tet, "v", elem_restriction_u_tet, basis_u_tet, CEED_VECTOR_ACTIVE);
116
117 // Hex Elements
118 // -- Restrictions
119 for (CeedInt i = 0; i < num_elem_hex; i++) {
120 col = i % n_x_hex;
121 row = i / n_x_hex;
122 offset = (n_x_tet * 2 + 1) * (n_y_tet * 2) * (1 + row) + col * 2;
123 for (CeedInt j = 0; j < p_hex; j++) {
124 for (CeedInt k = 0; k < p_hex; k++) ind_x_hex[p_hex * (p_hex * i + k) + j] = offset + k * (n_x_hex * 2 + 1) + j;
125 }
126 }
127 CeedElemRestrictionCreate(ceed, num_elem_hex, p_hex * p_hex, dim, num_dofs, dim * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x_hex,
128 &elem_restriction_x_hex);
129 CeedElemRestrictionCreate(ceed, num_elem_hex, p_hex * p_hex, 1, 1, num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x_hex, &elem_restriction_u_hex);
130
131 CeedInt strides_q_data_hex[3] = {1, q_hex * q_hex, q_hex * q_hex * dim * (dim + 1) / 2};
132 CeedElemRestrictionCreateStrided(ceed, num_elem_hex, q_hex * q_hex, dim * (dim + 1) / 2, dim * (dim + 1) / 2 * num_qpts_hex, strides_q_data_hex,
133 &elem_restriction_q_data_hex);
134
135 // -- Bases
136 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, p_hex, q_hex, CEED_GAUSS, &basis_x_hex);
137 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p_hex, q_hex, CEED_GAUSS, &basis_u_hex);
138
139 // -- QFunctions
140 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup_hex);
141 CeedQFunctionAddInput(qf_setup_hex, "weight", 1, CEED_EVAL_WEIGHT);
142 CeedQFunctionAddInput(qf_setup_hex, "dx", dim * dim, CEED_EVAL_GRAD);
143 CeedQFunctionAddOutput(qf_setup_hex, "rho", dim * (dim + 1) / 2, CEED_EVAL_NONE);
144
145 CeedQFunctionCreateInterior(ceed, 1, diff, diff_loc, &qf_diff_hex);
146 CeedQFunctionAddInput(qf_diff_hex, "rho", dim * (dim + 1) / 2, CEED_EVAL_NONE);
147 CeedQFunctionAddInput(qf_diff_hex, "u", dim, CEED_EVAL_GRAD);
148 CeedQFunctionAddOutput(qf_diff_hex, "v", dim, CEED_EVAL_GRAD);
149
150 // -- Operators
151 CeedOperatorCreate(ceed, qf_setup_hex, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_hex);
152 CeedOperatorSetField(op_setup_hex, "weight", CEED_ELEMRESTRICTION_NONE, basis_x_hex, CEED_VECTOR_NONE);
153 CeedOperatorSetField(op_setup_hex, "dx", elem_restriction_x_hex, basis_x_hex, CEED_VECTOR_ACTIVE);
154 CeedOperatorSetField(op_setup_hex, "rho", elem_restriction_q_data_hex, CEED_BASIS_NONE, q_data_hex);
155
156 CeedOperatorCreate(ceed, qf_diff_hex, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_diff_hex);
157 CeedOperatorSetField(op_diff_hex, "rho", elem_restriction_q_data_hex, CEED_BASIS_NONE, q_data_hex);
158 CeedOperatorSetField(op_diff_hex, "u", elem_restriction_u_hex, basis_u_hex, CEED_VECTOR_ACTIVE);
159 CeedOperatorSetField(op_diff_hex, "v", elem_restriction_u_hex, basis_u_hex, CEED_VECTOR_ACTIVE);
160
161 // Composite Operators
162 CeedOperatorCreateComposite(ceed, &op_setup);
163 CeedOperatorCompositeAddSub(op_setup, op_setup_tet);
164 CeedOperatorCompositeAddSub(op_setup, op_setup_hex);
165
166 CeedOperatorCreateComposite(ceed, &op_diff);
167 CeedOperatorCompositeAddSub(op_diff, op_diff_tet);
168 CeedOperatorCompositeAddSub(op_diff, op_diff_hex);
169
170 // Apply Setup Operator
171 CeedOperatorApply(op_setup, x, CEED_VECTOR_NONE, CEED_REQUEST_IMMEDIATE);
172
173 // Apply diff Operator
174 CeedVectorSetValue(u, 1.0);
175 CeedOperatorApply(op_diff, u, v, CEED_REQUEST_IMMEDIATE);
176
177 // Check output
178 {
179 const CeedScalar *v_array;
180
181 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
182 for (CeedInt i = 0; i < num_dofs; i++) {
183 if (fabs(v_array[i]) > 100. * CEED_EPSILON) printf("Computed: %f != True: 0.0\n", v_array[i]);
184 }
185 CeedVectorRestoreArrayRead(v, &v_array);
186 }
187
188 // Cleanup
189 CeedVectorDestroy(&x);
190 CeedVectorDestroy(&u);
191 CeedVectorDestroy(&v);
192 CeedVectorDestroy(&q_data_tet);
193 CeedVectorDestroy(&q_data_hex);
194 CeedElemRestrictionDestroy(&elem_restriction_u_tet);
195 CeedElemRestrictionDestroy(&elem_restriction_x_tet);
196 CeedElemRestrictionDestroy(&elem_restriction_q_data_tet);
197 CeedElemRestrictionDestroy(&elem_restriction_u_hex);
198 CeedElemRestrictionDestroy(&elem_restriction_x_hex);
199 CeedElemRestrictionDestroy(&elem_restriction_q_data_hex);
200 CeedBasisDestroy(&basis_u_tet);
201 CeedBasisDestroy(&basis_x_tet);
202 CeedBasisDestroy(&basis_u_hex);
203 CeedBasisDestroy(&basis_x_hex);
204 CeedQFunctionDestroy(&qf_setup_tet);
205 CeedQFunctionDestroy(&qf_diff_tet);
206 CeedOperatorDestroy(&op_setup_tet);
207 CeedOperatorDestroy(&op_diff_tet);
208 CeedQFunctionDestroy(&qf_setup_hex);
209 CeedQFunctionDestroy(&qf_diff_hex);
210 CeedOperatorDestroy(&op_setup_hex);
211 CeedOperatorDestroy(&op_diff_hex);
212 CeedOperatorDestroy(&op_setup);
213 CeedOperatorDestroy(&op_diff);
214 CeedDestroy(&ceed);
215 return 0;
216 }
217