1 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc-private/matorderimpl.h> /*I "petscmat.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMPlexCreateOrderingClosure_Static" 6 PetscErrorCode DMPlexCreateOrderingClosure_Static(DM dm, PetscInt numPoints, const PetscInt pperm[], PetscInt **clperm, PetscInt **invclperm) 7 { 8 PetscInt *perm, *iperm; 9 PetscInt depth, d, pStart, pEnd, fStart, fMax, fEnd, p; 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 14 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 15 ierr = PetscMalloc1(pEnd-pStart,&perm);CHKERRQ(ierr); 16 ierr = PetscMalloc1(pEnd-pStart,&iperm);CHKERRQ(ierr); 17 for (p = pStart; p < pEnd; ++p) iperm[p] = -1; 18 for (d = depth; d > 0; --d) { 19 ierr = DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);CHKERRQ(ierr); 20 ierr = DMPlexGetDepthStratum(dm, d-1, &fStart, &fEnd);CHKERRQ(ierr); 21 fMax = fStart; 22 for (p = pStart; p < pEnd; ++p) { 23 const PetscInt *cone; 24 PetscInt point, coneSize, c; 25 26 if (d == depth) { 27 perm[p] = pperm[p]; 28 iperm[pperm[p]] = p; 29 } 30 point = perm[p]; 31 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 32 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 33 for (c = 0; c < coneSize; ++c) { 34 const PetscInt oldc = cone[c]; 35 const PetscInt newc = iperm[oldc]; 36 37 if (newc < 0) { 38 perm[fMax] = oldc; 39 iperm[oldc] = fMax++; 40 } 41 } 42 } 43 if (fMax != fEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of depth %d faces %d does not match permuted nubmer %d", d, fEnd-fStart, fMax-fStart); 44 } 45 *clperm = perm; 46 *invclperm = iperm; 47 PetscFunctionReturn(0); 48 } 49 50 #undef __FUNCT__ 51 #define __FUNCT__ "DMPlexGetOrdering" 52 /*@ 53 DMPlexGetOrdering - Calculate a reordering of the mesh 54 55 Collective on DM 56 57 Input Parameter: 58 + dm - The DMPlex object 59 - otype - type of reordering, one of the following: 60 $ MATORDERINGNATURAL - Natural 61 $ MATORDERINGND - Nested Dissection 62 $ MATORDERING1WD - One-way Dissection 63 $ MATORDERINGRCM - Reverse Cuthill-McKee 64 $ MATORDERINGQMD - Quotient Minimum Degree 65 66 67 Output Parameter: 68 . perm - The point permutation as an IS 69 70 Level: intermediate 71 72 .keywords: mesh 73 .seealso: MatGetOrdering() 74 @*/ 75 PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, IS *perm) 76 { 77 PetscInt numCells = 0; 78 PetscInt *start = NULL, *adjacency = NULL, *cperm, *clperm, *invclperm, *mask, *xls, pStart, pEnd, c, i; 79 PetscErrorCode ierr; 80 81 PetscFunctionBegin; 82 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 83 PetscValidPointer(perm, 3); 84 ierr = DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency);CHKERRQ(ierr); 85 ierr = PetscMalloc3(numCells,&cperm,numCells,&mask,numCells*2,&xls);CHKERRQ(ierr); 86 if (numCells) { 87 /* Shift for Fortran numbering */ 88 for (i = 0; i < start[numCells]; ++i) ++adjacency[i]; 89 for (i = 0; i <= numCells; ++i) ++start[i]; 90 ierr = SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls);CHKERRQ(ierr); 91 } 92 ierr = PetscFree(start);CHKERRQ(ierr); 93 ierr = PetscFree(adjacency);CHKERRQ(ierr); 94 /* Shift for Fortran numbering */ 95 for (c = 0; c < numCells; ++c) --cperm[c]; 96 /* Construct closure */ 97 ierr = DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm); 98 ierr = PetscFree3(cperm,mask,xls);CHKERRQ(ierr); 99 ierr = PetscFree(clperm);CHKERRQ(ierr); 100 /* Invert permutation */ 101 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 102 ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm);CHKERRQ(ierr); 103 PetscFunctionReturn(0); 104 } 105 106 #undef __FUNCT__ 107 #define __FUNCT__ "DMPlexPermute" 108 /*@ 109 DMPlexPermute - Reorder the mesh according to the input permutation 110 111 Collective on DM 112 113 Input Parameter: 114 + dm - The DMPlex object 115 - perm - The point permutation 116 117 Output Parameter: 118 . pdm - The permuted DM 119 120 Level: intermediate 121 122 .keywords: mesh 123 .seealso: MatPermute() 124 @*/ 125 PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm) 126 { 127 DM_Plex *plex = (DM_Plex *) dm->data, *plexNew; 128 PetscSection section, sectionNew; 129 PetscInt dim; 130 PetscErrorCode ierr; 131 132 PetscFunctionBegin; 133 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 134 PetscValidHeaderSpecific(perm, IS_CLASSID, 2); 135 PetscValidPointer(pdm, 3); 136 ierr = DMCreate(PetscObjectComm((PetscObject) dm), pdm);CHKERRQ(ierr); 137 ierr = DMSetType(*pdm, DMPLEX);CHKERRQ(ierr); 138 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 139 ierr = DMSetDimension(*pdm, dim);CHKERRQ(ierr); 140 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 141 if (section) { 142 ierr = PetscSectionPermute(section, perm, §ionNew);CHKERRQ(ierr); 143 ierr = DMSetDefaultSection(*pdm, sectionNew);CHKERRQ(ierr); 144 ierr = PetscSectionDestroy(§ionNew);CHKERRQ(ierr); 145 } 146 plexNew = (DM_Plex *) (*pdm)->data; 147 /* Ignore ltogmap, ltogmapb */ 148 /* Ignore sf, defaultSF */ 149 /* Ignore globalVertexNumbers, globalCellNumbers */ 150 /* Remap coordinates */ 151 { 152 DM cdm, cdmNew; 153 PetscSection csection, csectionNew; 154 Vec coordinates, coordinatesNew; 155 PetscScalar *coords, *coordsNew; 156 const PetscInt *pperm; 157 PetscInt pStart, pEnd, p; 158 159 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 160 ierr = DMGetDefaultSection(cdm, &csection);CHKERRQ(ierr); 161 ierr = PetscSectionPermute(csection, perm, &csectionNew);CHKERRQ(ierr); 162 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 163 ierr = VecDuplicate(coordinates, &coordinatesNew);CHKERRQ(ierr); 164 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 165 ierr = VecGetArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 166 ierr = PetscSectionGetChart(csectionNew, &pStart, &pEnd);CHKERRQ(ierr); 167 ierr = ISGetIndices(perm, &pperm);CHKERRQ(ierr); 168 for (p = pStart; p < pEnd; ++p) { 169 PetscInt dof, off, offNew, d; 170 171 ierr = PetscSectionGetDof(csectionNew, p, &dof);CHKERRQ(ierr); 172 ierr = PetscSectionGetOffset(csection, p, &off);CHKERRQ(ierr); 173 ierr = PetscSectionGetOffset(csectionNew, pperm[p], &offNew);CHKERRQ(ierr); 174 for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d]; 175 } 176 ierr = ISRestoreIndices(perm, &pperm);CHKERRQ(ierr); 177 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 178 ierr = VecRestoreArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 179 ierr = DMGetCoordinateDM(*pdm, &cdmNew);CHKERRQ(ierr); 180 ierr = DMSetDefaultSection(cdmNew, csectionNew);CHKERRQ(ierr); 181 ierr = DMSetCoordinatesLocal(*pdm, coordinatesNew);CHKERRQ(ierr); 182 ierr = PetscSectionDestroy(&csectionNew);CHKERRQ(ierr); 183 ierr = VecDestroy(&coordinatesNew);CHKERRQ(ierr); 184 } 185 /* Reorder labels */ 186 { 187 PetscInt numLabels, l; 188 DMLabel label, labelNew; 189 190 ierr = DMPlexGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 191 for (l = numLabels-1; l >= 0; --l) { 192 ierr = DMPlexGetLabelByNum(dm, l, &label);CHKERRQ(ierr); 193 ierr = DMLabelPermute(label, perm, &labelNew);CHKERRQ(ierr); 194 ierr = DMPlexAddLabel(*pdm, labelNew);CHKERRQ(ierr); 195 } 196 if (plex->subpointMap) {ierr = DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);CHKERRQ(ierr);} 197 } 198 /* Reorder topology */ 199 { 200 const PetscInt *pperm; 201 PetscInt maxConeSize, maxSupportSize, n, pStart, pEnd, p; 202 203 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 204 plexNew->maxConeSize = maxConeSize; 205 plexNew->maxSupportSize = maxSupportSize; 206 ierr = PetscSectionDestroy(&plexNew->coneSection);CHKERRQ(ierr); 207 ierr = PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);CHKERRQ(ierr); 208 ierr = PetscSectionGetStorageSize(plexNew->coneSection, &n);CHKERRQ(ierr); 209 ierr = PetscMalloc1(n, &plexNew->cones);CHKERRQ(ierr); 210 ierr = PetscMalloc1(n, &plexNew->coneOrientations);CHKERRQ(ierr); 211 ierr = ISGetIndices(perm, &pperm);CHKERRQ(ierr); 212 ierr = PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);CHKERRQ(ierr); 213 for (p = pStart; p < pEnd; ++p) { 214 PetscInt dof, off, offNew, d; 215 216 ierr = PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);CHKERRQ(ierr); 217 ierr = PetscSectionGetOffset(plex->coneSection, p, &off);CHKERRQ(ierr); 218 ierr = PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);CHKERRQ(ierr); 219 for (d = 0; d < dof; ++d) { 220 plexNew->cones[offNew+d] = pperm[plex->cones[off+d]]; 221 plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d]; 222 } 223 } 224 ierr = PetscSectionDestroy(&plexNew->supportSection);CHKERRQ(ierr); 225 ierr = PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);CHKERRQ(ierr); 226 ierr = PetscSectionGetStorageSize(plexNew->supportSection, &n);CHKERRQ(ierr); 227 ierr = PetscMalloc1(n, &plexNew->supports);CHKERRQ(ierr); 228 ierr = PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);CHKERRQ(ierr); 229 for (p = pStart; p < pEnd; ++p) { 230 PetscInt dof, off, offNew, d; 231 232 ierr = PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);CHKERRQ(ierr); 233 ierr = PetscSectionGetOffset(plex->supportSection, p, &off);CHKERRQ(ierr); 234 ierr = PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);CHKERRQ(ierr); 235 for (d = 0; d < dof; ++d) { 236 plexNew->supports[offNew+d] = pperm[plex->supports[off+d]]; 237 } 238 } 239 ierr = ISRestoreIndices(perm, &pperm);CHKERRQ(ierr); 240 } 241 PetscFunctionReturn(0); 242 } 243