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 <petsc-private/petscfeimpl.h> 5 #include <petscsf.h> 6 #include <petscds.h> 7 8 /** hierarchy routines */ 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "DMPlexSetReferenceTree" 12 /*@ 13 DMPlexSetReferenceTree - set the reference tree for hierarchically non-conforming meshes. 14 15 Not collective 16 17 Input Parameters: 18 + dm - The DMPlex object 19 - ref - The reference tree DMPlex object 20 21 Level: intermediate 22 23 .seealso: DMPlexGetReferenceTree(), DMPlexCreateDefaultReferenceTree() 24 @*/ 25 PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref) 26 { 27 DM_Plex *mesh = (DM_Plex *)dm->data; 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 32 PetscValidHeaderSpecific(ref, DM_CLASSID, 2); 33 ierr = PetscObjectReference((PetscObject)ref);CHKERRQ(ierr); 34 ierr = DMDestroy(&mesh->referenceTree);CHKERRQ(ierr); 35 mesh->referenceTree = ref; 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "DMPlexGetReferenceTree" 41 /*@ 42 DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes. 43 44 Not collective 45 46 Input Parameters: 47 . dm - The DMPlex object 48 49 Output Parameters 50 . ref - The reference tree DMPlex object 51 52 Level: intermediate 53 54 .seealso: DMPlexSetReferenceTree(), DMPlexCreateDefaultReferenceTree() 55 @*/ 56 PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref) 57 { 58 DM_Plex *mesh = (DM_Plex *)dm->data; 59 60 PetscFunctionBegin; 61 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 62 PetscValidPointer(ref,2); 63 *ref = mesh->referenceTree; 64 PetscFunctionReturn(0); 65 } 66 67 #undef __FUNCT__ 68 #define __FUNCT__ "DMPlexReferenceTreeGetChildSymmetry_Default" 69 static PetscErrorCode DMPlexReferenceTreeGetChildSymmetry_Default(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB) 70 { 71 PetscInt coneSize, dStart, dEnd, dim, ABswap, oAvert, oBvert, ABswapVert; 72 PetscErrorCode ierr; 73 74 PetscFunctionBegin; 75 if (parentOrientA == parentOrientB) { 76 if (childOrientB) *childOrientB = childOrientA; 77 if (childB) *childB = childA; 78 PetscFunctionReturn(0); 79 } 80 for (dim = 0; dim < 3; dim++) { 81 ierr = DMPlexGetDepthStratum(dm,dim,&dStart,&dEnd);CHKERRQ(ierr); 82 if (parent >= dStart && parent <= dEnd) { 83 break; 84 } 85 } 86 if (dim > 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot perform child symmetry for %d-cells",dim); 87 if (!dim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A vertex has no children"); 88 if (childA < dStart || childA >= dEnd) { 89 /* this is a lower-dimensional child: bootstrap */ 90 PetscInt size, i, sA = -1, sB, sOrientB, sConeSize; 91 const PetscInt *supp, *coneA, *coneB, *oA, *oB; 92 93 ierr = DMPlexGetSupportSize(dm,childA,&size);CHKERRQ(ierr); 94 ierr = DMPlexGetSupport(dm,childA,&supp);CHKERRQ(ierr); 95 96 /* find a point sA in supp(childA) that has the same parent */ 97 for (i = 0; i < size; i++) { 98 PetscInt sParent; 99 100 sA = supp[i]; 101 if (sA == parent) continue; 102 ierr = DMPlexGetTreeParent(dm,sA,&sParent,NULL);CHKERRQ(ierr); 103 if (sParent == parent) { 104 break; 105 } 106 } 107 if (i == size) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"could not find support in children"); 108 /* find out which point sB is in an equivalent position to sA under 109 * parentOrientB */ 110 ierr = DMPlexReferenceTreeGetChildSymmetry_Default(dm,parent,parentOrientA,0,sA,parentOrientB,&sOrientB,&sB);CHKERRQ(ierr); 111 ierr = DMPlexGetConeSize(dm,sA,&sConeSize);CHKERRQ(ierr); 112 ierr = DMPlexGetCone(dm,sA,&coneA);CHKERRQ(ierr); 113 ierr = DMPlexGetCone(dm,sB,&coneB);CHKERRQ(ierr); 114 ierr = DMPlexGetConeOrientation(dm,sA,&oA);CHKERRQ(ierr); 115 ierr = DMPlexGetConeOrientation(dm,sB,&oB);CHKERRQ(ierr); 116 /* step through the cone of sA in natural order */ 117 for (i = 0; i < sConeSize; i++) { 118 if (coneA[i] == childA) { 119 /* if childA is at position i in coneA, 120 * then we want the point that is at sOrientB*i in coneB */ 121 PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize -(sOrientB+1) - i) % sConeSize); 122 if (childB) *childB = coneB[j]; 123 if (childOrientB) { 124 PetscInt oBtrue; 125 126 ierr = DMPlexGetConeSize(dm,childA,&coneSize);CHKERRQ(ierr); 127 /* compose sOrientB and oB[j] */ 128 if (coneSize != 0 && coneSize != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected a vertex or an edge"); 129 /* we may have to flip an edge */ 130 oBtrue = coneSize ? ((sOrientB >= 0) ? oB[j] : -(oB[j] + 2)) : 0; 131 ABswap = DihedralSwap(coneSize,oA[i],oBtrue);CHKERRQ(ierr); 132 *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap); 133 } 134 break; 135 } 136 } 137 if (i == sConeSize) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"support cone mismatch"); 138 PetscFunctionReturn(0); 139 } 140 /* get the cone size and symmetry swap */ 141 ierr = DMPlexGetConeSize(dm,parent,&coneSize);CHKERRQ(ierr); 142 ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB); 143 if (dim == 2) { 144 /* orientations refer to cones: we want them to refer to vertices: 145 * if it's a rotation, they are the same, but if the order is reversed, a 146 * permutation that puts side i first does *not* put vertex i first */ 147 oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1); 148 oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1); 149 ABswapVert = DihedralSwap(coneSize, oAvert, oBvert); 150 } 151 else { 152 oAvert = parentOrientA; 153 oBvert = parentOrientB; 154 ABswapVert = ABswap; 155 } 156 if (childB) { 157 /* assume that each child corresponds to a vertex, in the same order */ 158 PetscInt p, posA = -1, numChildren, i; 159 const PetscInt *children; 160 161 /* count which position the child is in */ 162 ierr = DMPlexGetTreeChildren(dm,parent,&numChildren,&children);CHKERRQ(ierr); 163 for (i = 0; i < numChildren; i++) { 164 p = children[i]; 165 if (p == childA) { 166 posA = i; 167 break; 168 } 169 } 170 if (posA >= coneSize) { 171 /* this is the triangle in the middle of a uniformly refined triangle: it 172 * is invariant */ 173 if (dim != 2 || posA != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Expected a middle triangle, got something else"); 174 *childB = childA; 175 } 176 else { 177 /* figure out position B by applying ABswapVert */ 178 PetscInt posB; 179 180 posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize -(ABswapVert + 1) - posA) % coneSize); 181 if (childB) *childB = children[posB]; 182 } 183 } 184 if (childOrientB) *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap); 185 PetscFunctionReturn(0); 186 } 187 188 #undef __FUNCT__ 189 #define __FUNCT__ "DMPlexReferenceTreeGetChildSymmetry" 190 /*@ 191 DMPlexReferenceTreeGetChildSymmetry - Given a reference tree, transform a childid and orientation from one parent frame to another 192 193 Input Parameters: 194 + dm - the reference tree DMPlex object 195 . parent - the parent point 196 . parentOrientA - the reference orientation for describing the parent 197 . childOrientA - the reference orientation for describing the child 198 . childA - the reference childID for describing the child 199 - parentOrientB - the new orientation for describing the parent 200 201 Output Parameters: 202 + childOrientB - if not NULL, set to the new oreintation for describing the child 203 . childB - if not NULL, the new childID for describing the child 204 205 Level: developer 206 207 .seealso: DMPlexGetReferenceTree(), DMPlexSetReferenceTree(), DMPlexSetTree() 208 @*/ 209 PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB) 210 { 211 DM_Plex *mesh = (DM_Plex *)dm->data; 212 PetscErrorCode ierr; 213 214 PetscFunctionBegin; 215 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 216 if (!mesh->getchildsymmetry) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlexReferenceTreeGetChildSymmetry not implemented"); 217 ierr = mesh->getchildsymmetry(dm,parent,parentOrientA,childOrientA,childA,parentOrientB,childOrientB,childB);CHKERRQ(ierr); 218 PetscFunctionReturn(0); 219 } 220 221 static PetscErrorCode DMPlexSetTree_Internal(DM,PetscSection,PetscInt*,PetscInt*,PetscBool,PetscBool); 222 223 #undef __FUNCT__ 224 #define __FUNCT__ "DMPlexCreateDefaultReferenceTree" 225 /*@ 226 DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement. 227 228 Collective on comm 229 230 Input Parameters: 231 + comm - the MPI communicator 232 . dim - the spatial dimension 233 - simplex - Flag for simplex, otherwise use a tensor-product cell 234 235 Output Parameters: 236 . ref - the reference tree DMPlex object 237 238 Level: intermediate 239 240 .keywords: reference cell 241 .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree() 242 @*/ 243 PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref) 244 { 245 DM_Plex *mesh; 246 DM K, Kref; 247 PetscInt p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs; 248 PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts; 249 DMLabel identity, identityRef; 250 PetscSection unionSection, unionConeSection, parentSection; 251 PetscScalar *unionCoords; 252 IS perm; 253 PetscErrorCode ierr; 254 255 PetscFunctionBegin; 256 /* create a reference element */ 257 ierr = DMPlexCreateReferenceCell(comm, dim, simplex, &K);CHKERRQ(ierr); 258 ierr = DMPlexCreateLabel(K, "identity");CHKERRQ(ierr); 259 ierr = DMPlexGetLabel(K, "identity", &identity);CHKERRQ(ierr); 260 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr); 261 for (p = pStart; p < pEnd; p++) { 262 ierr = DMLabelSetValue(identity, p, p);CHKERRQ(ierr); 263 } 264 /* refine it */ 265 ierr = DMRefine(K,comm,&Kref);CHKERRQ(ierr); 266 267 /* the reference tree is the union of these two, without duplicating 268 * points that appear in both */ 269 ierr = DMPlexGetLabel(Kref, "identity", &identityRef);CHKERRQ(ierr); 270 ierr = DMPlexGetChart(Kref, &pRefStart, &pRefEnd);CHKERRQ(ierr); 271 ierr = PetscSectionCreate(comm, &unionSection);CHKERRQ(ierr); 272 ierr = PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart));CHKERRQ(ierr); 273 /* count points that will go in the union */ 274 for (p = pStart; p < pEnd; p++) { 275 ierr = PetscSectionSetDof(unionSection, p - pStart, 1);CHKERRQ(ierr); 276 } 277 for (p = pRefStart; p < pRefEnd; p++) { 278 PetscInt q, qSize; 279 ierr = DMLabelGetValue(identityRef, p, &q);CHKERRQ(ierr); 280 ierr = DMLabelGetStratumSize(identityRef, q, &qSize);CHKERRQ(ierr); 281 if (qSize > 1) { 282 ierr = PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1);CHKERRQ(ierr); 283 } 284 } 285 ierr = PetscMalloc1((pEnd - pStart) + (pRefEnd - pRefStart),&permvals);CHKERRQ(ierr); 286 offset = 0; 287 /* stratify points in the union by topological dimension */ 288 for (d = 0; d <= dim; d++) { 289 PetscInt cStart, cEnd, c; 290 291 ierr = DMPlexGetHeightStratum(K, d, &cStart, &cEnd);CHKERRQ(ierr); 292 for (c = cStart; c < cEnd; c++) { 293 permvals[offset++] = c; 294 } 295 296 ierr = DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);CHKERRQ(ierr); 297 for (c = cStart; c < cEnd; c++) { 298 permvals[offset++] = c + (pEnd - pStart); 299 } 300 } 301 ierr = ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);CHKERRQ(ierr); 302 ierr = PetscSectionSetPermutation(unionSection,perm);CHKERRQ(ierr); 303 ierr = PetscSectionSetUp(unionSection);CHKERRQ(ierr); 304 ierr = PetscSectionGetStorageSize(unionSection,&numUnionPoints);CHKERRQ(ierr); 305 ierr = PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints);CHKERRQ(ierr); 306 /* count dimension points */ 307 for (d = 0; d <= dim; d++) { 308 PetscInt cStart, cOff, cOff2; 309 ierr = DMPlexGetHeightStratum(K,d,&cStart,NULL);CHKERRQ(ierr); 310 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff);CHKERRQ(ierr); 311 if (d < dim) { 312 ierr = DMPlexGetHeightStratum(K,d+1,&cStart,NULL);CHKERRQ(ierr); 313 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2);CHKERRQ(ierr); 314 } 315 else { 316 cOff2 = numUnionPoints; 317 } 318 numDimPoints[dim - d] = cOff2 - cOff; 319 } 320 ierr = PetscSectionCreate(comm, &unionConeSection);CHKERRQ(ierr); 321 ierr = PetscSectionSetChart(unionConeSection, 0, numUnionPoints);CHKERRQ(ierr); 322 /* count the cones in the union */ 323 for (p = pStart; p < pEnd; p++) { 324 PetscInt dof, uOff; 325 326 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 327 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 328 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 329 coneSizes[uOff] = dof; 330 } 331 for (p = pRefStart; p < pRefEnd; p++) { 332 PetscInt dof, uDof, uOff; 333 334 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 335 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 336 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 337 if (uDof) { 338 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 339 coneSizes[uOff] = dof; 340 } 341 } 342 ierr = PetscSectionSetUp(unionConeSection);CHKERRQ(ierr); 343 ierr = PetscSectionGetStorageSize(unionConeSection,&numCones);CHKERRQ(ierr); 344 ierr = PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations);CHKERRQ(ierr); 345 /* write the cones in the union */ 346 for (p = pStart; p < pEnd; p++) { 347 PetscInt dof, uOff, c, cOff; 348 const PetscInt *cone, *orientation; 349 350 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 351 ierr = DMPlexGetCone(K, p, &cone);CHKERRQ(ierr); 352 ierr = DMPlexGetConeOrientation(K, p, &orientation);CHKERRQ(ierr); 353 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 354 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 355 for (c = 0; c < dof; c++) { 356 PetscInt e, eOff; 357 e = cone[c]; 358 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 359 unionCones[cOff + c] = eOff; 360 unionOrientations[cOff + c] = orientation[c]; 361 } 362 } 363 for (p = pRefStart; p < pRefEnd; p++) { 364 PetscInt dof, uDof, uOff, c, cOff; 365 const PetscInt *cone, *orientation; 366 367 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 368 ierr = DMPlexGetCone(Kref, p, &cone);CHKERRQ(ierr); 369 ierr = DMPlexGetConeOrientation(Kref, p, &orientation);CHKERRQ(ierr); 370 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 371 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 372 if (uDof) { 373 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 374 for (c = 0; c < dof; c++) { 375 PetscInt e, eOff, eDof; 376 377 e = cone[c]; 378 ierr = PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof);CHKERRQ(ierr); 379 if (eDof) { 380 ierr = PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);CHKERRQ(ierr); 381 } 382 else { 383 ierr = DMLabelGetValue(identityRef, e, &e);CHKERRQ(ierr); 384 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 385 } 386 unionCones[cOff + c] = eOff; 387 unionOrientations[cOff + c] = orientation[c]; 388 } 389 } 390 } 391 /* get the coordinates */ 392 { 393 PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff; 394 PetscSection KcoordsSec, KrefCoordsSec; 395 Vec KcoordsVec, KrefCoordsVec; 396 PetscScalar *Kcoords; 397 398 DMGetCoordinateSection(K, &KcoordsSec);CHKERRQ(ierr); 399 DMGetCoordinatesLocal(K, &KcoordsVec);CHKERRQ(ierr); 400 DMGetCoordinateSection(Kref, &KrefCoordsSec);CHKERRQ(ierr); 401 DMGetCoordinatesLocal(Kref, &KrefCoordsVec);CHKERRQ(ierr); 402 403 numVerts = numDimPoints[0]; 404 ierr = PetscMalloc1(numVerts * dim,&unionCoords);CHKERRQ(ierr); 405 ierr = DMPlexGetDepthStratum(K,0,&vStart,&vEnd);CHKERRQ(ierr); 406 407 offset = 0; 408 for (v = vStart; v < vEnd; v++) { 409 ierr = PetscSectionGetOffset(unionSection,v - pStart,&vOff);CHKERRQ(ierr); 410 ierr = VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);CHKERRQ(ierr); 411 for (d = 0; d < dim; d++) { 412 unionCoords[offset * dim + d] = Kcoords[d]; 413 } 414 offset++; 415 } 416 ierr = DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd);CHKERRQ(ierr); 417 for (v = vRefStart; v < vRefEnd; v++) { 418 ierr = PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof);CHKERRQ(ierr); 419 ierr = PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff);CHKERRQ(ierr); 420 ierr = VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);CHKERRQ(ierr); 421 if (vDof) { 422 for (d = 0; d < dim; d++) { 423 unionCoords[offset * dim + d] = Kcoords[d]; 424 } 425 offset++; 426 } 427 } 428 } 429 ierr = DMCreate(comm,ref);CHKERRQ(ierr); 430 ierr = DMSetType(*ref,DMPLEX);CHKERRQ(ierr); 431 ierr = DMPlexSetDimension(*ref,dim);CHKERRQ(ierr); 432 ierr = DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords);CHKERRQ(ierr); 433 /* set the tree */ 434 ierr = PetscSectionCreate(comm,&parentSection);CHKERRQ(ierr); 435 ierr = PetscSectionSetChart(parentSection,0,numUnionPoints);CHKERRQ(ierr); 436 for (p = pRefStart; p < pRefEnd; p++) { 437 PetscInt uDof, uOff; 438 439 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 440 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 441 if (uDof) { 442 PetscSectionSetDof(parentSection,uOff,1);CHKERRQ(ierr); 443 } 444 } 445 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 446 ierr = PetscSectionGetStorageSize(parentSection,&parentSize);CHKERRQ(ierr); 447 ierr = PetscMalloc2(parentSize,&parents,parentSize,&childIDs);CHKERRQ(ierr); 448 for (p = pRefStart; p < pRefEnd; p++) { 449 PetscInt uDof, uOff; 450 451 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 452 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 453 if (uDof) { 454 PetscInt pOff, parent, parentU; 455 PetscSectionGetOffset(parentSection,uOff,&pOff);CHKERRQ(ierr); 456 DMLabelGetValue(identityRef,p,&parent);CHKERRQ(ierr); 457 ierr = PetscSectionGetOffset(unionSection, parent - pStart,&parentU);CHKERRQ(ierr); 458 parents[pOff] = parentU; 459 childIDs[pOff] = uOff; 460 } 461 } 462 ierr = DMPlexSetTree_Internal(*ref,parentSection,parents,childIDs,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 463 mesh = (DM_Plex *) (*ref)->data; 464 mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default; 465 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 466 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 467 468 /* clean up */ 469 ierr = PetscSectionDestroy(&unionSection);CHKERRQ(ierr); 470 ierr = PetscSectionDestroy(&unionConeSection);CHKERRQ(ierr); 471 ierr = ISDestroy(&perm);CHKERRQ(ierr); 472 ierr = PetscFree(unionCoords);CHKERRQ(ierr); 473 ierr = PetscFree2(unionCones,unionOrientations);CHKERRQ(ierr); 474 ierr = PetscFree2(coneSizes,numDimPoints);CHKERRQ(ierr); 475 ierr = DMDestroy(&K);CHKERRQ(ierr); 476 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 477 PetscFunctionReturn(0); 478 } 479 480 #undef __FUNCT__ 481 #define __FUNCT__ "DMPlexTreeSymmetrize" 482 static PetscErrorCode DMPlexTreeSymmetrize(DM dm) 483 { 484 DM_Plex *mesh = (DM_Plex *)dm->data; 485 PetscSection childSec, pSec; 486 PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT; 487 PetscInt *offsets, *children, pStart, pEnd; 488 PetscErrorCode ierr; 489 490 PetscFunctionBegin; 491 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 492 ierr = PetscSectionDestroy(&mesh->childSection);CHKERRQ(ierr); 493 ierr = PetscFree(mesh->children);CHKERRQ(ierr); 494 pSec = mesh->parentSection; 495 if (!pSec) PetscFunctionReturn(0); 496 ierr = PetscSectionGetStorageSize(pSec,&pSize);CHKERRQ(ierr); 497 for (p = 0; p < pSize; p++) { 498 PetscInt par = mesh->parents[p]; 499 500 parMax = PetscMax(parMax,par+1); 501 parMin = PetscMin(parMin,par); 502 } 503 if (parMin > parMax) { 504 parMin = -1; 505 parMax = -1; 506 } 507 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);CHKERRQ(ierr); 508 ierr = PetscSectionSetChart(childSec,parMin,parMax);CHKERRQ(ierr); 509 for (p = 0; p < pSize; p++) { 510 PetscInt par = mesh->parents[p]; 511 512 ierr = PetscSectionAddDof(childSec,par,1);CHKERRQ(ierr); 513 } 514 ierr = PetscSectionSetUp(childSec);CHKERRQ(ierr); 515 ierr = PetscSectionGetStorageSize(childSec,&cSize);CHKERRQ(ierr); 516 ierr = PetscMalloc1(cSize,&children);CHKERRQ(ierr); 517 ierr = PetscCalloc1(parMax-parMin,&offsets);CHKERRQ(ierr); 518 ierr = PetscSectionGetChart(pSec,&pStart,&pEnd);CHKERRQ(ierr); 519 for (p = pStart; p < pEnd; p++) { 520 PetscInt dof, off, i; 521 522 ierr = PetscSectionGetDof(pSec,p,&dof);CHKERRQ(ierr); 523 ierr = PetscSectionGetOffset(pSec,p,&off);CHKERRQ(ierr); 524 for (i = 0; i < dof; i++) { 525 PetscInt par = mesh->parents[off + i], cOff; 526 527 ierr = PetscSectionGetOffset(childSec,par,&cOff);CHKERRQ(ierr); 528 children[cOff + offsets[par-parMin]++] = p; 529 } 530 } 531 mesh->childSection = childSec; 532 mesh->children = children; 533 ierr = PetscFree(offsets);CHKERRQ(ierr); 534 PetscFunctionReturn(0); 535 } 536 537 #undef __FUNCT__ 538 #define __FUNCT__ "AnchorsFlatten" 539 static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew) 540 { 541 PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL; 542 const PetscInt *vals; 543 PetscSection secNew; 544 PetscBool anyNew, globalAnyNew; 545 PetscBool compress; 546 PetscErrorCode ierr; 547 548 PetscFunctionBegin; 549 ierr = PetscSectionGetChart(section,&pStart,&pEnd);CHKERRQ(ierr); 550 ierr = ISGetLocalSize(is,&size);CHKERRQ(ierr); 551 ierr = ISGetIndices(is,&vals);CHKERRQ(ierr); 552 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);CHKERRQ(ierr); 553 ierr = PetscSectionSetChart(secNew,pStart,pEnd);CHKERRQ(ierr); 554 for (i = 0; i < size; i++) { 555 PetscInt dof; 556 557 p = vals[i]; 558 if (p < pStart || p >= pEnd) continue; 559 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 560 if (dof) break; 561 } 562 if (i == size) { 563 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 564 anyNew = PETSC_FALSE; 565 compress = PETSC_FALSE; 566 sizeNew = 0; 567 } 568 else { 569 anyNew = PETSC_TRUE; 570 for (p = pStart; p < pEnd; p++) { 571 PetscInt dof, off; 572 573 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 574 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 575 for (i = 0; i < dof; i++) { 576 PetscInt q = vals[off + i], qDof = 0; 577 578 if (q >= pStart && q < pEnd) { 579 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 580 } 581 if (qDof) { 582 ierr = PetscSectionAddDof(secNew, p, qDof);CHKERRQ(ierr); 583 } 584 else { 585 ierr = PetscSectionAddDof(secNew, p, 1);CHKERRQ(ierr); 586 } 587 } 588 } 589 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 590 ierr = PetscSectionGetStorageSize(secNew,&sizeNew);CHKERRQ(ierr); 591 ierr = PetscMalloc1(sizeNew,&valsNew);CHKERRQ(ierr); 592 compress = PETSC_FALSE; 593 for (p = pStart; p < pEnd; p++) { 594 PetscInt dof, off, count, offNew, dofNew; 595 596 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 597 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 598 ierr = PetscSectionGetDof(secNew, p, &dofNew);CHKERRQ(ierr); 599 ierr = PetscSectionGetOffset(secNew, p, &offNew);CHKERRQ(ierr); 600 count = 0; 601 for (i = 0; i < dof; i++) { 602 PetscInt q = vals[off + i], qDof = 0, qOff = 0, j; 603 604 if (q >= pStart && q < pEnd) { 605 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 606 ierr = PetscSectionGetOffset(section, q, &qOff);CHKERRQ(ierr); 607 } 608 if (qDof) { 609 PetscInt oldCount = count; 610 611 for (j = 0; j < qDof; j++) { 612 PetscInt k, r = vals[qOff + j]; 613 614 for (k = 0; k < oldCount; k++) { 615 if (valsNew[offNew + k] == r) { 616 break; 617 } 618 } 619 if (k == oldCount) { 620 valsNew[offNew + count++] = r; 621 } 622 } 623 } 624 else { 625 PetscInt k, oldCount = count; 626 627 for (k = 0; k < oldCount; k++) { 628 if (valsNew[offNew + k] == q) { 629 break; 630 } 631 } 632 if (k == oldCount) { 633 valsNew[offNew + count++] = q; 634 } 635 } 636 } 637 if (count < dofNew) { 638 ierr = PetscSectionSetDof(secNew, p, count);CHKERRQ(ierr); 639 compress = PETSC_TRUE; 640 } 641 } 642 } 643 ierr = ISRestoreIndices(is,&vals);CHKERRQ(ierr); 644 ierr = MPI_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 645 if (!globalAnyNew) { 646 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 647 *sectionNew = NULL; 648 *isNew = NULL; 649 } 650 else { 651 PetscBool globalCompress; 652 653 ierr = MPI_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 654 if (compress) { 655 PetscSection secComp; 656 PetscInt *valsComp = NULL; 657 658 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);CHKERRQ(ierr); 659 ierr = PetscSectionSetChart(secComp,pStart,pEnd);CHKERRQ(ierr); 660 for (p = pStart; p < pEnd; p++) { 661 PetscInt dof; 662 663 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 664 ierr = PetscSectionSetDof(secComp, p, dof);CHKERRQ(ierr); 665 } 666 ierr = PetscSectionSetUp(secComp);CHKERRQ(ierr); 667 ierr = PetscSectionGetStorageSize(secComp,&sizeNew);CHKERRQ(ierr); 668 ierr = PetscMalloc1(sizeNew,&valsComp);CHKERRQ(ierr); 669 for (p = pStart; p < pEnd; p++) { 670 PetscInt dof, off, offNew, j; 671 672 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 673 ierr = PetscSectionGetOffset(secNew, p, &off);CHKERRQ(ierr); 674 ierr = PetscSectionGetOffset(secComp, p, &offNew);CHKERRQ(ierr); 675 for (j = 0; j < dof; j++) { 676 valsComp[offNew + j] = valsNew[off + j]; 677 } 678 } 679 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 680 secNew = secComp; 681 ierr = PetscFree(valsNew);CHKERRQ(ierr); 682 valsNew = valsComp; 683 } 684 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);CHKERRQ(ierr); 685 } 686 PetscFunctionReturn(0); 687 } 688 689 #undef __FUNCT__ 690 #define __FUNCT__ "DMPlexComputeConstraints_Tree" 691 static PetscErrorCode DMPlexComputeConstraints_Tree(DM dm) 692 { 693 PetscInt p, pStart, pEnd, *anchors, size; 694 PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT; 695 PetscSection aSec; 696 DMLabel canonLabel; 697 IS aIS; 698 PetscErrorCode ierr; 699 700 PetscFunctionBegin; 701 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 702 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 703 ierr = DMPlexGetLabel(dm,"canonical",&canonLabel);CHKERRQ(ierr); 704 for (p = pStart; p < pEnd; p++) { 705 PetscInt parent; 706 707 if (canonLabel) { 708 PetscInt canon; 709 710 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 711 if (p != canon) continue; 712 } 713 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 714 if (parent != p) { 715 aMin = PetscMin(aMin,p); 716 aMax = PetscMax(aMax,p+1); 717 } 718 } 719 if (aMin > aMax) { 720 aMin = -1; 721 aMax = -1; 722 } 723 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&aSec);CHKERRQ(ierr); 724 ierr = PetscSectionSetChart(aSec,aMin,aMax);CHKERRQ(ierr); 725 for (p = aMin; p < aMax; p++) { 726 PetscInt parent, ancestor = p; 727 728 if (canonLabel) { 729 PetscInt canon; 730 731 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 732 if (p != canon) continue; 733 } 734 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 735 while (parent != ancestor) { 736 ancestor = parent; 737 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 738 } 739 if (ancestor != p) { 740 PetscInt closureSize, *closure = NULL; 741 742 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 743 ierr = PetscSectionSetDof(aSec,p,closureSize);CHKERRQ(ierr); 744 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 745 } 746 } 747 ierr = PetscSectionSetUp(aSec);CHKERRQ(ierr); 748 ierr = PetscSectionGetStorageSize(aSec,&size);CHKERRQ(ierr); 749 ierr = PetscMalloc1(size,&anchors);CHKERRQ(ierr); 750 for (p = aMin; p < aMax; p++) { 751 PetscInt parent, ancestor = p; 752 753 if (canonLabel) { 754 PetscInt canon; 755 756 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 757 if (p != canon) continue; 758 } 759 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 760 while (parent != ancestor) { 761 ancestor = parent; 762 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 763 } 764 if (ancestor != p) { 765 PetscInt j, closureSize, *closure = NULL, aOff; 766 767 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 768 769 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 770 for (j = 0; j < closureSize; j++) { 771 anchors[aOff + j] = closure[2*j]; 772 } 773 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 774 } 775 } 776 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm),size,anchors,PETSC_OWN_POINTER,&aIS);CHKERRQ(ierr); 777 { 778 PetscSection aSecNew = aSec; 779 IS aISNew = aIS; 780 781 ierr = PetscObjectReference((PetscObject)aSec);CHKERRQ(ierr); 782 ierr = PetscObjectReference((PetscObject)aIS);CHKERRQ(ierr); 783 while (aSecNew) { 784 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 785 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 786 aSec = aSecNew; 787 aIS = aISNew; 788 aSecNew = NULL; 789 aISNew = NULL; 790 ierr = AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);CHKERRQ(ierr); 791 } 792 } 793 ierr = DMPlexSetConstraints(dm,aSec,aIS);CHKERRQ(ierr); 794 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 795 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 796 PetscFunctionReturn(0); 797 } 798 799 #undef __FUNCT__ 800 #define __FUNCT__ "DMPlexTreeExchangeSupports" 801 static PetscErrorCode DMPlexTreeExchangeSupports(DM dm) 802 { 803 DM_Plex *mesh = (DM_Plex *)dm->data; 804 PetscSection newSupportSection; 805 PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth; 806 PetscInt *offsets; 807 PetscErrorCode ierr; 808 809 PetscFunctionBegin; 810 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 811 /* symmetrize the hierarchy */ 812 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 813 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&newSupportSection);CHKERRQ(ierr); 814 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 815 ierr = PetscSectionSetChart(newSupportSection,pStart,pEnd);CHKERRQ(ierr); 816 ierr = PetscCalloc1(pEnd,&offsets);CHKERRQ(ierr); 817 /* if a point is in the support of q, it should be in the support of 818 * parent(q) */ 819 for (d = 0; d <= depth; d++) { 820 ierr = DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);CHKERRQ(ierr); 821 for (p = pStart; p < pEnd; ++p) { 822 PetscInt dof, q, qdof, parent; 823 824 ierr = PetscSectionGetDof(mesh->supportSection, p, &dof);CHKERRQ(ierr); 825 ierr = PetscSectionAddDof(newSupportSection, p, dof);CHKERRQ(ierr); 826 q = p; 827 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 828 while (parent != q && parent >= pStart && parent < pEnd) { 829 q = parent; 830 831 ierr = PetscSectionGetDof(mesh->supportSection, q, &qdof);CHKERRQ(ierr); 832 ierr = PetscSectionAddDof(newSupportSection,p,qdof);CHKERRQ(ierr); 833 ierr = PetscSectionAddDof(newSupportSection,q,dof);CHKERRQ(ierr); 834 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 835 } 836 } 837 } 838 ierr = PetscSectionSetUp(newSupportSection);CHKERRQ(ierr); 839 ierr = PetscSectionGetStorageSize(newSupportSection,&newSize);CHKERRQ(ierr); 840 ierr = PetscMalloc1(newSize,&newSupports);CHKERRQ(ierr); 841 for (d = 0; d <= depth; d++) { 842 ierr = DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);CHKERRQ(ierr); 843 for (p = pStart; p < pEnd; p++) { 844 PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent; 845 846 ierr = PetscSectionGetDof(mesh->supportSection, p, &dof);CHKERRQ(ierr); 847 ierr = PetscSectionGetOffset(mesh->supportSection, p, &off);CHKERRQ(ierr); 848 ierr = PetscSectionGetDof(newSupportSection, p, &newDof);CHKERRQ(ierr); 849 ierr = PetscSectionGetOffset(newSupportSection, p, &newOff);CHKERRQ(ierr); 850 for (i = 0; i < dof; i++) { 851 newSupports[newOff+offsets[p]++] = mesh->supports[off + i]; 852 } 853 mesh->maxSupportSize = PetscMax(mesh->maxSupportSize,newDof); 854 855 q = p; 856 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 857 while (parent != q && parent >= pStart && parent < pEnd) { 858 q = parent; 859 ierr = PetscSectionGetDof(mesh->supportSection, q, &qdof);CHKERRQ(ierr); 860 ierr = PetscSectionGetOffset(mesh->supportSection, q, &qoff);CHKERRQ(ierr); 861 ierr = PetscSectionGetOffset(newSupportSection, q, &newqOff);CHKERRQ(ierr); 862 for (i = 0; i < qdof; i++) { 863 newSupports[newOff+offsets[p]++] = mesh->supports[qoff + i]; 864 } 865 for (i = 0; i < dof; i++) { 866 newSupports[newqOff+offsets[q]++] = mesh->supports[off + i]; 867 } 868 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 869 } 870 } 871 } 872 ierr = PetscSectionDestroy(&mesh->supportSection);CHKERRQ(ierr); 873 mesh->supportSection = newSupportSection; 874 ierr = PetscFree(mesh->supports);CHKERRQ(ierr); 875 mesh->supports = newSupports; 876 ierr = PetscFree(offsets);CHKERRQ(ierr); 877 878 PetscFunctionReturn(0); 879 } 880 881 #undef __FUNCT__ 882 #define __FUNCT__ "DMPlexSetTree_Internal" 883 static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports) 884 { 885 DM_Plex *mesh = (DM_Plex *)dm->data; 886 DM refTree; 887 PetscInt size; 888 PetscErrorCode ierr; 889 890 PetscFunctionBegin; 891 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 892 PetscValidHeaderSpecific(parentSection, PETSC_SECTION_CLASSID, 2); 893 ierr = PetscObjectReference((PetscObject)parentSection);CHKERRQ(ierr); 894 ierr = PetscSectionDestroy(&mesh->parentSection);CHKERRQ(ierr); 895 mesh->parentSection = parentSection; 896 ierr = PetscSectionGetStorageSize(parentSection,&size);CHKERRQ(ierr); 897 if (parents != mesh->parents) { 898 ierr = PetscFree(mesh->parents);CHKERRQ(ierr); 899 ierr = PetscMalloc1(size,&mesh->parents);CHKERRQ(ierr); 900 ierr = PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));CHKERRQ(ierr); 901 } 902 if (childIDs != mesh->childIDs) { 903 ierr = PetscFree(mesh->childIDs);CHKERRQ(ierr); 904 ierr = PetscMalloc1(size,&mesh->childIDs);CHKERRQ(ierr); 905 ierr = PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));CHKERRQ(ierr); 906 } 907 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr); 908 if (refTree) { 909 DMLabel canonLabel; 910 911 ierr = DMPlexGetLabel(refTree,"canonical",&canonLabel);CHKERRQ(ierr); 912 if (canonLabel) { 913 PetscInt i; 914 915 for (i = 0; i < size; i++) { 916 PetscInt canon; 917 ierr = DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);CHKERRQ(ierr); 918 if (canon >= 0) { 919 mesh->childIDs[i] = canon; 920 } 921 } 922 } 923 } 924 ierr = DMPlexTreeSymmetrize(dm);CHKERRQ(ierr); 925 if (computeCanonical) { 926 PetscInt d, dim; 927 928 /* add the canonical label */ 929 ierr = DMPlexGetDimension(dm,&dim);CHKERRQ(ierr); 930 ierr = DMPlexCreateLabel(dm,"canonical");CHKERRQ(ierr); 931 for (d = 0; d <= dim; d++) { 932 PetscInt p, dStart, dEnd, canon = -1, cNumChildren; 933 const PetscInt *cChildren; 934 935 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr); 936 for (p = dStart; p < dEnd; p++) { 937 ierr = DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);CHKERRQ(ierr); 938 if (cNumChildren) { 939 canon = p; 940 break; 941 } 942 } 943 if (canon == -1) continue; 944 for (p = dStart; p < dEnd; p++) { 945 PetscInt numChildren, i; 946 const PetscInt *children; 947 948 ierr = DMPlexGetTreeChildren(dm,p,&numChildren,&children);CHKERRQ(ierr); 949 if (numChildren) { 950 if (numChildren != cNumChildren) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"All parent points in a stratum should have the same number of children: %d != %d", numChildren, cNumChildren); 951 ierr = DMPlexSetLabelValue(dm,"canonical",p,canon);CHKERRQ(ierr); 952 for (i = 0; i < numChildren; i++) { 953 ierr = DMPlexSetLabelValue(dm,"canonical",children[i],cChildren[i]);CHKERRQ(ierr); 954 } 955 } 956 } 957 } 958 } 959 if (exchangeSupports) { 960 ierr = DMPlexTreeExchangeSupports(dm);CHKERRQ(ierr); 961 } 962 ierr = DMPlexComputeConstraints_Tree(dm);CHKERRQ(ierr); 963 PetscFunctionReturn(0); 964 } 965 966 #undef __FUNCT__ 967 #define __FUNCT__ "DMPlexSetTree" 968 /*@ 969 DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates 970 the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its 971 tree root. 972 973 Collective on dm 974 975 Input Parameters: 976 + dm - the DMPlex object 977 . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 978 offset indexes the parent and childID list; the reference count of parentSection is incremented 979 . parents - a list of the point parents; copied, can be destroyed 980 - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 981 the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed 982 983 Level: intermediate 984 985 .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetConstraints(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 986 @*/ 987 PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[]) 988 { 989 PetscErrorCode ierr; 990 991 PetscFunctionBegin; 992 ierr = DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE);CHKERRQ(ierr); 993 PetscFunctionReturn(0); 994 } 995 996 #undef __FUNCT__ 997 #define __FUNCT__ "DMPlexGetTree" 998 /*@ 999 DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points. 1000 Collective on dm 1001 1002 Input Parameters: 1003 . dm - the DMPlex object 1004 1005 Output Parameters: 1006 + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 1007 offset indexes the parent and childID list 1008 . parents - a list of the point parents 1009 . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 1010 the child corresponds to the point in the reference tree with index childID 1011 . childSection - the inverse of the parent section 1012 - children - a list of the point children 1013 1014 Level: intermediate 1015 1016 .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetConstraints(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 1017 @*/ 1018 PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[]) 1019 { 1020 DM_Plex *mesh = (DM_Plex *)dm->data; 1021 1022 PetscFunctionBegin; 1023 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1024 if (parentSection) *parentSection = mesh->parentSection; 1025 if (parents) *parents = mesh->parents; 1026 if (childIDs) *childIDs = mesh->childIDs; 1027 if (childSection) *childSection = mesh->childSection; 1028 if (children) *children = mesh->children; 1029 PetscFunctionReturn(0); 1030 } 1031 1032 #undef __FUNCT__ 1033 #define __FUNCT__ "DMPlexGetTreeParent" 1034 /*@ 1035 DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the Sieve DAG) 1036 1037 Input Parameters: 1038 + dm - the DMPlex object 1039 - point - the query point 1040 1041 Output Parameters: 1042 + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent 1043 - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point 1044 does not have a parent 1045 1046 Level: intermediate 1047 1048 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren() 1049 @*/ 1050 PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID) 1051 { 1052 DM_Plex *mesh = (DM_Plex *)dm->data; 1053 PetscSection pSec; 1054 PetscErrorCode ierr; 1055 1056 PetscFunctionBegin; 1057 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1058 pSec = mesh->parentSection; 1059 if (pSec && point >= pSec->pStart && point < pSec->pEnd) { 1060 PetscInt dof; 1061 1062 ierr = PetscSectionGetDof (pSec, point, &dof);CHKERRQ(ierr); 1063 if (dof) { 1064 PetscInt off; 1065 1066 ierr = PetscSectionGetOffset (pSec, point, &off);CHKERRQ(ierr); 1067 if (parent) *parent = mesh->parents[off]; 1068 if (childID) *childID = mesh->childIDs[off]; 1069 PetscFunctionReturn(0); 1070 } 1071 } 1072 if (parent) { 1073 *parent = point; 1074 } 1075 if (childID) { 1076 *childID = 0; 1077 } 1078 PetscFunctionReturn(0); 1079 } 1080 1081 #undef __FUNCT__ 1082 #define __FUNCT__ "DMPlexGetTreeChildren" 1083 /*@C 1084 DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the Sieve DAG) 1085 1086 Input Parameters: 1087 + dm - the DMPlex object 1088 - point - the query point 1089 1090 Output Parameters: 1091 + numChildren - if not NULL, set to the number of children 1092 - children - if not NULL, set to a list children, or set to NULL if the point has no children 1093 1094 Level: intermediate 1095 1096 Fortran Notes: 1097 Since it returns an array, this routine is only available in Fortran 90, and you must 1098 include petsc.h90 in your code. 1099 1100 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent() 1101 @*/ 1102 PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[]) 1103 { 1104 DM_Plex *mesh = (DM_Plex *)dm->data; 1105 PetscSection childSec; 1106 PetscInt dof = 0; 1107 PetscErrorCode ierr; 1108 1109 PetscFunctionBegin; 1110 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1111 childSec = mesh->childSection; 1112 if (childSec && point >= childSec->pStart && point < childSec->pEnd) { 1113 ierr = PetscSectionGetDof (childSec, point, &dof);CHKERRQ(ierr); 1114 } 1115 if (numChildren) *numChildren = dof; 1116 if (children) { 1117 if (dof) { 1118 PetscInt off; 1119 1120 ierr = PetscSectionGetOffset (childSec, point, &off);CHKERRQ(ierr); 1121 *children = &mesh->children[off]; 1122 } 1123 else { 1124 *children = NULL; 1125 } 1126 } 1127 PetscFunctionReturn(0); 1128 } 1129 1130 #undef __FUNCT__ 1131 #define __FUNCT__ "DMPlexComputeConstraintMatrix_ReferenceTree" 1132 static PetscErrorCode DMPlexComputeConstraintMatrix_ReferenceTree(DM dm) 1133 { 1134 PetscDS ds; 1135 PetscInt spdim; 1136 PetscInt numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd; 1137 const PetscInt *anchors; 1138 PetscSection section, cSec, aSec; 1139 Mat cMat; 1140 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent; 1141 IS aIS; 1142 PetscErrorCode ierr; 1143 1144 PetscFunctionBegin; 1145 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 1146 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 1147 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1148 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 1149 ierr = DMGetDefaultSection(dm,§ion);CHKERRQ(ierr); 1150 ierr = DMPlexGetConstraintSection(dm,&cSec);CHKERRQ(ierr); 1151 ierr = DMPlexGetConstraints(dm,&aSec,&aIS);CHKERRQ(ierr); 1152 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 1153 ierr = PetscSectionGetChart(cSec,&conStart,&conEnd);CHKERRQ(ierr); 1154 ierr = DMPlexGetConstraintMatrix(dm,&cMat);CHKERRQ(ierr); 1155 ierr = DMPlexGetDimension(dm,&spdim);CHKERRQ(ierr); 1156 ierr = PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent);CHKERRQ(ierr); 1157 1158 for (f = 0; f < numFields; f++) { 1159 PetscFE fe; 1160 PetscDualSpace space; 1161 PetscInt i, j, k, nPoints, offset; 1162 PetscInt fSize, fComp; 1163 PetscScalar *B = NULL; 1164 PetscReal *weights, *pointsRef, *pointsReal; 1165 Mat Amat, Bmat, Xmat; 1166 1167 ierr = PetscDSGetDiscretization(ds,f,(PetscObject *)(&fe));CHKERRQ(ierr); 1168 ierr = PetscFEGetDualSpace(fe,&space);CHKERRQ(ierr); 1169 ierr = PetscDualSpaceGetDimension(space,&fSize);CHKERRQ(ierr); 1170 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 1171 ierr = MatCreate(PETSC_COMM_SELF,&Amat);CHKERRQ(ierr); 1172 ierr = MatSetSizes(Amat,fSize,fSize,fSize,fSize);CHKERRQ(ierr); 1173 ierr = MatSetType(Amat,MATSEQDENSE);CHKERRQ(ierr); 1174 ierr = MatSetUp(Amat);CHKERRQ(ierr); 1175 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);CHKERRQ(ierr); 1176 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);CHKERRQ(ierr); 1177 nPoints = 0; 1178 for (i = 0; i < fSize; i++) { 1179 PetscInt qPoints; 1180 PetscQuadrature quad; 1181 1182 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 1183 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 1184 nPoints += qPoints; 1185 } 1186 ierr = PetscMalloc3(nPoints,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal);CHKERRQ(ierr); 1187 offset = 0; 1188 for (i = 0; i < fSize; i++) { 1189 PetscInt qPoints; 1190 const PetscReal *p, *w; 1191 PetscQuadrature quad; 1192 1193 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 1194 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,&p,&w);CHKERRQ(ierr); 1195 ierr = PetscMemcpy(weights+offset,w,qPoints*sizeof(*w));CHKERRQ(ierr); 1196 ierr = PetscMemcpy(pointsRef+spdim*offset,p,spdim*qPoints*sizeof(*p));CHKERRQ(ierr); 1197 offset += qPoints; 1198 } 1199 ierr = PetscFEGetTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);CHKERRQ(ierr); 1200 offset = 0; 1201 for (i = 0; i < fSize; i++) { 1202 PetscInt qPoints; 1203 PetscQuadrature quad; 1204 1205 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 1206 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 1207 for (j = 0; j < fSize; j++) { 1208 PetscScalar val = 0.; 1209 1210 for (k = 0; k < qPoints; k++) { 1211 val += B[((offset + k) * fSize + j) * fComp] * weights[k]; 1212 } 1213 ierr = MatSetValue(Amat,i,j,val,INSERT_VALUES);CHKERRQ(ierr); 1214 } 1215 offset += qPoints; 1216 } 1217 ierr = PetscFERestoreTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);CHKERRQ(ierr); 1218 ierr = MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1219 ierr = MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1220 ierr = MatLUFactor(Amat,NULL,NULL,NULL);CHKERRQ(ierr); 1221 for (c = cStart; c < cEnd; c++) { 1222 PetscInt parent; 1223 PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL; 1224 PetscInt *childOffsets, *parentOffsets; 1225 1226 ierr = DMPlexGetTreeParent(dm,c,&parent,NULL);CHKERRQ(ierr); 1227 if (parent == c) continue; 1228 ierr = DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1229 for (i = 0; i < closureSize; i++) { 1230 PetscInt p = closure[2*i]; 1231 PetscInt conDof; 1232 1233 if (p < conStart || p >= conEnd) continue; 1234 if (numFields > 1) { 1235 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 1236 } 1237 else { 1238 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 1239 } 1240 if (conDof) break; 1241 } 1242 if (i == closureSize) { 1243 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1244 continue; 1245 } 1246 1247 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr); 1248 ierr = DMPlexComputeCellGeometry(dm, parent, v0parent, Jparent, invJparent, &detJparent);CHKERRQ(ierr); 1249 for (i = 0; i < nPoints; i++) { 1250 CoordinatesRefToReal(spdim, spdim, v0, J, &pointsRef[i*spdim],vtmp);CHKERRQ(ierr); 1251 CoordinatesRealToRef(spdim, spdim, v0parent, invJparent, vtmp, &pointsReal[i*spdim]);CHKERRQ(ierr); 1252 } 1253 ierr = PetscFEGetTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);CHKERRQ(ierr); 1254 offset = 0; 1255 for (i = 0; i < fSize; i++) { 1256 PetscInt qPoints; 1257 PetscQuadrature quad; 1258 1259 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 1260 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 1261 for (j = 0; j < fSize; j++) { 1262 PetscScalar val = 0.; 1263 1264 for (k = 0; k < qPoints; k++) { 1265 val += B[((offset + k) * fSize + j) * fComp] * weights[k]; 1266 } 1267 MatSetValue(Bmat,i,j,val,INSERT_VALUES);CHKERRQ(ierr); 1268 } 1269 offset += qPoints; 1270 } 1271 ierr = PetscFERestoreTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);CHKERRQ(ierr); 1272 ierr = MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1273 ierr = MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1274 ierr = MatMatSolve(Amat,Bmat,Xmat);CHKERRQ(ierr); 1275 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 1276 ierr = PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);CHKERRQ(ierr); 1277 childOffsets[0] = 0; 1278 for (i = 0; i < closureSize; i++) { 1279 PetscInt p = closure[2*i]; 1280 PetscInt dof; 1281 1282 if (numFields > 1) { 1283 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 1284 } 1285 else { 1286 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 1287 } 1288 childOffsets[i+1]=childOffsets[i]+dof / fComp; 1289 } 1290 parentOffsets[0] = 0; 1291 for (i = 0; i < closureSizeP; i++) { 1292 PetscInt p = closureP[2*i]; 1293 PetscInt dof; 1294 1295 if (numFields > 1) { 1296 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 1297 } 1298 else { 1299 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 1300 } 1301 parentOffsets[i+1]=parentOffsets[i]+dof / fComp; 1302 } 1303 for (i = 0; i < closureSize; i++) { 1304 PetscInt conDof, conOff, aDof, aOff; 1305 PetscInt p = closure[2*i]; 1306 PetscInt o = closure[2*i+1]; 1307 1308 if (p < conStart || p >= conEnd) continue; 1309 if (numFields > 1) { 1310 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 1311 ierr = PetscSectionGetFieldOffset(cSec,p,f,&conOff);CHKERRQ(ierr); 1312 } 1313 else { 1314 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 1315 ierr = PetscSectionGetOffset(cSec,p,&conOff);CHKERRQ(ierr); 1316 } 1317 if (!conDof) continue; 1318 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 1319 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 1320 for (k = 0; k < aDof; k++) { 1321 PetscInt a = anchors[aOff + k]; 1322 PetscInt aSecDof, aSecOff; 1323 1324 if (numFields > 1) { 1325 ierr = PetscSectionGetFieldDof(section,a,f,&aSecDof);CHKERRQ(ierr); 1326 ierr = PetscSectionGetFieldOffset(section,a,f,&aSecOff);CHKERRQ(ierr); 1327 } 1328 else { 1329 ierr = PetscSectionGetDof(section,a,&aSecDof);CHKERRQ(ierr); 1330 ierr = PetscSectionGetOffset(section,a,&aSecOff);CHKERRQ(ierr); 1331 } 1332 if (!aSecDof) continue; 1333 1334 for (j = 0; j < closureSizeP; j++) { 1335 PetscInt q = closureP[2*j]; 1336 PetscInt oq = closureP[2*j+1]; 1337 1338 if (q == a) { 1339 PetscInt r, s, t; 1340 1341 for (r = childOffsets[i]; r < childOffsets[i+1]; r++) { 1342 for (s = parentOffsets[j]; s < parentOffsets[j+1]; s++) { 1343 PetscScalar val; 1344 PetscInt insertCol, insertRow; 1345 1346 ierr = MatGetValue(Xmat,r,s,&val);CHKERRQ(ierr); 1347 if (o >= 0) { 1348 insertRow = conOff + fComp * (r - childOffsets[i]); 1349 } 1350 else { 1351 insertRow = conOff + fComp * (childOffsets[i + 1] - 1 - r); 1352 } 1353 if (oq >= 0) { 1354 insertCol = aSecOff + fComp * (s - parentOffsets[j]); 1355 } 1356 else { 1357 insertCol = aSecOff + fComp * (parentOffsets[j + 1] - 1 - s); 1358 } 1359 for (t = 0; t < fComp; t++) { 1360 ierr = MatSetValue(cMat,insertRow + t,insertCol + t,val,INSERT_VALUES);CHKERRQ(ierr); 1361 } 1362 } 1363 } 1364 } 1365 } 1366 } 1367 } 1368 ierr = PetscFree2(childOffsets,parentOffsets);CHKERRQ(ierr); 1369 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1370 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 1371 } 1372 ierr = MatDestroy(&Amat);CHKERRQ(ierr); 1373 ierr = MatDestroy(&Bmat);CHKERRQ(ierr); 1374 ierr = MatDestroy(&Xmat);CHKERRQ(ierr); 1375 ierr = PetscFree3(weights,pointsRef,pointsReal);CHKERRQ(ierr); 1376 } 1377 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1378 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1379 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);CHKERRQ(ierr); 1380 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 1381 1382 PetscFunctionReturn(0); 1383 } 1384 1385 #undef __FUNCT__ 1386 #define __FUNCT__ "DMPlexComputeConstraintMatrix_Tree" 1387 PetscErrorCode DMPlexComputeConstraintMatrix_Tree(DM dm) 1388 { 1389 DM refTree; 1390 PetscDS ds; 1391 Mat refCmat, cMat; 1392 PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN; 1393 PetscScalar ***refPointFieldMats, *pointWork; 1394 PetscSection refConSec, conSec, refAnSec, anSec, section, refSection; 1395 IS refAnIS, anIS; 1396 const PetscInt *refAnchors, *anchors; 1397 PetscErrorCode ierr; 1398 1399 PetscFunctionBegin; 1400 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1401 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 1402 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1403 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr); 1404 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 1405 ierr = DMPlexComputeConstraintMatrix_ReferenceTree(refTree);CHKERRQ(ierr); 1406 ierr = DMPlexGetConstraintMatrix(refTree,&refCmat);CHKERRQ(ierr); 1407 ierr = DMPlexGetChart(refTree,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1408 ierr = DMPlexGetConstraints(refTree,&refAnSec,&refAnIS);CHKERRQ(ierr); 1409 ierr = DMPlexGetConstraints(dm,&anSec,&anIS);CHKERRQ(ierr); 1410 ierr = ISGetIndices(refAnIS,&refAnchors);CHKERRQ(ierr); 1411 ierr = ISGetIndices(anIS,&anchors);CHKERRQ(ierr); 1412 ierr = DMPlexGetConstraintSection(refTree,&refConSec);CHKERRQ(ierr); 1413 ierr = DMPlexGetConstraintSection(dm,&conSec);CHKERRQ(ierr); 1414 ierr = DMGetDefaultSection(refTree,&refSection);CHKERRQ(ierr); 1415 ierr = DMGetDefaultSection(dm,§ion);CHKERRQ(ierr); 1416 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1417 ierr = PetscSectionGetChart(conSec,&conStart,&conEnd);CHKERRQ(ierr); 1418 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);CHKERRQ(ierr); 1419 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);CHKERRQ(ierr); 1420 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 1421 ierr = PetscSectionGetMaxDof(refAnSec,&maxAnDof);CHKERRQ(ierr); 1422 ierr = PetscMalloc1(maxDof,&rows);CHKERRQ(ierr); 1423 ierr = PetscMalloc1(maxDof*maxAnDof,&cols);CHKERRQ(ierr); 1424 ierr = PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);CHKERRQ(ierr); 1425 1426 /* step 1: get submats for every constrained point in the reference tree */ 1427 for (p = pRefStart; p < pRefEnd; p++) { 1428 PetscInt parent, closureSize, *closure = NULL, pDof; 1429 1430 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 1431 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 1432 if (!pDof || parent == p) continue; 1433 1434 ierr = PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);CHKERRQ(ierr); 1435 ierr = PetscCalloc1(numFields,&refPointFieldN[p-pRefStart]);CHKERRQ(ierr); 1436 ierr = DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1437 for (f = 0; f < numFields; f++) { 1438 PetscInt cDof, cOff, numCols, r, i, fComp; 1439 PetscFE fe; 1440 1441 ierr = PetscDSGetDiscretization(ds,f,(PetscObject *) &fe);CHKERRQ(ierr); 1442 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 1443 1444 if (numFields > 1) { 1445 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 1446 ierr = PetscSectionGetFieldOffset(refConSec,p,f,&cOff);CHKERRQ(ierr); 1447 } 1448 else { 1449 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 1450 ierr = PetscSectionGetOffset(refConSec,p,&cOff);CHKERRQ(ierr); 1451 } 1452 1453 if (!cDof) continue; 1454 for (r = 0; r < cDof; r++) { 1455 rows[r] = cOff + r; 1456 } 1457 numCols = 0; 1458 for (i = 0; i < closureSize; i++) { 1459 PetscInt q = closure[2*i]; 1460 PetscInt o = closure[2*i+1]; 1461 PetscInt aDof, aOff, j; 1462 1463 if (numFields > 1) { 1464 ierr = PetscSectionGetFieldDof(refSection,q,f,&aDof);CHKERRQ(ierr); 1465 ierr = PetscSectionGetFieldOffset(refSection,q,f,&aOff);CHKERRQ(ierr); 1466 } 1467 else { 1468 ierr = PetscSectionGetDof(refSection,q,&aDof);CHKERRQ(ierr); 1469 ierr = PetscSectionGetOffset(refSection,q,&aOff);CHKERRQ(ierr); 1470 } 1471 1472 for (j = 0; j < aDof; j++) { 1473 PetscInt node = (o >= 0) ? (j / fComp) : ((aDof - 1 - j) / fComp); 1474 PetscInt comp = (j % fComp); 1475 1476 cols[numCols++] = aOff + node * fComp + comp; 1477 } 1478 } 1479 refPointFieldN[p-pRefStart][f] = numCols; 1480 ierr = PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 1481 ierr = MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 1482 } 1483 ierr = DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1484 } 1485 1486 /* step 2: compute the preorder */ 1487 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 1488 ierr = PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);CHKERRQ(ierr); 1489 for (p = pStart; p < pEnd; p++) { 1490 perm[p - pStart] = p; 1491 iperm[p - pStart] = p-pStart; 1492 } 1493 for (p = 0; p < pEnd - pStart;) { 1494 PetscInt point = perm[p]; 1495 PetscInt parent; 1496 1497 ierr = DMPlexGetTreeParent(dm,point,&parent,NULL);CHKERRQ(ierr); 1498 if (parent == point) { 1499 p++; 1500 } 1501 else { 1502 PetscInt size, closureSize, *closure = NULL, i; 1503 1504 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1505 for (i = 0; i < closureSize; i++) { 1506 PetscInt q = closure[2*i]; 1507 if (iperm[q-pStart] > iperm[point-pStart]) { 1508 /* swap */ 1509 perm[p] = q; 1510 perm[iperm[q-pStart]] = point; 1511 iperm[point-pStart] = iperm[q-pStart]; 1512 iperm[q-pStart] = p; 1513 break; 1514 } 1515 } 1516 size = closureSize; 1517 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1518 if (i == size) { 1519 p++; 1520 } 1521 } 1522 } 1523 1524 /* step 3: fill the constraint matrix */ 1525 /* we are going to use a preorder progressive fill strategy. Mat doesn't 1526 * allow progressive fill without assembly, so we are going to set up the 1527 * values outside of the Mat first. 1528 */ 1529 { 1530 PetscInt nRows, row, nnz; 1531 PetscBool done; 1532 const PetscInt *ia, *ja; 1533 PetscScalar *vals; 1534 1535 ierr = DMPlexGetConstraintMatrix(dm,&cMat);CHKERRQ(ierr); 1536 ierr = MatGetRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);CHKERRQ(ierr); 1537 if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not get RowIJ of constraint matrix"); 1538 nnz = ia[nRows]; 1539 /* malloc and then zero rows right before we fill them: this way valgrind 1540 * can tell if we are doing progressive fill in the wrong order */ 1541 ierr = PetscMalloc1(nnz,&vals);CHKERRQ(ierr); 1542 for (p = 0; p < pEnd - pStart; p++) { 1543 PetscInt parent, childid, closureSize, *closure = NULL; 1544 PetscInt point = perm[p], pointDof; 1545 1546 ierr = DMPlexGetTreeParent(dm,point,&parent,&childid);CHKERRQ(ierr); 1547 if ((point < conStart) || (point >= conEnd) || (parent == point)) continue; 1548 ierr = PetscSectionGetDof(conSec,point,&pointDof);CHKERRQ(ierr); 1549 if (!pointDof) continue; 1550 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1551 for (f = 0; f < numFields; f++) { 1552 PetscInt cDof, cOff, numCols, numFillCols, i, r, fComp, matOffset, offset; 1553 PetscScalar *pointMat; 1554 PetscFE fe; 1555 1556 ierr = PetscDSGetDiscretization(ds,f,(PetscObject *) &fe);CHKERRQ(ierr); 1557 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 1558 1559 if (numFields > 1) { 1560 ierr = PetscSectionGetFieldDof(conSec,point,f,&cDof);CHKERRQ(ierr); 1561 ierr = PetscSectionGetFieldOffset(conSec,point,f,&cOff);CHKERRQ(ierr); 1562 } 1563 else { 1564 ierr = PetscSectionGetDof(conSec,point,&cDof);CHKERRQ(ierr); 1565 ierr = PetscSectionGetOffset(conSec,point,&cOff);CHKERRQ(ierr); 1566 } 1567 if (!cDof) continue; 1568 1569 /* make sure that every row for this point is the same size */ 1570 #if defined(PETSC_USE_DEBUG) 1571 for (r = 0; r < cDof; r++) { 1572 if (cDof > 1 && r) { 1573 if ((ia[cOff+r+1]-ia[cOff+r]) != (ia[cOff+r]-ia[cOff+r-1])) { 1574 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Two point rows have different nnz: %d vs. %d", (ia[rows[r]+1]-ia[rows[r]]), (ia[rows[r]]-ia[rows[r]-1])); 1575 } 1576 } 1577 } 1578 #endif 1579 /* zero rows */ 1580 for (i = ia[cOff] ; i< ia[cOff+cDof];i++) { 1581 vals[i] = 0.; 1582 } 1583 matOffset = ia[cOff]; 1584 numFillCols = ia[cOff+1] - matOffset; 1585 pointMat = refPointFieldMats[childid-pRefStart][f]; 1586 numCols = refPointFieldN[childid-pRefStart][f]; 1587 offset = 0; 1588 for (i = 0; i < closureSize; i++) { 1589 PetscInt q = closure[2*i]; 1590 PetscInt o = closure[2*i+1]; 1591 PetscInt aDof, aOff, j, k, qConDof, qConOff; 1592 1593 qConDof = qConOff = 0; 1594 if (numFields > 1) { 1595 ierr = PetscSectionGetFieldDof(section,q,f,&aDof);CHKERRQ(ierr); 1596 ierr = PetscSectionGetFieldOffset(section,q,f,&aOff);CHKERRQ(ierr); 1597 if (q >= conStart && q < conEnd) { 1598 ierr = PetscSectionGetFieldDof(conSec,q,f,&qConDof);CHKERRQ(ierr); 1599 ierr = PetscSectionGetFieldOffset(conSec,q,f,&qConOff);CHKERRQ(ierr); 1600 } 1601 } 1602 else { 1603 ierr = PetscSectionGetDof(section,q,&aDof);CHKERRQ(ierr); 1604 ierr = PetscSectionGetOffset(section,q,&aOff);CHKERRQ(ierr); 1605 if (q >= conStart && q < conEnd) { 1606 ierr = PetscSectionGetDof(conSec,q,&qConDof);CHKERRQ(ierr); 1607 ierr = PetscSectionGetOffset(conSec,q,&qConOff);CHKERRQ(ierr); 1608 } 1609 } 1610 if (!aDof) continue; 1611 if (qConDof) { 1612 /* this point has anchors: its rows of the matrix should already 1613 * be filled, thanks to preordering */ 1614 /* first multiply into pointWork, then set in matrix */ 1615 PetscInt aMatOffset = ia[qConOff]; 1616 PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset; 1617 for (r = 0; r < cDof; r++) { 1618 for (j = 0; j < aNumFillCols; j++) { 1619 PetscScalar inVal = 0; 1620 for (k = 0; k < aDof; k++) { 1621 PetscInt node = (o >= 0) ? (k / fComp) : ((aDof - 1 - k) / fComp); 1622 PetscInt comp = (k % fComp); 1623 PetscInt col = node * fComp + comp; 1624 1625 inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j]; 1626 } 1627 pointWork[r * aNumFillCols + j] = inVal; 1628 } 1629 } 1630 /* assume that the columns are sorted, spend less time searching */ 1631 for (j = 0, k = 0; j < aNumFillCols; j++) { 1632 PetscInt col = ja[aMatOffset + j]; 1633 for (;k < numFillCols; k++) { 1634 if (ja[matOffset + k] == col) { 1635 break; 1636 } 1637 } 1638 if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, col); 1639 for (r = 0; r < cDof; r++) { 1640 vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j]; 1641 } 1642 } 1643 } 1644 else { 1645 /* find where to put this portion of pointMat into the matrix */ 1646 for (k = 0; k < numFillCols; k++) { 1647 if (ja[matOffset + k] == aOff) { 1648 break; 1649 } 1650 } 1651 if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, aOff); 1652 for (j = 0; j < aDof; j++) { 1653 PetscInt node = (o >= 0) ? (j / fComp) : ((aDof - 1 - j) / fComp); 1654 PetscInt comp = (j % fComp); 1655 PetscInt col = node * fComp + comp; 1656 for (r = 0; r < cDof; r++) { 1657 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + col]; 1658 } 1659 } 1660 } 1661 offset += aDof; 1662 } 1663 } 1664 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1665 } 1666 for (row = 0; row < nRows; row++) { 1667 ierr = MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);CHKERRQ(ierr); 1668 } 1669 ierr = MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);CHKERRQ(ierr); 1670 if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix"); 1671 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1672 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1673 ierr = PetscFree(vals);CHKERRQ(ierr); 1674 } 1675 1676 /* clean up */ 1677 ierr = ISRestoreIndices(refAnIS,&refAnchors);CHKERRQ(ierr); 1678 ierr = ISRestoreIndices(anIS,&anchors);CHKERRQ(ierr); 1679 ierr = PetscFree2(perm,iperm);CHKERRQ(ierr); 1680 ierr = PetscFree(rows);CHKERRQ(ierr); 1681 ierr = PetscFree(cols);CHKERRQ(ierr); 1682 ierr = PetscFree(pointWork);CHKERRQ(ierr); 1683 for (p = pRefStart; p < pRefEnd; p++) { 1684 PetscInt parent, pDof; 1685 1686 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 1687 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 1688 if (!pDof || parent == p) continue; 1689 1690 for (f = 0; f < numFields; f++) { 1691 PetscInt cDof; 1692 1693 if (numFields > 1) { 1694 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 1695 } 1696 else { 1697 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 1698 } 1699 1700 if (!cDof) continue; 1701 ierr = PetscFree(refPointFieldMats[p - pRefStart][f]);CHKERRQ(ierr); 1702 } 1703 ierr = PetscFree(refPointFieldMats[p - pRefStart]);CHKERRQ(ierr); 1704 ierr = PetscFree(refPointFieldN[p - pRefStart]);CHKERRQ(ierr); 1705 } 1706 ierr = PetscFree(refPointFieldMats);CHKERRQ(ierr); 1707 ierr = PetscFree(refPointFieldN);CHKERRQ(ierr); 1708 PetscFunctionReturn(0); 1709 } 1710 1711 #undef __FUNCT__ 1712 #define __FUNCT__ "DMPlexTreeRefineCell" 1713 /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of 1714 * a non-conforming mesh. Local refinement comes later */ 1715 PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm) 1716 { 1717 DM K; 1718 PetscInt rank; 1719 PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd; 1720 PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations; 1721 PetscInt *Kembedding; 1722 PetscInt *cellClosure=NULL, nc; 1723 PetscScalar *newVertexCoords; 1724 PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset; 1725 PetscSection parentSection; 1726 PetscErrorCode ierr; 1727 1728 PetscFunctionBegin; 1729 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 1730 ierr = DMPlexGetDimension(dm,&dim);CHKERRQ(ierr); 1731 ierr = DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);CHKERRQ(ierr); 1732 ierr = DMPlexSetDimension(*ncdm,dim);CHKERRQ(ierr); 1733 1734 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1735 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection);CHKERRQ(ierr); 1736 ierr = DMPlexGetReferenceTree(dm,&K);CHKERRQ(ierr); 1737 if (!rank) { 1738 /* compute the new charts */ 1739 ierr = PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd);CHKERRQ(ierr); 1740 offset = 0; 1741 for (d = 0; d <= dim; d++) { 1742 PetscInt pOldCount, kStart, kEnd, k; 1743 1744 pNewStart[d] = offset; 1745 ierr = DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d]);CHKERRQ(ierr); 1746 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr); 1747 pOldCount = pOldEnd[d] - pOldStart[d]; 1748 /* adding the new points */ 1749 pNewCount[d] = pOldCount + kEnd - kStart; 1750 if (!d) { 1751 /* removing the cell */ 1752 pNewCount[d]--; 1753 } 1754 for (k = kStart; k < kEnd; k++) { 1755 PetscInt parent; 1756 ierr = DMPlexGetTreeParent(K,k,&parent,NULL);CHKERRQ(ierr); 1757 if (parent == k) { 1758 /* avoid double counting points that won't actually be new */ 1759 pNewCount[d]--; 1760 } 1761 } 1762 pNewEnd[d] = pNewStart[d] + pNewCount[d]; 1763 offset = pNewEnd[d]; 1764 1765 } 1766 if (cell < pOldStart[0] || cell >= pOldEnd[0]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"%d not in cell range [%d, %d)", cell, pOldStart[0], pOldEnd[0]); 1767 /* get the current closure of the cell that we are removing */ 1768 ierr = DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);CHKERRQ(ierr); 1769 1770 ierr = PetscMalloc1(pNewEnd[dim],&newConeSizes);CHKERRQ(ierr); 1771 { 1772 PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j; 1773 1774 ierr = DMPlexGetChart(K,&kStart,&kEnd);CHKERRQ(ierr); 1775 ierr = PetscMalloc4(kEnd-kStart,&Kembedding,kEnd-kStart,&perm,kEnd-kStart,&iperm,kEnd-kStart,&preOrient);CHKERRQ(ierr); 1776 1777 for (k = kStart; k < kEnd; k++) { 1778 perm[k - kStart] = k; 1779 iperm [k - kStart] = k - kStart; 1780 preOrient[k - kStart] = 0; 1781 } 1782 1783 ierr = DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);CHKERRQ(ierr); 1784 for (j = 1; j < closureSizeK; j++) { 1785 PetscInt parentOrientA = closureK[2*j+1]; 1786 PetscInt parentOrientB = cellClosure[2*j+1]; 1787 PetscInt p, q; 1788 1789 p = closureK[2*j]; 1790 q = cellClosure[2*j]; 1791 for (d = 0; d <= dim; d++) { 1792 if (q >= pOldStart[d] && q < pOldEnd[d]) { 1793 Kembedding[p] = (q - pOldStart[d]) + pNewStart[d]; 1794 } 1795 } 1796 if (parentOrientA != parentOrientB) { 1797 PetscInt numChildren, i; 1798 const PetscInt *children; 1799 1800 ierr = DMPlexGetTreeChildren(K,p,&numChildren,&children);CHKERRQ(ierr); 1801 for (i = 0; i < numChildren; i++) { 1802 PetscInt kPerm, oPerm; 1803 1804 k = children[i]; 1805 ierr = DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm);CHKERRQ(ierr); 1806 /* perm = what refTree position I'm in */ 1807 perm[kPerm-kStart] = k; 1808 /* iperm = who is at this position */ 1809 iperm[k-kStart] = kPerm-kStart; 1810 preOrient[kPerm-kStart] = oPerm; 1811 } 1812 } 1813 } 1814 ierr = DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);CHKERRQ(ierr); 1815 } 1816 ierr = PetscSectionSetChart(parentSection,0,pNewEnd[dim]);CHKERRQ(ierr); 1817 offset = 0; 1818 numNewCones = 0; 1819 for (d = 0; d <= dim; d++) { 1820 PetscInt kStart, kEnd, k; 1821 PetscInt p; 1822 PetscInt size; 1823 1824 for (p = pOldStart[d]; p < pOldEnd[d]; p++) { 1825 /* skip cell 0 */ 1826 if (p == cell) continue; 1827 /* old cones to new cones */ 1828 ierr = DMPlexGetConeSize(dm,p,&size);CHKERRQ(ierr); 1829 newConeSizes[offset++] = size; 1830 numNewCones += size; 1831 } 1832 1833 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr); 1834 for (k = kStart; k < kEnd; k++) { 1835 PetscInt kParent; 1836 1837 ierr = DMPlexGetTreeParent(K,k,&kParent,NULL);CHKERRQ(ierr); 1838 if (kParent != k) { 1839 Kembedding[k] = offset; 1840 ierr = DMPlexGetConeSize(K,k,&size);CHKERRQ(ierr); 1841 newConeSizes[offset++] = size; 1842 numNewCones += size; 1843 if (kParent != 0) { 1844 ierr = PetscSectionSetDof(parentSection,Kembedding[k],1);CHKERRQ(ierr); 1845 } 1846 } 1847 } 1848 } 1849 1850 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 1851 ierr = PetscSectionGetStorageSize(parentSection,&numPointsWithParents);CHKERRQ(ierr); 1852 ierr = PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations);CHKERRQ(ierr); 1853 ierr = PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs);CHKERRQ(ierr); 1854 1855 /* fill new cones */ 1856 offset = 0; 1857 for (d = 0; d <= dim; d++) { 1858 PetscInt kStart, kEnd, k, l; 1859 PetscInt p; 1860 PetscInt size; 1861 const PetscInt *cone, *orientation; 1862 1863 for (p = pOldStart[d]; p < pOldEnd[d]; p++) { 1864 /* skip cell 0 */ 1865 if (p == cell) continue; 1866 /* old cones to new cones */ 1867 ierr = DMPlexGetConeSize(dm,p,&size);CHKERRQ(ierr); 1868 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 1869 ierr = DMPlexGetConeOrientation(dm,p,&orientation);CHKERRQ(ierr); 1870 for (l = 0; l < size; l++) { 1871 newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1]; 1872 newOrientations[offset++] = orientation[l]; 1873 } 1874 } 1875 1876 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr); 1877 for (k = kStart; k < kEnd; k++) { 1878 PetscInt kPerm = perm[k], kParent; 1879 PetscInt preO = preOrient[k]; 1880 1881 ierr = DMPlexGetTreeParent(K,k,&kParent,NULL);CHKERRQ(ierr); 1882 if (kParent != k) { 1883 /* embed new cones */ 1884 ierr = DMPlexGetConeSize(K,k,&size);CHKERRQ(ierr); 1885 ierr = DMPlexGetCone(K,kPerm,&cone);CHKERRQ(ierr); 1886 ierr = DMPlexGetConeOrientation(K,kPerm,&orientation);CHKERRQ(ierr); 1887 for (l = 0; l < size; l++) { 1888 PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size); 1889 PetscInt newO, lSize, oTrue; 1890 1891 q = iperm[cone[m]]; 1892 newCones[offset] = Kembedding[q]; 1893 ierr = DMPlexGetConeSize(K,q,&lSize);CHKERRQ(ierr); 1894 oTrue = orientation[m]; 1895 oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2); 1896 newO = DihedralCompose(lSize,oTrue,preOrient[q]); 1897 newOrientations[offset++] = newO; 1898 } 1899 if (kParent != 0) { 1900 PetscInt newPoint = Kembedding[kParent]; 1901 ierr = PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset);CHKERRQ(ierr); 1902 parents[pOffset] = newPoint; 1903 childIDs[pOffset] = k; 1904 } 1905 } 1906 } 1907 } 1908 1909 ierr = PetscMalloc1(dim*(pNewEnd[dim]-pNewStart[dim]),&newVertexCoords);CHKERRQ(ierr); 1910 1911 /* fill coordinates */ 1912 offset = 0; 1913 { 1914 PetscInt k, kStart, kEnd, l; 1915 PetscSection vSection; 1916 PetscInt v; 1917 Vec coords; 1918 PetscScalar *coordvals; 1919 PetscInt dof, off; 1920 PetscScalar v0[3], J[9], detJ; 1921 1922 #if defined(PETSC_USE_DEBUG) 1923 ierr = DMPlexGetHeightStratum(K,0,&kStart,&kEnd);CHKERRQ(ierr); 1924 for (k = kStart; k < kEnd; k++) { 1925 ierr = DMPlexComputeCellGeometry(K, k, v0, J, NULL, &detJ);CHKERRQ(ierr); 1926 if (detJ <= 0.) SETERRQ1 (PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %d has bad determinant",k); 1927 } 1928 #endif 1929 ierr = DMPlexComputeCellGeometry(dm, cell, v0, J, NULL, &detJ);CHKERRQ(ierr); 1930 ierr = DMGetCoordinateSection(dm,&vSection);CHKERRQ(ierr); 1931 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 1932 ierr = VecGetArray(coords,&coordvals);CHKERRQ(ierr); 1933 for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) { 1934 1935 ierr = PetscSectionGetDof(vSection,v,&dof);CHKERRQ(ierr); 1936 ierr = PetscSectionGetOffset(vSection,v,&off);CHKERRQ(ierr); 1937 for (l = 0; l < dof; l++) { 1938 newVertexCoords[offset++] = coordvals[off + l]; 1939 } 1940 } 1941 ierr = VecRestoreArray(coords,&coordvals);CHKERRQ(ierr); 1942 1943 ierr = DMGetCoordinateSection(K,&vSection);CHKERRQ(ierr); 1944 ierr = DMGetCoordinatesLocal(K,&coords);CHKERRQ(ierr); 1945 ierr = VecGetArray(coords,&coordvals);CHKERRQ(ierr); 1946 ierr = DMPlexGetDepthStratum(K,0,&kStart,&kEnd);CHKERRQ(ierr); 1947 for (v = kStart; v < kEnd; v++) { 1948 PetscInt vPerm = perm[v]; 1949 PetscInt kParent; 1950 1951 ierr = DMPlexGetTreeParent(K,v,&kParent,NULL);CHKERRQ(ierr); 1952 if (kParent != v) { 1953 /* this is a new vertex */ 1954 ierr = PetscSectionGetOffset(vSection,vPerm,&off);CHKERRQ(ierr); 1955 CoordinatesRefToReal(dim, dim, v0, J, &coordvals[off],&newVertexCoords[offset]);CHKERRQ(ierr); 1956 offset += dim; 1957 } 1958 } 1959 ierr = VecRestoreArray(coords,&coordvals);CHKERRQ(ierr); 1960 } 1961 1962 /* need to reverse the order of pNewCount: vertices first, cells last */ 1963 for (d = 0; d < (dim + 1) / 2; d++) { 1964 PetscInt tmp; 1965 1966 tmp = pNewCount[d]; 1967 pNewCount[d] = pNewCount[dim - d]; 1968 pNewCount[dim - d] = tmp; 1969 } 1970 1971 ierr = DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);CHKERRQ(ierr); 1972 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr); 1973 ierr = DMPlexSetTree(*ncdm,parentSection,parents,childIDs);CHKERRQ(ierr); 1974 1975 /* clean up */ 1976 ierr = DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);CHKERRQ(ierr); 1977 ierr = PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);CHKERRQ(ierr); 1978 ierr = PetscFree(newConeSizes);CHKERRQ(ierr); 1979 ierr = PetscFree2(newCones,newOrientations);CHKERRQ(ierr); 1980 ierr = PetscFree(newVertexCoords);CHKERRQ(ierr); 1981 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 1982 ierr = PetscFree4(Kembedding,perm,iperm,preOrient);CHKERRQ(ierr); 1983 } 1984 else { 1985 PetscInt p, counts[4]; 1986 PetscInt *coneSizes, *cones, *orientations; 1987 Vec coordVec; 1988 PetscScalar *coords; 1989 1990 for (d = 0; d <= dim; d++) { 1991 PetscInt dStart, dEnd; 1992 1993 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr); 1994 counts[d] = dEnd - dStart; 1995 } 1996 ierr = PetscMalloc1(pEnd-pStart,&coneSizes);CHKERRQ(ierr); 1997 for (p = pStart; p < pEnd; p++) { 1998 ierr = DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);CHKERRQ(ierr); 1999 } 2000 ierr = DMPlexGetCones(dm, &cones);CHKERRQ(ierr); 2001 ierr = DMPlexGetConeOrientations(dm, &orientations);CHKERRQ(ierr); 2002 ierr = DMGetCoordinatesLocal(dm,&coordVec);CHKERRQ(ierr); 2003 ierr = VecGetArray(coordVec,&coords);CHKERRQ(ierr); 2004 2005 ierr = PetscSectionSetChart(parentSection,pStart,pEnd);CHKERRQ(ierr); 2006 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 2007 ierr = DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);CHKERRQ(ierr); 2008 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr); 2009 ierr = DMPlexSetTree(*ncdm,parentSection,NULL,NULL);CHKERRQ(ierr); 2010 ierr = VecRestoreArray(coordVec,&coords);CHKERRQ(ierr); 2011 } 2012 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 2013 2014 PetscFunctionReturn(0); 2015 } 2016