18a820b8eSToby Isaac 28a820b8eSToby Isaac #include <petscfe.h> 38a820b8eSToby Isaac #include <petscdmplex.h> 48a820b8eSToby Isaac #include <petsc/private/hashmap.h> 58a820b8eSToby Isaac #include <petsc/private/dmpleximpl.h> 68a820b8eSToby Isaac #include <petsc/private/petscfeimpl.h> 78a820b8eSToby Isaac 88a820b8eSToby Isaac const char help[] = "Test PETSCDUALSPACELAGRANGE\n"; 98a820b8eSToby Isaac 108a820b8eSToby Isaac typedef struct _PetscHashLagKey 118a820b8eSToby Isaac { 128a820b8eSToby Isaac PetscInt dim; 138a820b8eSToby Isaac PetscInt order; 148a820b8eSToby Isaac PetscInt formDegree; 158a820b8eSToby Isaac PetscBool trimmed; 166ff15688SToby Isaac PetscInt tensor; 178a820b8eSToby Isaac PetscBool continuous; 188a820b8eSToby Isaac } PetscHashLagKey; 198a820b8eSToby Isaac 208a820b8eSToby Isaac #define PetscHashLagKeyHash(key) \ 218a820b8eSToby Isaac PetscHashCombine(PetscHashCombine(PetscHashCombine(PetscHashInt((key).dim), \ 228a820b8eSToby Isaac PetscHashInt((key).order)), \ 238a820b8eSToby Isaac PetscHashInt((key).formDegree)), \ 248a820b8eSToby Isaac PetscHashCombine(PetscHashCombine(PetscHashInt((key).trimmed), \ 258a820b8eSToby Isaac PetscHashInt((key).tensor)), \ 268a820b8eSToby Isaac PetscHashInt((key).continuous))) 278a820b8eSToby Isaac 288a820b8eSToby Isaac #define PetscHashLagKeyEqual(k1,k2) \ 298a820b8eSToby Isaac (((k1).dim == (k2).dim) ? \ 308a820b8eSToby Isaac ((k1).order == (k2).order) ? \ 318a820b8eSToby Isaac ((k1).formDegree == (k2).formDegree) ? \ 328a820b8eSToby Isaac ((k1).trimmed == (k2).trimmed) ? \ 338a820b8eSToby Isaac ((k1).tensor == (k2).tensor) ? \ 348a820b8eSToby Isaac ((k1).continuous == (k2).continuous) : 0 : 0 : 0 : 0 : 0) 358a820b8eSToby Isaac 368a820b8eSToby Isaac PETSC_HASH_MAP(HashLag, PetscHashLagKey, PetscInt, PetscHashLagKeyHash, PetscHashLagKeyEqual, 0) 378a820b8eSToby Isaac 388a820b8eSToby Isaac static PetscErrorCode ExpectedNumDofs_Total(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs); 398a820b8eSToby Isaac static PetscErrorCode ExpectedNumDofs_Interior(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs); 408a820b8eSToby Isaac 418a820b8eSToby Isaac static PetscErrorCode ExpectedNumDofs_Total(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs) 428a820b8eSToby Isaac { 438a820b8eSToby Isaac PetscFunctionBegin; 448a820b8eSToby Isaac formDegree = PetscAbsInt(formDegree); 458a820b8eSToby Isaac /* see femtable.org for the source of most of these values */ 468a820b8eSToby Isaac *nDofs = -1; 478a820b8eSToby Isaac if (tensor == 0) { /* simplex spaces */ 488a820b8eSToby Isaac if (trimmed) { 498a820b8eSToby Isaac PetscInt rnchooserk; 508a820b8eSToby Isaac PetscInt rkm1choosek; 518a820b8eSToby Isaac 525f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDTBinomialInt(order + dim, order + formDegree, &rnchooserk)); 535f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDTBinomialInt(order + formDegree - 1, formDegree, &rkm1choosek)); 548a820b8eSToby Isaac *nDofs = rnchooserk * rkm1choosek * nCopies; 558a820b8eSToby Isaac } else { 568a820b8eSToby Isaac PetscInt rnchooserk; 578a820b8eSToby Isaac PetscInt rkchoosek; 588a820b8eSToby Isaac 595f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDTBinomialInt(order + dim, order + formDegree, &rnchooserk)); 605f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDTBinomialInt(order + formDegree, formDegree, &rkchoosek)); 618a820b8eSToby Isaac *nDofs = rnchooserk * rkchoosek * nCopies; 628a820b8eSToby Isaac } 638a820b8eSToby Isaac } else if (tensor == 1) { /* hypercubes */ 648a820b8eSToby Isaac if (trimmed) { 658a820b8eSToby Isaac PetscInt nchoosek; 668a820b8eSToby Isaac PetscInt rpowk, rp1pownmk; 678a820b8eSToby Isaac 685f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDTBinomialInt(dim, formDegree, &nchoosek)); 692f613bf5SBarry Smith rpowk = PetscPowInt(order, formDegree); 702f613bf5SBarry Smith rp1pownmk = PetscPowInt(order + 1, dim - formDegree); 718a820b8eSToby Isaac *nDofs = nchoosek * rpowk * rp1pownmk * nCopies; 728a820b8eSToby Isaac } else { 738a820b8eSToby Isaac PetscInt nchoosek; 748a820b8eSToby Isaac PetscInt rp1pown; 758a820b8eSToby Isaac 765f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDTBinomialInt(dim, formDegree, &nchoosek)); 772f613bf5SBarry Smith rp1pown = PetscPowInt(order + 1, dim); 788a820b8eSToby Isaac *nDofs = nchoosek * rp1pown * nCopies; 798a820b8eSToby Isaac } 808a820b8eSToby Isaac } else { /* prism */ 818a820b8eSToby Isaac PetscInt tracek = 0; 828a820b8eSToby Isaac PetscInt tracekm1 = 0; 838a820b8eSToby Isaac PetscInt fiber0 = 0; 848a820b8eSToby Isaac PetscInt fiber1 = 0; 858a820b8eSToby Isaac 868a820b8eSToby Isaac if (formDegree < dim) { 875f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Total(dim - 1, order, formDegree, trimmed, 0, 1, &tracek)); 885f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Total(1, order, 0, trimmed, 0, 1, &fiber0)); 898a820b8eSToby Isaac } 908a820b8eSToby Isaac if (formDegree > 0) { 915f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Total(dim - 1, order, formDegree - 1, trimmed, 0, 1, &tracekm1)); 925f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Total(1, order, 1, trimmed, 0, 1, &fiber1)); 938a820b8eSToby Isaac } 948a820b8eSToby Isaac *nDofs = (tracek * fiber0 + tracekm1 * fiber1) * nCopies; 958a820b8eSToby Isaac } 968a820b8eSToby Isaac PetscFunctionReturn(0); 978a820b8eSToby Isaac } 988a820b8eSToby Isaac 998a820b8eSToby Isaac static PetscErrorCode ExpectedNumDofs_Interior(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, 1008a820b8eSToby Isaac PetscInt tensor, PetscInt nCopies, PetscInt *nDofs) 1018a820b8eSToby Isaac { 1028a820b8eSToby Isaac PetscFunctionBegin; 1038a820b8eSToby Isaac formDegree = PetscAbsInt(formDegree); 1048a820b8eSToby Isaac /* see femtable.org for the source of most of these values */ 1058a820b8eSToby Isaac *nDofs = -1; 1068a820b8eSToby Isaac if (tensor == 0) { /* simplex spaces */ 1078a820b8eSToby Isaac if (trimmed) { 1088a820b8eSToby Isaac if (order + formDegree > dim) { 1098a820b8eSToby Isaac PetscInt eorder = order + formDegree - dim - 1; 1108a820b8eSToby Isaac PetscInt eformDegree = dim - formDegree; 1118a820b8eSToby Isaac PetscInt rnchooserk; 1128a820b8eSToby Isaac PetscInt rkchoosek; 1138a820b8eSToby Isaac 1145f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDTBinomialInt(eorder + dim, eorder + eformDegree, &rnchooserk)); 1155f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDTBinomialInt(eorder + eformDegree, eformDegree, &rkchoosek)); 1168a820b8eSToby Isaac *nDofs = rnchooserk * rkchoosek * nCopies; 1178a820b8eSToby Isaac } else { 1188a820b8eSToby Isaac *nDofs = 0; 1198a820b8eSToby Isaac } 1208a820b8eSToby Isaac 1218a820b8eSToby Isaac } else { 1228a820b8eSToby Isaac if (order + formDegree > dim) { 1238a820b8eSToby Isaac PetscInt eorder = order + formDegree - dim; 1248a820b8eSToby Isaac PetscInt eformDegree = dim - formDegree; 1258a820b8eSToby Isaac PetscInt rnchooserk; 1268a820b8eSToby Isaac PetscInt rkm1choosek; 1278a820b8eSToby Isaac 1285f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDTBinomialInt(eorder + dim, eorder + eformDegree, &rnchooserk)); 1295f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDTBinomialInt(eorder + eformDegree - 1, eformDegree, &rkm1choosek)); 1308a820b8eSToby Isaac *nDofs = rnchooserk * rkm1choosek * nCopies; 1318a820b8eSToby Isaac } else { 1328a820b8eSToby Isaac *nDofs = 0; 1338a820b8eSToby Isaac } 1348a820b8eSToby Isaac } 1358a820b8eSToby Isaac } else if (tensor == 1) { /* hypercubes */ 1368a820b8eSToby Isaac if (dim < 2) { 1375f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Interior(dim, order, formDegree, trimmed, 0, nCopies, nDofs)); 1388a820b8eSToby Isaac } else { 1398a820b8eSToby Isaac PetscInt tracek = 0; 1408a820b8eSToby Isaac PetscInt tracekm1 = 0; 1418a820b8eSToby Isaac PetscInt fiber0 = 0; 1428a820b8eSToby Isaac PetscInt fiber1 = 0; 1438a820b8eSToby Isaac 1448a820b8eSToby Isaac if (formDegree < dim) { 1455f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Interior(dim - 1, order, formDegree, trimmed, dim > 2 ? 1 : 0, 1, &tracek)); 1465f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Interior(1, order, 0, trimmed, 0, 1, &fiber0)); 1478a820b8eSToby Isaac } 1488a820b8eSToby Isaac if (formDegree > 0) { 1495f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Interior(dim - 1, order, formDegree - 1, trimmed, dim > 2 ? 1 : 0, 1, &tracekm1)); 1505f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Interior(1, order, 1, trimmed, 0, 1, &fiber1)); 1518a820b8eSToby Isaac } 1528a820b8eSToby Isaac *nDofs = (tracek * fiber0 + tracekm1 * fiber1) * nCopies; 1538a820b8eSToby Isaac } 1548a820b8eSToby Isaac } else { /* prism */ 1558a820b8eSToby Isaac PetscInt tracek = 0; 1568a820b8eSToby Isaac PetscInt tracekm1 = 0; 1578a820b8eSToby Isaac PetscInt fiber0 = 0; 1588a820b8eSToby Isaac PetscInt fiber1 = 0; 1598a820b8eSToby Isaac 1608a820b8eSToby Isaac if (formDegree < dim) { 1615f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Interior(dim - 1, order, formDegree, trimmed, 0, 1, &tracek)); 1625f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Interior(1, order, 0, trimmed, 0, 1, &fiber0)); 1638a820b8eSToby Isaac } 1648a820b8eSToby Isaac if (formDegree > 0) { 1655f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Interior(dim - 1, order, formDegree - 1, trimmed, 0, 1, &tracekm1)); 1665f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Interior(1, order, 1, trimmed, 0, 1, &fiber1)); 1678a820b8eSToby Isaac } 1688a820b8eSToby Isaac *nDofs = (tracek * fiber0 + tracekm1 * fiber1) * nCopies; 1698a820b8eSToby Isaac } 1708a820b8eSToby Isaac PetscFunctionReturn(0); 1718a820b8eSToby Isaac } 1728a820b8eSToby Isaac 1738a820b8eSToby Isaac PetscErrorCode testLagrange(PetscHashLag lagTable, DM K, PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscBool continuous, PetscInt nCopies) 1748a820b8eSToby Isaac { 1758a820b8eSToby Isaac PetscDualSpace sp; 1768a820b8eSToby Isaac MPI_Comm comm = PETSC_COMM_SELF; 1778a820b8eSToby Isaac PetscInt Nk; 1788a820b8eSToby Isaac PetscHashLagKey key; 1798a820b8eSToby Isaac PetscHashIter iter; 1808a820b8eSToby Isaac PetscBool missing; 1818a820b8eSToby Isaac PetscInt spdim, spintdim, exspdim, exspintdim; 1828a820b8eSToby Isaac 1838a820b8eSToby Isaac PetscFunctionBegin; 1845f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk)); 1855f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceCreate(comm, &sp)); 1865f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceSetType(sp, PETSCDUALSPACELAGRANGE)); 1875f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceSetDM(sp, K)); 1885f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceSetOrder(sp, order)); 1895f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceSetFormDegree(sp, formDegree)); 1905f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceSetNumComponents(sp, nCopies * Nk)); 1915f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceLagrangeSetContinuity(sp, continuous)); 1925f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceLagrangeSetTensor(sp, (PetscBool) tensor)); 1935f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceLagrangeSetTrimmed(sp, trimmed)); 1945f80ce2aSJacob Faibussowitsch CHKERRQ(PetscInfo(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)); 1955f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Total(dim, order, formDegree, trimmed, tensor, nCopies, &exspdim)); 1968a820b8eSToby Isaac if (continuous && dim > 0 && order > 0) { 1975f80ce2aSJacob Faibussowitsch CHKERRQ(ExpectedNumDofs_Interior(dim, order, formDegree, trimmed, tensor, nCopies, &exspintdim)); 1988a820b8eSToby Isaac } else { 1998a820b8eSToby Isaac exspintdim = exspdim; 2008a820b8eSToby Isaac } 2015f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceSetUp(sp)); 2025f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceGetDimension(sp, &spdim)); 2035f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceGetInteriorDimension(sp, &spintdim)); 2042c71b3e2SJacob Faibussowitsch PetscCheckFalse(spdim != exspdim,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected dual space dimension %D, got %D", exspdim, spdim); 2052c71b3e2SJacob Faibussowitsch PetscCheckFalse(spintdim != exspintdim,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected dual space interior dimension %D, got %D", exspintdim, spintdim); 2068a820b8eSToby Isaac key.dim = dim; 2078a820b8eSToby Isaac key.formDegree = formDegree; 2085f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceGetOrder(sp, &key.order)); 2095f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceLagrangeGetContinuity(sp, &key.continuous)); 2108a820b8eSToby Isaac if (tensor == 2) { 2118a820b8eSToby Isaac key.tensor = 2; 2128a820b8eSToby Isaac } else { 2136ff15688SToby Isaac PetscBool bTensor; 2146ff15688SToby Isaac 2155f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceLagrangeGetTensor(sp, &bTensor)); 2166ff15688SToby Isaac key.tensor = bTensor; 2178a820b8eSToby Isaac } 2185f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceLagrangeGetTrimmed(sp, &key.trimmed)); 2195f80ce2aSJacob Faibussowitsch CHKERRQ(PetscInfo(NULL, "After setup: order %D, trimmed %D, tensor %D, continuous %D\n", key.order, (PetscInt) key.trimmed, key.tensor, (PetscInt) key.continuous)); 2205f80ce2aSJacob Faibussowitsch CHKERRQ(PetscHashLagPut(lagTable, key, &iter, &missing)); 2218a820b8eSToby Isaac if (missing) { 2228a820b8eSToby Isaac DMPolytopeType type; 2238a820b8eSToby Isaac 2245f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetCellType(K, 0, &type)); 2255f80ce2aSJacob Faibussowitsch CHKERRQ(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)); 2265f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF)); 2278a820b8eSToby Isaac { 2288a820b8eSToby Isaac PetscQuadrature intNodes, allNodes; 2298a820b8eSToby Isaac Mat intMat, allMat; 2308a820b8eSToby Isaac MatInfo info; 2318a820b8eSToby Isaac PetscInt i, j, nodeIdxDim, nodeVecDim, nNodes; 2328a820b8eSToby Isaac const PetscInt *nodeIdx; 2338a820b8eSToby Isaac const PetscReal *nodeVec; 2348a820b8eSToby Isaac 2358a820b8eSToby Isaac PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 2368a820b8eSToby Isaac 2375f80ce2aSJacob Faibussowitsch CHKERRQ(PetscLagNodeIndicesGetData_Internal(lag->allNodeIndices, &nodeIdxDim, &nodeVecDim, &nNodes, &nodeIdx, &nodeVec)); 2382c71b3e2SJacob Faibussowitsch PetscCheckFalse(nodeVecDim != Nk,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect nodeVecDim"); 2392c71b3e2SJacob Faibussowitsch PetscCheckFalse(nNodes != spdim,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect nNodes"); 2408a820b8eSToby Isaac 2415f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceGetAllData(sp, &allNodes, &allMat)); 2428a820b8eSToby Isaac 2435f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All nodes:\n")); 2445f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF)); 2455f80ce2aSJacob Faibussowitsch CHKERRQ(PetscQuadratureView(allNodes, PETSC_VIEWER_STDOUT_SELF)); 2465f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF)); 2475f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All node indices:\n")); 2488a820b8eSToby Isaac for (i = 0; i < spdim; i++) { 2495f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "(")); 2508a820b8eSToby Isaac for (j = 0; j < nodeIdxDim; j++) { 2515f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " %D,", nodeIdx[i * nodeIdxDim + j])); 2528a820b8eSToby Isaac } 2535f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "): [")); 2548a820b8eSToby Isaac for (j = 0; j < nodeVecDim; j++) { 2555f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " %g,", (double) nodeVec[i * nodeVecDim + j])); 2568a820b8eSToby Isaac } 2575f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "]\n")); 2588a820b8eSToby Isaac } 2598a820b8eSToby Isaac 2605f80ce2aSJacob Faibussowitsch CHKERRQ(MatGetInfo(allMat, MAT_LOCAL, &info)); 2615f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All matrix: %D nonzeros\n", (PetscInt) info.nz_used)); 2628a820b8eSToby Isaac 2635f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceGetInteriorData(sp, &intNodes, &intMat)); 2648a820b8eSToby Isaac if (intMat && intMat != allMat) { 2658a820b8eSToby Isaac PetscInt intNodeIdxDim, intNodeVecDim, intNnodes; 2668a820b8eSToby Isaac const PetscInt *intNodeIdx; 2678a820b8eSToby Isaac const PetscReal *intNodeVec; 2688a820b8eSToby Isaac PetscBool same; 2698a820b8eSToby Isaac 2705f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior nodes:\n")); 2715f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF)); 2725f80ce2aSJacob Faibussowitsch CHKERRQ(PetscQuadratureView(intNodes, PETSC_VIEWER_STDOUT_SELF)); 2735f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF)); 2748a820b8eSToby Isaac 2755f80ce2aSJacob Faibussowitsch CHKERRQ(MatGetInfo(intMat, MAT_LOCAL, &info)); 2765f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior matrix: %D nonzeros\n", (PetscInt) info.nz_used)); 2775f80ce2aSJacob Faibussowitsch CHKERRQ(PetscLagNodeIndicesGetData_Internal(lag->intNodeIndices, &intNodeIdxDim, &intNodeVecDim, &intNnodes, &intNodeIdx, &intNodeVec)); 2782c71b3e2SJacob Faibussowitsch PetscCheckFalse(intNodeIdxDim != nodeIdxDim || intNodeVecDim != nodeVecDim,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node indices not the same shale as all node indices"); 2792c71b3e2SJacob Faibussowitsch PetscCheckFalse(intNnodes != spintdim,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect interior nNodes"); 2805f80ce2aSJacob Faibussowitsch CHKERRQ(PetscArraycmp(intNodeIdx, nodeIdx, nodeIdxDim * intNnodes, &same)); 281*28b400f6SJacob Faibussowitsch PetscCheck(same,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node indices not the same as start of all node indices"); 2825f80ce2aSJacob Faibussowitsch CHKERRQ(PetscArraycmp(intNodeVec, nodeVec, nodeVecDim * intNnodes, &same)); 283*28b400f6SJacob Faibussowitsch PetscCheck(same,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node vectors not the same as start of all node vectors"); 2848a820b8eSToby Isaac } else if (intMat) { 2855f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior data is the same as all data\n")); 2862c71b3e2SJacob Faibussowitsch PetscCheckFalse(intNodes != allNodes,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior nodes should be the same as all nodes"); 2872c71b3e2SJacob Faibussowitsch PetscCheckFalse(lag->intNodeIndices != lag->allNodeIndices,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node indices should be the same as all node indices"); 2888a820b8eSToby Isaac } 2898a820b8eSToby Isaac } 2908a820b8eSToby Isaac if (dim <= 2 && spintdim) { 291b5a892a1SMatthew G. Knepley PetscInt numFaces, o; 2928a820b8eSToby Isaac 293b5a892a1SMatthew G. Knepley { 294b5a892a1SMatthew G. Knepley DMPolytopeType ct; 295b5a892a1SMatthew G. Knepley /* The number of arrangements is no longer based on the number of faces */ 2965f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexGetCellType(K, 0, &ct)); 297b5a892a1SMatthew G. Knepley numFaces = DMPolytopeTypeGetNumArrangments(ct) / 2; 298b5a892a1SMatthew G. Knepley } 299b5a892a1SMatthew G. Knepley for (o = -numFaces; o < numFaces; ++o) { 3008a820b8eSToby Isaac Mat symMat; 3018a820b8eSToby Isaac 3025f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, o, &symMat)); 3035f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior node symmetry matrix for orientation %D:\n", o)); 3045f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF)); 3055f80ce2aSJacob Faibussowitsch CHKERRQ(MatView(symMat, PETSC_VIEWER_STDOUT_SELF)); 3065f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF)); 3075f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&symMat)); 3088a820b8eSToby Isaac } 3098a820b8eSToby Isaac } 3105f80ce2aSJacob Faibussowitsch CHKERRQ(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF)); 3118a820b8eSToby Isaac } 3125f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDualSpaceDestroy(&sp)); 3138a820b8eSToby Isaac PetscFunctionReturn(0); 3148a820b8eSToby Isaac } 3158a820b8eSToby Isaac 3168a820b8eSToby Isaac int main (int argc, char **argv) 3178a820b8eSToby Isaac { 3188a820b8eSToby Isaac PetscInt dim; 3198a820b8eSToby Isaac PetscHashLag lagTable; 3208a820b8eSToby Isaac PetscInt tensorCell; 3218a820b8eSToby Isaac PetscInt order, ordermin, ordermax; 3228a820b8eSToby Isaac PetscBool continuous; 3238a820b8eSToby Isaac PetscBool trimmed; 3248a820b8eSToby Isaac DM dm; 3258a820b8eSToby Isaac PetscErrorCode ierr; 3268a820b8eSToby Isaac 3278a820b8eSToby Isaac ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 3288a820b8eSToby Isaac dim = 3; 3298a820b8eSToby Isaac tensorCell = 0; 3308a820b8eSToby Isaac continuous = PETSC_FALSE; 3318a820b8eSToby Isaac trimmed = PETSC_FALSE; 3328a820b8eSToby Isaac ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Options for PETSCDUALSPACELAGRANGE test","none");CHKERRQ(ierr); 3335f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsRangeInt("-dim", "The spatial dimension","ex1.c",dim,&dim,NULL,0,3)); 3345f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsRangeInt("-tensor", "(0) simplex (1) hypercube (2) wedge","ex1.c",tensorCell,&tensorCell,NULL,0,2)); 3355f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBool("-continuous", "Whether the dual space has continuity","ex1.c",continuous,&continuous,NULL)); 3365f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsBool("-trimmed", "Whether the dual space matches a trimmed polynomial space","ex1.c",trimmed,&trimmed,NULL)); 3371e1ea65dSPierre Jolivet ierr = PetscOptionsEnd();CHKERRQ(ierr); 3385f80ce2aSJacob Faibussowitsch CHKERRQ(PetscHashLagCreate(&lagTable)); 3398a820b8eSToby Isaac 3408a820b8eSToby Isaac if (tensorCell < 2) { 3415f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, (PetscBool) !tensorCell), &dm)); 3428a820b8eSToby Isaac } else { 3435f80ce2aSJacob Faibussowitsch CHKERRQ(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DM_POLYTOPE_TRI_PRISM, &dm)); 3448a820b8eSToby Isaac } 3458a820b8eSToby Isaac ordermin = trimmed ? 1 : 0; 3468a820b8eSToby Isaac ordermax = tensorCell == 2 ? 4 : tensorCell == 1 ? 3 : dim + 2; 3478a820b8eSToby Isaac for (order = ordermin; order <= ordermax; order++) { 3488a820b8eSToby Isaac PetscInt formDegree; 3498a820b8eSToby Isaac 3508a820b8eSToby Isaac for (formDegree = PetscMin(0,-dim+1); formDegree <= dim; formDegree++) { 3518a820b8eSToby Isaac PetscInt nCopies; 3528a820b8eSToby Isaac 3538a820b8eSToby Isaac for (nCopies = 1; nCopies <= 3; nCopies++) { 3545f80ce2aSJacob Faibussowitsch CHKERRQ(testLagrange(lagTable, dm, dim, order, formDegree, trimmed, (PetscBool) tensorCell, continuous, nCopies)); 3558a820b8eSToby Isaac } 3568a820b8eSToby Isaac } 3578a820b8eSToby Isaac } 3585f80ce2aSJacob Faibussowitsch CHKERRQ(DMDestroy(&dm)); 3595f80ce2aSJacob Faibussowitsch CHKERRQ(PetscHashLagDestroy(&lagTable)); 3608a820b8eSToby Isaac ierr = PetscFinalize(); 3618a820b8eSToby Isaac return ierr; 3628a820b8eSToby Isaac } 3638a820b8eSToby Isaac 3648a820b8eSToby Isaac /*TEST 3658a820b8eSToby Isaac 3668a820b8eSToby Isaac test: 3678a820b8eSToby Isaac suffix: 0 3688a820b8eSToby Isaac args: -dim 0 3698a820b8eSToby Isaac 3708a820b8eSToby Isaac test: 3718a820b8eSToby Isaac suffix: 1_discontinuous_full 3728a820b8eSToby Isaac args: -dim 1 -continuous 0 -trimmed 0 3738a820b8eSToby Isaac 3748a820b8eSToby Isaac test: 3758a820b8eSToby Isaac suffix: 1_continuous_full 3768a820b8eSToby Isaac args: -dim 1 -continuous 1 -trimmed 0 3778a820b8eSToby Isaac 3788a820b8eSToby Isaac test: 3798a820b8eSToby Isaac suffix: 2_simplex_discontinuous_full 3808a820b8eSToby Isaac args: -dim 2 -tensor 0 -continuous 0 -trimmed 0 3818a820b8eSToby Isaac 3828a820b8eSToby Isaac test: 3838a820b8eSToby Isaac suffix: 2_simplex_continuous_full 3848a820b8eSToby Isaac args: -dim 2 -tensor 0 -continuous 1 -trimmed 0 3858a820b8eSToby Isaac 3868a820b8eSToby Isaac test: 3878a820b8eSToby Isaac suffix: 2_tensor_discontinuous_full 3888a820b8eSToby Isaac args: -dim 2 -tensor 1 -continuous 0 -trimmed 0 3898a820b8eSToby Isaac 3908a820b8eSToby Isaac test: 3918a820b8eSToby Isaac suffix: 2_tensor_continuous_full 3928a820b8eSToby Isaac args: -dim 2 -tensor 1 -continuous 1 -trimmed 0 3938a820b8eSToby Isaac 3948a820b8eSToby Isaac test: 3958a820b8eSToby Isaac suffix: 3_simplex_discontinuous_full 3968a820b8eSToby Isaac args: -dim 3 -tensor 0 -continuous 0 -trimmed 0 3978a820b8eSToby Isaac 3988a820b8eSToby Isaac test: 3998a820b8eSToby Isaac suffix: 3_simplex_continuous_full 4008a820b8eSToby Isaac args: -dim 3 -tensor 0 -continuous 1 -trimmed 0 4018a820b8eSToby Isaac 4028a820b8eSToby Isaac test: 4038a820b8eSToby Isaac suffix: 3_tensor_discontinuous_full 4048a820b8eSToby Isaac args: -dim 3 -tensor 1 -continuous 0 -trimmed 0 4058a820b8eSToby Isaac 4068a820b8eSToby Isaac test: 4078a820b8eSToby Isaac suffix: 3_tensor_continuous_full 4088a820b8eSToby Isaac args: -dim 3 -tensor 1 -continuous 1 -trimmed 0 4098a820b8eSToby Isaac 4108a820b8eSToby Isaac test: 4118a820b8eSToby Isaac suffix: 3_wedge_discontinuous_full 4128a820b8eSToby Isaac args: -dim 3 -tensor 2 -continuous 0 -trimmed 0 4138a820b8eSToby Isaac 4148a820b8eSToby Isaac test: 4158a820b8eSToby Isaac suffix: 3_wedge_continuous_full 4168a820b8eSToby Isaac args: -dim 3 -tensor 2 -continuous 1 -trimmed 0 4178a820b8eSToby Isaac 4188a820b8eSToby Isaac test: 4198a820b8eSToby Isaac suffix: 1_discontinuous_trimmed 4208a820b8eSToby Isaac args: -dim 1 -continuous 0 -trimmed 1 4218a820b8eSToby Isaac 4228a820b8eSToby Isaac test: 4238a820b8eSToby Isaac suffix: 1_continuous_trimmed 4248a820b8eSToby Isaac args: -dim 1 -continuous 1 -trimmed 1 4258a820b8eSToby Isaac 4268a820b8eSToby Isaac test: 4278a820b8eSToby Isaac suffix: 2_simplex_discontinuous_trimmed 4288a820b8eSToby Isaac args: -dim 2 -tensor 0 -continuous 0 -trimmed 1 4298a820b8eSToby Isaac 4308a820b8eSToby Isaac test: 4318a820b8eSToby Isaac suffix: 2_simplex_continuous_trimmed 4328a820b8eSToby Isaac args: -dim 2 -tensor 0 -continuous 1 -trimmed 1 4338a820b8eSToby Isaac 4348a820b8eSToby Isaac test: 4358a820b8eSToby Isaac suffix: 2_tensor_discontinuous_trimmed 4368a820b8eSToby Isaac args: -dim 2 -tensor 1 -continuous 0 -trimmed 1 4378a820b8eSToby Isaac 4388a820b8eSToby Isaac test: 4398a820b8eSToby Isaac suffix: 2_tensor_continuous_trimmed 4408a820b8eSToby Isaac args: -dim 2 -tensor 1 -continuous 1 -trimmed 1 4418a820b8eSToby Isaac 4428a820b8eSToby Isaac test: 4438a820b8eSToby Isaac suffix: 3_simplex_discontinuous_trimmed 4448a820b8eSToby Isaac args: -dim 3 -tensor 0 -continuous 0 -trimmed 1 4458a820b8eSToby Isaac 4468a820b8eSToby Isaac test: 4478a820b8eSToby Isaac suffix: 3_simplex_continuous_trimmed 4488a820b8eSToby Isaac args: -dim 3 -tensor 0 -continuous 1 -trimmed 1 4498a820b8eSToby Isaac 4508a820b8eSToby Isaac test: 4518a820b8eSToby Isaac suffix: 3_tensor_discontinuous_trimmed 4528a820b8eSToby Isaac args: -dim 3 -tensor 1 -continuous 0 -trimmed 1 4538a820b8eSToby Isaac 4548a820b8eSToby Isaac test: 4558a820b8eSToby Isaac suffix: 3_tensor_continuous_trimmed 4568a820b8eSToby Isaac args: -dim 3 -tensor 1 -continuous 1 -trimmed 1 4578a820b8eSToby Isaac 4588a820b8eSToby Isaac test: 4598a820b8eSToby Isaac suffix: 3_wedge_discontinuous_trimmed 4608a820b8eSToby Isaac args: -dim 3 -tensor 2 -continuous 0 -trimmed 1 4618a820b8eSToby Isaac 4628a820b8eSToby Isaac test: 4638a820b8eSToby Isaac suffix: 3_wedge_continuous_trimmed 4648a820b8eSToby Isaac args: -dim 3 -tensor 2 -continuous 1 -trimmed 1 4658a820b8eSToby Isaac 4668a820b8eSToby Isaac TEST*/ 467