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; 89 PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts; 90 DMLabel identity, identityRef; 91 PetscSection unionSection, unionConeSection; 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 /* clean up */ 275 ierr = PetscSectionDestroy(&unionSection);CHKERRQ(ierr); 276 ierr = PetscSectionDestroy(&unionConeSection);CHKERRQ(ierr); 277 ierr = ISDestroy(&perm);CHKERRQ(ierr); 278 ierr = PetscFree(unionCoords);CHKERRQ(ierr); 279 ierr = PetscFree2(unionCones,unionOrientations);CHKERRQ(ierr); 280 ierr = PetscFree2(coneSizes,numDimPoints);CHKERRQ(ierr); 281 ierr = DMDestroy(&K);CHKERRQ(ierr); 282 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 283 PetscFunctionReturn(0); 284 } 285 286 #undef __FUNCT__ 287 #define __FUNCT__ "DMPlexSetTree" 288 /*@ 289 DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates 290 the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its 291 tree root. 292 293 Collective on dm 294 295 Input Parameters: 296 + dm - the DMPlex object 297 . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 298 offset indexes the parent and childID list; the reference count of parentSection is incremented 299 . parents - a list of the point parents; copied, can be destroyed 300 - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 301 the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed 302 303 Level: intermediate 304 305 .seealso: DMPlexSetReferenceTree(), DMPlexSetConstraints(), DMPlexGetTreeParent() 306 @*/ 307 PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs) 308 { 309 DM_Plex *mesh = (DM_Plex *)dm->data; 310 PetscInt size; 311 PetscErrorCode ierr; 312 313 PetscFunctionBegin; 314 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 315 PetscValidHeaderSpecific(parentSection, PETSC_SECTION_CLASSID, 2); 316 ierr = PetscObjectReference((PetscObject)parentSection);CHKERRQ(ierr); 317 ierr = PetscSectionDestroy(&mesh->parentSection);CHKERRQ(ierr); 318 mesh->parentSection = parentSection; 319 ierr = PetscSectionGetStorageSize(parentSection,&size);CHKERRQ(ierr); 320 if (parents != mesh->parents) { 321 ierr = PetscFree(mesh->parents);CHKERRQ(ierr); 322 ierr = PetscMalloc1(size,&mesh->parents);CHKERRQ(ierr); 323 ierr = PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));CHKERRQ(ierr); 324 } 325 if (childIDs != mesh->childIDs) { 326 ierr = PetscFree(mesh->childIDs);CHKERRQ(ierr); 327 ierr = PetscMalloc1(size,&mesh->childIDs);CHKERRQ(ierr); 328 ierr = PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));CHKERRQ(ierr); 329 } 330 PetscFunctionReturn(0); 331 } 332 333 #undef __FUNCT__ 334 #define __FUNCT__ "DMPlexGetTreeParent" 335 /*@ 336 DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the Sieve DAG) 337 338 Input Parameters: 339 + dm - the DMPlex object 340 - point - the query point 341 342 Output Parameters: 343 + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent 344 - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point 345 does not have a parent 346 347 Level: intermediate 348 349 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren() 350 @*/ 351 PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID) 352 { 353 DM_Plex *mesh = (DM_Plex *)dm->data; 354 PetscSection pSec; 355 PetscErrorCode ierr; 356 357 PetscFunctionBegin; 358 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 359 pSec = mesh->parentSection; 360 if (pSec && point >= pSec->pStart && point < pSec->pEnd) { 361 PetscInt dof; 362 363 ierr = PetscSectionGetDof (pSec, point, &dof);CHKERRQ(ierr); 364 if (dof) { 365 PetscInt off; 366 367 ierr = PetscSectionGetOffset (pSec, point, &off);CHKERRQ(ierr); 368 if (parent) *parent = mesh->parents[off]; 369 if (childID) *childID = mesh->childIDs[off]; 370 PetscFunctionReturn(0); 371 } 372 } 373 if (parent) { 374 *parent = point; 375 } 376 if (childID) { 377 *childID = 0; 378 } 379 PetscFunctionReturn(0); 380 } 381 382 #undef __FUNCT__ 383 #define __FUNCT__ "DMPlexGetTreeChildren" 384 /*@C 385 DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the Sieve DAG) 386 387 Input Parameters: 388 + dm - the DMPlex object 389 - point - the query point 390 391 Output Parameters: 392 + numChildren - if not NULL, set to the number of children 393 - children - if not NULL, set to a list children, or set to NULL if the point has no children 394 395 Level: intermediate 396 397 Fortran Notes: 398 Since it returns an array, this routine is only available in Fortran 90, and you must 399 include petsc.h90 in your code. 400 401 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent() 402 @*/ 403 PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[]) 404 { 405 DM_Plex *mesh = (DM_Plex *)dm->data; 406 PetscSection childSec; 407 PetscInt dof = 0; 408 PetscErrorCode ierr; 409 410 PetscFunctionBegin; 411 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 412 childSec = mesh->childSection; 413 if (childSec && point >= childSec->pStart && point < childSec->pEnd) { 414 ierr = PetscSectionGetDof (childSec, point, &dof);CHKERRQ(ierr); 415 } 416 if (numChildren) *numChildren = dof; 417 if (children) { 418 if (dof) { 419 PetscInt off; 420 421 ierr = PetscSectionGetOffset (childSec, point, &off);CHKERRQ(ierr); 422 *children = &mesh->children[off]; 423 } 424 else { 425 *children = NULL; 426 } 427 } 428 PetscFunctionReturn(0); 429 } 430