1 /// @file
2 /// Test assembly of H(div) mass matrix operator diagonal
3 /// \test Test assembly of H(div) mass matrix operator diagonal
4 #include "t580-operator.h"
5
6 #include <ceed.h>
7 #include <math.h>
8 #include <stdio.h>
9 #include <stdlib.h>
10
11 #include "t330-basis.h"
12
main(int argc,char ** argv)13 int main(int argc, char **argv) {
14 Ceed ceed;
15 CeedElemRestriction elem_restriction_x, elem_restriction_u;
16 CeedBasis basis_x, basis_u;
17 CeedQFunction qf_mass;
18 CeedOperator op_mass;
19 CeedVector x, assembled, u, v;
20 CeedInt dim = 2, p = 8, q = 3, px = 2;
21 CeedInt n_x = 1, n_y = 1; // Currently only implemented for single element
22 CeedInt row, column, offset;
23 CeedInt num_elem = n_x * n_y, num_faces = (n_x + 1) * n_y + (n_y + 1) * n_x;
24 CeedInt num_dofs_x = (n_x + 1) * (n_y + 1), num_dofs_u = num_faces * 2, num_qpts = q * q;
25 CeedInt ind_x[num_elem * px * px], ind_u[num_elem * p];
26 bool orient_u[num_elem * p];
27 CeedScalar assembled_true[num_dofs_u];
28 CeedScalar q_ref[dim * num_qpts], q_weight[num_qpts];
29 CeedScalar interp[dim * p * num_qpts], div[p * num_qpts];
30
31 CeedInit(argv[1], &ceed);
32
33 // Vectors
34 CeedVectorCreate(ceed, dim * num_dofs_x, &x);
35 {
36 CeedScalar x_array[dim * num_dofs_x];
37
38 for (CeedInt i = 0; i < n_x + 1; i++) {
39 for (CeedInt j = 0; j < n_y + 1; j++) {
40 x_array[i + j * (n_x + 1) + 0 * num_dofs_x] = i / (CeedScalar)n_x;
41 x_array[i + j * (n_x + 1) + 1 * num_dofs_x] = j / (CeedScalar)n_y;
42 }
43 }
44 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
45 }
46 CeedVectorCreate(ceed, num_dofs_u, &u);
47 CeedVectorCreate(ceed, num_dofs_u, &v);
48
49 // Restrictions
50 for (CeedInt i = 0; i < num_elem; i++) {
51 column = i % n_x;
52 row = i / n_x;
53 offset = column * (px - 1) + row * (n_x + 1) * (px - 1);
54
55 for (CeedInt j = 0; j < px; j++) {
56 for (CeedInt k = 0; k < px; k++) {
57 ind_x[px * (px * i + k) + j] = offset + k * (n_x + 1) + j;
58 }
59 }
60 }
61 bool orient_u_local[8] = {false, false, false, false, true, true, true, true};
62 CeedInt ind_u_local[8] = {0, 1, 6, 7, 2, 3, 4, 5};
63 for (CeedInt j = 0; j < n_y; j++) {
64 for (CeedInt i = 0; i < n_x; i++) {
65 for (CeedInt k = 0; k < p; k++) {
66 ind_u[k] = ind_u_local[k];
67 orient_u[k] = orient_u_local[k];
68 }
69 }
70 }
71 CeedElemRestrictionCreate(ceed, num_elem, px * px, dim, num_dofs_x, dim * num_dofs_x, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x);
72 CeedElemRestrictionCreateOriented(ceed, num_elem, p, 1, 1, num_dofs_u, CEED_MEM_HOST, CEED_COPY_VALUES, ind_u, orient_u, &elem_restriction_u);
73
74 // Bases
75 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, px, q, CEED_GAUSS, &basis_x);
76
77 BuildHdivQuadrilateral(q, q_ref, q_weight, interp, div, CEED_GAUSS);
78 CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, 1, p, num_qpts, interp, div, q_ref, q_weight, &basis_u);
79
80 // QFunctions
81 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass);
82 CeedQFunctionAddInput(qf_mass, "weight", 1, CEED_EVAL_WEIGHT);
83 CeedQFunctionAddInput(qf_mass, "dx", dim * dim, CEED_EVAL_GRAD);
84 CeedQFunctionAddInput(qf_mass, "u", dim, CEED_EVAL_INTERP);
85 CeedQFunctionAddOutput(qf_mass, "v", dim, CEED_EVAL_INTERP);
86
87 // Operators
88 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass);
89 CeedOperatorSetField(op_mass, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
90 CeedOperatorSetField(op_mass, "dx", elem_restriction_x, basis_x, x);
91 CeedOperatorSetField(op_mass, "u", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
92 CeedOperatorSetField(op_mass, "v", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
93
94 // Assemble diagonal
95 CeedVectorCreate(ceed, num_dofs_u, &assembled);
96 CeedOperatorLinearAssembleDiagonal(op_mass, assembled, CEED_REQUEST_IMMEDIATE);
97
98 // Manually assemble diagonal
99 CeedVectorSetValue(u, 0.0);
100 for (int i = 0; i < num_dofs_u; i++) {
101 CeedScalar *u_array;
102 const CeedScalar *v_array;
103
104 // Set input
105 CeedVectorGetArray(u, CEED_MEM_HOST, &u_array);
106 u_array[i] = 1.0;
107 if (i) u_array[i - 1] = 0.0;
108 CeedVectorRestoreArray(u, &u_array);
109
110 // Compute diag entry for DoF i
111 CeedOperatorApply(op_mass, u, v, CEED_REQUEST_IMMEDIATE);
112
113 // Retrieve entry
114 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
115 assembled_true[i] = v_array[i];
116 CeedVectorRestoreArrayRead(v, &v_array);
117 }
118
119 // Check output
120 {
121 const CeedScalar *assembled_array;
122
123 CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array);
124 for (int i = 0; i < num_dofs_u; i++) {
125 if (fabs(assembled_array[i] - assembled_true[i]) > 100. * CEED_EPSILON) {
126 // LCOV_EXCL_START
127 printf("[%" CeedInt_FMT "] Error in assembly: %f != %f\n", i, assembled_array[i], assembled_true[i]);
128 // LCOV_EXCL_STOP
129 }
130 }
131 CeedVectorRestoreArrayRead(assembled, &assembled_array);
132 }
133
134 // Cleanup
135 CeedVectorDestroy(&x);
136 CeedVectorDestroy(&assembled);
137 CeedVectorDestroy(&u);
138 CeedVectorDestroy(&v);
139 CeedElemRestrictionDestroy(&elem_restriction_u);
140 CeedElemRestrictionDestroy(&elem_restriction_x);
141 CeedBasisDestroy(&basis_u);
142 CeedBasisDestroy(&basis_x);
143 CeedQFunctionDestroy(&qf_mass);
144 CeedOperatorDestroy(&op_mass);
145 CeedDestroy(&ceed);
146 return 0;
147 }
148