1af0996ceSBarry Smith #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2af0996ceSBarry Smith #include <petsc/private/matorderimpl.h> /*I "petscmat.h" I*/ 3695799ffSMatthew G. Knepley #include <petsc/private/dmlabelimpl.h> 48ed5f475SMatthew G. Knepley 5d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMPlexCreateOrderingClosure_Static(DM dm, PetscInt numPoints, const PetscInt pperm[], PetscInt **clperm, PetscInt **invclperm) 6d71ae5a4SJacob Faibussowitsch { 78ed5f475SMatthew G. Knepley PetscInt *perm, *iperm; 88ed5f475SMatthew G. Knepley PetscInt depth, d, pStart, pEnd, fStart, fMax, fEnd, p; 98ed5f475SMatthew G. Knepley 108ed5f475SMatthew G. Knepley PetscFunctionBegin; 119566063dSJacob Faibussowitsch PetscCall(DMPlexGetDepth(dm, &depth)); 129566063dSJacob Faibussowitsch PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 139566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(pEnd - pStart, &perm)); 149566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(pEnd - pStart, &iperm)); 158ed5f475SMatthew G. Knepley for (p = pStart; p < pEnd; ++p) iperm[p] = -1; 168ed5f475SMatthew G. Knepley for (d = depth; d > 0; --d) { 179566063dSJacob Faibussowitsch PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd)); 189566063dSJacob Faibussowitsch PetscCall(DMPlexGetDepthStratum(dm, d - 1, &fStart, &fEnd)); 198ed5f475SMatthew G. Knepley fMax = fStart; 208ed5f475SMatthew G. Knepley for (p = pStart; p < pEnd; ++p) { 218ed5f475SMatthew G. Knepley const PetscInt *cone; 228ed5f475SMatthew G. Knepley PetscInt point, coneSize, c; 238ed5f475SMatthew G. Knepley 248ed5f475SMatthew G. Knepley if (d == depth) { 258ed5f475SMatthew G. Knepley perm[p] = pperm[p]; 268ed5f475SMatthew G. Knepley iperm[pperm[p]] = p; 278ed5f475SMatthew G. Knepley } 288ed5f475SMatthew G. Knepley point = perm[p]; 299566063dSJacob Faibussowitsch PetscCall(DMPlexGetConeSize(dm, point, &coneSize)); 309566063dSJacob Faibussowitsch PetscCall(DMPlexGetCone(dm, point, &cone)); 318ed5f475SMatthew G. Knepley for (c = 0; c < coneSize; ++c) { 328ed5f475SMatthew G. Knepley const PetscInt oldc = cone[c]; 338ed5f475SMatthew G. Knepley const PetscInt newc = iperm[oldc]; 348ed5f475SMatthew G. Knepley 358ed5f475SMatthew G. Knepley if (newc < 0) { 368ed5f475SMatthew G. Knepley perm[fMax] = oldc; 378ed5f475SMatthew G. Knepley iperm[oldc] = fMax++; 388ed5f475SMatthew G. Knepley } 398ed5f475SMatthew G. Knepley } 408ed5f475SMatthew G. Knepley } 4163a3b9bcSJacob Faibussowitsch PetscCheck(fMax == fEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of depth %" PetscInt_FMT " faces %" PetscInt_FMT " does not match permuted number %" PetscInt_FMT, d, fEnd - fStart, fMax - fStart); 428ed5f475SMatthew G. Knepley } 438ed5f475SMatthew G. Knepley *clperm = perm; 448ed5f475SMatthew G. Knepley *invclperm = iperm; 453ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 468ed5f475SMatthew G. Knepley } 478ed5f475SMatthew G. Knepley 4860225df5SJacob Faibussowitsch /*@C 498ed5f475SMatthew G. Knepley DMPlexGetOrdering - Calculate a reordering of the mesh 508ed5f475SMatthew G. Knepley 5120f4b53cSBarry Smith Collective 528ed5f475SMatthew G. Knepley 53d8d19677SJose E. Roman Input Parameters: 548ed5f475SMatthew G. Knepley + dm - The DMPlex object 5520f4b53cSBarry Smith . otype - type of reordering, see `MatOrderingType` 5620f4b53cSBarry Smith - label - [Optional] Label used to segregate ordering into sets, or `NULL` 578ed5f475SMatthew G. Knepley 588ed5f475SMatthew G. Knepley Output Parameter: 5920f4b53cSBarry Smith . perm - The point permutation as an `IS`, `perm`[old point number] = new point number 608ed5f475SMatthew G. Knepley 618ed5f475SMatthew G. Knepley Level: intermediate 628ed5f475SMatthew G. Knepley 6320f4b53cSBarry Smith Note: 6420f4b53cSBarry Smith The label is used to group sets of points together by label value. This makes it easy to reorder a mesh which 6520f4b53cSBarry Smith has different types of cells, and then loop over each set of reordered cells for assembly. 6620f4b53cSBarry Smith 6720f4b53cSBarry Smith .seealso: `DMPLEX`, `DMPlexPermute()`, `MatOrderingType`, `MatGetOrdering()` 688ed5f475SMatthew G. Knepley @*/ 69d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm) 70d71ae5a4SJacob Faibussowitsch { 718ed5f475SMatthew G. Knepley PetscInt numCells = 0; 723c67d5adSSatish Balay PetscInt *start = NULL, *adjacency = NULL, *cperm, *clperm = NULL, *invclperm = NULL, *mask, *xls, pStart, pEnd, c, i; 738ed5f475SMatthew G. Knepley 748ed5f475SMatthew G. Knepley PetscFunctionBegin; 758ed5f475SMatthew G. Knepley PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 764f572ea9SToby Isaac PetscAssertPointer(perm, 4); 779566063dSJacob Faibussowitsch PetscCall(DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency)); 789566063dSJacob Faibussowitsch PetscCall(PetscMalloc3(numCells, &cperm, numCells, &mask, numCells * 2, &xls)); 79d80c2872SMatthew G. Knepley if (numCells) { 808ed5f475SMatthew G. Knepley /* Shift for Fortran numbering */ 818ed5f475SMatthew G. Knepley for (i = 0; i < start[numCells]; ++i) ++adjacency[i]; 828ed5f475SMatthew G. Knepley for (i = 0; i <= numCells; ++i) ++start[i]; 839566063dSJacob Faibussowitsch PetscCall(SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls)); 84d80c2872SMatthew G. Knepley } 859566063dSJacob Faibussowitsch PetscCall(PetscFree(start)); 869566063dSJacob Faibussowitsch PetscCall(PetscFree(adjacency)); 878ed5f475SMatthew G. Knepley /* Shift for Fortran numbering */ 888ed5f475SMatthew G. Knepley for (c = 0; c < numCells; ++c) --cperm[c]; 89c99f2573SMatthew G. Knepley /* Segregate */ 90c99f2573SMatthew G. Knepley if (label) { 91c99f2573SMatthew G. Knepley IS valueIS; 9250be1cf5SMatthew G. Knepley const PetscInt *valuesTmp; 9350be1cf5SMatthew G. Knepley PetscInt *values; 94c99f2573SMatthew G. Knepley PetscInt numValues, numPoints = 0; 95c99f2573SMatthew G. Knepley PetscInt *sperm, *vsize, *voff, v; 96c99f2573SMatthew G. Knepley 9750be1cf5SMatthew G. Knepley // Can't directly sort the valueIS, since it is a view into the DMLabel 989566063dSJacob Faibussowitsch PetscCall(DMLabelGetValueIS(label, &valueIS)); 999566063dSJacob Faibussowitsch PetscCall(ISGetLocalSize(valueIS, &numValues)); 10050be1cf5SMatthew G. Knepley PetscCall(ISGetIndices(valueIS, &valuesTmp)); 10150be1cf5SMatthew G. Knepley PetscCall(PetscCalloc4(numCells, &sperm, numValues, &values, numValues, &vsize, numValues + 1, &voff)); 10250be1cf5SMatthew G. Knepley PetscCall(PetscArraycpy(values, valuesTmp, numValues)); 10350be1cf5SMatthew G. Knepley PetscCall(PetscSortInt(numValues, values)); 10450be1cf5SMatthew G. Knepley PetscCall(ISRestoreIndices(valueIS, &valuesTmp)); 10550be1cf5SMatthew G. Knepley PetscCall(ISDestroy(&valueIS)); 106c99f2573SMatthew G. Knepley for (v = 0; v < numValues; ++v) { 1079566063dSJacob Faibussowitsch PetscCall(DMLabelGetStratumSize(label, values[v], &vsize[v])); 108c99f2573SMatthew G. Knepley if (v < numValues - 1) voff[v + 2] += vsize[v] + voff[v + 1]; 109c99f2573SMatthew G. Knepley numPoints += vsize[v]; 110c99f2573SMatthew G. Knepley } 1116469bafaSMatthew G. Knepley PetscCheck(numPoints == numCells, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label only covers %" PetscInt_FMT " cells != %" PetscInt_FMT " total cells", numPoints, numCells); 112c99f2573SMatthew G. Knepley for (c = 0; c < numCells; ++c) { 113c99f2573SMatthew G. Knepley const PetscInt oldc = cperm[c]; 114c99f2573SMatthew G. Knepley PetscInt val, vloc; 115c99f2573SMatthew G. Knepley 1169566063dSJacob Faibussowitsch PetscCall(DMLabelGetValue(label, oldc, &val)); 11763a3b9bcSJacob Faibussowitsch PetscCheck(val != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %" PetscInt_FMT " not present in label", oldc); 1189566063dSJacob Faibussowitsch PetscCall(PetscFindInt(val, numValues, values, &vloc)); 11963a3b9bcSJacob Faibussowitsch PetscCheck(vloc >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %" PetscInt_FMT " not present label", val); 120c99f2573SMatthew G. Knepley sperm[voff[vloc + 1]++] = oldc; 121c99f2573SMatthew G. Knepley } 122ad540459SPierre Jolivet for (v = 0; v < numValues; ++v) PetscCheck(voff[v + 1] - voff[v] == vsize[v], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of %" PetscInt_FMT " values found is %" PetscInt_FMT " != %" PetscInt_FMT, values[v], voff[v + 1] - voff[v], vsize[v]); 1239566063dSJacob Faibussowitsch PetscCall(PetscArraycpy(cperm, sperm, numCells)); 12450be1cf5SMatthew G. Knepley PetscCall(PetscFree4(sperm, values, vsize, voff)); 125c99f2573SMatthew G. Knepley } 1268ed5f475SMatthew G. Knepley /* Construct closure */ 1279566063dSJacob Faibussowitsch PetscCall(DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm)); 1289566063dSJacob Faibussowitsch PetscCall(PetscFree3(cperm, mask, xls)); 1299566063dSJacob Faibussowitsch PetscCall(PetscFree(clperm)); 1308ed5f475SMatthew G. Knepley /* Invert permutation */ 1319566063dSJacob Faibussowitsch PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1329566063dSJacob Faibussowitsch PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, invclperm, PETSC_OWN_POINTER, perm)); 1333ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1348ed5f475SMatthew G. Knepley } 1358ed5f475SMatthew G. Knepley 1368ed5f475SMatthew G. Knepley /*@ 1376913077dSMatthew G. Knepley DMPlexGetOrdering1D - Reorder the vertices so that the mesh is in a line 1386913077dSMatthew G. Knepley 13920f4b53cSBarry Smith Collective 1406913077dSMatthew G. Knepley 1416913077dSMatthew G. Knepley Input Parameter: 14220f4b53cSBarry Smith . dm - The `DMPLEX` object 1436913077dSMatthew G. Knepley 1446913077dSMatthew G. Knepley Output Parameter: 14520f4b53cSBarry Smith . perm - The point permutation as an `IS`, `perm`[old point number] = new point number 1466913077dSMatthew G. Knepley 1476913077dSMatthew G. Knepley Level: intermediate 1486913077dSMatthew G. Knepley 14920f4b53cSBarry Smith .seealso: `DMPLEX`, `DMPlexGetOrdering()`, `DMPlexPermute()`, `MatGetOrdering()` 1506913077dSMatthew G. Knepley @*/ 151d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexGetOrdering1D(DM dm, IS *perm) 152d71ae5a4SJacob Faibussowitsch { 1536913077dSMatthew G. Knepley PetscInt *points; 1546913077dSMatthew G. Knepley const PetscInt *support, *cone; 1556913077dSMatthew G. Knepley PetscInt dim, pStart, pEnd, cStart, cEnd, c, vStart, vEnd, v, suppSize, lastCell = 0; 1566913077dSMatthew G. Knepley 1576913077dSMatthew G. Knepley PetscFunctionBegin; 1589566063dSJacob Faibussowitsch PetscCall(DMGetDimension(dm, &dim)); 1596913077dSMatthew G. Knepley PetscCheck(dim == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Input mesh must be one dimensional, not %" PetscInt_FMT, dim); 1609566063dSJacob Faibussowitsch PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1619566063dSJacob Faibussowitsch PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1629566063dSJacob Faibussowitsch PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1639566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(pEnd - pStart, &points)); 1646913077dSMatthew G. Knepley for (c = cStart; c < cEnd; ++c) points[c] = c; 1656913077dSMatthew G. Knepley for (v = vStart; v < vEnd; ++v) points[v] = v; 1666913077dSMatthew G. Knepley for (v = vStart; v < vEnd; ++v) { 1679566063dSJacob Faibussowitsch PetscCall(DMPlexGetSupportSize(dm, v, &suppSize)); 1689566063dSJacob Faibussowitsch PetscCall(DMPlexGetSupport(dm, v, &support)); 1699371c9d4SSatish Balay if (suppSize == 1) { 1709371c9d4SSatish Balay lastCell = support[0]; 1719371c9d4SSatish Balay break; 1729371c9d4SSatish Balay } 1736913077dSMatthew G. Knepley } 1746913077dSMatthew G. Knepley if (v < vEnd) { 1756913077dSMatthew G. Knepley PetscInt pos = cEnd; 1766913077dSMatthew G. Knepley 1776913077dSMatthew G. Knepley points[v] = pos++; 1786913077dSMatthew G. Knepley while (lastCell >= cStart) { 1799566063dSJacob Faibussowitsch PetscCall(DMPlexGetCone(dm, lastCell, &cone)); 1806913077dSMatthew G. Knepley if (cone[0] == v) v = cone[1]; 1816913077dSMatthew G. Knepley else v = cone[0]; 1829566063dSJacob Faibussowitsch PetscCall(DMPlexGetSupport(dm, v, &support)); 1839566063dSJacob Faibussowitsch PetscCall(DMPlexGetSupportSize(dm, v, &suppSize)); 1849371c9d4SSatish Balay if (suppSize == 1) { 1859371c9d4SSatish Balay lastCell = -1; 1869371c9d4SSatish Balay } else { 1876913077dSMatthew G. Knepley if (support[0] == lastCell) lastCell = support[1]; 1886913077dSMatthew G. Knepley else lastCell = support[0]; 1896913077dSMatthew G. Knepley } 1906913077dSMatthew G. Knepley points[v] = pos++; 1916913077dSMatthew G. Knepley } 1926913077dSMatthew G. Knepley PetscCheck(pos == pEnd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Last vertex was %" PetscInt_FMT ", not %" PetscInt_FMT, pos, pEnd); 1936913077dSMatthew G. Knepley } 1949566063dSJacob Faibussowitsch PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), pEnd - pStart, points, PETSC_OWN_POINTER, perm)); 1953ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 1966913077dSMatthew G. Knepley } 1976913077dSMatthew G. Knepley 198d71ae5a4SJacob Faibussowitsch static PetscErrorCode DMPlexRemapCoordinates_Private(IS perm, PetscSection cs, Vec coordinates, PetscSection *csNew, Vec *coordinatesNew) 199d71ae5a4SJacob Faibussowitsch { 2006858538eSMatthew G. Knepley PetscScalar *coords, *coordsNew; 2016858538eSMatthew G. Knepley const PetscInt *pperm; 2026858538eSMatthew G. Knepley PetscInt pStart, pEnd, p; 2036858538eSMatthew G. Knepley const char *name; 2046858538eSMatthew G. Knepley 2056858538eSMatthew G. Knepley PetscFunctionBegin; 2066858538eSMatthew G. Knepley PetscCall(PetscSectionPermute(cs, perm, csNew)); 2076858538eSMatthew G. Knepley PetscCall(VecDuplicate(coordinates, coordinatesNew)); 2086858538eSMatthew G. Knepley PetscCall(PetscObjectGetName((PetscObject)coordinates, &name)); 2096858538eSMatthew G. Knepley PetscCall(PetscObjectSetName((PetscObject)*coordinatesNew, name)); 2106858538eSMatthew G. Knepley PetscCall(VecGetArray(coordinates, &coords)); 2116858538eSMatthew G. Knepley PetscCall(VecGetArray(*coordinatesNew, &coordsNew)); 2126858538eSMatthew G. Knepley PetscCall(PetscSectionGetChart(*csNew, &pStart, &pEnd)); 2136858538eSMatthew G. Knepley PetscCall(ISGetIndices(perm, &pperm)); 2146858538eSMatthew G. Knepley for (p = pStart; p < pEnd; ++p) { 2156858538eSMatthew G. Knepley PetscInt dof, off, offNew, d; 2166858538eSMatthew G. Knepley 2176858538eSMatthew G. Knepley PetscCall(PetscSectionGetDof(*csNew, p, &dof)); 2186858538eSMatthew G. Knepley PetscCall(PetscSectionGetOffset(cs, p, &off)); 2196858538eSMatthew G. Knepley PetscCall(PetscSectionGetOffset(*csNew, pperm[p], &offNew)); 2206858538eSMatthew G. Knepley for (d = 0; d < dof; ++d) coordsNew[offNew + d] = coords[off + d]; 2216858538eSMatthew G. Knepley } 2226858538eSMatthew G. Knepley PetscCall(ISRestoreIndices(perm, &pperm)); 2236858538eSMatthew G. Knepley PetscCall(VecRestoreArray(coordinates, &coords)); 2246858538eSMatthew G. Knepley PetscCall(VecRestoreArray(*coordinatesNew, &coordsNew)); 2253ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 2266858538eSMatthew G. Knepley } 2276858538eSMatthew G. Knepley 2286913077dSMatthew G. Knepley /*@ 2298ed5f475SMatthew G. Knepley DMPlexPermute - Reorder the mesh according to the input permutation 2308ed5f475SMatthew G. Knepley 23120f4b53cSBarry Smith Collective 2328ed5f475SMatthew G. Knepley 233d8d19677SJose E. Roman Input Parameters: 23420f4b53cSBarry Smith + dm - The `DMPLEX` object 23520f4b53cSBarry Smith - perm - The point permutation, `perm`[old point number] = new point number 2368ed5f475SMatthew G. Knepley 2378ed5f475SMatthew G. Knepley Output Parameter: 23820f4b53cSBarry Smith . pdm - The permuted `DM` 2398ed5f475SMatthew G. Knepley 2408ed5f475SMatthew G. Knepley Level: intermediate 2418ed5f475SMatthew G. Knepley 24220f4b53cSBarry Smith .seealso: `DMPLEX`, `MatPermute()` 2438ed5f475SMatthew G. Knepley @*/ 244d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm) 245d71ae5a4SJacob Faibussowitsch { 2468ed5f475SMatthew G. Knepley DM_Plex *plex = (DM_Plex *)dm->data, *plexNew; 2472f300c92SMatthew Knepley PetscInt dim, cdim; 2486b74b800SMatthew G. Knepley const char *name; 2498ed5f475SMatthew G. Knepley 2508ed5f475SMatthew G. Knepley PetscFunctionBegin; 2518ed5f475SMatthew G. Knepley PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2528ed5f475SMatthew G. Knepley PetscValidHeaderSpecific(perm, IS_CLASSID, 2); 2534f572ea9SToby Isaac PetscAssertPointer(pdm, 3); 2549566063dSJacob Faibussowitsch PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), pdm)); 2559566063dSJacob Faibussowitsch PetscCall(DMSetType(*pdm, DMPLEX)); 2569566063dSJacob Faibussowitsch PetscCall(PetscObjectGetName((PetscObject)dm, &name)); 2579566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)*pdm, name)); 2589566063dSJacob Faibussowitsch PetscCall(DMGetDimension(dm, &dim)); 2599566063dSJacob Faibussowitsch PetscCall(DMSetDimension(*pdm, dim)); 2609566063dSJacob Faibussowitsch PetscCall(DMGetCoordinateDim(dm, &cdim)); 2619566063dSJacob Faibussowitsch PetscCall(DMSetCoordinateDim(*pdm, cdim)); 2629566063dSJacob Faibussowitsch PetscCall(DMCopyDisc(dm, *pdm)); 2636b74b800SMatthew G. Knepley if (dm->localSection) { 2646b74b800SMatthew G. Knepley PetscSection section, sectionNew; 2656b74b800SMatthew G. Knepley 2669566063dSJacob Faibussowitsch PetscCall(DMGetLocalSection(dm, §ion)); 2679566063dSJacob Faibussowitsch PetscCall(PetscSectionPermute(section, perm, §ionNew)); 2689566063dSJacob Faibussowitsch PetscCall(DMSetLocalSection(*pdm, sectionNew)); 2699566063dSJacob Faibussowitsch PetscCall(PetscSectionDestroy(§ionNew)); 270dd742cf1SMatthew G. Knepley } 2718ed5f475SMatthew G. Knepley plexNew = (DM_Plex *)(*pdm)->data; 2728ed5f475SMatthew G. Knepley /* Ignore ltogmap, ltogmapb */ 2731bb6d2a8SBarry Smith /* Ignore sf, sectionSF */ 2748ed5f475SMatthew G. Knepley /* Ignore globalVertexNumbers, globalCellNumbers */ 2758ed5f475SMatthew G. Knepley /* Reorder labels */ 2768ed5f475SMatthew G. Knepley { 277a85475f2SMatthew G. Knepley PetscInt numLabels, l; 278a85475f2SMatthew G. Knepley DMLabel label, labelNew; 2798ed5f475SMatthew G. Knepley 2809566063dSJacob Faibussowitsch PetscCall(DMGetNumLabels(dm, &numLabels)); 2815d80c0bfSVaclav Hapla for (l = 0; l < numLabels; ++l) { 2829566063dSJacob Faibussowitsch PetscCall(DMGetLabelByNum(dm, l, &label)); 2839566063dSJacob Faibussowitsch PetscCall(DMLabelPermute(label, perm, &labelNew)); 2849566063dSJacob Faibussowitsch PetscCall(DMAddLabel(*pdm, labelNew)); 2859566063dSJacob Faibussowitsch PetscCall(DMLabelDestroy(&labelNew)); 2868ed5f475SMatthew G. Knepley } 2879566063dSJacob Faibussowitsch PetscCall(DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel)); 2889566063dSJacob Faibussowitsch if (plex->subpointMap) PetscCall(DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap)); 2898ed5f475SMatthew G. Knepley } 290695799ffSMatthew G. Knepley if ((*pdm)->celltypeLabel) { 291695799ffSMatthew G. Knepley DMLabel ctLabel; 292695799ffSMatthew G. Knepley 293695799ffSMatthew G. Knepley // Reset label for fast lookup 294695799ffSMatthew G. Knepley PetscCall(DMPlexGetCellTypeLabel(*pdm, &ctLabel)); 295695799ffSMatthew G. Knepley PetscCall(DMLabelMakeAllInvalid_Internal(ctLabel)); 296695799ffSMatthew G. Knepley } 2978ed5f475SMatthew G. Knepley /* Reorder topology */ 2988ed5f475SMatthew G. Knepley { 2998ed5f475SMatthew G. Knepley const PetscInt *pperm; 3006302a7fbSVaclav Hapla PetscInt n, pStart, pEnd, p; 3018ed5f475SMatthew G. Knepley 3029566063dSJacob Faibussowitsch PetscCall(PetscSectionDestroy(&plexNew->coneSection)); 3039566063dSJacob Faibussowitsch PetscCall(PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection)); 3049566063dSJacob Faibussowitsch PetscCall(PetscSectionGetStorageSize(plexNew->coneSection, &n)); 3059566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(n, &plexNew->cones)); 3069566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(n, &plexNew->coneOrientations)); 3079566063dSJacob Faibussowitsch PetscCall(ISGetIndices(perm, &pperm)); 3089566063dSJacob Faibussowitsch PetscCall(PetscSectionGetChart(plex->coneSection, &pStart, &pEnd)); 3098ed5f475SMatthew G. Knepley for (p = pStart; p < pEnd; ++p) { 3108ed5f475SMatthew G. Knepley PetscInt dof, off, offNew, d; 3118ed5f475SMatthew G. Knepley 3129566063dSJacob Faibussowitsch PetscCall(PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof)); 3139566063dSJacob Faibussowitsch PetscCall(PetscSectionGetOffset(plex->coneSection, p, &off)); 3149566063dSJacob Faibussowitsch PetscCall(PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew)); 3158ed5f475SMatthew G. Knepley for (d = 0; d < dof; ++d) { 3168ed5f475SMatthew G. Knepley plexNew->cones[offNew + d] = pperm[plex->cones[off + d]]; 3178ed5f475SMatthew G. Knepley plexNew->coneOrientations[offNew + d] = plex->coneOrientations[off + d]; 3188ed5f475SMatthew G. Knepley } 3198ed5f475SMatthew G. Knepley } 3209566063dSJacob Faibussowitsch PetscCall(PetscSectionDestroy(&plexNew->supportSection)); 3219566063dSJacob Faibussowitsch PetscCall(PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection)); 3229566063dSJacob Faibussowitsch PetscCall(PetscSectionGetStorageSize(plexNew->supportSection, &n)); 3239566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(n, &plexNew->supports)); 3249566063dSJacob Faibussowitsch PetscCall(PetscSectionGetChart(plex->supportSection, &pStart, &pEnd)); 3258ed5f475SMatthew G. Knepley for (p = pStart; p < pEnd; ++p) { 3268ed5f475SMatthew G. Knepley PetscInt dof, off, offNew, d; 3278ed5f475SMatthew G. Knepley 3289566063dSJacob Faibussowitsch PetscCall(PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof)); 3299566063dSJacob Faibussowitsch PetscCall(PetscSectionGetOffset(plex->supportSection, p, &off)); 3309566063dSJacob Faibussowitsch PetscCall(PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew)); 331ad540459SPierre Jolivet for (d = 0; d < dof; ++d) plexNew->supports[offNew + d] = pperm[plex->supports[off + d]]; 3328ed5f475SMatthew G. Knepley } 3339566063dSJacob Faibussowitsch PetscCall(ISRestoreIndices(perm, &pperm)); 3348ed5f475SMatthew G. Knepley } 335a6e0b375SMatthew G. Knepley /* Remap coordinates */ 336a6e0b375SMatthew G. Knepley { 337a6e0b375SMatthew G. Knepley DM cdm, cdmNew; 3386858538eSMatthew G. Knepley PetscSection cs, csNew; 339a6e0b375SMatthew G. Knepley Vec coordinates, coordinatesNew; 340a6e0b375SMatthew G. Knepley 3416858538eSMatthew G. Knepley PetscCall(DMGetCoordinateSection(dm, &cs)); 3429566063dSJacob Faibussowitsch PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 3436858538eSMatthew G. Knepley PetscCall(DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew)); 3446858538eSMatthew G. Knepley PetscCall(DMSetCoordinateSection(*pdm, PETSC_DETERMINE, csNew)); 3459566063dSJacob Faibussowitsch PetscCall(DMSetCoordinatesLocal(*pdm, coordinatesNew)); 3466858538eSMatthew G. Knepley PetscCall(PetscSectionDestroy(&csNew)); 3479566063dSJacob Faibussowitsch PetscCall(VecDestroy(&coordinatesNew)); 3486858538eSMatthew G. Knepley 3496858538eSMatthew G. Knepley PetscCall(DMGetCellCoordinateDM(dm, &cdm)); 3506858538eSMatthew G. Knepley if (cdm) { 3516858538eSMatthew G. Knepley PetscCall(DMGetCoordinateDM(*pdm, &cdm)); 3526858538eSMatthew G. Knepley PetscCall(DMClone(cdm, &cdmNew)); 3536858538eSMatthew G. Knepley PetscCall(DMSetCellCoordinateDM(*pdm, cdmNew)); 3546858538eSMatthew G. Knepley PetscCall(DMDestroy(&cdmNew)); 3556858538eSMatthew G. Knepley PetscCall(DMGetCellCoordinateSection(dm, &cs)); 3566858538eSMatthew G. Knepley PetscCall(DMGetCellCoordinatesLocal(dm, &coordinates)); 3576858538eSMatthew G. Knepley PetscCall(DMPlexRemapCoordinates_Private(perm, cs, coordinates, &csNew, &coordinatesNew)); 3586858538eSMatthew G. Knepley PetscCall(DMSetCellCoordinateSection(*pdm, PETSC_DETERMINE, csNew)); 3596858538eSMatthew G. Knepley PetscCall(DMSetCellCoordinatesLocal(*pdm, coordinatesNew)); 3606858538eSMatthew G. Knepley PetscCall(PetscSectionDestroy(&csNew)); 3616858538eSMatthew G. Knepley PetscCall(VecDestroy(&coordinatesNew)); 3626858538eSMatthew G. Knepley } 363a6e0b375SMatthew G. Knepley } 3645de52c6dSVaclav Hapla PetscCall(DMPlexCopy_Internal(dm, PETSC_TRUE, PETSC_TRUE, *pdm)); 365e66f2fa0SMatthew G. Knepley (*pdm)->setupcalled = PETSC_TRUE; 3663ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3678ed5f475SMatthew G. Knepley } 3686bc1bd01Sksagiyam 369d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexReorderSetDefault_Plex(DM dm, DMPlexReorderDefaultFlag reorder) 370d71ae5a4SJacob Faibussowitsch { 3716bc1bd01Sksagiyam DM_Plex *mesh = (DM_Plex *)dm->data; 3726bc1bd01Sksagiyam 3736bc1bd01Sksagiyam PetscFunctionBegin; 3746bc1bd01Sksagiyam mesh->reorderDefault = reorder; 3753ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3766bc1bd01Sksagiyam } 3776bc1bd01Sksagiyam 3786bc1bd01Sksagiyam /*@ 3796bc1bd01Sksagiyam DMPlexReorderSetDefault - Set flag indicating whether the DM should be reordered by default 3806bc1bd01Sksagiyam 38120f4b53cSBarry Smith Logically Collective 3826bc1bd01Sksagiyam 3836bc1bd01Sksagiyam Input Parameters: 38420f4b53cSBarry Smith + dm - The `DM` 3856bc1bd01Sksagiyam - reorder - Flag for reordering 3866bc1bd01Sksagiyam 3876bc1bd01Sksagiyam Level: intermediate 3886bc1bd01Sksagiyam 3896bc1bd01Sksagiyam .seealso: `DMPlexReorderGetDefault()` 3906bc1bd01Sksagiyam @*/ 391d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexReorderSetDefault(DM dm, DMPlexReorderDefaultFlag reorder) 392d71ae5a4SJacob Faibussowitsch { 3936bc1bd01Sksagiyam PetscFunctionBegin; 3946bc1bd01Sksagiyam PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 3956bc1bd01Sksagiyam PetscTryMethod(dm, "DMPlexReorderSetDefault_C", (DM, DMPlexReorderDefaultFlag), (dm, reorder)); 3963ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 3976bc1bd01Sksagiyam } 3986bc1bd01Sksagiyam 399d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexReorderGetDefault_Plex(DM dm, DMPlexReorderDefaultFlag *reorder) 400d71ae5a4SJacob Faibussowitsch { 4016bc1bd01Sksagiyam DM_Plex *mesh = (DM_Plex *)dm->data; 4026bc1bd01Sksagiyam 4036bc1bd01Sksagiyam PetscFunctionBegin; 4046bc1bd01Sksagiyam *reorder = mesh->reorderDefault; 4053ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 4066bc1bd01Sksagiyam } 4076bc1bd01Sksagiyam 4086bc1bd01Sksagiyam /*@ 4096bc1bd01Sksagiyam DMPlexReorderGetDefault - Get flag indicating whether the DM should be reordered by default 4106bc1bd01Sksagiyam 41120f4b53cSBarry Smith Not Collective 4126bc1bd01Sksagiyam 4136bc1bd01Sksagiyam Input Parameter: 41420f4b53cSBarry Smith . dm - The `DM` 4156bc1bd01Sksagiyam 4166bc1bd01Sksagiyam Output Parameter: 4176bc1bd01Sksagiyam . reorder - Flag for reordering 4186bc1bd01Sksagiyam 4196bc1bd01Sksagiyam Level: intermediate 4206bc1bd01Sksagiyam 4216bc1bd01Sksagiyam .seealso: `DMPlexReorderSetDefault()` 4226bc1bd01Sksagiyam @*/ 423d71ae5a4SJacob Faibussowitsch PetscErrorCode DMPlexReorderGetDefault(DM dm, DMPlexReorderDefaultFlag *reorder) 424d71ae5a4SJacob Faibussowitsch { 4256bc1bd01Sksagiyam PetscFunctionBegin; 4266bc1bd01Sksagiyam PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 4274f572ea9SToby Isaac PetscAssertPointer(reorder, 2); 4286bc1bd01Sksagiyam PetscUseMethod(dm, "DMPlexReorderGetDefault_C", (DM, DMPlexReorderDefaultFlag *), (dm, reorder)); 4293ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 4306bc1bd01Sksagiyam } 431*d02c7345SMatthew G. Knepley 432*d02c7345SMatthew G. Knepley // Reorder to group split nodes 433*d02c7345SMatthew G. Knepley PetscErrorCode DMPlexCreateSectionPermutation_Internal(DM dm, IS *permutation, PetscBT *blockStarts) 434*d02c7345SMatthew G. Knepley { 435*d02c7345SMatthew G. Knepley IS permIS; 436*d02c7345SMatthew G. Knepley PetscBT bt, blst; 437*d02c7345SMatthew G. Knepley PetscInt *perm; 438*d02c7345SMatthew G. Knepley PetscInt pStart, pEnd, i = 0; 439*d02c7345SMatthew G. Knepley 440*d02c7345SMatthew G. Knepley PetscFunctionBegin; 441*d02c7345SMatthew G. Knepley *permutation = NULL; 442*d02c7345SMatthew G. Knepley *blockStarts = NULL; 443*d02c7345SMatthew G. Knepley { 444*d02c7345SMatthew G. Knepley DMPlexReorderDefaultFlag reorder; 445*d02c7345SMatthew G. Knepley PetscCall(DMPlexReorderSectionGetDefault(dm, &reorder)); 446*d02c7345SMatthew G. Knepley if (reorder != DMPLEX_REORDER_DEFAULT_TRUE) PetscFunctionReturn(PETSC_SUCCESS); 447*d02c7345SMatthew G. Knepley } 448*d02c7345SMatthew G. Knepley PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 449*d02c7345SMatthew G. Knepley PetscCall(PetscMalloc1(pEnd - pStart, &perm)); 450*d02c7345SMatthew G. Knepley PetscCall(PetscBTCreate(pEnd - pStart, &bt)); 451*d02c7345SMatthew G. Knepley PetscCall(PetscBTCreate(pEnd - pStart, &blst)); 452*d02c7345SMatthew G. Knepley for (PetscInt p = pStart; p < pEnd; ++p) { 453*d02c7345SMatthew G. Knepley const PetscInt *supp, *cone; 454*d02c7345SMatthew G. Knepley PetscInt suppSize; 455*d02c7345SMatthew G. Knepley 456*d02c7345SMatthew G. Knepley if (PetscBTLookupSet(bt, p)) continue; 457*d02c7345SMatthew G. Knepley PetscCall(PetscBTSet(blst, p)); 458*d02c7345SMatthew G. Knepley perm[i++] = p; 459*d02c7345SMatthew G. Knepley // Check for tensor cells in the support 460*d02c7345SMatthew G. Knepley PetscCall(DMPlexGetSupport(dm, p, &supp)); 461*d02c7345SMatthew G. Knepley PetscCall(DMPlexGetSupportSize(dm, p, &suppSize)); 462*d02c7345SMatthew G. Knepley for (PetscInt s = 0; s < suppSize; ++s) { 463*d02c7345SMatthew G. Knepley DMPolytopeType ct; 464*d02c7345SMatthew G. Knepley PetscInt q, qq; 465*d02c7345SMatthew G. Knepley 466*d02c7345SMatthew G. Knepley PetscCall(DMPlexGetCellType(dm, supp[s], &ct)); 467*d02c7345SMatthew G. Knepley switch (ct) { 468*d02c7345SMatthew G. Knepley case DM_POLYTOPE_POINT_PRISM_TENSOR: 469*d02c7345SMatthew G. Knepley case DM_POLYTOPE_SEG_PRISM_TENSOR: 470*d02c7345SMatthew G. Knepley case DM_POLYTOPE_TRI_PRISM_TENSOR: 471*d02c7345SMatthew G. Knepley case DM_POLYTOPE_QUAD_PRISM_TENSOR: 472*d02c7345SMatthew G. Knepley // If found, move up the split partner of the tensor cell, and the cell itself 473*d02c7345SMatthew G. Knepley PetscCall(DMPlexGetCone(dm, supp[s], &cone)); 474*d02c7345SMatthew G. Knepley qq = supp[s]; 475*d02c7345SMatthew G. Knepley q = (cone[0] == p) ? cone[1] : cone[0]; 476*d02c7345SMatthew G. Knepley if (!PetscBTLookupSet(bt, q)) { 477*d02c7345SMatthew G. Knepley perm[i++] = q; 478*d02c7345SMatthew G. Knepley s = suppSize; 479*d02c7345SMatthew G. Knepley } 480*d02c7345SMatthew G. Knepley if (!PetscBTLookupSet(bt, qq)) { 481*d02c7345SMatthew G. Knepley perm[i++] = qq; 482*d02c7345SMatthew G. Knepley s = suppSize; 483*d02c7345SMatthew G. Knepley } 484*d02c7345SMatthew G. Knepley break; 485*d02c7345SMatthew G. Knepley default: 486*d02c7345SMatthew G. Knepley break; 487*d02c7345SMatthew G. Knepley } 488*d02c7345SMatthew G. Knepley } 489*d02c7345SMatthew G. Knepley } 490*d02c7345SMatthew G. Knepley PetscCall(PetscBTDestroy(&bt)); 491*d02c7345SMatthew G. Knepley PetscCheck(i == pEnd - pStart, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of points in permutation %" PetscInt_FMT " does not match chart size %" PetscInt_FMT, i, pEnd - pStart); 492*d02c7345SMatthew G. Knepley PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, perm, PETSC_OWN_POINTER, &permIS)); 493*d02c7345SMatthew G. Knepley PetscCall(ISSetPermutation(permIS)); 494*d02c7345SMatthew G. Knepley *permutation = permIS; 495*d02c7345SMatthew G. Knepley *blockStarts = blst; 496*d02c7345SMatthew G. Knepley PetscFunctionReturn(PETSC_SUCCESS); 497*d02c7345SMatthew G. Knepley } 498*d02c7345SMatthew G. Knepley 499*d02c7345SMatthew G. Knepley PetscErrorCode DMPlexReorderSectionSetDefault_Plex(DM dm, DMPlexReorderDefaultFlag reorder) 500*d02c7345SMatthew G. Knepley { 501*d02c7345SMatthew G. Knepley DM_Plex *mesh = (DM_Plex *)dm->data; 502*d02c7345SMatthew G. Knepley 503*d02c7345SMatthew G. Knepley PetscFunctionBegin; 504*d02c7345SMatthew G. Knepley mesh->reorderSection = reorder; 505*d02c7345SMatthew G. Knepley PetscFunctionReturn(PETSC_SUCCESS); 506*d02c7345SMatthew G. Knepley } 507*d02c7345SMatthew G. Knepley 508*d02c7345SMatthew G. Knepley /*@ 509*d02c7345SMatthew G. Knepley DMPlexReorderSectionSetDefault - Set flag indicating whether the local section should be reordered by default 510*d02c7345SMatthew G. Knepley 511*d02c7345SMatthew G. Knepley Logically collective 512*d02c7345SMatthew G. Knepley 513*d02c7345SMatthew G. Knepley Input Parameters: 514*d02c7345SMatthew G. Knepley + dm - The DM 515*d02c7345SMatthew G. Knepley - reorder - Flag for reordering 516*d02c7345SMatthew G. Knepley 517*d02c7345SMatthew G. Knepley Level: intermediate 518*d02c7345SMatthew G. Knepley 519*d02c7345SMatthew G. Knepley .seealso: `DMPlexReorderSectionGetDefault()` 520*d02c7345SMatthew G. Knepley @*/ 521*d02c7345SMatthew G. Knepley PetscErrorCode DMPlexReorderSectionSetDefault(DM dm, DMPlexReorderDefaultFlag reorder) 522*d02c7345SMatthew G. Knepley { 523*d02c7345SMatthew G. Knepley PetscFunctionBegin; 524*d02c7345SMatthew G. Knepley PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 525*d02c7345SMatthew G. Knepley PetscTryMethod(dm, "DMPlexReorderSectionSetDefault_C", (DM, DMPlexReorderDefaultFlag), (dm, reorder)); 526*d02c7345SMatthew G. Knepley PetscFunctionReturn(PETSC_SUCCESS); 527*d02c7345SMatthew G. Knepley } 528*d02c7345SMatthew G. Knepley 529*d02c7345SMatthew G. Knepley PetscErrorCode DMPlexReorderSectionGetDefault_Plex(DM dm, DMPlexReorderDefaultFlag *reorder) 530*d02c7345SMatthew G. Knepley { 531*d02c7345SMatthew G. Knepley DM_Plex *mesh = (DM_Plex *)dm->data; 532*d02c7345SMatthew G. Knepley 533*d02c7345SMatthew G. Knepley PetscFunctionBegin; 534*d02c7345SMatthew G. Knepley *reorder = mesh->reorderSection; 535*d02c7345SMatthew G. Knepley PetscFunctionReturn(PETSC_SUCCESS); 536*d02c7345SMatthew G. Knepley } 537*d02c7345SMatthew G. Knepley 538*d02c7345SMatthew G. Knepley /*@ 539*d02c7345SMatthew G. Knepley DMPlexReorderSectionGetDefault - Get flag indicating whether the local section should be reordered by default 540*d02c7345SMatthew G. Knepley 541*d02c7345SMatthew G. Knepley Not collective 542*d02c7345SMatthew G. Knepley 543*d02c7345SMatthew G. Knepley Input Parameter: 544*d02c7345SMatthew G. Knepley . dm - The DM 545*d02c7345SMatthew G. Knepley 546*d02c7345SMatthew G. Knepley Output Parameter: 547*d02c7345SMatthew G. Knepley . reorder - Flag for reordering 548*d02c7345SMatthew G. Knepley 549*d02c7345SMatthew G. Knepley Level: intermediate 550*d02c7345SMatthew G. Knepley 551*d02c7345SMatthew G. Knepley .seealso: `DMPlexReorderSetDefault()` 552*d02c7345SMatthew G. Knepley @*/ 553*d02c7345SMatthew G. Knepley PetscErrorCode DMPlexReorderSectionGetDefault(DM dm, DMPlexReorderDefaultFlag *reorder) 554*d02c7345SMatthew G. Knepley { 555*d02c7345SMatthew G. Knepley PetscFunctionBegin; 556*d02c7345SMatthew G. Knepley PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 557*d02c7345SMatthew G. Knepley PetscAssertPointer(reorder, 2); 558*d02c7345SMatthew G. Knepley PetscUseMethod(dm, "DMPlexReorderSectionGetDefault_C", (DM, DMPlexReorderDefaultFlag *), (dm, reorder)); 559*d02c7345SMatthew G. Knepley PetscFunctionReturn(PETSC_SUCCESS); 560*d02c7345SMatthew G. Knepley } 561