1 2 #include <petscfe.h> 3 #include <petscdmplex.h> 4 #include <petsc/private/hashmap.h> 5 #include <petsc/private/dmpleximpl.h> 6 #include <petsc/private/petscfeimpl.h> 7 8 const char help[] = "Test PETSCDUALSPACELAGRANGE\n"; 9 10 typedef struct _PetscHashLagKey 11 { 12 PetscInt dim; 13 PetscInt order; 14 PetscInt formDegree; 15 PetscBool trimmed; 16 PetscInt tensor; 17 PetscBool continuous; 18 } PetscHashLagKey; 19 20 #define PetscHashLagKeyHash(key) \ 21 PetscHashCombine(PetscHashCombine(PetscHashCombine(PetscHashInt((key).dim), \ 22 PetscHashInt((key).order)), \ 23 PetscHashInt((key).formDegree)), \ 24 PetscHashCombine(PetscHashCombine(PetscHashInt((key).trimmed), \ 25 PetscHashInt((key).tensor)), \ 26 PetscHashInt((key).continuous))) 27 28 #define PetscHashLagKeyEqual(k1,k2) \ 29 (((k1).dim == (k2).dim) ? \ 30 ((k1).order == (k2).order) ? \ 31 ((k1).formDegree == (k2).formDegree) ? \ 32 ((k1).trimmed == (k2).trimmed) ? \ 33 ((k1).tensor == (k2).tensor) ? \ 34 ((k1).continuous == (k2).continuous) : 0 : 0 : 0 : 0 : 0) 35 36 PETSC_HASH_MAP(HashLag, PetscHashLagKey, PetscInt, PetscHashLagKeyHash, PetscHashLagKeyEqual, 0) 37 38 static PetscErrorCode ExpectedNumDofs_Total(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs); 39 static PetscErrorCode ExpectedNumDofs_Interior(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs); 40 41 static PetscErrorCode ExpectedNumDofs_Total(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs) 42 { 43 PetscErrorCode ierr; 44 45 PetscFunctionBegin; 46 formDegree = PetscAbsInt(formDegree); 47 /* see femtable.org for the source of most of these values */ 48 *nDofs = -1; 49 if (tensor == 0) { /* simplex spaces */ 50 if (trimmed) { 51 PetscInt rnchooserk; 52 PetscInt rkm1choosek; 53 54 ierr = PetscDTBinomialInt(order + dim, order + formDegree, &rnchooserk);CHKERRQ(ierr); 55 ierr = PetscDTBinomialInt(order + formDegree - 1, formDegree, &rkm1choosek);CHKERRQ(ierr); 56 *nDofs = rnchooserk * rkm1choosek * nCopies; 57 } else { 58 PetscInt rnchooserk; 59 PetscInt rkchoosek; 60 61 ierr = PetscDTBinomialInt(order + dim, order + formDegree, &rnchooserk);CHKERRQ(ierr); 62 ierr = PetscDTBinomialInt(order + formDegree, formDegree, &rkchoosek);CHKERRQ(ierr); 63 *nDofs = rnchooserk * rkchoosek * nCopies; 64 } 65 } else if (tensor == 1) { /* hypercubes */ 66 if (trimmed) { 67 PetscInt nchoosek; 68 PetscInt rpowk, rp1pownmk; 69 70 ierr = PetscDTBinomialInt(dim, formDegree, &nchoosek);CHKERRQ(ierr); 71 rpowk = PetscPowInt(order, formDegree);CHKERRQ(ierr); 72 rp1pownmk = PetscPowInt(order + 1, dim - formDegree);CHKERRQ(ierr); 73 *nDofs = nchoosek * rpowk * rp1pownmk * nCopies; 74 } else { 75 PetscInt nchoosek; 76 PetscInt rp1pown; 77 78 ierr = PetscDTBinomialInt(dim, formDegree, &nchoosek);CHKERRQ(ierr); 79 rp1pown = PetscPowInt(order + 1, dim);CHKERRQ(ierr); 80 *nDofs = nchoosek * rp1pown * nCopies; 81 } 82 } else { /* prism */ 83 PetscInt tracek = 0; 84 PetscInt tracekm1 = 0; 85 PetscInt fiber0 = 0; 86 PetscInt fiber1 = 0; 87 88 if (formDegree < dim) { 89 ierr = ExpectedNumDofs_Total(dim - 1, order, formDegree, trimmed, 0, 1, &tracek);CHKERRQ(ierr); 90 ierr = ExpectedNumDofs_Total(1, order, 0, trimmed, 0, 1, &fiber0);CHKERRQ(ierr); 91 } 92 if (formDegree > 0) { 93 ierr = ExpectedNumDofs_Total(dim - 1, order, formDegree - 1, trimmed, 0, 1, &tracekm1);CHKERRQ(ierr); 94 ierr = ExpectedNumDofs_Total(1, order, 1, trimmed, 0, 1, &fiber1);CHKERRQ(ierr); 95 } 96 *nDofs = (tracek * fiber0 + tracekm1 * fiber1) * nCopies; 97 } 98 PetscFunctionReturn(0); 99 } 100 101 static PetscErrorCode ExpectedNumDofs_Interior(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, 102 PetscInt tensor, PetscInt nCopies, PetscInt *nDofs) 103 { 104 PetscErrorCode ierr; 105 106 PetscFunctionBegin; 107 formDegree = PetscAbsInt(formDegree); 108 /* see femtable.org for the source of most of these values */ 109 *nDofs = -1; 110 if (tensor == 0) { /* simplex spaces */ 111 if (trimmed) { 112 if (order + formDegree > dim) { 113 PetscInt eorder = order + formDegree - dim - 1; 114 PetscInt eformDegree = dim - formDegree; 115 PetscInt rnchooserk; 116 PetscInt rkchoosek; 117 118 ierr = PetscDTBinomialInt(eorder + dim, eorder + eformDegree, &rnchooserk);CHKERRQ(ierr); 119 ierr = PetscDTBinomialInt(eorder + eformDegree, eformDegree, &rkchoosek);CHKERRQ(ierr); 120 *nDofs = rnchooserk * rkchoosek * nCopies; 121 } else { 122 *nDofs = 0; 123 } 124 125 } else { 126 if (order + formDegree > dim) { 127 PetscInt eorder = order + formDegree - dim; 128 PetscInt eformDegree = dim - formDegree; 129 PetscInt rnchooserk; 130 PetscInt rkm1choosek; 131 132 ierr = PetscDTBinomialInt(eorder + dim, eorder + eformDegree, &rnchooserk);CHKERRQ(ierr); 133 ierr = PetscDTBinomialInt(eorder + eformDegree - 1, eformDegree, &rkm1choosek);CHKERRQ(ierr); 134 *nDofs = rnchooserk * rkm1choosek * nCopies; 135 } else { 136 *nDofs = 0; 137 } 138 } 139 } else if (tensor == 1) { /* hypercubes */ 140 if (dim < 2) { 141 ierr = ExpectedNumDofs_Interior(dim, order, formDegree, trimmed, 0, nCopies, nDofs);CHKERRQ(ierr); 142 } else { 143 PetscInt tracek = 0; 144 PetscInt tracekm1 = 0; 145 PetscInt fiber0 = 0; 146 PetscInt fiber1 = 0; 147 148 if (formDegree < dim) { 149 ierr = ExpectedNumDofs_Interior(dim - 1, order, formDegree, trimmed, dim > 2 ? 1 : 0, 1, &tracek);CHKERRQ(ierr); 150 ierr = ExpectedNumDofs_Interior(1, order, 0, trimmed, 0, 1, &fiber0);CHKERRQ(ierr); 151 } 152 if (formDegree > 0) { 153 ierr = ExpectedNumDofs_Interior(dim - 1, order, formDegree - 1, trimmed, dim > 2 ? 1 : 0, 1, &tracekm1);CHKERRQ(ierr); 154 ierr = ExpectedNumDofs_Interior(1, order, 1, trimmed, 0, 1, &fiber1);CHKERRQ(ierr); 155 } 156 *nDofs = (tracek * fiber0 + tracekm1 * fiber1) * nCopies; 157 } 158 } else { /* prism */ 159 PetscInt tracek = 0; 160 PetscInt tracekm1 = 0; 161 PetscInt fiber0 = 0; 162 PetscInt fiber1 = 0; 163 164 if (formDegree < dim) { 165 ierr = ExpectedNumDofs_Interior(dim - 1, order, formDegree, trimmed, 0, 1, &tracek);CHKERRQ(ierr); 166 ierr = ExpectedNumDofs_Interior(1, order, 0, trimmed, 0, 1, &fiber0);CHKERRQ(ierr); 167 } 168 if (formDegree > 0) { 169 ierr = ExpectedNumDofs_Interior(dim - 1, order, formDegree - 1, trimmed, 0, 1, &tracekm1);CHKERRQ(ierr); 170 ierr = ExpectedNumDofs_Interior(1, order, 1, trimmed, 0, 1, &fiber1);CHKERRQ(ierr); 171 } 172 *nDofs = (tracek * fiber0 + tracekm1 * fiber1) * nCopies; 173 } 174 PetscFunctionReturn(0); 175 } 176 177 PetscErrorCode testLagrange(PetscHashLag lagTable, DM K, PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscBool continuous, PetscInt nCopies) 178 { 179 PetscDualSpace sp; 180 MPI_Comm comm = PETSC_COMM_SELF; 181 PetscInt Nk; 182 PetscHashLagKey key; 183 PetscHashIter iter; 184 PetscBool missing; 185 PetscInt spdim, spintdim, exspdim, exspintdim; 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);CHKERRQ(ierr); 190 ierr = PetscDualSpaceCreate(comm, &sp);CHKERRQ(ierr); 191 ierr = PetscDualSpaceSetType(sp, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 192 ierr = PetscDualSpaceSetDM(sp, K);CHKERRQ(ierr); 193 ierr = PetscDualSpaceSetOrder(sp, order);CHKERRQ(ierr); 194 ierr = PetscDualSpaceSetFormDegree(sp, formDegree);CHKERRQ(ierr); 195 ierr = PetscDualSpaceSetNumComponents(sp, nCopies * Nk);CHKERRQ(ierr); 196 ierr = PetscDualSpaceLagrangeSetContinuity(sp, continuous);CHKERRQ(ierr); 197 ierr = PetscDualSpaceLagrangeSetTensor(sp, (PetscBool) tensor);CHKERRQ(ierr); 198 ierr = PetscDualSpaceLagrangeSetTrimmed(sp, trimmed);CHKERRQ(ierr); 199 ierr = PetscInfo7(NULL, "Input: dim %D, order %D, trimmed %D, tensor %D, continuous %D, formDegree %D, nCopies %D\n", dim, order, (PetscInt) trimmed, tensor, (PetscInt) continuous, formDegree, nCopies);CHKERRQ(ierr); 200 ierr = ExpectedNumDofs_Total(dim, order, formDegree, trimmed, tensor, nCopies, &exspdim);CHKERRQ(ierr); 201 if (continuous && dim > 0 && order > 0) { 202 ierr = ExpectedNumDofs_Interior(dim, order, formDegree, trimmed, tensor, nCopies, &exspintdim);CHKERRQ(ierr); 203 } else { 204 exspintdim = exspdim; 205 } 206 ierr = PetscDualSpaceSetUp(sp);CHKERRQ(ierr); 207 ierr = PetscDualSpaceGetDimension(sp, &spdim);CHKERRQ(ierr); 208 ierr = PetscDualSpaceGetInteriorDimension(sp, &spintdim);CHKERRQ(ierr); 209 if (spdim != exspdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected dual space dimension %D, got %D\n", exspdim, spdim);CHKERRQ(ierr); 210 if (spintdim != exspintdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected dual space interior dimension %D, got %D\n", exspintdim, spintdim);CHKERRQ(ierr); 211 key.dim = dim; 212 key.formDegree = formDegree; 213 ierr = PetscDualSpaceGetOrder(sp, &key.order);CHKERRQ(ierr); 214 ierr = PetscDualSpaceLagrangeGetContinuity(sp, &key.continuous);CHKERRQ(ierr); 215 if (tensor == 2) { 216 key.tensor = 2; 217 } else { 218 PetscBool bTensor; 219 220 ierr = PetscDualSpaceLagrangeGetTensor(sp, &bTensor);CHKERRQ(ierr); 221 key.tensor = bTensor; 222 } 223 ierr = PetscDualSpaceLagrangeGetTrimmed(sp, &key.trimmed);CHKERRQ(ierr); 224 ierr = PetscInfo4(NULL, "After setup: order %D, trimmed %D, tensor %D, continuous %D\n", key.order, (PetscInt) key.trimmed, key.tensor, (PetscInt) key.continuous);CHKERRQ(ierr); 225 ierr = PetscHashLagPut(lagTable, key, &iter, &missing);CHKERRQ(ierr); 226 if (missing) { 227 DMPolytopeType type; 228 229 ierr = DMPlexGetCellType(K, 0, &type);CHKERRQ(ierr); 230 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "New space: %s, order %D, trimmed %D, tensor %D, continuous %D, form degree %D\n", DMPolytopeTypes[type], order, (PetscInt) trimmed, tensor, (PetscInt) continuous, formDegree);CHKERRQ(ierr); 231 ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 232 { 233 PetscQuadrature intNodes, allNodes; 234 Mat intMat, allMat; 235 MatInfo info; 236 PetscInt i, j, nodeIdxDim, nodeVecDim, nNodes; 237 const PetscInt *nodeIdx; 238 const PetscReal *nodeVec; 239 240 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 241 242 ierr = PetscLagNodeIndicesGetData_Internal(lag->allNodeIndices, &nodeIdxDim, &nodeVecDim, &nNodes, &nodeIdx, &nodeVec);CHKERRQ(ierr); 243 if (nodeVecDim != Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect nodeVecDim");CHKERRQ(ierr); 244 if (nNodes != spdim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect nNodes");CHKERRQ(ierr); 245 246 ierr = PetscDualSpaceGetAllData(sp, &allNodes, &allMat);CHKERRQ(ierr); 247 248 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All nodes:\n");CHKERRQ(ierr); 249 ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 250 ierr = PetscQuadratureView(allNodes, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 251 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 252 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All node indices:\n");CHKERRQ(ierr); 253 for (i = 0; i < spdim; i++) { 254 ierr = PetscPrintf(PETSC_COMM_SELF, "(");CHKERRQ(ierr); 255 for (j = 0; j < nodeIdxDim; j++) { 256 ierr = PetscPrintf(PETSC_COMM_SELF, " %D,", nodeIdx[i * nodeIdxDim + j]);CHKERRQ(ierr); 257 } 258 ierr = PetscPrintf(PETSC_COMM_SELF, "): [");CHKERRQ(ierr); 259 for (j = 0; j < nodeVecDim; j++) { 260 ierr = PetscPrintf(PETSC_COMM_SELF, " %g,", (double) nodeVec[i * nodeVecDim + j]);CHKERRQ(ierr); 261 } 262 ierr = PetscPrintf(PETSC_COMM_SELF, "]\n");CHKERRQ(ierr); 263 } 264 265 ierr = MatGetInfo(allMat, MAT_LOCAL, &info);CHKERRQ(ierr); 266 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All matrix: %D nonzeros\n", (PetscInt) info.nz_used);CHKERRQ(ierr); 267 268 ierr = PetscDualSpaceGetInteriorData(sp, &intNodes, &intMat);CHKERRQ(ierr); 269 if (intMat && intMat != allMat) { 270 PetscInt intNodeIdxDim, intNodeVecDim, intNnodes; 271 const PetscInt *intNodeIdx; 272 const PetscReal *intNodeVec; 273 PetscBool same; 274 275 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior nodes:\n");CHKERRQ(ierr); 276 ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 277 ierr = PetscQuadratureView(intNodes, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 278 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 279 280 ierr = MatGetInfo(intMat, MAT_LOCAL, &info);CHKERRQ(ierr); 281 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior matrix: %D nonzeros\n", (PetscInt) info.nz_used);CHKERRQ(ierr); 282 ierr = PetscLagNodeIndicesGetData_Internal(lag->intNodeIndices, &intNodeIdxDim, &intNodeVecDim, &intNnodes, &intNodeIdx, &intNodeVec);CHKERRQ(ierr); 283 if (intNodeIdxDim != nodeIdxDim || intNodeVecDim != nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node indices not the same shale as all node indices"); 284 if (intNnodes != spintdim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect interior nNodes");CHKERRQ(ierr); 285 ierr = PetscArraycmp(intNodeIdx, nodeIdx, nodeIdxDim * intNnodes, &same);CHKERRQ(ierr); 286 if (!same) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node indices not the same as start of all node indices"); 287 ierr = PetscArraycmp(intNodeVec, nodeVec, nodeVecDim * intNnodes, &same);CHKERRQ(ierr); 288 if (!same) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node vectors not the same as start of all node vectors"); 289 } else if (intMat) { 290 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior data is the same as all data\n");CHKERRQ(ierr); 291 if (intNodes != allNodes) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior nodes should be the same as all nodes"); 292 if (lag->intNodeIndices != lag->allNodeIndices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node indices should be the same as all node indices"); 293 } 294 } 295 if (dim <= 2 && spintdim) { 296 PetscInt coneSize, o; 297 298 ierr = DMPlexGetConeSize(K, 0, &coneSize);CHKERRQ(ierr); 299 for (o = -coneSize; o < coneSize; o++) { 300 Mat symMat; 301 302 ierr = PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, o, &symMat);CHKERRQ(ierr); 303 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior node symmetry matrix for orientation %D:\n", o);CHKERRQ(ierr); 304 ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 305 ierr = MatView(symMat, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 306 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 307 ierr = MatDestroy(&symMat);CHKERRQ(ierr); 308 } 309 } 310 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 311 } 312 ierr = PetscDualSpaceDestroy(&sp);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 static PetscErrorCode DMPlexCreateReferenceWedge(MPI_Comm comm, DM *refdm) 317 { 318 PetscInt dim = 3; 319 DM rdm; 320 PetscErrorCode ierr; 321 322 PetscFunctionBeginUser; 323 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 324 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 325 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 326 { 327 PetscInt numPoints[4] = {6, 9, 5, 1}; 328 PetscInt coneSize[21] = {5, 329 3, 3, 330 4, 4, 4, 331 2, 2, 2, 2, 2, 2, 2, 2, 2, 332 0, 0, 0, 0, 0, 0}; 333 PetscInt cones[41] = {1, 2, 3, 4, 5, 334 6, 7, 8, 335 9, 10, 11, 336 8, 12, 9, 13, 337 7, 14, 10, 12, 338 6, 13, 11, 14, 339 15, 16, 16, 17, 17, 15, 340 18, 19, 19, 20, 20, 18, 341 17, 19, 18, 15, 16, 20}; 342 PetscInt coneOrientations[41] = {0, 0, 0, 0, 0, 343 0, 0, 0, 344 0, 0, 0, 345 -2, 0, -2, 0, 346 -2, 0, -2, -2, 347 -2, -2, -2, -2, 348 0, 0, 0, 0, 0, 0, 349 0, 0, 0, 0, 0, 0, 350 0, 0, 0, 0, 0, 0}; 351 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 352 -1.0, 1.0, -1.0, 353 1.0, -1.0, -1.0, 354 -1.0, -1.0, 1.0, 355 1.0, -1.0, 1.0, 356 -1.0, 1.0, 1.0}; 357 358 ierr = DMPlexCreateFromDAG(rdm, 3, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 359 } 360 *refdm = rdm; 361 PetscFunctionReturn(0); 362 } 363 364 int main (int argc, char **argv) 365 { 366 PetscInt dim; 367 PetscHashLag lagTable; 368 PetscInt tensorCell; 369 PetscInt order, ordermin, ordermax; 370 PetscBool continuous; 371 PetscBool trimmed; 372 DM dm; 373 PetscErrorCode ierr; 374 375 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 376 dim = 3; 377 tensorCell = 0; 378 continuous = PETSC_FALSE; 379 trimmed = PETSC_FALSE; 380 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Options for PETSCDUALSPACELAGRANGE test","none");CHKERRQ(ierr); 381 ierr = PetscOptionsRangeInt("-dim", "The spatial dimension","ex1.c",dim,&dim,NULL,0,3);CHKERRQ(ierr); 382 ierr = PetscOptionsRangeInt("-tensor", "(0) simplex (1) hypercube (2) wedge","ex1.c",tensorCell,&tensorCell,NULL,0,2);CHKERRQ(ierr); 383 ierr = PetscOptionsBool("-continuous", "Whether the dual space has continuity","ex1.c",continuous,&continuous,NULL);CHKERRQ(ierr); 384 ierr = PetscOptionsBool("-trimmed", "Whether the dual space matches a trimmed polynomial space","ex1.c",trimmed,&trimmed,NULL);CHKERRQ(ierr); 385 ierr = PetscOptionsEnd(); 386 ierr = PetscHashLagCreate(&lagTable);CHKERRQ(ierr); 387 388 if (tensorCell < 2) { 389 ierr = DMPlexCreateReferenceCell(PETSC_COMM_SELF, dim, (PetscBool) !tensorCell, &dm);CHKERRQ(ierr); 390 } else { 391 ierr = DMPlexCreateReferenceWedge(PETSC_COMM_SELF, &dm);CHKERRQ(ierr); 392 } 393 ordermin = trimmed ? 1 : 0; 394 ordermax = tensorCell == 2 ? 4 : tensorCell == 1 ? 3 : dim + 2; 395 for (order = ordermin; order <= ordermax; order++) { 396 PetscInt formDegree; 397 398 for (formDegree = PetscMin(0,-dim+1); formDegree <= dim; formDegree++) { 399 PetscInt nCopies; 400 401 for (nCopies = 1; nCopies <= 3; nCopies++) { 402 ierr = testLagrange(lagTable, dm, dim, order, formDegree, trimmed, (PetscBool) tensorCell, continuous, nCopies);CHKERRQ(ierr); 403 } 404 } 405 } 406 ierr = DMDestroy(&dm);CHKERRQ(ierr); 407 ierr = PetscHashLagDestroy(&lagTable);CHKERRQ(ierr); 408 ierr = PetscFinalize(); 409 return ierr; 410 } 411 412 /*TEST 413 414 test: 415 suffix: 0 416 args: -dim 0 417 418 test: 419 suffix: 1_discontinuous_full 420 args: -dim 1 -continuous 0 -trimmed 0 421 422 test: 423 suffix: 1_continuous_full 424 args: -dim 1 -continuous 1 -trimmed 0 425 426 test: 427 suffix: 2_simplex_discontinuous_full 428 args: -dim 2 -tensor 0 -continuous 0 -trimmed 0 429 430 test: 431 suffix: 2_simplex_continuous_full 432 args: -dim 2 -tensor 0 -continuous 1 -trimmed 0 433 434 test: 435 suffix: 2_tensor_discontinuous_full 436 args: -dim 2 -tensor 1 -continuous 0 -trimmed 0 437 438 test: 439 suffix: 2_tensor_continuous_full 440 args: -dim 2 -tensor 1 -continuous 1 -trimmed 0 441 442 test: 443 suffix: 3_simplex_discontinuous_full 444 args: -dim 3 -tensor 0 -continuous 0 -trimmed 0 445 446 test: 447 suffix: 3_simplex_continuous_full 448 args: -dim 3 -tensor 0 -continuous 1 -trimmed 0 449 450 test: 451 suffix: 3_tensor_discontinuous_full 452 args: -dim 3 -tensor 1 -continuous 0 -trimmed 0 453 454 test: 455 suffix: 3_tensor_continuous_full 456 args: -dim 3 -tensor 1 -continuous 1 -trimmed 0 457 458 test: 459 suffix: 3_wedge_discontinuous_full 460 args: -dim 3 -tensor 2 -continuous 0 -trimmed 0 461 462 test: 463 suffix: 3_wedge_continuous_full 464 args: -dim 3 -tensor 2 -continuous 1 -trimmed 0 465 466 test: 467 suffix: 1_discontinuous_trimmed 468 args: -dim 1 -continuous 0 -trimmed 1 469 470 test: 471 suffix: 1_continuous_trimmed 472 args: -dim 1 -continuous 1 -trimmed 1 473 474 test: 475 suffix: 2_simplex_discontinuous_trimmed 476 args: -dim 2 -tensor 0 -continuous 0 -trimmed 1 477 478 test: 479 suffix: 2_simplex_continuous_trimmed 480 args: -dim 2 -tensor 0 -continuous 1 -trimmed 1 481 482 test: 483 suffix: 2_tensor_discontinuous_trimmed 484 args: -dim 2 -tensor 1 -continuous 0 -trimmed 1 485 486 test: 487 suffix: 2_tensor_continuous_trimmed 488 args: -dim 2 -tensor 1 -continuous 1 -trimmed 1 489 490 test: 491 suffix: 3_simplex_discontinuous_trimmed 492 args: -dim 3 -tensor 0 -continuous 0 -trimmed 1 493 494 test: 495 suffix: 3_simplex_continuous_trimmed 496 args: -dim 3 -tensor 0 -continuous 1 -trimmed 1 497 498 test: 499 suffix: 3_tensor_discontinuous_trimmed 500 args: -dim 3 -tensor 1 -continuous 0 -trimmed 1 501 502 test: 503 suffix: 3_tensor_continuous_trimmed 504 args: -dim 3 -tensor 1 -continuous 1 -trimmed 1 505 506 test: 507 suffix: 3_wedge_discontinuous_trimmed 508 args: -dim 3 -tensor 2 -continuous 0 -trimmed 1 509 510 test: 511 suffix: 3_wedge_continuous_trimmed 512 args: -dim 3 -tensor 2 -continuous 1 -trimmed 1 513 514 TEST*/ 515