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__ "DMPlexSetTree" 321 /*@ 322 DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates 323 the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its 324 tree root. 325 326 Collective on dm 327 328 Input Parameters: 329 + dm - the DMPlex object 330 . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 331 offset indexes the parent and childID list; the reference count of parentSection is incremented 332 . parents - a list of the point parents; copied, can be destroyed 333 - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 334 the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed 335 336 Level: intermediate 337 338 .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetConstraints(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 339 @*/ 340 PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[]) 341 { 342 DM_Plex *mesh = (DM_Plex *)dm->data; 343 PetscInt size; 344 PetscErrorCode ierr; 345 346 PetscFunctionBegin; 347 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 348 PetscValidHeaderSpecific(parentSection, PETSC_SECTION_CLASSID, 2); 349 ierr = PetscObjectReference((PetscObject)parentSection);CHKERRQ(ierr); 350 ierr = PetscSectionDestroy(&mesh->parentSection);CHKERRQ(ierr); 351 mesh->parentSection = parentSection; 352 ierr = PetscSectionGetStorageSize(parentSection,&size);CHKERRQ(ierr); 353 if (parents != mesh->parents) { 354 ierr = PetscFree(mesh->parents);CHKERRQ(ierr); 355 ierr = PetscMalloc1(size,&mesh->parents);CHKERRQ(ierr); 356 ierr = PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));CHKERRQ(ierr); 357 } 358 if (childIDs != mesh->childIDs) { 359 ierr = PetscFree(mesh->childIDs);CHKERRQ(ierr); 360 ierr = PetscMalloc1(size,&mesh->childIDs);CHKERRQ(ierr); 361 ierr = PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));CHKERRQ(ierr); 362 } 363 PetscFunctionReturn(0); 364 } 365 366 #undef __FUNCT__ 367 #define __FUNCT__ "DMPlexGetTree" 368 /*@ 369 DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points. 370 Collective on dm 371 372 Input Parameters: 373 . dm - the DMPlex object 374 375 Output Parameters: 376 + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 377 offset indexes the parent and childID list 378 . parents - a list of the point parents 379 . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 380 the child corresponds to the point in the reference tree with index childID 381 . childSection - the inverse of the parent section 382 - children - a list of the point children 383 384 Level: intermediate 385 386 .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetConstraints(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 387 @*/ 388 PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[]) 389 { 390 DM_Plex *mesh = (DM_Plex *)dm->data; 391 392 PetscFunctionBegin; 393 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 394 if (parentSection) *parentSection = mesh->parentSection; 395 if (parents) *parents = mesh->parents; 396 if (childIDs) *childIDs = mesh->childIDs; 397 if (childSection) *childSection = mesh->childSection; 398 if (children) *children = mesh->children; 399 PetscFunctionReturn(0); 400 } 401 402 #undef __FUNCT__ 403 #define __FUNCT__ "DMPlexGetTreeParent" 404 /*@ 405 DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the Sieve DAG) 406 407 Input Parameters: 408 + dm - the DMPlex object 409 - point - the query point 410 411 Output Parameters: 412 + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent 413 - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point 414 does not have a parent 415 416 Level: intermediate 417 418 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren() 419 @*/ 420 PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID) 421 { 422 DM_Plex *mesh = (DM_Plex *)dm->data; 423 PetscSection pSec; 424 PetscErrorCode ierr; 425 426 PetscFunctionBegin; 427 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 428 pSec = mesh->parentSection; 429 if (pSec && point >= pSec->pStart && point < pSec->pEnd) { 430 PetscInt dof; 431 432 ierr = PetscSectionGetDof (pSec, point, &dof);CHKERRQ(ierr); 433 if (dof) { 434 PetscInt off; 435 436 ierr = PetscSectionGetOffset (pSec, point, &off);CHKERRQ(ierr); 437 if (parent) *parent = mesh->parents[off]; 438 if (childID) *childID = mesh->childIDs[off]; 439 PetscFunctionReturn(0); 440 } 441 } 442 if (parent) { 443 *parent = point; 444 } 445 if (childID) { 446 *childID = 0; 447 } 448 PetscFunctionReturn(0); 449 } 450 451 #undef __FUNCT__ 452 #define __FUNCT__ "DMPlexGetTreeChildren" 453 /*@C 454 DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the Sieve DAG) 455 456 Input Parameters: 457 + dm - the DMPlex object 458 - point - the query point 459 460 Output Parameters: 461 + numChildren - if not NULL, set to the number of children 462 - children - if not NULL, set to a list children, or set to NULL if the point has no children 463 464 Level: intermediate 465 466 Fortran Notes: 467 Since it returns an array, this routine is only available in Fortran 90, and you must 468 include petsc.h90 in your code. 469 470 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent() 471 @*/ 472 PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[]) 473 { 474 DM_Plex *mesh = (DM_Plex *)dm->data; 475 PetscSection childSec; 476 PetscInt dof = 0; 477 PetscErrorCode ierr; 478 479 PetscFunctionBegin; 480 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 481 childSec = mesh->childSection; 482 if (childSec && point >= childSec->pStart && point < childSec->pEnd) { 483 ierr = PetscSectionGetDof (childSec, point, &dof);CHKERRQ(ierr); 484 } 485 if (numChildren) *numChildren = dof; 486 if (children) { 487 if (dof) { 488 PetscInt off; 489 490 ierr = PetscSectionGetOffset (childSec, point, &off);CHKERRQ(ierr); 491 *children = &mesh->children[off]; 492 } 493 else { 494 *children = NULL; 495 } 496 } 497 PetscFunctionReturn(0); 498 } 499