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