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