1 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <../src/sys/utils/hash.h> 3 #include <petsc-private/isimpl.h> 4 #include <petscsf.h> 5 6 /** hierarchy routines */ 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "DMPlexSetReferenceTree" 10 /*@ 11 DMPlexSetReferenceTree - set the reference tree for hierarchically non-conforming meshes. 12 13 Not collective 14 15 Input Parameters: 16 + dm - The DMPlex object 17 - ref - The reference tree DMPlex object 18 19 Level: intermediate 20 21 .seealso: DMPlexGetReferenceTree(), DMPlexCreateDefaultReferenceTree() 22 @*/ 23 PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref) 24 { 25 DM_Plex *mesh = (DM_Plex *)dm->data; 26 PetscErrorCode ierr; 27 28 PetscFunctionBegin; 29 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 30 PetscValidHeaderSpecific(ref, DM_CLASSID, 2); 31 ierr = PetscObjectReference((PetscObject)ref);CHKERRQ(ierr); 32 ierr = DMDestroy(&mesh->referenceTree);CHKERRQ(ierr); 33 mesh->referenceTree = ref; 34 PetscFunctionReturn(0); 35 } 36 37 #undef __FUNCT__ 38 #define __FUNCT__ "DMPlexGetReferenceTree" 39 /*@ 40 DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes. 41 42 Not collective 43 44 Input Parameters: 45 . dm - The DMPlex object 46 47 Output Parameters 48 . ref - The reference tree DMPlex object 49 50 Level: intermediate 51 52 .seealso: DMPlexSetReferenceTree(), DMPlexCreateDefaultReferenceTree() 53 @*/ 54 PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref) 55 { 56 DM_Plex *mesh = (DM_Plex *)dm->data; 57 58 PetscFunctionBegin; 59 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 60 PetscValidPointer(ref,2); 61 *ref = mesh->referenceTree; 62 PetscFunctionReturn(0); 63 } 64 65 #undef __FUNCT__ 66 #define __FUNCT__ "DMPlexCreateDefaultReferenceTree" 67 /*@ 68 DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement. 69 70 Collective on comm 71 72 Input Parameters: 73 + comm - the MPI communicator 74 . dim - the spatial dimension 75 - simplex - Flag for simplex, otherwise use a tensor-product cell 76 77 Output Parameters: 78 . ref - the reference tree DMPlex object 79 80 Level: intermediate 81 82 .keywords: reference cell 83 .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree() 84 @*/ 85 PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref) 86 { 87 DM K, Kref; 88 PetscInt p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs; 89 PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts; 90 DMLabel identity, identityRef; 91 PetscSection unionSection, unionConeSection, parentSection; 92 PetscScalar *unionCoords; 93 IS perm; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 /* create a reference element */ 98 ierr = DMPlexCreateReferenceCell(comm, dim, simplex, &K);CHKERRQ(ierr); 99 ierr = DMPlexCreateLabel(K, "identity");CHKERRQ(ierr); 100 ierr = DMPlexGetLabel(K, "identity", &identity);CHKERRQ(ierr); 101 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr); 102 for (p = pStart; p < pEnd; p++) { 103 ierr = DMLabelSetValue(identity, p, p);CHKERRQ(ierr); 104 } 105 /* refine it */ 106 ierr = DMRefine(K,comm,&Kref);CHKERRQ(ierr); 107 108 /* the reference tree is the union of these two, without duplicating 109 * points that appear in both */ 110 ierr = DMPlexGetLabel(Kref, "identity", &identityRef);CHKERRQ(ierr); 111 ierr = DMPlexGetChart(Kref, &pRefStart, &pRefEnd);CHKERRQ(ierr); 112 ierr = PetscSectionCreate(comm, &unionSection);CHKERRQ(ierr); 113 ierr = PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart));CHKERRQ(ierr); 114 /* count points that will go in the union */ 115 for (p = pStart; p < pEnd; p++) { 116 ierr = PetscSectionSetDof(unionSection, p - pStart, 1);CHKERRQ(ierr); 117 } 118 for (p = pRefStart; p < pRefEnd; p++) { 119 PetscInt q, qSize; 120 ierr = DMLabelGetValue(identityRef, p, &q);CHKERRQ(ierr); 121 ierr = DMLabelGetStratumSize(identityRef, q, &qSize);CHKERRQ(ierr); 122 if (qSize > 1) { 123 ierr = PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1);CHKERRQ(ierr); 124 } 125 } 126 ierr = PetscMalloc1((pEnd - pStart) + (pRefEnd - pRefStart),&permvals);CHKERRQ(ierr); 127 offset = 0; 128 /* stratify points in the union by topological dimension */ 129 for (d = 0; d <= dim; d++) { 130 PetscInt cStart, cEnd, c; 131 132 ierr = DMPlexGetHeightStratum(K, d, &cStart, &cEnd);CHKERRQ(ierr); 133 for (c = cStart; c < cEnd; c++) { 134 permvals[offset++] = c; 135 } 136 137 ierr = DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);CHKERRQ(ierr); 138 for (c = cStart; c < cEnd; c++) { 139 permvals[offset++] = c + (pEnd - pStart); 140 } 141 } 142 ierr = ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);CHKERRQ(ierr); 143 ierr = PetscSectionSetPermutation(unionSection,perm);CHKERRQ(ierr); 144 ierr = PetscSectionSetUp(unionSection);CHKERRQ(ierr); 145 ierr = PetscSectionGetStorageSize(unionSection,&numUnionPoints);CHKERRQ(ierr); 146 ierr = PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints);CHKERRQ(ierr); 147 /* count dimension points */ 148 for (d = 0; d <= dim; d++) { 149 PetscInt cStart, cOff, cOff2; 150 ierr = DMPlexGetHeightStratum(K,d,&cStart,NULL);CHKERRQ(ierr); 151 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff);CHKERRQ(ierr); 152 if (d < dim) { 153 ierr = DMPlexGetHeightStratum(K,d+1,&cStart,NULL);CHKERRQ(ierr); 154 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2);CHKERRQ(ierr); 155 } 156 else { 157 cOff2 = numUnionPoints; 158 } 159 numDimPoints[dim - d] = cOff2 - cOff; 160 } 161 ierr = PetscSectionCreate(comm, &unionConeSection);CHKERRQ(ierr); 162 ierr = PetscSectionSetChart(unionConeSection, 0, numUnionPoints);CHKERRQ(ierr); 163 /* count the cones in the union */ 164 for (p = pStart; p < pEnd; p++) { 165 PetscInt dof, uOff; 166 167 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 168 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 169 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 170 coneSizes[uOff] = dof; 171 } 172 for (p = pRefStart; p < pRefEnd; p++) { 173 PetscInt dof, uDof, uOff; 174 175 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 176 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 177 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 178 if (uDof) { 179 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 180 coneSizes[uOff] = dof; 181 } 182 } 183 ierr = PetscSectionSetUp(unionConeSection);CHKERRQ(ierr); 184 ierr = PetscSectionGetStorageSize(unionConeSection,&numCones);CHKERRQ(ierr); 185 ierr = PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations);CHKERRQ(ierr); 186 /* write the cones in the union */ 187 for (p = pStart; p < pEnd; p++) { 188 PetscInt dof, uOff, c, cOff; 189 const PetscInt *cone, *orientation; 190 191 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 192 ierr = DMPlexGetCone(K, p, &cone);CHKERRQ(ierr); 193 ierr = DMPlexGetConeOrientation(K, p, &orientation);CHKERRQ(ierr); 194 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 195 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 196 for (c = 0; c < dof; c++) { 197 PetscInt e, eOff; 198 e = cone[c]; 199 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 200 unionCones[cOff + c] = eOff; 201 unionOrientations[cOff + c] = orientation[c]; 202 } 203 } 204 for (p = pRefStart; p < pRefEnd; p++) { 205 PetscInt dof, uDof, uOff, c, cOff; 206 const PetscInt *cone, *orientation; 207 208 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 209 ierr = DMPlexGetCone(Kref, p, &cone);CHKERRQ(ierr); 210 ierr = DMPlexGetConeOrientation(Kref, p, &orientation);CHKERRQ(ierr); 211 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 212 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 213 if (uDof) { 214 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 215 for (c = 0; c < dof; c++) { 216 PetscInt e, eOff, eDof; 217 218 e = cone[c]; 219 ierr = PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof);CHKERRQ(ierr); 220 if (eDof) { 221 ierr = PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);CHKERRQ(ierr); 222 } 223 else { 224 ierr = DMLabelGetValue(identityRef, e, &e);CHKERRQ(ierr); 225 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 226 } 227 unionCones[cOff + c] = eOff; 228 unionOrientations[cOff + c] = orientation[c]; 229 } 230 } 231 } 232 /* get the coordinates */ 233 { 234 PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff; 235 PetscSection KcoordsSec, KrefCoordsSec; 236 Vec KcoordsVec, KrefCoordsVec; 237 PetscScalar *Kcoords; 238 239 DMGetCoordinateSection(K, &KcoordsSec);CHKERRQ(ierr); 240 DMGetCoordinatesLocal(K, &KcoordsVec);CHKERRQ(ierr); 241 DMGetCoordinateSection(Kref, &KrefCoordsSec);CHKERRQ(ierr); 242 DMGetCoordinatesLocal(Kref, &KrefCoordsVec);CHKERRQ(ierr); 243 244 numVerts = numDimPoints[0]; 245 ierr = PetscMalloc1(numVerts * dim,&unionCoords);CHKERRQ(ierr); 246 ierr = DMPlexGetDepthStratum(K,0,&vStart,&vEnd);CHKERRQ(ierr); 247 248 offset = 0; 249 for (v = vStart; v < vEnd; v++) { 250 ierr = PetscSectionGetOffset(unionSection,v - pStart,&vOff);CHKERRQ(ierr); 251 ierr = VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);CHKERRQ(ierr); 252 for (d = 0; d < dim; d++) { 253 unionCoords[offset * dim + d] = Kcoords[d]; 254 } 255 offset++; 256 } 257 ierr = DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd);CHKERRQ(ierr); 258 for (v = vRefStart; v < vRefEnd; v++) { 259 ierr = PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof);CHKERRQ(ierr); 260 ierr = PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff);CHKERRQ(ierr); 261 ierr = VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);CHKERRQ(ierr); 262 if (vDof) { 263 for (d = 0; d < dim; d++) { 264 unionCoords[offset * dim + d] = Kcoords[d]; 265 } 266 offset++; 267 } 268 } 269 } 270 ierr = DMCreate(comm,ref);CHKERRQ(ierr); 271 ierr = DMSetType(*ref,DMPLEX);CHKERRQ(ierr); 272 ierr = DMPlexSetDimension(*ref,dim);CHKERRQ(ierr); 273 ierr = DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords);CHKERRQ(ierr); 274 /* set the tree */ 275 ierr = PetscSectionCreate(comm,&parentSection);CHKERRQ(ierr); 276 ierr = PetscSectionSetChart(parentSection,0,numUnionPoints);CHKERRQ(ierr); 277 for (p = pRefStart; p < pRefEnd; p++) { 278 PetscInt uDof, uOff; 279 280 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 281 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 282 if (uDof) { 283 PetscSectionSetDof(parentSection,uOff,1);CHKERRQ(ierr); 284 } 285 } 286 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 287 ierr = PetscSectionGetStorageSize(parentSection,&parentSize);CHKERRQ(ierr); 288 ierr = PetscMalloc2(parentSize,&parents,parentSize,&childIDs);CHKERRQ(ierr); 289 for (p = pRefStart; p < pRefEnd; p++) { 290 PetscInt uDof, uOff; 291 292 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 293 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 294 if (uDof) { 295 PetscInt pOff, parent, parentU; 296 PetscSectionGetOffset(parentSection,uOff,&pOff);CHKERRQ(ierr); 297 DMLabelGetValue(identityRef,p,&parent);CHKERRQ(ierr); 298 ierr = PetscSectionGetOffset(unionSection, parent - pStart,&parentU);CHKERRQ(ierr); 299 parents[pOff] = parentU; 300 childIDs[pOff] = uOff; 301 } 302 } 303 ierr = DMPlexSetTree(*ref,parentSection,parents,childIDs);CHKERRQ(ierr); 304 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 305 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 306 307 /* clean up */ 308 ierr = PetscSectionDestroy(&unionSection);CHKERRQ(ierr); 309 ierr = PetscSectionDestroy(&unionConeSection);CHKERRQ(ierr); 310 ierr = ISDestroy(&perm);CHKERRQ(ierr); 311 ierr = PetscFree(unionCoords);CHKERRQ(ierr); 312 ierr = PetscFree2(unionCones,unionOrientations);CHKERRQ(ierr); 313 ierr = PetscFree2(coneSizes,numDimPoints);CHKERRQ(ierr); 314 ierr = DMDestroy(&K);CHKERRQ(ierr); 315 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 316 PetscFunctionReturn(0); 317 } 318 319 #undef __FUNCT__ 320 #define __FUNCT__ "DMPlexTreeSymmetrize" 321 static PetscErrorCode DMPlexTreeSymmetrize(DM dm) 322 { 323 DM_Plex *mesh = (DM_Plex *)dm->data; 324 PetscSection childSec, pSec; 325 PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT; 326 PetscInt *offsets, *children, pStart, pEnd; 327 PetscErrorCode ierr; 328 329 PetscFunctionBegin; 330 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 331 ierr = PetscSectionDestroy(&mesh->childSection);CHKERRQ(ierr); 332 ierr = PetscFree(mesh->children);CHKERRQ(ierr); 333 pSec = mesh->parentSection; 334 if (!pSec) PetscFunctionReturn(0); 335 ierr = PetscSectionGetStorageSize(pSec,&pSize);CHKERRQ(ierr); 336 for (p = 0; p < pSize; p++) { 337 PetscInt par = mesh->parents[p]; 338 339 parMax = PetscMax(parMax,par+1); 340 parMin = PetscMin(parMin,par); 341 } 342 if (parMin > parMax) { 343 parMin = -1; 344 parMax = -1; 345 } 346 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);CHKERRQ(ierr); 347 ierr = PetscSectionSetChart(childSec,parMin,parMax);CHKERRQ(ierr); 348 for (p = 0; p < pSize; p++) { 349 PetscInt par = mesh->parents[p]; 350 351 ierr = PetscSectionAddDof(childSec,par,1);CHKERRQ(ierr); 352 } 353 ierr = PetscSectionSetUp(childSec);CHKERRQ(ierr); 354 ierr = PetscSectionGetStorageSize(childSec,&cSize);CHKERRQ(ierr); 355 ierr = PetscMalloc1(cSize,&children);CHKERRQ(ierr); 356 ierr = PetscCalloc1(parMax-parMin,&offsets);CHKERRQ(ierr); 357 ierr = PetscSectionGetChart(pSec,&pStart,&pEnd);CHKERRQ(ierr); 358 for (p = pStart; p < pEnd; p++) { 359 PetscInt dof, off, i; 360 361 ierr = PetscSectionGetDof(pSec,p,&dof);CHKERRQ(ierr); 362 ierr = PetscSectionGetOffset(pSec,p,&off);CHKERRQ(ierr); 363 for (i = 0; i < dof; i++) { 364 PetscInt par = mesh->parents[off + i], cOff; 365 366 ierr = PetscSectionGetOffset(childSec,par,&cOff);CHKERRQ(ierr); 367 children[cOff + offsets[par-parMin]++] = p; 368 } 369 } 370 mesh->childSection = childSec; 371 mesh->children = children; 372 ierr = PetscFree(offsets);CHKERRQ(ierr); 373 PetscFunctionReturn(0); 374 } 375 376 #undef __FUNCT__ 377 #define __FUNCT__ "AnchorsFlatten" 378 static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew) 379 { 380 PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL; 381 const PetscInt *vals; 382 PetscSection secNew; 383 PetscBool anyNew, globalAnyNew; 384 PetscBool compress; 385 PetscErrorCode ierr; 386 387 PetscFunctionBegin; 388 ierr = PetscSectionGetChart(section,&pStart,&pEnd);CHKERRQ(ierr); 389 ierr = ISGetLocalSize(is,&size);CHKERRQ(ierr); 390 ierr = ISGetIndices(is,&vals);CHKERRQ(ierr); 391 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);CHKERRQ(ierr); 392 ierr = PetscSectionSetChart(secNew,pStart,pEnd);CHKERRQ(ierr); 393 for (i = 0; i < size; i++) { 394 PetscInt dof; 395 396 p = vals[i]; 397 if (p < pStart || p >= pEnd) continue; 398 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 399 if (dof) break; 400 } 401 if (i == size) { 402 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 403 anyNew = PETSC_FALSE; 404 compress = PETSC_FALSE; 405 sizeNew = 0; 406 } 407 else { 408 anyNew = PETSC_TRUE; 409 for (p = pStart; p < pEnd; p++) { 410 PetscInt dof, off; 411 412 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 413 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 414 for (i = 0; i < dof; i++) { 415 PetscInt q = vals[off + i], qDof = 0; 416 417 if (q >= pStart && q < pEnd) { 418 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 419 } 420 if (qDof) { 421 ierr = PetscSectionAddDof(secNew, p, qDof);CHKERRQ(ierr); 422 } 423 else { 424 ierr = PetscSectionAddDof(secNew, p, 1);CHKERRQ(ierr); 425 } 426 } 427 } 428 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 429 ierr = PetscSectionGetStorageSize(secNew,&sizeNew);CHKERRQ(ierr); 430 ierr = PetscMalloc1(sizeNew,&valsNew);CHKERRQ(ierr); 431 compress = PETSC_FALSE; 432 for (p = pStart; p < pEnd; p++) { 433 PetscInt dof, off, count, offNew, dofNew; 434 435 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 436 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 437 ierr = PetscSectionGetDof(secNew, p, &dofNew);CHKERRQ(ierr); 438 ierr = PetscSectionGetOffset(secNew, p, &offNew);CHKERRQ(ierr); 439 count = 0; 440 for (i = 0; i < dof; i++) { 441 PetscInt q = vals[off + i], qDof = 0, qOff = 0, j; 442 443 if (q >= pStart && q < pEnd) { 444 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 445 ierr = PetscSectionGetOffset(section, q, &qOff);CHKERRQ(ierr); 446 } 447 if (qDof) { 448 PetscInt oldCount = count; 449 450 for (j = 0; j < qDof; j++) { 451 PetscInt k, r = vals[qOff + j]; 452 453 for (k = 0; k < oldCount; k++) { 454 if (valsNew[offNew + k] == r) { 455 break; 456 } 457 } 458 if (k == oldCount) { 459 valsNew[offNew + count++] = r; 460 } 461 } 462 } 463 else { 464 PetscInt k, oldCount = count; 465 466 for (k = 0; k < oldCount; k++) { 467 if (valsNew[offNew + k] == q) { 468 break; 469 } 470 } 471 if (k == oldCount) { 472 valsNew[offNew + count++] = q; 473 } 474 } 475 } 476 if (count < dofNew) { 477 ierr = PetscSectionSetDof(secNew, p, count);CHKERRQ(ierr); 478 compress = PETSC_TRUE; 479 } 480 } 481 } 482 ierr = ISRestoreIndices(is,&vals);CHKERRQ(ierr); 483 ierr = MPI_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 484 if (!globalAnyNew) { 485 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 486 *sectionNew = NULL; 487 *isNew = NULL; 488 } 489 else { 490 PetscBool globalCompress; 491 492 ierr = MPI_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 493 if (compress) { 494 PetscSection secComp; 495 PetscInt *valsComp = NULL; 496 497 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);CHKERRQ(ierr); 498 ierr = PetscSectionSetChart(secComp,pStart,pEnd);CHKERRQ(ierr); 499 for (p = pStart; p < pEnd; p++) { 500 PetscInt dof; 501 502 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 503 ierr = PetscSectionSetDof(secComp, p, dof);CHKERRQ(ierr); 504 } 505 ierr = PetscSectionSetUp(secComp);CHKERRQ(ierr); 506 ierr = PetscSectionGetStorageSize(secComp,&sizeNew);CHKERRQ(ierr); 507 ierr = PetscMalloc1(sizeNew,&valsComp);CHKERRQ(ierr); 508 for (p = pStart; p < pEnd; p++) { 509 PetscInt dof, off, offNew, j; 510 511 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 512 ierr = PetscSectionGetOffset(secNew, p, &off);CHKERRQ(ierr); 513 ierr = PetscSectionGetOffset(secComp, p, &offNew);CHKERRQ(ierr); 514 for (j = 0; j < dof; j++) { 515 valsComp[offNew + j] = valsNew[off + j]; 516 } 517 } 518 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 519 secNew = secComp; 520 ierr = PetscFree(valsNew);CHKERRQ(ierr); 521 valsNew = valsComp; 522 } 523 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);CHKERRQ(ierr); 524 } 525 PetscFunctionReturn(0); 526 } 527 528 #undef __FUNCT__ 529 #define __FUNCT__ "DMPlexComputeConstraints_Tree" 530 static PetscErrorCode DMPlexComputeConstraints_Tree(DM dm) 531 { 532 PetscInt p, pStart, pEnd, *anchors, size; 533 PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT; 534 PetscSection aSec; 535 IS aIS; 536 PetscErrorCode ierr; 537 538 PetscFunctionBegin; 539 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 540 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 541 for (p = pStart; p < pEnd; p++) { 542 PetscInt parent; 543 544 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 545 if (parent != p) { 546 aMin = PetscMin(aMin,p); 547 aMax = PetscMax(aMax,p+1); 548 } 549 } 550 if (aMin > aMax) { 551 aMin = -1; 552 aMax = -1; 553 } 554 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&aSec);CHKERRQ(ierr); 555 ierr = PetscSectionSetChart(aSec,aMin,aMax);CHKERRQ(ierr); 556 for (p = aMin; p < aMax; p++) { 557 PetscInt parent, ancestor = p; 558 559 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 560 while (parent != ancestor) { 561 ancestor = parent; 562 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 563 } 564 if (ancestor != p) { 565 PetscInt closureSize, *closure = NULL; 566 567 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 568 ierr = PetscSectionSetDof(aSec,p,closureSize);CHKERRQ(ierr); 569 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 570 } 571 } 572 ierr = PetscSectionSetUp(aSec);CHKERRQ(ierr); 573 ierr = PetscSectionGetStorageSize(aSec,&size);CHKERRQ(ierr); 574 ierr = PetscMalloc1(size,&anchors);CHKERRQ(ierr); 575 for (p = aMin; p < aMax; p++) { 576 PetscInt parent, ancestor = p; 577 578 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 579 while (parent != ancestor) { 580 ancestor = parent; 581 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 582 } 583 if (ancestor != p) { 584 PetscInt j, closureSize, *closure = NULL, aOff; 585 586 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 587 588 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 589 for (j = 0; j < closureSize; j++) { 590 anchors[aOff + j] = closure[2*j]; 591 } 592 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 593 } 594 } 595 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm),size,anchors,PETSC_OWN_POINTER,&aIS);CHKERRQ(ierr); 596 { 597 PetscSection aSecNew = aSec; 598 IS aISNew = aIS; 599 600 ierr = PetscObjectReference((PetscObject)aSec);CHKERRQ(ierr); 601 ierr = PetscObjectReference((PetscObject)aIS);CHKERRQ(ierr); 602 while (aSecNew) { 603 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 604 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 605 aSec = aSecNew; 606 aIS = aISNew; 607 aSecNew = NULL; 608 aISNew = NULL; 609 ierr = AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);CHKERRQ(ierr); 610 } 611 } 612 ierr = DMPlexSetConstraints(dm,aSec,aIS);CHKERRQ(ierr); 613 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 614 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 615 PetscFunctionReturn(0); 616 } 617 618 #undef __FUNCT__ 619 #define __FUNCT__ "DMPlexSetTree" 620 /*@ 621 DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates 622 the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its 623 tree root. 624 625 Collective on dm 626 627 Input Parameters: 628 + dm - the DMPlex object 629 . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 630 offset indexes the parent and childID list; the reference count of parentSection is incremented 631 . parents - a list of the point parents; copied, can be destroyed 632 - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 633 the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed 634 635 Level: intermediate 636 637 .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetConstraints(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 638 @*/ 639 PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[]) 640 { 641 DM_Plex *mesh = (DM_Plex *)dm->data; 642 PetscInt size; 643 PetscErrorCode ierr; 644 645 PetscFunctionBegin; 646 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 647 PetscValidHeaderSpecific(parentSection, PETSC_SECTION_CLASSID, 2); 648 ierr = PetscObjectReference((PetscObject)parentSection);CHKERRQ(ierr); 649 ierr = PetscSectionDestroy(&mesh->parentSection);CHKERRQ(ierr); 650 mesh->parentSection = parentSection; 651 ierr = PetscSectionGetStorageSize(parentSection,&size);CHKERRQ(ierr); 652 if (parents != mesh->parents) { 653 ierr = PetscFree(mesh->parents);CHKERRQ(ierr); 654 ierr = PetscMalloc1(size,&mesh->parents);CHKERRQ(ierr); 655 ierr = PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));CHKERRQ(ierr); 656 } 657 if (childIDs != mesh->childIDs) { 658 ierr = PetscFree(mesh->childIDs);CHKERRQ(ierr); 659 ierr = PetscMalloc1(size,&mesh->childIDs);CHKERRQ(ierr); 660 ierr = PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));CHKERRQ(ierr); 661 } 662 ierr = DMPlexTreeSymmetrize(dm);CHKERRQ(ierr); 663 ierr = DMPlexComputeConstraints_Tree(dm);CHKERRQ(ierr); 664 PetscFunctionReturn(0); 665 } 666 667 #undef __FUNCT__ 668 #define __FUNCT__ "DMPlexGetTree" 669 /*@ 670 DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points. 671 Collective on dm 672 673 Input Parameters: 674 . dm - the DMPlex object 675 676 Output Parameters: 677 + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 678 offset indexes the parent and childID list 679 . parents - a list of the point parents 680 . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 681 the child corresponds to the point in the reference tree with index childID 682 . childSection - the inverse of the parent section 683 - children - a list of the point children 684 685 Level: intermediate 686 687 .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetConstraints(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 688 @*/ 689 PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[]) 690 { 691 DM_Plex *mesh = (DM_Plex *)dm->data; 692 693 PetscFunctionBegin; 694 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 695 if (parentSection) *parentSection = mesh->parentSection; 696 if (parents) *parents = mesh->parents; 697 if (childIDs) *childIDs = mesh->childIDs; 698 if (childSection) *childSection = mesh->childSection; 699 if (children) *children = mesh->children; 700 PetscFunctionReturn(0); 701 } 702 703 #undef __FUNCT__ 704 #define __FUNCT__ "DMPlexGetTreeParent" 705 /*@ 706 DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the Sieve DAG) 707 708 Input Parameters: 709 + dm - the DMPlex object 710 - point - the query point 711 712 Output Parameters: 713 + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent 714 - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point 715 does not have a parent 716 717 Level: intermediate 718 719 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren() 720 @*/ 721 PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID) 722 { 723 DM_Plex *mesh = (DM_Plex *)dm->data; 724 PetscSection pSec; 725 PetscErrorCode ierr; 726 727 PetscFunctionBegin; 728 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 729 pSec = mesh->parentSection; 730 if (pSec && point >= pSec->pStart && point < pSec->pEnd) { 731 PetscInt dof; 732 733 ierr = PetscSectionGetDof (pSec, point, &dof);CHKERRQ(ierr); 734 if (dof) { 735 PetscInt off; 736 737 ierr = PetscSectionGetOffset (pSec, point, &off);CHKERRQ(ierr); 738 if (parent) *parent = mesh->parents[off]; 739 if (childID) *childID = mesh->childIDs[off]; 740 PetscFunctionReturn(0); 741 } 742 } 743 if (parent) { 744 *parent = point; 745 } 746 if (childID) { 747 *childID = 0; 748 } 749 PetscFunctionReturn(0); 750 } 751 752 #undef __FUNCT__ 753 #define __FUNCT__ "DMPlexGetTreeChildren" 754 /*@C 755 DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the Sieve DAG) 756 757 Input Parameters: 758 + dm - the DMPlex object 759 - point - the query point 760 761 Output Parameters: 762 + numChildren - if not NULL, set to the number of children 763 - children - if not NULL, set to a list children, or set to NULL if the point has no children 764 765 Level: intermediate 766 767 Fortran Notes: 768 Since it returns an array, this routine is only available in Fortran 90, and you must 769 include petsc.h90 in your code. 770 771 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent() 772 @*/ 773 PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[]) 774 { 775 DM_Plex *mesh = (DM_Plex *)dm->data; 776 PetscSection childSec; 777 PetscInt dof = 0; 778 PetscErrorCode ierr; 779 780 PetscFunctionBegin; 781 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 782 childSec = mesh->childSection; 783 if (childSec && point >= childSec->pStart && point < childSec->pEnd) { 784 ierr = PetscSectionGetDof (childSec, point, &dof);CHKERRQ(ierr); 785 } 786 if (numChildren) *numChildren = dof; 787 if (children) { 788 if (dof) { 789 PetscInt off; 790 791 ierr = PetscSectionGetOffset (childSec, point, &off);CHKERRQ(ierr); 792 *children = &mesh->children[off]; 793 } 794 else { 795 *children = NULL; 796 } 797 } 798 PetscFunctionReturn(0); 799 } 800