18a820b8eSToby Isaac #include <petscfe.h> 28a820b8eSToby Isaac #include <petscdmplex.h> 38a820b8eSToby Isaac #include <petsc/private/hashmap.h> 48a820b8eSToby Isaac #include <petsc/private/dmpleximpl.h> 58a820b8eSToby Isaac #include <petsc/private/petscfeimpl.h> 68a820b8eSToby Isaac 78a820b8eSToby Isaac const char help[] = "Test PETSCDUALSPACELAGRANGE\n"; 88a820b8eSToby Isaac 99371c9d4SSatish Balay typedef struct _PetscHashLagKey { 108a820b8eSToby Isaac PetscInt dim; 118a820b8eSToby Isaac PetscInt order; 128a820b8eSToby Isaac PetscInt formDegree; 138a820b8eSToby Isaac PetscBool trimmed; 146ff15688SToby Isaac PetscInt tensor; 158a820b8eSToby Isaac PetscBool continuous; 168a820b8eSToby Isaac } PetscHashLagKey; 178a820b8eSToby Isaac 188a820b8eSToby Isaac #define PetscHashLagKeyHash(key) \ 199371c9d4SSatish Balay PetscHashCombine(PetscHashCombine(PetscHashCombine(PetscHashInt((key).dim), PetscHashInt((key).order)), PetscHashInt((key).formDegree)), PetscHashCombine(PetscHashCombine(PetscHashInt((key).trimmed), PetscHashInt((key).tensor)), PetscHashInt((key).continuous))) 208a820b8eSToby Isaac 218a820b8eSToby Isaac #define PetscHashLagKeyEqual(k1, k2) \ 229371c9d4SSatish Balay (((k1).dim == (k2).dim) ? ((k1).order == (k2).order) ? ((k1).formDegree == (k2).formDegree) ? ((k1).trimmed == (k2).trimmed) ? ((k1).tensor == (k2).tensor) ? ((k1).continuous == (k2).continuous) : 0 : 0 : 0 : 0 : 0) 238a820b8eSToby Isaac 248a820b8eSToby Isaac PETSC_HASH_MAP(HashLag, PetscHashLagKey, PetscInt, PetscHashLagKeyHash, PetscHashLagKeyEqual, 0) 258a820b8eSToby Isaac 268a820b8eSToby Isaac static PetscErrorCode ExpectedNumDofs_Total(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs); 278a820b8eSToby Isaac static PetscErrorCode ExpectedNumDofs_Interior(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs); 288a820b8eSToby Isaac 29d71ae5a4SJacob Faibussowitsch static PetscErrorCode ExpectedNumDofs_Total(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs) 30d71ae5a4SJacob Faibussowitsch { 318a820b8eSToby Isaac PetscFunctionBegin; 328a820b8eSToby Isaac formDegree = PetscAbsInt(formDegree); 338a820b8eSToby Isaac /* see femtable.org for the source of most of these values */ 348a820b8eSToby Isaac *nDofs = -1; 358a820b8eSToby Isaac if (tensor == 0) { /* simplex spaces */ 368a820b8eSToby Isaac if (trimmed) { 378a820b8eSToby Isaac PetscInt rnchooserk; 388a820b8eSToby Isaac PetscInt rkm1choosek; 398a820b8eSToby Isaac 409566063dSJacob Faibussowitsch PetscCall(PetscDTBinomialInt(order + dim, order + formDegree, &rnchooserk)); 419566063dSJacob Faibussowitsch PetscCall(PetscDTBinomialInt(order + formDegree - 1, formDegree, &rkm1choosek)); 428a820b8eSToby Isaac *nDofs = rnchooserk * rkm1choosek * nCopies; 438a820b8eSToby Isaac } else { 448a820b8eSToby Isaac PetscInt rnchooserk; 458a820b8eSToby Isaac PetscInt rkchoosek; 468a820b8eSToby Isaac 479566063dSJacob Faibussowitsch PetscCall(PetscDTBinomialInt(order + dim, order + formDegree, &rnchooserk)); 489566063dSJacob Faibussowitsch PetscCall(PetscDTBinomialInt(order + formDegree, formDegree, &rkchoosek)); 498a820b8eSToby Isaac *nDofs = rnchooserk * rkchoosek * nCopies; 508a820b8eSToby Isaac } 518a820b8eSToby Isaac } else if (tensor == 1) { /* hypercubes */ 528a820b8eSToby Isaac if (trimmed) { 538a820b8eSToby Isaac PetscInt nchoosek; 548a820b8eSToby Isaac PetscInt rpowk, rp1pownmk; 558a820b8eSToby Isaac 569566063dSJacob Faibussowitsch PetscCall(PetscDTBinomialInt(dim, formDegree, &nchoosek)); 572f613bf5SBarry Smith rpowk = PetscPowInt(order, formDegree); 582f613bf5SBarry Smith rp1pownmk = PetscPowInt(order + 1, dim - formDegree); 598a820b8eSToby Isaac *nDofs = nchoosek * rpowk * rp1pownmk * nCopies; 608a820b8eSToby Isaac } else { 618a820b8eSToby Isaac PetscInt nchoosek; 628a820b8eSToby Isaac PetscInt rp1pown; 638a820b8eSToby Isaac 649566063dSJacob Faibussowitsch PetscCall(PetscDTBinomialInt(dim, formDegree, &nchoosek)); 652f613bf5SBarry Smith rp1pown = PetscPowInt(order + 1, dim); 668a820b8eSToby Isaac *nDofs = nchoosek * rp1pown * nCopies; 678a820b8eSToby Isaac } 688a820b8eSToby Isaac } else { /* prism */ 698a820b8eSToby Isaac PetscInt tracek = 0; 708a820b8eSToby Isaac PetscInt tracekm1 = 0; 718a820b8eSToby Isaac PetscInt fiber0 = 0; 728a820b8eSToby Isaac PetscInt fiber1 = 0; 738a820b8eSToby Isaac 748a820b8eSToby Isaac if (formDegree < dim) { 759566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Total(dim - 1, order, formDegree, trimmed, 0, 1, &tracek)); 769566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Total(1, order, 0, trimmed, 0, 1, &fiber0)); 778a820b8eSToby Isaac } 788a820b8eSToby Isaac if (formDegree > 0) { 799566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Total(dim - 1, order, formDegree - 1, trimmed, 0, 1, &tracekm1)); 809566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Total(1, order, 1, trimmed, 0, 1, &fiber1)); 818a820b8eSToby Isaac } 828a820b8eSToby Isaac *nDofs = (tracek * fiber0 + tracekm1 * fiber1) * nCopies; 838a820b8eSToby Isaac } 843ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 858a820b8eSToby Isaac } 868a820b8eSToby Isaac 87d71ae5a4SJacob Faibussowitsch static PetscErrorCode ExpectedNumDofs_Interior(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs) 88d71ae5a4SJacob Faibussowitsch { 898a820b8eSToby Isaac PetscFunctionBegin; 908a820b8eSToby Isaac formDegree = PetscAbsInt(formDegree); 918a820b8eSToby Isaac /* see femtable.org for the source of most of these values */ 928a820b8eSToby Isaac *nDofs = -1; 938a820b8eSToby Isaac if (tensor == 0) { /* simplex spaces */ 948a820b8eSToby Isaac if (trimmed) { 958a820b8eSToby Isaac if (order + formDegree > dim) { 968a820b8eSToby Isaac PetscInt eorder = order + formDegree - dim - 1; 978a820b8eSToby Isaac PetscInt eformDegree = dim - formDegree; 988a820b8eSToby Isaac PetscInt rnchooserk; 998a820b8eSToby Isaac PetscInt rkchoosek; 1008a820b8eSToby Isaac 1019566063dSJacob Faibussowitsch PetscCall(PetscDTBinomialInt(eorder + dim, eorder + eformDegree, &rnchooserk)); 1029566063dSJacob Faibussowitsch PetscCall(PetscDTBinomialInt(eorder + eformDegree, eformDegree, &rkchoosek)); 1038a820b8eSToby Isaac *nDofs = rnchooserk * rkchoosek * nCopies; 1048a820b8eSToby Isaac } else { 1058a820b8eSToby Isaac *nDofs = 0; 1068a820b8eSToby Isaac } 1078a820b8eSToby Isaac 1088a820b8eSToby Isaac } else { 1098a820b8eSToby Isaac if (order + formDegree > dim) { 1108a820b8eSToby Isaac PetscInt eorder = order + formDegree - dim; 1118a820b8eSToby Isaac PetscInt eformDegree = dim - formDegree; 1128a820b8eSToby Isaac PetscInt rnchooserk; 1138a820b8eSToby Isaac PetscInt rkm1choosek; 1148a820b8eSToby Isaac 1159566063dSJacob Faibussowitsch PetscCall(PetscDTBinomialInt(eorder + dim, eorder + eformDegree, &rnchooserk)); 1169566063dSJacob Faibussowitsch PetscCall(PetscDTBinomialInt(eorder + eformDegree - 1, eformDegree, &rkm1choosek)); 1178a820b8eSToby Isaac *nDofs = rnchooserk * rkm1choosek * nCopies; 1188a820b8eSToby Isaac } else { 1198a820b8eSToby Isaac *nDofs = 0; 1208a820b8eSToby Isaac } 1218a820b8eSToby Isaac } 1228a820b8eSToby Isaac } else if (tensor == 1) { /* hypercubes */ 1238a820b8eSToby Isaac if (dim < 2) { 1249566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Interior(dim, order, formDegree, trimmed, 0, nCopies, nDofs)); 1258a820b8eSToby Isaac } else { 1268a820b8eSToby Isaac PetscInt tracek = 0; 1278a820b8eSToby Isaac PetscInt tracekm1 = 0; 1288a820b8eSToby Isaac PetscInt fiber0 = 0; 1298a820b8eSToby Isaac PetscInt fiber1 = 0; 1308a820b8eSToby Isaac 1318a820b8eSToby Isaac if (formDegree < dim) { 1329566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Interior(dim - 1, order, formDegree, trimmed, dim > 2 ? 1 : 0, 1, &tracek)); 1339566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Interior(1, order, 0, trimmed, 0, 1, &fiber0)); 1348a820b8eSToby Isaac } 1358a820b8eSToby Isaac if (formDegree > 0) { 1369566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Interior(dim - 1, order, formDegree - 1, trimmed, dim > 2 ? 1 : 0, 1, &tracekm1)); 1379566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Interior(1, order, 1, trimmed, 0, 1, &fiber1)); 1388a820b8eSToby Isaac } 1398a820b8eSToby Isaac *nDofs = (tracek * fiber0 + tracekm1 * fiber1) * nCopies; 1408a820b8eSToby Isaac } 1418a820b8eSToby Isaac } else { /* prism */ 1428a820b8eSToby Isaac PetscInt tracek = 0; 1438a820b8eSToby Isaac PetscInt tracekm1 = 0; 1448a820b8eSToby Isaac PetscInt fiber0 = 0; 1458a820b8eSToby Isaac PetscInt fiber1 = 0; 1468a820b8eSToby Isaac 1478a820b8eSToby Isaac if (formDegree < dim) { 1489566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Interior(dim - 1, order, formDegree, trimmed, 0, 1, &tracek)); 1499566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Interior(1, order, 0, trimmed, 0, 1, &fiber0)); 1508a820b8eSToby Isaac } 1518a820b8eSToby Isaac if (formDegree > 0) { 1529566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Interior(dim - 1, order, formDegree - 1, trimmed, 0, 1, &tracekm1)); 1539566063dSJacob Faibussowitsch PetscCall(ExpectedNumDofs_Interior(1, order, 1, trimmed, 0, 1, &fiber1)); 1548a820b8eSToby Isaac } 1558a820b8eSToby Isaac *nDofs = (tracek * fiber0 + tracekm1 * fiber1) * nCopies; 1568a820b8eSToby Isaac } 1573ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1588a820b8eSToby Isaac } 1598a820b8eSToby Isaac 160*b8b5be36SMartin Diehl PetscErrorCode testLagrange(PetscHashLag lagTable, DM K, PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensorCell, PetscBool continuous, PetscInt nCopies) 161d71ae5a4SJacob Faibussowitsch { 1628a820b8eSToby Isaac PetscDualSpace sp; 1638a820b8eSToby Isaac MPI_Comm comm = PETSC_COMM_SELF; 1648a820b8eSToby Isaac PetscInt Nk; 1658a820b8eSToby Isaac PetscHashLagKey key; 1668a820b8eSToby Isaac PetscHashIter iter; 1678a820b8eSToby Isaac PetscBool missing; 1688a820b8eSToby Isaac PetscInt spdim, spintdim, exspdim, exspintdim; 1698a820b8eSToby Isaac 1708a820b8eSToby Isaac PetscFunctionBegin; 1719566063dSJacob Faibussowitsch PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk)); 1729566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceCreate(comm, &sp)); 1739566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceSetType(sp, PETSCDUALSPACELAGRANGE)); 1749566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceSetDM(sp, K)); 1759566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceSetOrder(sp, order)); 1769566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceSetFormDegree(sp, formDegree)); 1779566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceSetNumComponents(sp, nCopies * Nk)); 1789566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceLagrangeSetContinuity(sp, continuous)); 179*b8b5be36SMartin Diehl PetscCall(PetscDualSpaceLagrangeSetTensor(sp, (PetscBool)(tensorCell > 0))); 1809566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceLagrangeSetTrimmed(sp, trimmed)); 181*b8b5be36SMartin Diehl PetscCall(PetscInfo(NULL, "Input: dim %" PetscInt_FMT ", order %" PetscInt_FMT ", trimmed %" PetscInt_FMT ", tensorCell %" PetscInt_FMT ", continuous %" PetscInt_FMT ", formDegree %" PetscInt_FMT ", nCopies %" PetscInt_FMT "\n", dim, order, (PetscInt)trimmed, tensorCell, (PetscInt)continuous, formDegree, nCopies)); 182*b8b5be36SMartin Diehl PetscCall(ExpectedNumDofs_Total(dim, order, formDegree, trimmed, tensorCell, nCopies, &exspdim)); 1838a820b8eSToby Isaac if (continuous && dim > 0 && order > 0) { 184*b8b5be36SMartin Diehl PetscCall(ExpectedNumDofs_Interior(dim, order, formDegree, trimmed, tensorCell, nCopies, &exspintdim)); 1858a820b8eSToby Isaac } else { 1868a820b8eSToby Isaac exspintdim = exspdim; 1878a820b8eSToby Isaac } 1889566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceSetUp(sp)); 1899566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceGetDimension(sp, &spdim)); 1909566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceGetInteriorDimension(sp, &spintdim)); 19163a3b9bcSJacob Faibussowitsch PetscCheck(spdim == exspdim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected dual space dimension %" PetscInt_FMT ", got %" PetscInt_FMT, exspdim, spdim); 19263a3b9bcSJacob Faibussowitsch PetscCheck(spintdim == exspintdim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected dual space interior dimension %" PetscInt_FMT ", got %" PetscInt_FMT, exspintdim, spintdim); 1938a820b8eSToby Isaac key.dim = dim; 1948a820b8eSToby Isaac key.formDegree = formDegree; 1959566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceGetOrder(sp, &key.order)); 1969566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceLagrangeGetContinuity(sp, &key.continuous)); 197*b8b5be36SMartin Diehl if (tensorCell == 2) { 1988a820b8eSToby Isaac key.tensor = 2; 1998a820b8eSToby Isaac } else { 2006ff15688SToby Isaac PetscBool bTensor; 2016ff15688SToby Isaac 2029566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceLagrangeGetTensor(sp, &bTensor)); 2036ff15688SToby Isaac key.tensor = bTensor; 2048a820b8eSToby Isaac } 2059566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceLagrangeGetTrimmed(sp, &key.trimmed)); 20663a3b9bcSJacob Faibussowitsch PetscCall(PetscInfo(NULL, "After setup: order %" PetscInt_FMT ", trimmed %" PetscInt_FMT ", tensor %" PetscInt_FMT ", continuous %" PetscInt_FMT "\n", key.order, (PetscInt)key.trimmed, key.tensor, (PetscInt)key.continuous)); 2079566063dSJacob Faibussowitsch PetscCall(PetscHashLagPut(lagTable, key, &iter, &missing)); 2088a820b8eSToby Isaac if (missing) { 2098a820b8eSToby Isaac DMPolytopeType type; 2108a820b8eSToby Isaac 2119566063dSJacob Faibussowitsch PetscCall(DMPlexGetCellType(K, 0, &type)); 212*b8b5be36SMartin Diehl PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "New space: %s, order %" PetscInt_FMT ", trimmed %" PetscInt_FMT ", tensor %" PetscInt_FMT ", continuous %" PetscInt_FMT ", form degree %" PetscInt_FMT "\n", DMPolytopeTypes[type], order, (PetscInt)trimmed, tensorCell, (PetscInt)continuous, formDegree)); 2139566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF)); 2148a820b8eSToby Isaac { 2158a820b8eSToby Isaac PetscQuadrature intNodes, allNodes; 2168a820b8eSToby Isaac Mat intMat, allMat; 2178a820b8eSToby Isaac MatInfo info; 2188a820b8eSToby Isaac PetscInt i, j, nodeIdxDim, nodeVecDim, nNodes; 2198a820b8eSToby Isaac const PetscInt *nodeIdx; 2208a820b8eSToby Isaac const PetscReal *nodeVec; 2218a820b8eSToby Isaac 2228a820b8eSToby Isaac PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 2238a820b8eSToby Isaac 2249566063dSJacob Faibussowitsch PetscCall(PetscLagNodeIndicesGetData_Internal(lag->allNodeIndices, &nodeIdxDim, &nodeVecDim, &nNodes, &nodeIdx, &nodeVec)); 22508401ef6SPierre Jolivet PetscCheck(nodeVecDim == Nk, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect nodeVecDim"); 22608401ef6SPierre Jolivet PetscCheck(nNodes == spdim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect nNodes"); 2278a820b8eSToby Isaac 2289566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceGetAllData(sp, &allNodes, &allMat)); 2298a820b8eSToby Isaac 2309566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All nodes:\n")); 2319566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF)); 2329566063dSJacob Faibussowitsch PetscCall(PetscQuadratureView(allNodes, PETSC_VIEWER_STDOUT_SELF)); 2339566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF)); 2349566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All node indices:\n")); 2358a820b8eSToby Isaac for (i = 0; i < spdim; i++) { 2369566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_SELF, "(")); 23748a46eb9SPierre Jolivet for (j = 0; j < nodeIdxDim; j++) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %" PetscInt_FMT ",", nodeIdx[i * nodeIdxDim + j])); 2389566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_SELF, "): [")); 23948a46eb9SPierre Jolivet for (j = 0; j < nodeVecDim; j++) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g,", (double)nodeVec[i * nodeVecDim + j])); 2409566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_SELF, "]\n")); 2418a820b8eSToby Isaac } 2428a820b8eSToby Isaac 2439566063dSJacob Faibussowitsch PetscCall(MatGetInfo(allMat, MAT_LOCAL, &info)); 24463a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All matrix: %" PetscInt_FMT " nonzeros\n", (PetscInt)info.nz_used)); 2458a820b8eSToby Isaac 2469566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceGetInteriorData(sp, &intNodes, &intMat)); 2478a820b8eSToby Isaac if (intMat && intMat != allMat) { 2488a820b8eSToby Isaac PetscInt intNodeIdxDim, intNodeVecDim, intNnodes; 2498a820b8eSToby Isaac const PetscInt *intNodeIdx; 2508a820b8eSToby Isaac const PetscReal *intNodeVec; 2518a820b8eSToby Isaac PetscBool same; 2528a820b8eSToby Isaac 2539566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior nodes:\n")); 2549566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF)); 2559566063dSJacob Faibussowitsch PetscCall(PetscQuadratureView(intNodes, PETSC_VIEWER_STDOUT_SELF)); 2569566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF)); 2578a820b8eSToby Isaac 2589566063dSJacob Faibussowitsch PetscCall(MatGetInfo(intMat, MAT_LOCAL, &info)); 25963a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior matrix: %" PetscInt_FMT " nonzeros\n", (PetscInt)info.nz_used)); 2609566063dSJacob Faibussowitsch PetscCall(PetscLagNodeIndicesGetData_Internal(lag->intNodeIndices, &intNodeIdxDim, &intNodeVecDim, &intNnodes, &intNodeIdx, &intNodeVec)); 2611dca8a05SBarry Smith PetscCheck(intNodeIdxDim == nodeIdxDim && intNodeVecDim == nodeVecDim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node indices not the same shale as all node indices"); 26208401ef6SPierre Jolivet PetscCheck(intNnodes == spintdim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect interior nNodes"); 2639566063dSJacob Faibussowitsch PetscCall(PetscArraycmp(intNodeIdx, nodeIdx, nodeIdxDim * intNnodes, &same)); 26428b400f6SJacob Faibussowitsch PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node indices not the same as start of all node indices"); 2659566063dSJacob Faibussowitsch PetscCall(PetscArraycmp(intNodeVec, nodeVec, nodeVecDim * intNnodes, &same)); 26628b400f6SJacob Faibussowitsch PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node vectors not the same as start of all node vectors"); 2678a820b8eSToby Isaac } else if (intMat) { 2689566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior data is the same as all data\n")); 26908401ef6SPierre Jolivet PetscCheck(intNodes == allNodes, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior nodes should be the same as all nodes"); 27008401ef6SPierre Jolivet PetscCheck(lag->intNodeIndices == lag->allNodeIndices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Interior node indices should be the same as all node indices"); 2718a820b8eSToby Isaac } 2728a820b8eSToby Isaac } 2738a820b8eSToby Isaac if (dim <= 2 && spintdim) { 274b5a892a1SMatthew G. Knepley PetscInt numFaces, o; 2758a820b8eSToby Isaac 276b5a892a1SMatthew G. Knepley { 277b5a892a1SMatthew G. Knepley DMPolytopeType ct; 278b5a892a1SMatthew G. Knepley /* The number of arrangements is no longer based on the number of faces */ 2799566063dSJacob Faibussowitsch PetscCall(DMPlexGetCellType(K, 0, &ct)); 28085036b15SMatthew G. Knepley numFaces = DMPolytopeTypeGetNumArrangements(ct) / 2; 281b5a892a1SMatthew G. Knepley } 282b5a892a1SMatthew G. Knepley for (o = -numFaces; o < numFaces; ++o) { 2838a820b8eSToby Isaac Mat symMat; 2848a820b8eSToby Isaac 2859566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, o, &symMat)); 28663a3b9bcSJacob Faibussowitsch PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior node symmetry matrix for orientation %" PetscInt_FMT ":\n", o)); 2879566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF)); 2889566063dSJacob Faibussowitsch PetscCall(MatView(symMat, PETSC_VIEWER_STDOUT_SELF)); 2899566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF)); 2909566063dSJacob Faibussowitsch PetscCall(MatDestroy(&symMat)); 2918a820b8eSToby Isaac } 2928a820b8eSToby Isaac } 2939566063dSJacob Faibussowitsch PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF)); 2948a820b8eSToby Isaac } 2959566063dSJacob Faibussowitsch PetscCall(PetscDualSpaceDestroy(&sp)); 2963ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2978a820b8eSToby Isaac } 2988a820b8eSToby Isaac 299d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv) 300d71ae5a4SJacob Faibussowitsch { 3018a820b8eSToby Isaac PetscInt dim; 3028a820b8eSToby Isaac PetscHashLag lagTable; 3038a820b8eSToby Isaac PetscInt tensorCell; 3048a820b8eSToby Isaac PetscInt order, ordermin, ordermax; 3058a820b8eSToby Isaac PetscBool continuous; 3068a820b8eSToby Isaac PetscBool trimmed; 3078a820b8eSToby Isaac DM dm; 3088a820b8eSToby Isaac 309327415f7SBarry Smith PetscFunctionBeginUser; 3109566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 3118a820b8eSToby Isaac dim = 3; 3128a820b8eSToby Isaac tensorCell = 0; 3138a820b8eSToby Isaac continuous = PETSC_FALSE; 3148a820b8eSToby Isaac trimmed = PETSC_FALSE; 315d0609cedSBarry Smith PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for PETSCDUALSPACELAGRANGE test", "none"); 3169566063dSJacob Faibussowitsch PetscCall(PetscOptionsRangeInt("-dim", "The spatial dimension", "ex1.c", dim, &dim, NULL, 0, 3)); 3179566063dSJacob Faibussowitsch PetscCall(PetscOptionsRangeInt("-tensor", "(0) simplex (1) hypercube (2) wedge", "ex1.c", tensorCell, &tensorCell, NULL, 0, 2)); 3189566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-continuous", "Whether the dual space has continuity", "ex1.c", continuous, &continuous, NULL)); 3199566063dSJacob Faibussowitsch PetscCall(PetscOptionsBool("-trimmed", "Whether the dual space matches a trimmed polynomial space", "ex1.c", trimmed, &trimmed, NULL)); 320d0609cedSBarry Smith PetscOptionsEnd(); 3219566063dSJacob Faibussowitsch PetscCall(PetscHashLagCreate(&lagTable)); 3228a820b8eSToby Isaac 3238a820b8eSToby Isaac if (tensorCell < 2) { 324*b8b5be36SMartin Diehl PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, (PetscBool)(tensorCell == 0)), &dm)); 3258a820b8eSToby Isaac } else { 3269566063dSJacob Faibussowitsch PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DM_POLYTOPE_TRI_PRISM, &dm)); 3278a820b8eSToby Isaac } 3288a820b8eSToby Isaac ordermin = trimmed ? 1 : 0; 3298a820b8eSToby Isaac ordermax = tensorCell == 2 ? 4 : tensorCell == 1 ? 3 : dim + 2; 3308a820b8eSToby Isaac for (order = ordermin; order <= ordermax; order++) { 3318a820b8eSToby Isaac PetscInt formDegree; 3328a820b8eSToby Isaac 3338a820b8eSToby Isaac for (formDegree = PetscMin(0, -dim + 1); formDegree <= dim; formDegree++) { 3348a820b8eSToby Isaac PetscInt nCopies; 3358a820b8eSToby Isaac 336*b8b5be36SMartin Diehl for (nCopies = 1; nCopies <= 3; nCopies++) PetscCall(testLagrange(lagTable, dm, dim, order, formDegree, trimmed, tensorCell, continuous, nCopies)); 3378a820b8eSToby Isaac } 3388a820b8eSToby Isaac } 3399566063dSJacob Faibussowitsch PetscCall(DMDestroy(&dm)); 3409566063dSJacob Faibussowitsch PetscCall(PetscHashLagDestroy(&lagTable)); 3419566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 342b122ec5aSJacob Faibussowitsch return 0; 3438a820b8eSToby Isaac } 3448a820b8eSToby Isaac 3458a820b8eSToby Isaac /*TEST 3468a820b8eSToby Isaac 3478a820b8eSToby Isaac test: 3488a820b8eSToby Isaac suffix: 0 3498a820b8eSToby Isaac args: -dim 0 3508a820b8eSToby Isaac 3518a820b8eSToby Isaac test: 3528a820b8eSToby Isaac suffix: 1_discontinuous_full 3538a820b8eSToby Isaac args: -dim 1 -continuous 0 -trimmed 0 3548a820b8eSToby Isaac 3558a820b8eSToby Isaac test: 3568a820b8eSToby Isaac suffix: 1_continuous_full 3578a820b8eSToby Isaac args: -dim 1 -continuous 1 -trimmed 0 3588a820b8eSToby Isaac 3598a820b8eSToby Isaac test: 3608a820b8eSToby Isaac suffix: 2_simplex_discontinuous_full 3618a820b8eSToby Isaac args: -dim 2 -tensor 0 -continuous 0 -trimmed 0 3628a820b8eSToby Isaac 3638a820b8eSToby Isaac test: 3648a820b8eSToby Isaac suffix: 2_simplex_continuous_full 3658a820b8eSToby Isaac args: -dim 2 -tensor 0 -continuous 1 -trimmed 0 3668a820b8eSToby Isaac 3678a820b8eSToby Isaac test: 3688a820b8eSToby Isaac suffix: 2_tensor_discontinuous_full 3698a820b8eSToby Isaac args: -dim 2 -tensor 1 -continuous 0 -trimmed 0 3708a820b8eSToby Isaac 3718a820b8eSToby Isaac test: 3728a820b8eSToby Isaac suffix: 2_tensor_continuous_full 3738a820b8eSToby Isaac args: -dim 2 -tensor 1 -continuous 1 -trimmed 0 3748a820b8eSToby Isaac 3758a820b8eSToby Isaac test: 3768a820b8eSToby Isaac suffix: 3_simplex_discontinuous_full 3778a820b8eSToby Isaac args: -dim 3 -tensor 0 -continuous 0 -trimmed 0 3788a820b8eSToby Isaac 3798a820b8eSToby Isaac test: 3808a820b8eSToby Isaac suffix: 3_simplex_continuous_full 3818a820b8eSToby Isaac args: -dim 3 -tensor 0 -continuous 1 -trimmed 0 3828a820b8eSToby Isaac 3838a820b8eSToby Isaac test: 3848a820b8eSToby Isaac suffix: 3_tensor_discontinuous_full 3858a820b8eSToby Isaac args: -dim 3 -tensor 1 -continuous 0 -trimmed 0 3868a820b8eSToby Isaac 3878a820b8eSToby Isaac test: 3888a820b8eSToby Isaac suffix: 3_tensor_continuous_full 3898a820b8eSToby Isaac args: -dim 3 -tensor 1 -continuous 1 -trimmed 0 3908a820b8eSToby Isaac 3918a820b8eSToby Isaac test: 3928a820b8eSToby Isaac suffix: 3_wedge_discontinuous_full 3938a820b8eSToby Isaac args: -dim 3 -tensor 2 -continuous 0 -trimmed 0 3948a820b8eSToby Isaac 3958a820b8eSToby Isaac test: 3968a820b8eSToby Isaac suffix: 3_wedge_continuous_full 3978a820b8eSToby Isaac args: -dim 3 -tensor 2 -continuous 1 -trimmed 0 3988a820b8eSToby Isaac 3998a820b8eSToby Isaac test: 4008a820b8eSToby Isaac suffix: 1_discontinuous_trimmed 4018a820b8eSToby Isaac args: -dim 1 -continuous 0 -trimmed 1 4028a820b8eSToby Isaac 4038a820b8eSToby Isaac test: 4048a820b8eSToby Isaac suffix: 1_continuous_trimmed 4058a820b8eSToby Isaac args: -dim 1 -continuous 1 -trimmed 1 4068a820b8eSToby Isaac 4078a820b8eSToby Isaac test: 4088a820b8eSToby Isaac suffix: 2_simplex_discontinuous_trimmed 4098a820b8eSToby Isaac args: -dim 2 -tensor 0 -continuous 0 -trimmed 1 4108a820b8eSToby Isaac 4118a820b8eSToby Isaac test: 4128a820b8eSToby Isaac suffix: 2_simplex_continuous_trimmed 4138a820b8eSToby Isaac args: -dim 2 -tensor 0 -continuous 1 -trimmed 1 4148a820b8eSToby Isaac 4158a820b8eSToby Isaac test: 4168a820b8eSToby Isaac suffix: 2_tensor_discontinuous_trimmed 4178a820b8eSToby Isaac args: -dim 2 -tensor 1 -continuous 0 -trimmed 1 4188a820b8eSToby Isaac 4198a820b8eSToby Isaac test: 4208a820b8eSToby Isaac suffix: 2_tensor_continuous_trimmed 4218a820b8eSToby Isaac args: -dim 2 -tensor 1 -continuous 1 -trimmed 1 4228a820b8eSToby Isaac 4238a820b8eSToby Isaac test: 4248a820b8eSToby Isaac suffix: 3_simplex_discontinuous_trimmed 4258a820b8eSToby Isaac args: -dim 3 -tensor 0 -continuous 0 -trimmed 1 4268a820b8eSToby Isaac 4278a820b8eSToby Isaac test: 4288a820b8eSToby Isaac suffix: 3_simplex_continuous_trimmed 4298a820b8eSToby Isaac args: -dim 3 -tensor 0 -continuous 1 -trimmed 1 4308a820b8eSToby Isaac 4318a820b8eSToby Isaac test: 4328a820b8eSToby Isaac suffix: 3_tensor_discontinuous_trimmed 4338a820b8eSToby Isaac args: -dim 3 -tensor 1 -continuous 0 -trimmed 1 4348a820b8eSToby Isaac 4358a820b8eSToby Isaac test: 4368a820b8eSToby Isaac suffix: 3_tensor_continuous_trimmed 4378a820b8eSToby Isaac args: -dim 3 -tensor 1 -continuous 1 -trimmed 1 4388a820b8eSToby Isaac 4398a820b8eSToby Isaac test: 4408a820b8eSToby Isaac suffix: 3_wedge_discontinuous_trimmed 4418a820b8eSToby Isaac args: -dim 3 -tensor 2 -continuous 0 -trimmed 1 4428a820b8eSToby Isaac 4438a820b8eSToby Isaac test: 4448a820b8eSToby Isaac suffix: 3_wedge_continuous_trimmed 4458a820b8eSToby Isaac args: -dim 3 -tensor 2 -continuous 1 -trimmed 1 4468a820b8eSToby Isaac 4478a820b8eSToby Isaac TEST*/ 448