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 CeedTransposeMode lmode = CEED_NOTRANSPOSE; 22 CeedElemRestriction ErestrictxTet, ErestrictuTet, 23 ErestrictxiTet, ErestrictuiTet, 24 ErestrictxHex, ErestrictuHex, 25 ErestrictxiHex, ErestrictuiHex; 26 CeedBasis bxTet, buTet, 27 bxHex, buHex; 28 CeedQFunction qf_setupTet, qf_massTet, 29 qf_setupHex, qf_massHex; 30 CeedOperator op_setupTet, op_massTet, 31 op_setupHex, op_massHex, 32 op_setup, op_mass; 33 CeedVector qdataTet, qdataHex, X; 34 CeedInt nelemTet = 6, PTet = 6, QTet = 4, 35 nelemHex = 6, PHex = 3, QHex = 4, dim = 2; 36 CeedInt nx = 3, ny = 3, 37 nxTet = 3, nyTet = 1, nxHex = 3; 38 CeedInt row, col, offset; 39 CeedInt ndofs = (nx*2+1)*(ny*2+1), 40 nqptsTet = nelemTet*QTet, 41 nqptsHex = nelemHex*QHex*QHex; 42 CeedInt indxTet[nelemTet*PTet], 43 indxHex[nelemHex*PHex*PHex]; 44 CeedScalar qref[dim*QTet], qweight[QTet]; 45 CeedScalar interp[PTet*QTet], grad[dim*PTet*QTet]; 46 47 CeedInit(argv[1], &ceed); 48 49 // DoF Coordinates 50 CeedVectorCreate(ceed, dim*ndofs, &X); 51 52 // Qdata Vectors 53 CeedVectorCreate(ceed, nqptsTet, &qdataTet); 54 CeedVectorCreate(ceed, nqptsHex, &qdataHex); 55 56 // Set up Tet Elements 57 for (CeedInt i=0; i<nelemTet/2; i++) { 58 col = i % nxTet; 59 row = i / nxTet; 60 offset = col*2 + row*(nxTet*2+1)*2; 61 62 indxTet[i*2*PTet + 0] = 2 + offset; 63 indxTet[i*2*PTet + 1] = 9 + offset; 64 indxTet[i*2*PTet + 2] = 16 + offset; 65 indxTet[i*2*PTet + 3] = 1 + offset; 66 indxTet[i*2*PTet + 4] = 8 + offset; 67 indxTet[i*2*PTet + 5] = 0 + offset; 68 69 indxTet[i*2*PTet + 6] = 14 + offset; 70 indxTet[i*2*PTet + 7] = 7 + offset; 71 indxTet[i*2*PTet + 8] = 0 + offset; 72 indxTet[i*2*PTet + 9] = 15 + offset; 73 indxTet[i*2*PTet + 10] = 8 + offset; 74 indxTet[i*2*PTet + 11] = 16 + offset; 75 } 76 77 // -- Restrictions 78 CeedElemRestrictionCreate(ceed, lmode, nelemTet, PTet, ndofs, dim, 79 CEED_MEM_HOST, CEED_USE_POINTER, indxTet, 80 &ErestrictxTet); 81 CeedElemRestrictionCreateIdentity(ceed, lmode, nelemTet, PTet, nelemTet*PTet, 82 dim, &ErestrictxiTet); 83 84 CeedElemRestrictionCreate(ceed, lmode, nelemTet, PTet, ndofs, 1, 85 CEED_MEM_HOST, CEED_USE_POINTER, indxTet, 86 &ErestrictuTet); 87 CeedElemRestrictionCreateIdentity(ceed, lmode, nelemTet, QTet, nqptsTet, 1, 88 &ErestrictuiTet); 89 90 // -- Bases 91 buildmats(qref, qweight, interp, grad); 92 CeedBasisCreateH1(ceed, CEED_TRIANGLE, dim, PTet, QTet, interp, grad, qref, 93 qweight, &bxTet); 94 95 buildmats(qref, qweight, interp, grad); 96 CeedBasisCreateH1(ceed, CEED_TRIANGLE, 1, PTet, QTet, interp, grad, qref, 97 qweight, &buTet); 98 99 // -- QFunctions 100 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setupTet); 101 CeedQFunctionAddInput(qf_setupTet, "_weight", 1, CEED_EVAL_WEIGHT); 102 CeedQFunctionAddInput(qf_setupTet, "dx", dim*dim, CEED_EVAL_GRAD); 103 CeedQFunctionAddOutput(qf_setupTet, "rho", 1, CEED_EVAL_NONE); 104 105 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_massTet); 106 CeedQFunctionAddInput(qf_massTet, "rho", 1, CEED_EVAL_NONE); 107 CeedQFunctionAddInput(qf_massTet, "u", 1, CEED_EVAL_INTERP); 108 CeedQFunctionAddOutput(qf_massTet, "v", 1, CEED_EVAL_INTERP); 109 110 // -- Operators 111 // ---- Setup Tet 112 CeedOperatorCreate(ceed, qf_setupTet, CEED_QFUNCTION_NONE, 113 CEED_QFUNCTION_NONE, &op_setupTet); 114 CeedOperatorSetField(op_setupTet, "_weight", ErestrictxiTet, bxTet, 115 CEED_VECTOR_NONE); 116 CeedOperatorSetField(op_setupTet, "dx", ErestrictxTet, bxTet, 117 CEED_VECTOR_ACTIVE); 118 CeedOperatorSetField(op_setupTet, "rho", ErestrictuiTet, 119 CEED_BASIS_COLLOCATED, qdataTet); 120 // ---- Mass Tet 121 CeedOperatorCreate(ceed, qf_massTet, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 122 &op_massTet); 123 CeedOperatorSetField(op_massTet, "rho", ErestrictuiTet, CEED_BASIS_COLLOCATED, 124 qdataTet); 125 CeedOperatorSetField(op_massTet, "u", ErestrictuTet, buTet, 126 CEED_VECTOR_ACTIVE); 127 CeedOperatorSetField(op_massTet, "v", ErestrictuTet, buTet, 128 CEED_VECTOR_ACTIVE); 129 130 // Set up Hex Elements 131 for (CeedInt i=0; i<nelemHex; i++) { 132 col = i % nxHex; 133 row = i / nxHex; 134 offset = (nxTet*2+1)*(nyTet*2)*(1+row) + col*2; 135 for (CeedInt j=0; j<PHex; j++) 136 for (CeedInt k=0; k<PHex; k++) 137 indxHex[PHex*(PHex*i+k)+j] = offset + k*(nxHex*2+1) + j; 138 } 139 140 // -- Restrictions 141 CeedElemRestrictionCreate(ceed, lmode, nelemHex, PHex*PHex, ndofs, dim, 142 CEED_MEM_HOST, CEED_USE_POINTER, indxHex, 143 &ErestrictxHex); 144 CeedElemRestrictionCreateIdentity(ceed, lmode, nelemHex, PHex*PHex, 145 nelemHex*PHex*PHex, dim, &ErestrictxiHex); 146 147 CeedElemRestrictionCreate(ceed, lmode, nelemHex, PHex*PHex, ndofs, 1, 148 CEED_MEM_HOST, CEED_USE_POINTER, indxHex, 149 &ErestrictuHex); 150 CeedElemRestrictionCreateIdentity(ceed, lmode, nelemHex, QHex*QHex, nqptsHex, 151 1, &ErestrictuiHex); 152 153 // -- Bases 154 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, PHex, QHex, CEED_GAUSS, 155 &bxHex); 156 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, PHex, QHex, CEED_GAUSS, &buHex); 157 158 // -- QFunctions 159 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setupHex); 160 CeedQFunctionAddInput(qf_setupHex, "_weight", 1, CEED_EVAL_WEIGHT); 161 CeedQFunctionAddInput(qf_setupHex, "dx", dim*dim, CEED_EVAL_GRAD); 162 CeedQFunctionAddOutput(qf_setupHex, "rho", 1, CEED_EVAL_NONE); 163 164 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_massHex); 165 CeedQFunctionAddInput(qf_massHex, "rho", 1, CEED_EVAL_NONE); 166 CeedQFunctionAddInput(qf_massHex, "u", 1, CEED_EVAL_INTERP); 167 CeedQFunctionAddOutput(qf_massHex, "v", 1, CEED_EVAL_INTERP); 168 169 // -- Operators 170 CeedOperatorCreate(ceed, qf_setupHex, CEED_QFUNCTION_NONE, 171 CEED_QFUNCTION_NONE, &op_setupHex); 172 CeedOperatorSetField(op_setupHex, "_weight", ErestrictxiHex, bxHex, 173 CEED_VECTOR_NONE); 174 CeedOperatorSetField(op_setupHex, "dx", ErestrictxHex, bxHex, 175 CEED_VECTOR_ACTIVE); 176 CeedOperatorSetField(op_setupHex, "rho", ErestrictuiHex, 177 CEED_BASIS_COLLOCATED, qdataHex); 178 179 CeedOperatorCreate(ceed, qf_massHex, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 180 &op_massHex); 181 CeedOperatorSetField(op_massHex, "rho", ErestrictuiHex, CEED_BASIS_COLLOCATED, 182 qdataHex); 183 CeedOperatorSetField(op_massHex, "u", ErestrictuHex, buHex, 184 CEED_VECTOR_ACTIVE); 185 CeedOperatorSetField(op_massHex, "v", ErestrictuHex, buHex, 186 CEED_VECTOR_ACTIVE); 187 188 // Set up Composite Operators 189 // -- Create 190 CeedCompositeOperatorCreate(ceed, &op_setup); 191 // -- Add SubOperators 192 CeedCompositeOperatorAddSub(op_setup, op_setupTet); 193 CeedCompositeOperatorAddSub(op_setup, op_setupHex); 194 195 // -- Create 196 CeedCompositeOperatorCreate(ceed, &op_mass); 197 // -- Add SubOperators 198 CeedCompositeOperatorAddSub(op_mass, op_massTet); 199 CeedCompositeOperatorAddSub(op_mass, op_massHex); 200 201 // View 202 CeedOperatorView(op_setup, stdout); 203 CeedOperatorView(op_mass, stdout); 204 205 // Cleanup 206 CeedQFunctionDestroy(&qf_setupTet); 207 CeedQFunctionDestroy(&qf_massTet); 208 CeedOperatorDestroy(&op_setupTet); 209 CeedOperatorDestroy(&op_massTet); 210 CeedQFunctionDestroy(&qf_setupHex); 211 CeedQFunctionDestroy(&qf_massHex); 212 CeedOperatorDestroy(&op_setupHex); 213 CeedOperatorDestroy(&op_massHex); 214 CeedOperatorDestroy(&op_setup); 215 CeedOperatorDestroy(&op_mass); 216 CeedElemRestrictionDestroy(&ErestrictuTet); 217 CeedElemRestrictionDestroy(&ErestrictxTet); 218 CeedElemRestrictionDestroy(&ErestrictuiTet); 219 CeedElemRestrictionDestroy(&ErestrictxiTet); 220 CeedElemRestrictionDestroy(&ErestrictuHex); 221 CeedElemRestrictionDestroy(&ErestrictxHex); 222 CeedElemRestrictionDestroy(&ErestrictuiHex); 223 CeedElemRestrictionDestroy(&ErestrictxiHex); 224 CeedBasisDestroy(&buTet); 225 CeedBasisDestroy(&bxTet); 226 CeedBasisDestroy(&buHex); 227 CeedBasisDestroy(&bxHex); 228 CeedVectorDestroy(&X); 229 CeedVectorDestroy(&qdataTet); 230 CeedVectorDestroy(&qdataHex); 231 CeedDestroy(&ceed); 232 return 0; 233 } 234