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 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