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