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 PetscBool 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, trimmed, tensor, 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 ierr = PetscDualSpaceLagrangeGetTensor(sp, &key.tensor);CHKERRQ(ierr); 219 } 220 ierr = PetscDualSpaceLagrangeGetTrimmed(sp, &key.trimmed);CHKERRQ(ierr); 221 ierr = PetscInfo4(NULL, "After setup: order %D, trimmed %D, tensor %D, continuous %D\n", key.order, key.trimmed, key.tensor, key.continuous);CHKERRQ(ierr); 222 ierr = PetscHashLagPut(lagTable, key, &iter, &missing);CHKERRQ(ierr); 223 if (missing) { 224 DMPolytopeType type; 225 226 ierr = DMPlexGetCellType(K, 0, &type);CHKERRQ(ierr); 227 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "New space: %s, order %D, trimmed %D, tensor %D, continuous %D, form degree %D\n", DMPolytopeTypes[type], order, trimmed, tensor, continuous, formDegree);CHKERRQ(ierr); 228 ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 229 { 230 PetscQuadrature intNodes, allNodes; 231 Mat intMat, allMat; 232 MatInfo info; 233 PetscInt i, j, nodeIdxDim, nodeVecDim, nNodes; 234 const PetscInt *nodeIdx; 235 const PetscReal *nodeVec; 236 237 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 238 239 ierr = PetscLagNodeIndicesGetData_Internal(lag->allNodeIndices, &nodeIdxDim, &nodeVecDim, &nNodes, &nodeIdx, &nodeVec);CHKERRQ(ierr); 240 if (nodeVecDim != Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect nodeVecDim");CHKERRQ(ierr); 241 if (nNodes != spdim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect nNodes");CHKERRQ(ierr); 242 243 ierr = PetscDualSpaceGetAllData(sp, &allNodes, &allMat);CHKERRQ(ierr); 244 245 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All nodes:\n");CHKERRQ(ierr); 246 ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 247 ierr = PetscQuadratureView(allNodes, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 248 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 249 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All node indices:\n");CHKERRQ(ierr); 250 for (i = 0; i < spdim; i++) { 251 ierr = PetscPrintf(PETSC_COMM_SELF, "(");CHKERRQ(ierr); 252 for (j = 0; j < nodeIdxDim; j++) { 253 ierr = PetscPrintf(PETSC_COMM_SELF, " %D,", nodeIdx[i * nodeIdxDim + j]);CHKERRQ(ierr); 254 } 255 ierr = PetscPrintf(PETSC_COMM_SELF, "): [");CHKERRQ(ierr); 256 for (j = 0; j < nodeVecDim; j++) { 257 ierr = PetscPrintf(PETSC_COMM_SELF, " %g,", nodeVec[i * nodeVecDim + j]);CHKERRQ(ierr); 258 } 259 ierr = PetscPrintf(PETSC_COMM_SELF, "]\n");CHKERRQ(ierr); 260 } 261 262 ierr = MatGetInfo(allMat, MAT_LOCAL, &info);CHKERRQ(ierr); 263 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All matrix: %D nonzeros\n", (PetscInt) info.nz_used);CHKERRQ(ierr); 264 265 ierr = PetscDualSpaceGetInteriorData(sp, &intNodes, &intMat);CHKERRQ(ierr); 266 if (intMat && intMat != allMat) { 267 PetscInt intNodeIdxDim, intNodeVecDim, intNnodes; 268 const PetscInt *intNodeIdx; 269 const PetscReal *intNodeVec; 270 PetscBool same; 271 272 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior nodes:\n");CHKERRQ(ierr); 273 ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 274 ierr = PetscQuadratureView(intNodes, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 275 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 276 277 ierr = MatGetInfo(intMat, MAT_LOCAL, &info);CHKERRQ(ierr); 278 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior matrix: %D nonzeros\n", (PetscInt) info.nz_used);CHKERRQ(ierr); 279 ierr = PetscLagNodeIndicesGetData_Internal(lag->intNodeIndices, &intNodeIdxDim, &intNodeVecDim, &intNnodes, &intNodeIdx, &intNodeVec);CHKERRQ(ierr); 280 if (intNodeIdxDim != nodeIdxDim || intNodeVecDim != nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node indices not the same shale as all node indices"); 281 if (intNnodes != spintdim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect interior nNodes");CHKERRQ(ierr); 282 ierr = PetscArraycmp(intNodeIdx, nodeIdx, nodeIdxDim * intNnodes, &same);CHKERRQ(ierr); 283 if (!same) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node indices not the same as start of all node indices"); 284 ierr = PetscArraycmp(intNodeVec, nodeVec, nodeVecDim * intNnodes, &same);CHKERRQ(ierr); 285 if (!same) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node vectors not the same as start of all node vectors"); 286 } else if (intMat) { 287 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior data is the same as all data\n");CHKERRQ(ierr); 288 if (intNodes != allNodes) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior nodes should be the same as all nodes"); 289 if (lag->intNodeIndices != lag->allNodeIndices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node indices should be the same as all node indices"); 290 } 291 } 292 if (dim <= 2 && spintdim) { 293 PetscInt coneSize, o; 294 295 ierr = DMPlexGetConeSize(K, 0, &coneSize);CHKERRQ(ierr); 296 for (o = -coneSize; o < coneSize; o++) { 297 Mat symMat; 298 299 ierr = PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, o, &symMat);CHKERRQ(ierr); 300 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior node symmetry matrix for orientation %D:\n", o);CHKERRQ(ierr); 301 ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 302 ierr = MatView(symMat, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 303 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 304 ierr = MatDestroy(&symMat);CHKERRQ(ierr); 305 } 306 } 307 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 308 } 309 ierr = PetscDualSpaceDestroy(&sp);CHKERRQ(ierr); 310 PetscFunctionReturn(0); 311 } 312 313 static PetscErrorCode DMPlexCreateReferenceWedge(MPI_Comm comm, DM *refdm) 314 { 315 PetscInt dim = 3; 316 DM rdm; 317 PetscErrorCode ierr; 318 319 PetscFunctionBeginUser; 320 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr); 321 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr); 322 ierr = DMSetDimension(rdm, dim);CHKERRQ(ierr); 323 { 324 PetscInt numPoints[4] = {6, 9, 5, 1}; 325 PetscInt coneSize[21] = {5, 326 3, 3, 327 4, 4, 4, 328 2, 2, 2, 2, 2, 2, 2, 2, 2, 329 0, 0, 0, 0, 0, 0}; 330 PetscInt cones[41] = {1, 2, 3, 4, 5, 331 6, 7, 8, 332 9, 10, 11, 333 8, 12, 9, 13, 334 7, 14, 10, 12, 335 6, 13, 11, 14, 336 15, 16, 16, 17, 17, 15, 337 18, 19, 19, 20, 20, 18, 338 17, 19, 18, 15, 16, 20}; 339 PetscInt coneOrientations[41] = {0, 0, 0, 0, 0, 340 0, 0, 0, 341 0, 0, 0, 342 -2, 0, -2, 0, 343 -2, 0, -2, -2, 344 -2, -2, -2, -2, 345 0, 0, 0, 0, 0, 0, 346 0, 0, 0, 0, 0, 0, 347 0, 0, 0, 0, 0, 0}; 348 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 349 -1.0, 1.0, -1.0, 350 1.0, -1.0, -1.0, 351 -1.0, -1.0, 1.0, 352 1.0, -1.0, 1.0, 353 -1.0, 1.0, 1.0}; 354 355 ierr = DMPlexCreateFromDAG(rdm, 3, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr); 356 } 357 *refdm = rdm; 358 PetscFunctionReturn(0); 359 } 360 361 int main (int argc, char **argv) 362 { 363 PetscInt dim; 364 PetscHashLag lagTable; 365 PetscInt tensorCell; 366 PetscInt order, ordermin, ordermax; 367 PetscBool continuous; 368 PetscBool trimmed; 369 DM dm; 370 PetscErrorCode ierr; 371 372 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 373 dim = 3; 374 tensorCell = 0; 375 continuous = PETSC_FALSE; 376 trimmed = PETSC_FALSE; 377 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Options for PETSCDUALSPACELAGRANGE test","none");CHKERRQ(ierr); 378 ierr = PetscOptionsRangeInt("-dim", "The spatial dimension","ex1.c",dim,&dim,NULL,0,3);CHKERRQ(ierr); 379 ierr = PetscOptionsRangeInt("-tensor", "(0) simplex (1) hypercube (2) wedge","ex1.c",tensorCell,&tensorCell,NULL,0,2);CHKERRQ(ierr); 380 ierr = PetscOptionsBool("-continuous", "Whether the dual space has continuity","ex1.c",continuous,&continuous,NULL);CHKERRQ(ierr); 381 ierr = PetscOptionsBool("-trimmed", "Whether the dual space matches a trimmed polynomial space","ex1.c",trimmed,&trimmed,NULL);CHKERRQ(ierr); 382 ierr = PetscOptionsEnd(); 383 ierr = PetscHashLagCreate(&lagTable);CHKERRQ(ierr); 384 385 if (tensorCell < 2) { 386 ierr = DMPlexCreateReferenceCell(PETSC_COMM_SELF, dim, (PetscBool) !tensorCell, &dm);CHKERRQ(ierr); 387 } else { 388 ierr = DMPlexCreateReferenceWedge(PETSC_COMM_SELF, &dm);CHKERRQ(ierr); 389 } 390 ordermin = trimmed ? 1 : 0; 391 ordermax = tensorCell == 2 ? 4 : tensorCell == 1 ? 3 : dim + 2; 392 for (order = ordermin; order <= ordermax; order++) { 393 PetscInt formDegree; 394 395 for (formDegree = PetscMin(0,-dim+1); formDegree <= dim; formDegree++) { 396 PetscInt nCopies; 397 398 for (nCopies = 1; nCopies <= 3; nCopies++) { 399 ierr = testLagrange(lagTable, dm, dim, order, formDegree, trimmed, (PetscBool) tensorCell, continuous, nCopies);CHKERRQ(ierr); 400 } 401 } 402 } 403 ierr = DMDestroy(&dm);CHKERRQ(ierr); 404 ierr = PetscHashLagDestroy(&lagTable);CHKERRQ(ierr); 405 ierr = PetscFinalize(); 406 return ierr; 407 } 408 409 /*TEST 410 411 test: 412 suffix: 0 413 args: -dim 0 414 415 test: 416 suffix: 1_discontinuous_full 417 args: -dim 1 -continuous 0 -trimmed 0 418 419 test: 420 suffix: 1_continuous_full 421 args: -dim 1 -continuous 1 -trimmed 0 422 423 test: 424 suffix: 2_simplex_discontinuous_full 425 args: -dim 2 -tensor 0 -continuous 0 -trimmed 0 426 427 test: 428 suffix: 2_simplex_continuous_full 429 args: -dim 2 -tensor 0 -continuous 1 -trimmed 0 430 431 test: 432 suffix: 2_tensor_discontinuous_full 433 args: -dim 2 -tensor 1 -continuous 0 -trimmed 0 434 435 test: 436 suffix: 2_tensor_continuous_full 437 args: -dim 2 -tensor 1 -continuous 1 -trimmed 0 438 439 test: 440 suffix: 3_simplex_discontinuous_full 441 args: -dim 3 -tensor 0 -continuous 0 -trimmed 0 442 443 test: 444 suffix: 3_simplex_continuous_full 445 args: -dim 3 -tensor 0 -continuous 1 -trimmed 0 446 447 test: 448 suffix: 3_tensor_discontinuous_full 449 args: -dim 3 -tensor 1 -continuous 0 -trimmed 0 450 451 test: 452 suffix: 3_tensor_continuous_full 453 args: -dim 3 -tensor 1 -continuous 1 -trimmed 0 454 455 test: 456 suffix: 3_wedge_discontinuous_full 457 args: -dim 3 -tensor 2 -continuous 0 -trimmed 0 458 459 test: 460 suffix: 3_wedge_continuous_full 461 args: -dim 3 -tensor 2 -continuous 1 -trimmed 0 462 463 test: 464 suffix: 1_discontinuous_trimmed 465 args: -dim 1 -continuous 0 -trimmed 1 466 467 test: 468 suffix: 1_continuous_trimmed 469 args: -dim 1 -continuous 1 -trimmed 1 470 471 test: 472 suffix: 2_simplex_discontinuous_trimmed 473 args: -dim 2 -tensor 0 -continuous 0 -trimmed 1 474 475 test: 476 suffix: 2_simplex_continuous_trimmed 477 args: -dim 2 -tensor 0 -continuous 1 -trimmed 1 478 479 test: 480 suffix: 2_tensor_discontinuous_trimmed 481 args: -dim 2 -tensor 1 -continuous 0 -trimmed 1 482 483 test: 484 suffix: 2_tensor_continuous_trimmed 485 args: -dim 2 -tensor 1 -continuous 1 -trimmed 1 486 487 test: 488 suffix: 3_simplex_discontinuous_trimmed 489 args: -dim 3 -tensor 0 -continuous 0 -trimmed 1 490 491 test: 492 suffix: 3_simplex_continuous_trimmed 493 args: -dim 3 -tensor 0 -continuous 1 -trimmed 1 494 495 test: 496 suffix: 3_tensor_discontinuous_trimmed 497 args: -dim 3 -tensor 1 -continuous 0 -trimmed 1 498 499 test: 500 suffix: 3_tensor_continuous_trimmed 501 args: -dim 3 -tensor 1 -continuous 1 -trimmed 1 502 503 test: 504 suffix: 3_wedge_discontinuous_trimmed 505 args: -dim 3 -tensor 2 -continuous 0 -trimmed 1 506 507 test: 508 suffix: 3_wedge_continuous_trimmed 509 args: -dim 3 -tensor 2 -continuous 1 -trimmed 1 510 511 TEST*/ 512