1 /// @file 2 /// Test creation, action, and destruction for mass matrix operator 3 /// \test Test creation, action, and destruction for 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_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; 33 CeedInt nelem_tet = 6, P_tet = 6, Q_tet = 4, 34 nelem_hex = 6, P_hex = 3, Q_hex = 4, dim = 2; 35 CeedInt nx = 3, ny = 3, 36 nx_tet = 3, ny_tet = 1, nx_hex = 3; 37 CeedInt row, col, offset; 38 CeedInt num_dofs = (nx*2+1)*(ny*2+1), 39 num_qpts_tet = nelem_tet*Q_tet, 40 num_qpts_hex = nelem_hex*Q_hex*Q_hex; 41 CeedInt ind_x_tet[nelem_tet*P_tet], 42 ind_x_hex[nelem_hex*P_hex*P_hex]; 43 CeedScalar q_ref[dim*Q_tet], q_weight[Q_tet]; 44 CeedScalar interp[P_tet*Q_tet], grad[dim*P_tet*Q_tet]; 45 46 CeedInit(argv[1], &ceed); 47 48 // DoF Coordinates 49 CeedVectorCreate(ceed, dim*num_dofs, &X); 50 51 // Qdata Vectors 52 CeedVectorCreate(ceed, num_qpts_tet, &q_data_tet); 53 CeedVectorCreate(ceed, num_qpts_hex, &q_data_hex); 54 55 // Set up _tet Elements 56 for (CeedInt i=0; i<nelem_tet/2; i++) { 57 col = i % nx_tet; 58 row = i / nx_tet; 59 offset = col*2 + row*(nx_tet*2+1)*2; 60 61 ind_x_tet[i*2*P_tet + 0] = 2 + offset; 62 ind_x_tet[i*2*P_tet + 1] = 9 + offset; 63 ind_x_tet[i*2*P_tet + 2] = 16 + offset; 64 ind_x_tet[i*2*P_tet + 3] = 1 + offset; 65 ind_x_tet[i*2*P_tet + 4] = 8 + offset; 66 ind_x_tet[i*2*P_tet + 5] = 0 + offset; 67 68 ind_x_tet[i*2*P_tet + 6] = 14 + offset; 69 ind_x_tet[i*2*P_tet + 7] = 7 + offset; 70 ind_x_tet[i*2*P_tet + 8] = 0 + offset; 71 ind_x_tet[i*2*P_tet + 9] = 15 + offset; 72 ind_x_tet[i*2*P_tet + 10] = 8 + offset; 73 ind_x_tet[i*2*P_tet + 11] = 16 + offset; 74 } 75 76 // -- Restrictions 77 CeedElemRestrictionCreate(ceed, nelem_tet, P_tet, dim, num_dofs, dim*num_dofs, 78 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_tet, 79 &elem_restr_x_tet); 80 81 CeedElemRestrictionCreate(ceed, nelem_tet, P_tet, 1, 1, num_dofs, 82 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_tet, 83 &elem_restr_u_tet); 84 CeedInt strides_qd_tet[3] = {1, Q_tet, Q_tet}; 85 CeedElemRestrictionCreateStrided(ceed, nelem_tet, Q_tet, 1, num_qpts_tet, 86 strides_qd_tet, &elem_restr_qd_tet); 87 88 // -- Bases 89 buildmats(q_ref, q_weight, interp, grad); 90 CeedBasisCreateH1(ceed, CEED_TRIANGLE, dim, P_tet, Q_tet, interp, grad, q_ref, 91 q_weight, &basis_x_tet); 92 93 buildmats(q_ref, q_weight, interp, grad); 94 CeedBasisCreateH1(ceed, CEED_TRIANGLE, 1, P_tet, Q_tet, interp, grad, q_ref, 95 q_weight, &basis_u_tet); 96 97 // -- QFunctions 98 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup_tet); 99 CeedQFunctionAddInput(qf_setup_tet, "_weight", 1, CEED_EVAL_WEIGHT); 100 CeedQFunctionAddInput(qf_setup_tet, "dx", dim*dim, CEED_EVAL_GRAD); 101 CeedQFunctionAddOutput(qf_setup_tet, "rho", 1, CEED_EVAL_NONE); 102 103 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass_tet); 104 CeedQFunctionAddInput(qf_mass_tet, "rho", 1, CEED_EVAL_NONE); 105 CeedQFunctionAddInput(qf_mass_tet, "u", 1, CEED_EVAL_INTERP); 106 CeedQFunctionAddOutput(qf_mass_tet, "v", 1, CEED_EVAL_INTERP); 107 108 // -- Operators 109 // ---- Setup _tet 110 CeedOperatorCreate(ceed, qf_setup_tet, CEED_QFUNCTION_NONE, 111 CEED_QFUNCTION_NONE, &op_setup_tet); 112 CeedOperatorSetField(op_setup_tet, "_weight", CEED_ELEMRESTRICTION_NONE, 113 basis_x_tet, 114 CEED_VECTOR_NONE); 115 CeedOperatorSetField(op_setup_tet, "dx", elem_restr_x_tet, basis_x_tet, 116 CEED_VECTOR_ACTIVE); 117 CeedOperatorSetField(op_setup_tet, "rho", elem_restr_qd_tet, 118 CEED_BASIS_COLLOCATED, q_data_tet); 119 // ---- Mass _tet 120 CeedOperatorCreate(ceed, qf_mass_tet, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 121 &op_mass_tet); 122 CeedOperatorSetField(op_mass_tet, "rho", elem_restr_qd_tet, 123 CEED_BASIS_COLLOCATED, 124 q_data_tet); 125 CeedOperatorSetField(op_mass_tet, "u", elem_restr_u_tet, basis_u_tet, 126 CEED_VECTOR_ACTIVE); 127 CeedOperatorSetField(op_mass_tet, "v", elem_restr_u_tet, basis_u_tet, 128 CEED_VECTOR_ACTIVE); 129 130 // Set up _hex Elements 131 for (CeedInt i=0; i<nelem_hex; i++) { 132 col = i % nx_hex; 133 row = i / nx_hex; 134 offset = (nx_tet*2+1)*(ny_tet*2)*(1+row) + col*2; 135 for (CeedInt j=0; j<P_hex; j++) 136 for (CeedInt k=0; k<P_hex; k++) 137 ind_x_hex[P_hex*(P_hex*i+k)+j] = offset + k*(nx_hex*2+1) + j; 138 } 139 140 // -- Restrictions 141 CeedElemRestrictionCreate(ceed, nelem_hex, P_hex*P_hex, dim, num_dofs, 142 dim*num_dofs, 143 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_hex, 144 &elem_restr_x_hex); 145 146 CeedElemRestrictionCreate(ceed, nelem_hex, P_hex*P_hex, 1, 1, num_dofs, 147 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_hex, 148 &elem_restr_u_hex); 149 CeedInt strides_qd_hex[3] = {1, Q_hex*Q_hex, Q_hex*Q_hex}; 150 CeedElemRestrictionCreateStrided(ceed, nelem_hex, Q_hex*Q_hex, 1, num_qpts_hex, 151 strides_qd_hex, &elem_restr_qd_i_hex); 152 153 // -- Bases 154 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, P_hex, Q_hex, CEED_GAUSS, 155 &basis_x_hex); 156 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P_hex, Q_hex, CEED_GAUSS, 157 &basis_u_hex); 158 159 // -- QFunctions 160 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup_hex); 161 CeedQFunctionAddInput(qf_setup_hex, "_weight", 1, CEED_EVAL_WEIGHT); 162 CeedQFunctionAddInput(qf_setup_hex, "dx", dim*dim, CEED_EVAL_GRAD); 163 CeedQFunctionAddOutput(qf_setup_hex, "rho", 1, CEED_EVAL_NONE); 164 165 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass_hex); 166 CeedQFunctionAddInput(qf_mass_hex, "rho", 1, CEED_EVAL_NONE); 167 CeedQFunctionAddInput(qf_mass_hex, "u", 1, CEED_EVAL_INTERP); 168 CeedQFunctionAddOutput(qf_mass_hex, "v", 1, CEED_EVAL_INTERP); 169 170 // -- Operators 171 CeedOperatorCreate(ceed, qf_setup_hex, CEED_QFUNCTION_NONE, 172 CEED_QFUNCTION_NONE, &op_setup_hex); 173 CeedOperatorSetField(op_setup_hex, "_weight", CEED_ELEMRESTRICTION_NONE, 174 basis_x_hex, 175 CEED_VECTOR_NONE); 176 CeedOperatorSetField(op_setup_hex, "dx", elem_restr_x_hex, basis_x_hex, 177 CEED_VECTOR_ACTIVE); 178 CeedOperatorSetField(op_setup_hex, "rho", elem_restr_qd_i_hex, 179 CEED_BASIS_COLLOCATED, q_data_hex); 180 181 CeedOperatorCreate(ceed, qf_mass_hex, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 182 &op_mass_hex); 183 CeedOperatorSetField(op_mass_hex, "rho", elem_restr_qd_i_hex, 184 CEED_BASIS_COLLOCATED, 185 q_data_hex); 186 CeedOperatorSetField(op_mass_hex, "u", elem_restr_u_hex, basis_u_hex, 187 CEED_VECTOR_ACTIVE); 188 CeedOperatorSetField(op_mass_hex, "v", elem_restr_u_hex, basis_u_hex, 189 CEED_VECTOR_ACTIVE); 190 191 // Set up Composite Operators 192 // -- Create 193 CeedCompositeOperatorCreate(ceed, &op_setup); 194 // -- Add SubOperators 195 CeedCompositeOperatorAddSub(op_setup, op_setup_tet); 196 CeedCompositeOperatorAddSub(op_setup, op_setup_hex); 197 198 // -- Create 199 CeedCompositeOperatorCreate(ceed, &op_mass); 200 // -- Add SubOperators 201 CeedCompositeOperatorAddSub(op_mass, op_mass_tet); 202 CeedCompositeOperatorAddSub(op_mass, op_mass_hex); 203 204 // View 205 CeedOperatorView(op_setup, stdout); 206 CeedOperatorView(op_mass, stdout); 207 208 // Cleanup 209 CeedQFunctionDestroy(&qf_setup_tet); 210 CeedQFunctionDestroy(&qf_mass_tet); 211 CeedOperatorDestroy(&op_setup_tet); 212 CeedOperatorDestroy(&op_mass_tet); 213 CeedQFunctionDestroy(&qf_setup_hex); 214 CeedQFunctionDestroy(&qf_mass_hex); 215 CeedOperatorDestroy(&op_setup_hex); 216 CeedOperatorDestroy(&op_mass_hex); 217 CeedOperatorDestroy(&op_setup); 218 CeedOperatorDestroy(&op_mass); 219 CeedElemRestrictionDestroy(&elem_restr_u_tet); 220 CeedElemRestrictionDestroy(&elem_restr_x_tet); 221 CeedElemRestrictionDestroy(&elem_restr_qd_tet); 222 CeedElemRestrictionDestroy(&elem_restr_u_hex); 223 CeedElemRestrictionDestroy(&elem_restr_x_hex); 224 CeedElemRestrictionDestroy(&elem_restr_qd_i_hex); 225 CeedBasisDestroy(&basis_u_tet); 226 CeedBasisDestroy(&basis_x_tet); 227 CeedBasisDestroy(&basis_u_hex); 228 CeedBasisDestroy(&basis_x_hex); 229 CeedVectorDestroy(&X); 230 CeedVectorDestroy(&q_data_tet); 231 CeedVectorDestroy(&q_data_hex); 232 CeedDestroy(&ceed); 233 return 0; 234 } 235