1 /// @file
2 /// Test full assembly of a non-square mass matrix operator (see t553)
3 /// \test Test full assembly of a non-square mass matrix operator
4 #include <ceed.h>
5 #include <math.h>
6 #include <stdio.h>
7 #include <stdlib.h>
8
main(int argc,char ** argv)9 int main(int argc, char **argv) {
10 Ceed ceed;
11 CeedElemRestriction elem_restriction_x, elem_restriction_q_data, elem_restriction_u_coarse, elem_restriction_u_fine;
12 CeedBasis basis_x, basis_u_coarse, basis_u_fine;
13 CeedQFunction qf_setup, qf_mass;
14 CeedOperator op_setup, op_mass;
15 CeedVector q_data, x, u, v;
16 CeedInt num_elem = 15, p_coarse = 3, p_fine = 5, q = 8;
17 CeedInt num_dofs_x = num_elem + 1, num_dofs_u_coarse = num_elem * (p_coarse - 1) + 1, num_dofs_u_fine = num_elem * (p_fine - 1) + 1;
18 CeedInt ind_u_coarse[num_elem * p_coarse], ind_u_fine[num_elem * p_fine], ind_x[num_elem * 2];
19 CeedScalar assembled_values[num_dofs_u_coarse * num_dofs_u_fine];
20 CeedScalar assembled_true[num_dofs_u_coarse * num_dofs_u_fine];
21
22 CeedInit(argv[1], &ceed);
23
24 CeedVectorCreate(ceed, num_dofs_x, &x);
25 {
26 CeedScalar x_array[num_dofs_x];
27
28 for (CeedInt i = 0; i < num_dofs_x; i++) x_array[i] = (CeedScalar)i / (num_dofs_x - 1);
29 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
30 }
31 CeedVectorCreate(ceed, num_dofs_u_coarse, &u);
32 CeedVectorCreate(ceed, num_dofs_u_fine, &v);
33 CeedVectorCreate(ceed, num_elem * q, &q_data);
34
35 // Restrictions
36 for (CeedInt i = 0; i < num_elem; i++) {
37 ind_x[2 * i + 0] = i;
38 ind_x[2 * i + 1] = i + 1;
39 }
40 CeedElemRestrictionCreate(ceed, num_elem, 2, 1, 1, num_dofs_x, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x);
41
42 for (CeedInt i = 0; i < num_elem; i++) {
43 for (CeedInt j = 0; j < p_coarse; j++) {
44 ind_u_coarse[p_coarse * i + j] = i * (p_coarse - 1) + j;
45 }
46 }
47 CeedElemRestrictionCreate(ceed, num_elem, p_coarse, 1, 1, num_dofs_u_coarse, CEED_MEM_HOST, CEED_USE_POINTER, ind_u_coarse,
48 &elem_restriction_u_coarse);
49
50 for (CeedInt i = 0; i < num_elem; i++) {
51 for (CeedInt j = 0; j < p_fine; j++) {
52 ind_u_fine[p_fine * i + j] = i * (p_fine - 1) + j;
53 }
54 }
55 CeedElemRestrictionCreate(ceed, num_elem, p_fine, 1, 1, num_dofs_u_fine, CEED_MEM_HOST, CEED_USE_POINTER, ind_u_fine, &elem_restriction_u_fine);
56
57 CeedInt strides_q_data[3] = {1, q, q};
58 CeedElemRestrictionCreateStrided(ceed, num_elem, q, 1, q * num_elem, strides_q_data, &elem_restriction_q_data);
59
60 // Bases
61 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS, &basis_x);
62 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p_coarse, q, CEED_GAUSS, &basis_u_coarse);
63 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p_fine, q, CEED_GAUSS, &basis_u_fine);
64
65 // QFunctions
66 CeedQFunctionCreateInteriorByName(ceed, "Mass1DBuild", &qf_setup);
67 CeedQFunctionCreateInteriorByName(ceed, "MassApply", &qf_mass);
68
69 // Operators
70 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup);
71
72 CeedOperatorSetField(op_setup, "weights", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
73 CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE);
74 CeedOperatorSetField(op_setup, "qdata", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
75
76 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass);
77 CeedOperatorSetField(op_mass, "qdata", elem_restriction_q_data, CEED_BASIS_NONE, q_data);
78 CeedOperatorSetField(op_mass, "u", elem_restriction_u_coarse, basis_u_coarse, CEED_VECTOR_ACTIVE);
79 CeedOperatorSetField(op_mass, "v", elem_restriction_u_fine, basis_u_fine, CEED_VECTOR_ACTIVE);
80
81 CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE);
82
83 // Fully assemble operator
84 CeedSize num_entries;
85 CeedInt *rows;
86 CeedInt *cols;
87 CeedVector assembled;
88
89 for (CeedInt k = 0; k < num_dofs_u_coarse * num_dofs_u_fine; ++k) {
90 assembled_values[k] = 0.0;
91 assembled_true[k] = 0.0;
92 }
93 CeedOperatorLinearAssembleSymbolic(op_mass, &num_entries, &rows, &cols);
94 CeedVectorCreate(ceed, num_entries, &assembled);
95 CeedOperatorLinearAssemble(op_mass, assembled);
96 {
97 const CeedScalar *assembled_array;
98
99 CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array);
100 for (CeedInt k = 0; k < num_entries; ++k) {
101 assembled_values[rows[k] * num_dofs_u_coarse + cols[k]] += assembled_array[k];
102 }
103 CeedVectorRestoreArrayRead(assembled, &assembled_array);
104 }
105
106 // Manually assemble operator
107 CeedVectorSetValue(u, 0.0);
108 for (CeedInt j = 0; j < num_dofs_u_coarse; j++) {
109 CeedScalar *u_array;
110 const CeedScalar *v_array;
111
112 // Set input
113 CeedVectorGetArray(u, CEED_MEM_HOST, &u_array);
114 u_array[j] = 1.0;
115 if (j) u_array[j - 1] = 0.0;
116 CeedVectorRestoreArray(u, &u_array);
117
118 // Compute entries for column j
119 CeedOperatorApply(op_mass, u, v, CEED_REQUEST_IMMEDIATE);
120
121 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
122 for (CeedInt i = 0; i < num_dofs_u_fine; i++) assembled_true[i * num_dofs_u_coarse + j] = v_array[i];
123 CeedVectorRestoreArrayRead(v, &v_array);
124 }
125
126 // Check output
127 for (CeedInt i = 0; i < num_dofs_u_fine; i++) {
128 for (CeedInt j = 0; j < num_dofs_u_coarse; j++) {
129 if (fabs(assembled_values[i * num_dofs_u_coarse + j] - assembled_true[i * num_dofs_u_coarse + j]) > 100. * CEED_EPSILON) {
130 // LCOV_EXCL_START
131 printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] Error in assembly: %f != %f\n", i, j, assembled_values[i * num_dofs_u_coarse + j],
132 assembled_true[i * num_dofs_u_coarse + j]);
133 // LCOV_EXCL_STOP
134 }
135 }
136 }
137
138 // Cleanup
139 free(rows);
140 free(cols);
141 CeedVectorDestroy(&x);
142 CeedVectorDestroy(&q_data);
143 CeedVectorDestroy(&u);
144 CeedVectorDestroy(&v);
145 CeedVectorDestroy(&assembled);
146 CeedElemRestrictionDestroy(&elem_restriction_u_coarse);
147 CeedElemRestrictionDestroy(&elem_restriction_u_fine);
148 CeedElemRestrictionDestroy(&elem_restriction_x);
149 CeedElemRestrictionDestroy(&elem_restriction_q_data);
150 CeedBasisDestroy(&basis_u_coarse);
151 CeedBasisDestroy(&basis_u_fine);
152 CeedBasisDestroy(&basis_x);
153 CeedQFunctionDestroy(&qf_setup);
154 CeedQFunctionDestroy(&qf_mass);
155 CeedOperatorDestroy(&op_setup);
156 CeedOperatorDestroy(&op_mass);
157 CeedDestroy(&ceed);
158 return 0;
159 }
160