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