1 /// @file 2 /// Test creation, action, and destruction for composite mass matrix operator 3 /// \test Test creation, action, and destruction for composite mass matrix operator 4 #include <ceed.h> 5 #include <stdlib.h> 6 #include <math.h> 7 #include "t320-basis.h" 8 #include "t510-operator.h" 9 10 /* The mesh comprises of two rows of 3 quadralaterals followed by one row 11 of 6 triangles: 12 _ _ _ 13 |_|_|_| 14 |_|_|_| 15 |/|/|/| 16 17 */ 18 19 int main(int argc, char **argv) { 20 Ceed ceed; 21 CeedElemRestriction elem_restr_x_tet, elem_restr_u_tet, 22 elem_restr_qd_i_tet, 23 elem_restr_x_hex, elem_restr_u_hex, 24 elem_restr_qd_i_hex; 25 CeedBasis basis_x_tet, basis_u_tet, 26 basis_x_hex, basis_u_hex; 27 CeedQFunction qf_setup_tet, qf_mass_tet, 28 qf_setup_hex, qf_mass_hex; 29 CeedOperator op_setup_tet, op_mass_tet, 30 op_setup_hex, op_mass_hex, 31 op_setup, op_mass; 32 CeedVector q_data_tet, q_data_hex, X, U, V; 33 const CeedScalar *hv; 34 CeedInt num_elem_tet = 6, P_tet = 6, Q_tet = 4, 35 num_elem_hex = 6, P_hex = 3, Q_hex = 4, dim = 2; 36 CeedInt n_x = 3, n_y = 3, 37 n_x_tet = 3, n_y_tet = 1, n_x_hex = 3; 38 CeedInt row, col, offset; 39 CeedInt num_dofs = (n_x*2+1)*(n_y*2+1), 40 num_qpts_tet = num_elem_tet*Q_tet, 41 num_qpts_hex = num_elem_hex*Q_hex*Q_hex; 42 CeedInt ind_x_tet[num_elem_tet*P_tet], 43 ind_x_hex[num_elem_hex*P_hex*P_hex]; 44 CeedScalar x[dim*num_dofs]; 45 CeedScalar q_ref[dim*Q_tet], q_weight[Q_tet]; 46 CeedScalar interp[P_tet*Q_tet], grad[dim*P_tet*Q_tet]; 47 CeedScalar sum; 48 49 CeedInit(argv[1], &ceed); 50 51 // DoF Coordinates 52 for (CeedInt i=0; i<n_y*2+1; i++) 53 for (CeedInt j=0; j<n_x*2+1; j++) { 54 x[i+j*(n_y*2+1)+0*num_dofs] = (CeedScalar) i / (2*n_y); 55 x[i+j*(n_y*2+1)+1*num_dofs] = (CeedScalar) j / (2*n_x); 56 } 57 CeedVectorCreate(ceed, dim*num_dofs, &X); 58 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 59 60 // Qdata Vectors 61 CeedVectorCreate(ceed, num_qpts_tet, &q_data_tet); 62 CeedVectorCreate(ceed, num_qpts_hex, &q_data_hex); 63 64 // Tet Elements 65 for (CeedInt i=0; i<num_elem_tet/2; i++) { 66 col = i % n_x_tet; 67 row = i / n_x_tet; 68 offset = col*2 + row*(n_x_tet*2+1)*2; 69 70 ind_x_tet[i*2*P_tet + 0] = 2 + offset; 71 ind_x_tet[i*2*P_tet + 1] = 9 + offset; 72 ind_x_tet[i*2*P_tet + 2] = 16 + offset; 73 ind_x_tet[i*2*P_tet + 3] = 1 + offset; 74 ind_x_tet[i*2*P_tet + 4] = 8 + offset; 75 ind_x_tet[i*2*P_tet + 5] = 0 + offset; 76 77 ind_x_tet[i*2*P_tet + 6] = 14 + offset; 78 ind_x_tet[i*2*P_tet + 7] = 7 + offset; 79 ind_x_tet[i*2*P_tet + 8] = 0 + offset; 80 ind_x_tet[i*2*P_tet + 9] = 15 + offset; 81 ind_x_tet[i*2*P_tet + 10] = 8 + offset; 82 ind_x_tet[i*2*P_tet + 11] = 16 + offset; 83 } 84 85 // -- Restrictions 86 CeedElemRestrictionCreate(ceed, num_elem_tet, P_tet, dim, num_dofs, 87 dim*num_dofs, 88 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_tet, 89 &elem_restr_x_tet); 90 91 CeedElemRestrictionCreate(ceed, num_elem_tet, P_tet, 1, 1, num_dofs, 92 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_tet, 93 &elem_restr_u_tet); 94 CeedInt strides_qd_tet[3] = {1, Q_tet, Q_tet}; 95 CeedElemRestrictionCreateStrided(ceed, num_elem_tet, Q_tet, 1, num_qpts_tet, 96 strides_qd_tet, &elem_restr_qd_i_tet); 97 98 // -- Bases 99 buildmats(q_ref, q_weight, interp, grad); 100 CeedBasisCreateH1(ceed, CEED_TOPOLOGY_TRIANGLE, dim, P_tet, Q_tet, interp, grad, 101 q_ref, q_weight, &basis_x_tet); 102 103 buildmats(q_ref, q_weight, interp, grad); 104 CeedBasisCreateH1(ceed, CEED_TOPOLOGY_TRIANGLE, 1, P_tet, Q_tet, interp, grad, 105 q_ref, q_weight, &basis_u_tet); 106 107 // -- QFunctions 108 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup_tet); 109 CeedQFunctionAddInput(qf_setup_tet, "weight", 1, CEED_EVAL_WEIGHT); 110 CeedQFunctionAddInput(qf_setup_tet, "dx", dim*dim, CEED_EVAL_GRAD); 111 CeedQFunctionAddOutput(qf_setup_tet, "rho", 1, CEED_EVAL_NONE); 112 113 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass_tet); 114 CeedQFunctionAddInput(qf_mass_tet, "rho", 1, CEED_EVAL_NONE); 115 CeedQFunctionAddInput(qf_mass_tet, "u", 1, CEED_EVAL_INTERP); 116 CeedQFunctionAddOutput(qf_mass_tet, "v", 1, CEED_EVAL_INTERP); 117 118 // -- Operators 119 // ---- Setup Tet 120 CeedOperatorCreate(ceed, qf_setup_tet, CEED_QFUNCTION_NONE, 121 CEED_QFUNCTION_NONE, &op_setup_tet); 122 CeedOperatorSetField(op_setup_tet, "weight", CEED_ELEMRESTRICTION_NONE, 123 basis_x_tet, 124 CEED_VECTOR_NONE); 125 CeedOperatorSetField(op_setup_tet, "dx", elem_restr_x_tet, basis_x_tet, 126 CEED_VECTOR_ACTIVE); 127 CeedOperatorSetField(op_setup_tet, "rho", elem_restr_qd_i_tet, 128 CEED_BASIS_COLLOCATED, q_data_tet); 129 // ---- Mass Tet 130 CeedOperatorCreate(ceed, qf_mass_tet, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 131 &op_mass_tet); 132 CeedOperatorSetField(op_mass_tet, "rho", elem_restr_qd_i_tet, 133 CEED_BASIS_COLLOCATED, 134 q_data_tet); 135 CeedOperatorSetField(op_mass_tet, "u", elem_restr_u_tet, basis_u_tet, 136 CEED_VECTOR_ACTIVE); 137 CeedOperatorSetField(op_mass_tet, "v", elem_restr_u_tet, basis_u_tet, 138 CEED_VECTOR_ACTIVE); 139 140 // Hex Elements 141 for (CeedInt i=0; i<num_elem_hex; i++) { 142 col = i % n_x_hex; 143 row = i / n_x_hex; 144 offset = (n_x_tet*2+1)*(n_y_tet*2)*(1+row) + col*2; 145 for (CeedInt j=0; j<P_hex; j++) 146 for (CeedInt k=0; k<P_hex; k++) 147 ind_x_hex[P_hex*(P_hex*i+k)+j] = offset + k*(n_x_hex*2+1) + j; 148 } 149 150 // -- Restrictions 151 CeedElemRestrictionCreate(ceed, num_elem_hex, P_hex*P_hex, dim, num_dofs, 152 dim*num_dofs, 153 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_hex, 154 &elem_restr_x_hex); 155 156 CeedElemRestrictionCreate(ceed, num_elem_hex, P_hex*P_hex, 1, 1, num_dofs, 157 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_hex, 158 &elem_restr_u_hex); 159 CeedInt strides_qd_hex[3] = {1, Q_hex*Q_hex, Q_hex*Q_hex}; 160 CeedElemRestrictionCreateStrided(ceed, num_elem_hex, Q_hex*Q_hex, 1, 161 num_qpts_hex, 162 strides_qd_hex, &elem_restr_qd_i_hex); 163 164 // -- Bases 165 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, P_hex, Q_hex, CEED_GAUSS, 166 &basis_x_hex); 167 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P_hex, Q_hex, CEED_GAUSS, 168 &basis_u_hex); 169 170 // -- QFunctions 171 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup_hex); 172 CeedQFunctionAddInput(qf_setup_hex, "_weight", 1, CEED_EVAL_WEIGHT); 173 CeedQFunctionAddInput(qf_setup_hex, "dx", dim*dim, CEED_EVAL_GRAD); 174 CeedQFunctionAddOutput(qf_setup_hex, "rho", 1, CEED_EVAL_NONE); 175 176 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass_hex); 177 CeedQFunctionAddInput(qf_mass_hex, "rho", 1, CEED_EVAL_NONE); 178 CeedQFunctionAddInput(qf_mass_hex, "u", 1, CEED_EVAL_INTERP); 179 CeedQFunctionAddOutput(qf_mass_hex, "v", 1, CEED_EVAL_INTERP); 180 181 // -- Operators 182 CeedOperatorCreate(ceed, qf_setup_hex, CEED_QFUNCTION_NONE, 183 CEED_QFUNCTION_NONE, &op_setup_hex); 184 CeedOperatorSetField(op_setup_hex, "_weight", CEED_ELEMRESTRICTION_NONE, 185 basis_x_hex, 186 CEED_VECTOR_NONE); 187 CeedOperatorSetField(op_setup_hex, "dx", elem_restr_x_hex, basis_x_hex, 188 CEED_VECTOR_ACTIVE); 189 CeedOperatorSetField(op_setup_hex, "rho", elem_restr_qd_i_hex, 190 CEED_BASIS_COLLOCATED, q_data_hex); 191 192 CeedOperatorCreate(ceed, qf_mass_hex, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 193 &op_mass_hex); 194 CeedOperatorSetField(op_mass_hex, "rho", elem_restr_qd_i_hex, 195 CEED_BASIS_COLLOCATED, 196 q_data_hex); 197 CeedOperatorSetField(op_mass_hex, "u", elem_restr_u_hex, basis_u_hex, 198 CEED_VECTOR_ACTIVE); 199 CeedOperatorSetField(op_mass_hex, "v", elem_restr_u_hex, basis_u_hex, 200 CEED_VECTOR_ACTIVE); 201 202 // Composite Operators 203 CeedCompositeOperatorCreate(ceed, &op_setup); 204 CeedCompositeOperatorAddSub(op_setup, op_setup_tet); 205 CeedCompositeOperatorAddSub(op_setup, op_setup_hex); 206 207 CeedCompositeOperatorCreate(ceed, &op_mass); 208 CeedCompositeOperatorAddSub(op_mass, op_mass_tet); 209 CeedCompositeOperatorAddSub(op_mass, op_mass_hex); 210 211 // Apply Setup Operator 212 CeedOperatorApply(op_setup, X, CEED_VECTOR_NONE, CEED_REQUEST_IMMEDIATE); 213 214 // Apply Mass Operator 215 CeedVectorCreate(ceed, num_dofs, &U); 216 CeedVectorSetValue(U, 1.0); 217 CeedVectorCreate(ceed, num_dofs, &V); 218 219 CeedOperatorApply(op_mass, U, V, CEED_REQUEST_IMMEDIATE); 220 221 // Check output 222 CeedVectorGetArrayRead(V, CEED_MEM_HOST, &hv); 223 sum = 0.; 224 for (CeedInt i=0; i<num_dofs; i++) 225 sum += hv[i]; 226 if (fabs(sum-1.)>1000.*CEED_EPSILON) 227 // LCOV_EXCL_START 228 printf("Computed Area: %f != True Area: 1.0\n", sum); 229 // LCOV_EXCL_STOP 230 CeedVectorRestoreArrayRead(V, &hv); 231 232 // Cleanup 233 CeedQFunctionDestroy(&qf_setup_tet); 234 CeedQFunctionDestroy(&qf_mass_tet); 235 CeedOperatorDestroy(&op_setup_tet); 236 CeedOperatorDestroy(&op_mass_tet); 237 CeedQFunctionDestroy(&qf_setup_hex); 238 CeedQFunctionDestroy(&qf_mass_hex); 239 CeedOperatorDestroy(&op_setup_hex); 240 CeedOperatorDestroy(&op_mass_hex); 241 CeedOperatorDestroy(&op_setup); 242 CeedOperatorDestroy(&op_mass); 243 CeedElemRestrictionDestroy(&elem_restr_u_tet); 244 CeedElemRestrictionDestroy(&elem_restr_x_tet); 245 CeedElemRestrictionDestroy(&elem_restr_qd_i_tet); 246 CeedElemRestrictionDestroy(&elem_restr_u_hex); 247 CeedElemRestrictionDestroy(&elem_restr_x_hex); 248 CeedElemRestrictionDestroy(&elem_restr_qd_i_hex); 249 CeedBasisDestroy(&basis_u_tet); 250 CeedBasisDestroy(&basis_x_tet); 251 CeedBasisDestroy(&basis_u_hex); 252 CeedBasisDestroy(&basis_x_hex); 253 CeedVectorDestroy(&X); 254 CeedVectorDestroy(&U); 255 CeedVectorDestroy(&V); 256 CeedVectorDestroy(&q_data_tet); 257 CeedVectorDestroy(&q_data_hex); 258 CeedDestroy(&ceed); 259 return 0; 260 } 261