1 /// @file
2 /// Test full assembly of an identity operator (see t509)
3 /// \test Test full assembly of an identity 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_u, elem_restriction_u_i;
12 CeedBasis basis_u;
13 CeedQFunction qf_identity;
14 CeedOperator op_identity;
15 CeedVector u, v;
16 CeedInt num_elem = 15, p = 5, q = 8;
17 CeedInt num_nodes = num_elem * (p - 1) + 1;
18 CeedInt ind_u[num_elem * p];
19 CeedScalar assembled_values[num_nodes * q * num_elem];
20 CeedScalar assembled_true[num_nodes * q * num_elem];
21
22 CeedInit(argv[1], &ceed);
23
24 CeedVectorCreate(ceed, num_nodes, &u);
25 CeedVectorCreate(ceed, q * num_elem, &v);
26
27 // Restrictions
28 for (CeedInt i = 0; i < num_elem; i++) {
29 for (CeedInt j = 0; j < p; j++) {
30 ind_u[p * i + j] = i * (p - 1) + j;
31 }
32 }
33 CeedElemRestrictionCreate(ceed, num_elem, p, 1, 1, num_nodes, CEED_MEM_HOST, CEED_USE_POINTER, ind_u, &elem_restriction_u);
34
35 CeedInt strides_u_i[3] = {1, q, q};
36 CeedElemRestrictionCreateStrided(ceed, num_elem, q, 1, q * num_elem, strides_u_i, &elem_restriction_u_i);
37
38 // Bases
39 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, q, CEED_GAUSS, &basis_u);
40
41 // QFunction
42 CeedQFunctionCreateIdentity(ceed, 1, CEED_EVAL_INTERP, CEED_EVAL_NONE, &qf_identity);
43
44 // Operators
45 CeedOperatorCreate(ceed, qf_identity, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_identity);
46 CeedOperatorSetField(op_identity, "input", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
47 CeedOperatorSetField(op_identity, "output", elem_restriction_u_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
48
49 // Fully assemble operator
50 CeedSize num_entries;
51 CeedInt *rows;
52 CeedInt *cols;
53 CeedVector assembled;
54
55 for (CeedInt k = 0; k < num_nodes * q * num_elem; ++k) {
56 assembled_values[k] = 0.0;
57 assembled_true[k] = 0.0;
58 }
59 CeedOperatorLinearAssembleSymbolic(op_identity, &num_entries, &rows, &cols);
60 CeedVectorCreate(ceed, num_entries, &assembled);
61 CeedOperatorLinearAssemble(op_identity, assembled);
62 {
63 const CeedScalar *assembled_array;
64
65 CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array);
66 for (CeedInt k = 0; k < num_entries; ++k) {
67 assembled_values[rows[k] * num_nodes + cols[k]] += assembled_array[k];
68 }
69 CeedVectorRestoreArrayRead(assembled, &assembled_array);
70 }
71
72 // Manually assemble operator
73 CeedVectorSetValue(u, 0.0);
74 for (CeedInt j = 0; j < num_nodes; j++) {
75 CeedScalar *u_array;
76 const CeedScalar *v_array;
77
78 // Set input
79 CeedVectorGetArray(u, CEED_MEM_HOST, &u_array);
80 u_array[j] = 1.0;
81 if (j) u_array[j - 1] = 0.0;
82 CeedVectorRestoreArray(u, &u_array);
83
84 // Compute entries for column j
85 CeedOperatorApply(op_identity, u, v, CEED_REQUEST_IMMEDIATE);
86
87 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
88 for (CeedInt i = 0; i < q * num_elem; i++) assembled_true[i * num_nodes + j] = v_array[i];
89 CeedVectorRestoreArrayRead(v, &v_array);
90 }
91
92 // Check output
93 for (CeedInt i = 0; i < q * num_elem; i++) {
94 for (CeedInt j = 0; j < num_nodes; j++) {
95 if (fabs(assembled_values[i * num_nodes + j] - assembled_true[i * num_nodes + j]) > 100. * CEED_EPSILON) {
96 // LCOV_EXCL_START
97 printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] Error in assembly: %f != %f\n", i, j, assembled_values[i * num_nodes + j],
98 assembled_true[i * num_nodes + j]);
99 // LCOV_EXCL_STOP
100 }
101 }
102 }
103
104 // Cleanup
105 free(rows);
106 free(cols);
107 CeedVectorDestroy(&u);
108 CeedVectorDestroy(&v);
109 CeedVectorDestroy(&assembled);
110 CeedElemRestrictionDestroy(&elem_restriction_u);
111 CeedElemRestrictionDestroy(&elem_restriction_u_i);
112 CeedBasisDestroy(&basis_u);
113 CeedQFunctionDestroy(&qf_identity);
114 CeedOperatorDestroy(&op_identity);
115 CeedDestroy(&ceed);
116 return 0;
117 }
118