1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/matorderimpl.h> /*I "petscmat.h" I*/ 3 4 static PetscErrorCode DMPlexCreateOrderingClosure_Static(DM dm, PetscInt numPoints, const PetscInt pperm[], PetscInt **clperm, PetscInt **invclperm) 5 { 6 PetscInt *perm, *iperm; 7 PetscInt depth, d, pStart, pEnd, fStart, fMax, fEnd, p; 8 PetscErrorCode ierr; 9 10 PetscFunctionBegin; 11 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 12 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 13 ierr = PetscMalloc1(pEnd-pStart,&perm);CHKERRQ(ierr); 14 ierr = PetscMalloc1(pEnd-pStart,&iperm);CHKERRQ(ierr); 15 for (p = pStart; p < pEnd; ++p) iperm[p] = -1; 16 for (d = depth; d > 0; --d) { 17 ierr = DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);CHKERRQ(ierr); 18 ierr = DMPlexGetDepthStratum(dm, d-1, &fStart, &fEnd);CHKERRQ(ierr); 19 fMax = fStart; 20 for (p = pStart; p < pEnd; ++p) { 21 const PetscInt *cone; 22 PetscInt point, coneSize, c; 23 24 if (d == depth) { 25 perm[p] = pperm[p]; 26 iperm[pperm[p]] = p; 27 } 28 point = perm[p]; 29 ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr); 30 ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr); 31 for (c = 0; c < coneSize; ++c) { 32 const PetscInt oldc = cone[c]; 33 const PetscInt newc = iperm[oldc]; 34 35 if (newc < 0) { 36 perm[fMax] = oldc; 37 iperm[oldc] = fMax++; 38 } 39 } 40 } 41 if (fMax != fEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of depth %d faces %d does not match permuted number %d", d, fEnd-fStart, fMax-fStart); 42 } 43 *clperm = perm; 44 *invclperm = iperm; 45 PetscFunctionReturn(0); 46 } 47 48 /*@ 49 DMPlexGetOrdering - Calculate a reordering of the mesh 50 51 Collective on dm 52 53 Input Parameters: 54 + dm - The DMPlex object 55 . otype - type of reordering, one of the following: 56 $ MATORDERINGNATURAL - Natural 57 $ MATORDERINGND - Nested Dissection 58 $ MATORDERING1WD - One-way Dissection 59 $ MATORDERINGRCM - Reverse Cuthill-McKee 60 $ MATORDERINGQMD - Quotient Minimum Degree 61 - label - [Optional] Label used to segregate ordering into sets, or NULL 62 63 Output Parameter: 64 . perm - The point permutation as an IS, perm[old point number] = new point number 65 66 Note: The label is used to group sets of points together by label value. This makes it easy to reorder a mesh which 67 has different types of cells, and then loop over each set of reordered cells for assembly. 68 69 Level: intermediate 70 71 .seealso: MatGetOrdering() 72 @*/ 73 PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, DMLabel label, IS *perm) 74 { 75 PetscInt numCells = 0; 76 PetscInt *start = NULL, *adjacency = NULL, *cperm, *clperm = NULL, *invclperm = NULL, *mask, *xls, pStart, pEnd, c, i; 77 PetscErrorCode ierr; 78 79 PetscFunctionBegin; 80 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 81 PetscValidPointer(perm, 4); 82 ierr = DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency);CHKERRQ(ierr); 83 ierr = PetscMalloc3(numCells,&cperm,numCells,&mask,numCells*2,&xls);CHKERRQ(ierr); 84 if (numCells) { 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 } 90 ierr = PetscFree(start);CHKERRQ(ierr); 91 ierr = PetscFree(adjacency);CHKERRQ(ierr); 92 /* Shift for Fortran numbering */ 93 for (c = 0; c < numCells; ++c) --cperm[c]; 94 /* Segregate */ 95 if (label) { 96 IS valueIS; 97 const PetscInt *values; 98 PetscInt numValues, numPoints = 0; 99 PetscInt *sperm, *vsize, *voff, v; 100 101 ierr = DMLabelGetValueIS(label, &valueIS);CHKERRQ(ierr); 102 ierr = ISSort(valueIS);CHKERRQ(ierr); 103 ierr = ISGetLocalSize(valueIS, &numValues);CHKERRQ(ierr); 104 ierr = ISGetIndices(valueIS, &values);CHKERRQ(ierr); 105 ierr = PetscCalloc3(numCells,&sperm,numValues,&vsize,numValues+1,&voff);CHKERRQ(ierr); 106 for (v = 0; v < numValues; ++v) { 107 ierr = DMLabelGetStratumSize(label, values[v], &vsize[v]);CHKERRQ(ierr); 108 if (v < numValues-1) voff[v+2] += vsize[v] + voff[v+1]; 109 numPoints += vsize[v]; 110 } 111 if (numPoints != numCells) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label only covers %D cells < %D total", numPoints, numCells); 112 for (c = 0; c < numCells; ++c) { 113 const PetscInt oldc = cperm[c]; 114 PetscInt val, vloc; 115 116 ierr = DMLabelGetValue(label, oldc, &val);CHKERRQ(ierr); 117 if (val == -1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cell %D not present in label", oldc); 118 ierr = PetscFindInt(val, numValues, values, &vloc);CHKERRQ(ierr); 119 if (vloc < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %D not present label", val); 120 sperm[voff[vloc+1]++] = oldc; 121 } 122 for (v = 0; v < numValues; ++v) { 123 if (voff[v+1] - voff[v] != vsize[v]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of %D values found is %D != %D", values[v], voff[v+1] - voff[v], vsize[v]); 124 } 125 ierr = ISRestoreIndices(valueIS, &values);CHKERRQ(ierr); 126 ierr = ISDestroy(&valueIS);CHKERRQ(ierr); 127 ierr = PetscArraycpy(cperm, sperm, numCells);CHKERRQ(ierr); 128 ierr = PetscFree3(sperm, vsize, voff);CHKERRQ(ierr); 129 } 130 /* Construct closure */ 131 ierr = DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm);CHKERRQ(ierr); 132 ierr = PetscFree3(cperm,mask,xls);CHKERRQ(ierr); 133 ierr = PetscFree(clperm);CHKERRQ(ierr); 134 /* Invert permutation */ 135 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 136 ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm);CHKERRQ(ierr); 137 PetscFunctionReturn(0); 138 } 139 140 /*@ 141 DMPlexPermute - Reorder the mesh according to the input permutation 142 143 Collective on dm 144 145 Input Parameters: 146 + dm - The DMPlex object 147 - perm - The point permutation, perm[old point number] = new point number 148 149 Output Parameter: 150 . pdm - The permuted DM 151 152 Level: intermediate 153 154 .seealso: MatPermute() 155 @*/ 156 PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm) 157 { 158 DM_Plex *plex = (DM_Plex *) dm->data, *plexNew; 159 PetscSection section, sectionNew; 160 PetscInt dim, cdim; 161 PetscErrorCode ierr; 162 163 PetscFunctionBegin; 164 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 165 PetscValidHeaderSpecific(perm, IS_CLASSID, 2); 166 PetscValidPointer(pdm, 3); 167 ierr = DMCreate(PetscObjectComm((PetscObject) dm), pdm);CHKERRQ(ierr); 168 ierr = DMSetType(*pdm, DMPLEX);CHKERRQ(ierr); 169 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 170 ierr = DMSetDimension(*pdm, dim);CHKERRQ(ierr); 171 ierr = DMGetCoordinateDim(dm, &cdim);CHKERRQ(ierr); 172 ierr = DMSetCoordinateDim(*pdm, cdim);CHKERRQ(ierr); 173 ierr = DMCopyDisc(dm, *pdm);CHKERRQ(ierr); 174 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 175 if (section) { 176 ierr = PetscSectionPermute(section, perm, §ionNew);CHKERRQ(ierr); 177 ierr = DMSetLocalSection(*pdm, sectionNew);CHKERRQ(ierr); 178 ierr = PetscSectionDestroy(§ionNew);CHKERRQ(ierr); 179 } 180 plexNew = (DM_Plex *) (*pdm)->data; 181 /* Ignore ltogmap, ltogmapb */ 182 /* Ignore sf, sectionSF */ 183 /* Ignore globalVertexNumbers, globalCellNumbers */ 184 /* Reorder labels */ 185 { 186 PetscInt numLabels, l; 187 DMLabel label, labelNew; 188 189 ierr = DMGetNumLabels(dm, &numLabels);CHKERRQ(ierr); 190 for (l = 0; l < numLabels; ++l) { 191 ierr = DMGetLabelByNum(dm, l, &label);CHKERRQ(ierr); 192 ierr = DMLabelPermute(label, perm, &labelNew);CHKERRQ(ierr); 193 ierr = DMAddLabel(*pdm, labelNew);CHKERRQ(ierr); 194 ierr = DMLabelDestroy(&labelNew);CHKERRQ(ierr); 195 } 196 ierr = DMGetLabel(*pdm, "depth", &(*pdm)->depthLabel);CHKERRQ(ierr); 197 if (plex->subpointMap) {ierr = DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);CHKERRQ(ierr);} 198 } 199 /* Reorder topology */ 200 { 201 const PetscInt *pperm; 202 PetscInt maxConeSize, maxSupportSize, n, pStart, pEnd, p; 203 204 ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr); 205 plexNew->maxConeSize = maxConeSize; 206 plexNew->maxSupportSize = maxSupportSize; 207 ierr = PetscSectionDestroy(&plexNew->coneSection);CHKERRQ(ierr); 208 ierr = PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);CHKERRQ(ierr); 209 ierr = PetscSectionGetStorageSize(plexNew->coneSection, &n);CHKERRQ(ierr); 210 ierr = PetscMalloc1(n, &plexNew->cones);CHKERRQ(ierr); 211 ierr = PetscMalloc1(n, &plexNew->coneOrientations);CHKERRQ(ierr); 212 ierr = ISGetIndices(perm, &pperm);CHKERRQ(ierr); 213 ierr = PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);CHKERRQ(ierr); 214 for (p = pStart; p < pEnd; ++p) { 215 PetscInt dof, off, offNew, d; 216 217 ierr = PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);CHKERRQ(ierr); 218 ierr = PetscSectionGetOffset(plex->coneSection, p, &off);CHKERRQ(ierr); 219 ierr = PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);CHKERRQ(ierr); 220 for (d = 0; d < dof; ++d) { 221 plexNew->cones[offNew+d] = pperm[plex->cones[off+d]]; 222 plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d]; 223 } 224 } 225 ierr = PetscSectionDestroy(&plexNew->supportSection);CHKERRQ(ierr); 226 ierr = PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);CHKERRQ(ierr); 227 ierr = PetscSectionGetStorageSize(plexNew->supportSection, &n);CHKERRQ(ierr); 228 ierr = PetscMalloc1(n, &plexNew->supports);CHKERRQ(ierr); 229 ierr = PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);CHKERRQ(ierr); 230 for (p = pStart; p < pEnd; ++p) { 231 PetscInt dof, off, offNew, d; 232 233 ierr = PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);CHKERRQ(ierr); 234 ierr = PetscSectionGetOffset(plex->supportSection, p, &off);CHKERRQ(ierr); 235 ierr = PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);CHKERRQ(ierr); 236 for (d = 0; d < dof; ++d) { 237 plexNew->supports[offNew+d] = pperm[plex->supports[off+d]]; 238 } 239 } 240 ierr = ISRestoreIndices(perm, &pperm);CHKERRQ(ierr); 241 } 242 /* Remap coordinates */ 243 { 244 DM cdm, cdmNew; 245 PetscSection csection, csectionNew; 246 Vec coordinates, coordinatesNew; 247 PetscScalar *coords, *coordsNew; 248 const PetscInt *pperm; 249 PetscInt pStart, pEnd, p; 250 const char *name; 251 252 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 253 ierr = DMGetLocalSection(cdm, &csection);CHKERRQ(ierr); 254 ierr = PetscSectionPermute(csection, perm, &csectionNew);CHKERRQ(ierr); 255 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 256 ierr = VecDuplicate(coordinates, &coordinatesNew);CHKERRQ(ierr); 257 ierr = PetscObjectGetName((PetscObject)coordinates,&name);CHKERRQ(ierr); 258 ierr = PetscObjectSetName((PetscObject)coordinatesNew,name);CHKERRQ(ierr); 259 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 260 ierr = VecGetArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 261 ierr = PetscSectionGetChart(csectionNew, &pStart, &pEnd);CHKERRQ(ierr); 262 ierr = ISGetIndices(perm, &pperm);CHKERRQ(ierr); 263 for (p = pStart; p < pEnd; ++p) { 264 PetscInt dof, off, offNew, d; 265 266 ierr = PetscSectionGetDof(csectionNew, p, &dof);CHKERRQ(ierr); 267 ierr = PetscSectionGetOffset(csection, p, &off);CHKERRQ(ierr); 268 ierr = PetscSectionGetOffset(csectionNew, pperm[p], &offNew);CHKERRQ(ierr); 269 for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d]; 270 } 271 ierr = ISRestoreIndices(perm, &pperm);CHKERRQ(ierr); 272 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 273 ierr = VecRestoreArray(coordinatesNew, &coordsNew);CHKERRQ(ierr); 274 ierr = DMGetCoordinateDM(*pdm, &cdmNew);CHKERRQ(ierr); 275 ierr = DMSetLocalSection(cdmNew, csectionNew);CHKERRQ(ierr); 276 ierr = DMSetCoordinatesLocal(*pdm, coordinatesNew);CHKERRQ(ierr); 277 ierr = PetscSectionDestroy(&csectionNew);CHKERRQ(ierr); 278 ierr = VecDestroy(&coordinatesNew);CHKERRQ(ierr); 279 } 280 (*pdm)->setupcalled = PETSC_TRUE; 281 PetscFunctionReturn(0); 282 } 283