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_TOPOLOGY_TRIANGLE, dim, P_tet, Q_tet, interp, grad, 91 q_ref, q_weight, &basis_x_tet); 92 93 buildmats(q_ref, q_weight, interp, grad); 94 CeedBasisCreateH1(ceed, CEED_TOPOLOGY_TRIANGLE, 1, P_tet, Q_tet, interp, grad, 95 q_ref, 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 CeedOperatorSetName(op_setup_tet, "triangle elements"); 113 CeedOperatorSetField(op_setup_tet, "weight", CEED_ELEMRESTRICTION_NONE, 114 basis_x_tet, 115 CEED_VECTOR_NONE); 116 CeedOperatorSetField(op_setup_tet, "dx", elem_restr_x_tet, basis_x_tet, 117 CEED_VECTOR_ACTIVE); 118 CeedOperatorSetField(op_setup_tet, "rho", elem_restr_qd_tet, 119 CEED_BASIS_COLLOCATED, q_data_tet); 120 // ---- Mass _tet 121 CeedOperatorCreate(ceed, qf_mass_tet, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 122 &op_mass_tet); 123 CeedOperatorSetName(op_mass_tet, "triangle elements"); 124 CeedOperatorSetField(op_mass_tet, "rho", elem_restr_qd_tet, 125 CEED_BASIS_COLLOCATED, 126 q_data_tet); 127 CeedOperatorSetField(op_mass_tet, "u", elem_restr_u_tet, basis_u_tet, 128 CEED_VECTOR_ACTIVE); 129 CeedOperatorSetField(op_mass_tet, "v", elem_restr_u_tet, basis_u_tet, 130 CEED_VECTOR_ACTIVE); 131 132 // Set up _hex Elements 133 for (CeedInt i=0; i<nelem_hex; i++) { 134 col = i % nx_hex; 135 row = i / nx_hex; 136 offset = (nx_tet*2+1)*(ny_tet*2)*(1+row) + col*2; 137 for (CeedInt j=0; j<P_hex; j++) 138 for (CeedInt k=0; k<P_hex; k++) 139 ind_x_hex[P_hex*(P_hex*i+k)+j] = offset + k*(nx_hex*2+1) + j; 140 } 141 142 // -- Restrictions 143 CeedElemRestrictionCreate(ceed, nelem_hex, P_hex*P_hex, dim, num_dofs, 144 dim*num_dofs, 145 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_hex, 146 &elem_restr_x_hex); 147 148 CeedElemRestrictionCreate(ceed, nelem_hex, P_hex*P_hex, 1, 1, num_dofs, 149 CEED_MEM_HOST, CEED_USE_POINTER, ind_x_hex, 150 &elem_restr_u_hex); 151 CeedInt strides_qd_hex[3] = {1, Q_hex*Q_hex, Q_hex*Q_hex}; 152 CeedElemRestrictionCreateStrided(ceed, nelem_hex, Q_hex*Q_hex, 1, num_qpts_hex, 153 strides_qd_hex, &elem_restr_qd_i_hex); 154 155 // -- Bases 156 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, P_hex, Q_hex, CEED_GAUSS, 157 &basis_x_hex); 158 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P_hex, Q_hex, CEED_GAUSS, 159 &basis_u_hex); 160 161 // -- QFunctions 162 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup_hex); 163 CeedQFunctionAddInput(qf_setup_hex, "weight", 1, CEED_EVAL_WEIGHT); 164 CeedQFunctionAddInput(qf_setup_hex, "dx", dim*dim, CEED_EVAL_GRAD); 165 CeedQFunctionAddOutput(qf_setup_hex, "rho", 1, CEED_EVAL_NONE); 166 167 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass_hex); 168 CeedQFunctionAddInput(qf_mass_hex, "rho", 1, CEED_EVAL_NONE); 169 CeedQFunctionAddInput(qf_mass_hex, "u", 1, CEED_EVAL_INTERP); 170 CeedQFunctionAddOutput(qf_mass_hex, "v", 1, CEED_EVAL_INTERP); 171 172 // -- Operators 173 CeedOperatorCreate(ceed, qf_setup_hex, CEED_QFUNCTION_NONE, 174 CEED_QFUNCTION_NONE, &op_setup_hex); 175 CeedOperatorSetName(op_setup_hex, "quadralateral elements"); 176 CeedOperatorSetField(op_setup_hex, "weight", CEED_ELEMRESTRICTION_NONE, 177 basis_x_hex, 178 CEED_VECTOR_NONE); 179 CeedOperatorSetField(op_setup_hex, "dx", elem_restr_x_hex, basis_x_hex, 180 CEED_VECTOR_ACTIVE); 181 CeedOperatorSetField(op_setup_hex, "rho", elem_restr_qd_i_hex, 182 CEED_BASIS_COLLOCATED, q_data_hex); 183 184 CeedOperatorCreate(ceed, qf_mass_hex, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 185 &op_mass_hex); 186 CeedOperatorSetName(op_mass_hex, "quadralateral elements"); 187 CeedOperatorSetField(op_mass_hex, "rho", elem_restr_qd_i_hex, 188 CEED_BASIS_COLLOCATED, 189 q_data_hex); 190 CeedOperatorSetField(op_mass_hex, "u", elem_restr_u_hex, basis_u_hex, 191 CEED_VECTOR_ACTIVE); 192 CeedOperatorSetField(op_mass_hex, "v", elem_restr_u_hex, basis_u_hex, 193 CEED_VECTOR_ACTIVE); 194 195 // Set up Composite Operators 196 // -- Create 197 CeedCompositeOperatorCreate(ceed, &op_setup); 198 CeedOperatorSetName(op_setup, "setup"); 199 // -- Add SubOperators 200 CeedCompositeOperatorAddSub(op_setup, op_setup_tet); 201 CeedCompositeOperatorAddSub(op_setup, op_setup_hex); 202 203 // -- Create 204 CeedCompositeOperatorCreate(ceed, &op_mass); 205 CeedOperatorSetName(op_mass, "mass"); 206 // -- Add SubOperators 207 CeedCompositeOperatorAddSub(op_mass, op_mass_tet); 208 CeedCompositeOperatorAddSub(op_mass, op_mass_hex); 209 210 // View 211 CeedOperatorView(op_setup, stdout); 212 CeedOperatorView(op_mass, stdout); 213 214 // Cleanup 215 CeedQFunctionDestroy(&qf_setup_tet); 216 CeedQFunctionDestroy(&qf_mass_tet); 217 CeedOperatorDestroy(&op_setup_tet); 218 CeedOperatorDestroy(&op_mass_tet); 219 CeedQFunctionDestroy(&qf_setup_hex); 220 CeedQFunctionDestroy(&qf_mass_hex); 221 CeedOperatorDestroy(&op_setup_hex); 222 CeedOperatorDestroy(&op_mass_hex); 223 CeedOperatorDestroy(&op_setup); 224 CeedOperatorDestroy(&op_mass); 225 CeedElemRestrictionDestroy(&elem_restr_u_tet); 226 CeedElemRestrictionDestroy(&elem_restr_x_tet); 227 CeedElemRestrictionDestroy(&elem_restr_qd_tet); 228 CeedElemRestrictionDestroy(&elem_restr_u_hex); 229 CeedElemRestrictionDestroy(&elem_restr_x_hex); 230 CeedElemRestrictionDestroy(&elem_restr_qd_i_hex); 231 CeedBasisDestroy(&basis_u_tet); 232 CeedBasisDestroy(&basis_x_tet); 233 CeedBasisDestroy(&basis_u_hex); 234 CeedBasisDestroy(&basis_x_hex); 235 CeedVectorDestroy(&X); 236 CeedVectorDestroy(&q_data_tet); 237 CeedVectorDestroy(&q_data_hex); 238 CeedDestroy(&ceed); 239 return 0; 240 } 241