1 /// @file 2 /// Test full assembly of mass matrix operator with oriented and curl-oriented element restrictions (see t560) 3 /// \test Test full assembly of mass matrix operator with oriented and curl-oriented element restrictions 4 #include <ceed.h> 5 #include <math.h> 6 #include <stdio.h> 7 #include <stdlib.h> 8 9 #include "t510-operator.h" 10 11 int main(int argc, char **argv) { 12 Ceed ceed; 13 CeedElemRestriction elem_restriction_x, elem_restriction_q_data; 14 CeedElemRestriction elem_restriction_u, oriented_elem_restriction_u, curl_oriented_elem_restriction_u; 15 CeedBasis basis_x, basis_u; 16 CeedQFunction qf_setup, qf_mass; 17 CeedOperator op_setup, op_mass, op_mass_oriented, op_mass_curl_oriented; 18 CeedVector q_data, x; 19 CeedInt p = 3, q = 4, dim = 2; 20 CeedInt n_x = 3, n_y = 2; 21 CeedInt num_elem = n_x * n_y; 22 CeedInt num_dofs = (n_x * 2 + 1) * (n_y * 2 + 1), num_qpts = num_elem * q * q; 23 CeedInt ind_x[num_elem * p * p]; 24 bool orients_u[num_elem * p * p]; 25 CeedInt8 curl_orients_u[3 * num_elem * p * p]; 26 CeedScalar assembled_values[num_dofs * num_dofs]; 27 CeedScalar assembled_values_oriented[num_dofs * num_dofs]; 28 CeedScalar assembled_values_curl_oriented[num_dofs * num_dofs]; 29 30 CeedInit(argv[1], &ceed); 31 32 // Vectors 33 CeedVectorCreate(ceed, dim * num_dofs, &x); 34 { 35 CeedScalar x_array[dim * num_dofs]; 36 37 for (CeedInt i = 0; i < n_x * 2 + 1; i++) { 38 for (CeedInt j = 0; j < n_y * 2 + 1; j++) { 39 x_array[i + j * (n_x * 2 + 1) + 0 * num_dofs] = (CeedScalar)i / (2 * n_x); 40 x_array[i + j * (n_x * 2 + 1) + 1 * num_dofs] = (CeedScalar)j / (2 * n_y); 41 } 42 } 43 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); 44 } 45 CeedVectorCreate(ceed, num_qpts, &q_data); 46 47 // Restrictions 48 for (CeedInt i = 0; i < num_elem; i++) { 49 CeedInt col, row, offset; 50 col = i % n_x; 51 row = i / n_x; 52 offset = col * (p - 1) + row * (n_x * 2 + 1) * (p - 1); 53 for (CeedInt j = 0; j < p; j++) { 54 for (CeedInt k = 0; k < p; k++) { 55 ind_x[p * (p * i + k) + j] = offset + k * (n_x * 2 + 1) + j; 56 orients_u[p * (p * i + k) + j] = false; 57 curl_orients_u[3 * (p * (p * i + k) + j) + 0] = 0; 58 curl_orients_u[3 * (p * (p * i + k) + j) + 1] = 1; 59 curl_orients_u[3 * (p * (p * i + k) + j) + 2] = 0; 60 } 61 } 62 } 63 CeedElemRestrictionCreate(ceed, num_elem, p * p, dim, num_dofs, dim * num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x); 64 CeedElemRestrictionCreate(ceed, num_elem, p * p, 1, 1, num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_u); 65 CeedElemRestrictionCreateOriented(ceed, num_elem, p * p, 1, 1, num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, orients_u, 66 &oriented_elem_restriction_u); 67 CeedElemRestrictionCreateCurlOriented(ceed, num_elem, p * p, 1, 1, num_dofs, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, curl_orients_u, 68 &curl_oriented_elem_restriction_u); 69 70 CeedInt strides_q_data[3] = {1, q * q, q * q}; 71 CeedElemRestrictionCreateStrided(ceed, num_elem, q * q, 1, num_qpts, strides_q_data, &elem_restriction_q_data); 72 73 // Bases 74 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, p, q, CEED_GAUSS, &basis_x); 75 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u); 76 77 // QFunctions 78 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup); 79 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); 80 CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD); 81 CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE); 82 83 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass); 84 CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE); 85 CeedQFunctionAddInput(qf_mass, "u", 1, CEED_EVAL_INTERP); 86 CeedQFunctionAddOutput(qf_mass, "v", 1, CEED_EVAL_INTERP); 87 88 // Operators 89 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup); 90 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 91 CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE); 92 CeedOperatorSetField(op_setup, "rho", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 93 94 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass); 95 CeedOperatorSetField(op_mass, "rho", elem_restriction_q_data, CEED_BASIS_NONE, q_data); 96 CeedOperatorSetField(op_mass, "u", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); 97 CeedOperatorSetField(op_mass, "v", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); 98 99 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass_oriented); 100 CeedOperatorSetField(op_mass_oriented, "rho", elem_restriction_q_data, CEED_BASIS_NONE, q_data); 101 CeedOperatorSetField(op_mass_oriented, "u", oriented_elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); 102 CeedOperatorSetField(op_mass_oriented, "v", oriented_elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); 103 104 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass_curl_oriented); 105 CeedOperatorSetField(op_mass_curl_oriented, "rho", elem_restriction_q_data, CEED_BASIS_NONE, q_data); 106 CeedOperatorSetField(op_mass_curl_oriented, "u", curl_oriented_elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); 107 CeedOperatorSetField(op_mass_curl_oriented, "v", curl_oriented_elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE); 108 109 // Apply Setup Operator 110 CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE); 111 112 // Fully assemble operators 113 CeedSize num_entries, num_entries_oriented, num_entries_curl_oriented; 114 CeedInt *rows, *rows_oriented, *rows_curl_oriented; 115 CeedInt *cols, *cols_oriented, *cols_curl_oriented; 116 CeedVector assembled, assembled_oriented, assembled_curl_oriented; 117 118 for (CeedInt k = 0; k < num_dofs * num_dofs; ++k) { 119 assembled_values[k] = 0.0; 120 assembled_values_oriented[k] = 0.0; 121 assembled_values_curl_oriented[k] = 0.0; 122 } 123 CeedOperatorLinearAssembleSymbolic(op_mass, &num_entries, &rows, &cols); 124 CeedOperatorLinearAssembleSymbolic(op_mass_oriented, &num_entries_oriented, &rows_oriented, &cols_oriented); 125 CeedOperatorLinearAssembleSymbolic(op_mass_curl_oriented, &num_entries_curl_oriented, &rows_curl_oriented, &cols_curl_oriented); 126 CeedVectorCreate(ceed, num_entries, &assembled); 127 CeedVectorCreate(ceed, num_entries_oriented, &assembled_oriented); 128 CeedVectorCreate(ceed, num_entries_curl_oriented, &assembled_curl_oriented); 129 CeedOperatorLinearAssemble(op_mass, assembled); 130 CeedOperatorLinearAssemble(op_mass_oriented, assembled_oriented); 131 CeedOperatorLinearAssemble(op_mass_curl_oriented, assembled_curl_oriented); 132 { 133 const CeedScalar *assembled_array; 134 135 CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array); 136 for (CeedInt k = 0; k < num_entries; ++k) { 137 assembled_values[rows[k] * num_dofs + cols[k]] += assembled_array[k]; 138 } 139 CeedVectorRestoreArrayRead(assembled, &assembled_array); 140 } 141 { 142 const CeedScalar *assembled_array; 143 144 CeedVectorGetArrayRead(assembled_oriented, CEED_MEM_HOST, &assembled_array); 145 for (CeedInt k = 0; k < num_entries_oriented; ++k) { 146 assembled_values_oriented[rows_oriented[k] * num_dofs + cols_oriented[k]] += assembled_array[k]; 147 } 148 CeedVectorRestoreArrayRead(assembled_oriented, &assembled_array); 149 } 150 { 151 const CeedScalar *assembled_array; 152 153 CeedVectorGetArrayRead(assembled_curl_oriented, CEED_MEM_HOST, &assembled_array); 154 for (CeedInt k = 0; k < num_entries_curl_oriented; ++k) { 155 assembled_values_curl_oriented[rows_curl_oriented[k] * num_dofs + cols_curl_oriented[k]] += assembled_array[k]; 156 } 157 CeedVectorRestoreArrayRead(assembled_curl_oriented, &assembled_array); 158 } 159 160 // Check output 161 for (CeedInt i = 0; i < num_dofs; i++) { 162 for (CeedInt j = 0; j < num_dofs; j++) { 163 if (fabs(assembled_values_oriented[i * num_dofs + j] - assembled_values[i * num_dofs + j]) > 100. * CEED_EPSILON) { 164 // LCOV_EXCL_START 165 printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] Error in oriented assembly: %f != %f\n", i, j, assembled_values_oriented[i * num_dofs + j], 166 assembled_values[i * num_dofs + j]); 167 // LCOV_EXCL_STOP 168 } 169 if (fabs(assembled_values_curl_oriented[i * num_dofs + j] - assembled_values[i * num_dofs + j]) > 100. * CEED_EPSILON) { 170 // LCOV_EXCL_START 171 printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] Error in curl-oriented assembly: %f != %f\n", i, j, 172 assembled_values_curl_oriented[i * num_dofs + j], assembled_values[i * num_dofs + j]); 173 // LCOV_EXCL_STOP 174 } 175 } 176 } 177 178 // Cleanup 179 free(rows); 180 free(cols); 181 free(rows_oriented); 182 free(cols_oriented); 183 free(rows_curl_oriented); 184 free(cols_curl_oriented); 185 CeedVectorDestroy(&x); 186 CeedVectorDestroy(&q_data); 187 CeedVectorDestroy(&assembled); 188 CeedVectorDestroy(&assembled_oriented); 189 CeedVectorDestroy(&assembled_curl_oriented); 190 CeedElemRestrictionDestroy(&elem_restriction_u); 191 CeedElemRestrictionDestroy(&oriented_elem_restriction_u); 192 CeedElemRestrictionDestroy(&curl_oriented_elem_restriction_u); 193 CeedElemRestrictionDestroy(&elem_restriction_x); 194 CeedElemRestrictionDestroy(&elem_restriction_q_data); 195 CeedBasisDestroy(&basis_u); 196 CeedBasisDestroy(&basis_x); 197 CeedQFunctionDestroy(&qf_setup); 198 CeedQFunctionDestroy(&qf_mass); 199 CeedOperatorDestroy(&op_setup); 200 CeedOperatorDestroy(&op_mass); 201 CeedOperatorDestroy(&op_mass_oriented); 202 CeedOperatorDestroy(&op_mass_curl_oriented); 203 CeedDestroy(&ceed); 204 return 0; 205 } 206