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