1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/isimpl.h> 3 #include <petsc/private/petscfeimpl.h> 4 #include <petscsf.h> 5 #include <petscds.h> 6 7 /* hierarchy routines */ 8 9 /*@ 10 DMPlexSetReferenceTree - set the reference tree for hierarchically non-conforming meshes. 11 12 Not collective 13 14 Input Parameters: 15 + dm - The DMPlex object 16 - ref - The reference tree DMPlex object 17 18 Level: intermediate 19 20 .seealso: `DMPlexGetReferenceTree()`, `DMPlexCreateDefaultReferenceTree()` 21 @*/ 22 PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref) 23 { 24 DM_Plex *mesh = (DM_Plex *)dm->data; 25 26 PetscFunctionBegin; 27 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 28 if (ref) {PetscValidHeaderSpecific(ref, DM_CLASSID, 2);} 29 PetscCall(PetscObjectReference((PetscObject)ref)); 30 PetscCall(DMDestroy(&mesh->referenceTree)); 31 mesh->referenceTree = ref; 32 PetscFunctionReturn(0); 33 } 34 35 /*@ 36 DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes. 37 38 Not collective 39 40 Input Parameters: 41 . dm - The DMPlex object 42 43 Output Parameters: 44 . ref - The reference tree DMPlex object 45 46 Level: intermediate 47 48 .seealso: `DMPlexSetReferenceTree()`, `DMPlexCreateDefaultReferenceTree()` 49 @*/ 50 PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref) 51 { 52 DM_Plex *mesh = (DM_Plex *)dm->data; 53 54 PetscFunctionBegin; 55 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 56 PetscValidPointer(ref,2); 57 *ref = mesh->referenceTree; 58 PetscFunctionReturn(0); 59 } 60 61 static PetscErrorCode DMPlexReferenceTreeGetChildSymmetry_Default(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB) 62 { 63 PetscInt coneSize, dStart, dEnd, dim, ABswap, oAvert, oBvert, ABswapVert; 64 65 PetscFunctionBegin; 66 if (parentOrientA == parentOrientB) { 67 if (childOrientB) *childOrientB = childOrientA; 68 if (childB) *childB = childA; 69 PetscFunctionReturn(0); 70 } 71 for (dim = 0; dim < 3; dim++) { 72 PetscCall(DMPlexGetDepthStratum(dm,dim,&dStart,&dEnd)); 73 if (parent >= dStart && parent <= dEnd) { 74 break; 75 } 76 } 77 PetscCheck(dim <= 2,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot perform child symmetry for %" PetscInt_FMT "-cells",dim); 78 PetscCheck(dim,PETSC_COMM_SELF,PETSC_ERR_PLIB,"A vertex has no children"); 79 if (childA < dStart || childA >= dEnd) { 80 /* this is a lower-dimensional child: bootstrap */ 81 PetscInt size, i, sA = -1, sB, sOrientB, sConeSize; 82 const PetscInt *supp, *coneA, *coneB, *oA, *oB; 83 84 PetscCall(DMPlexGetSupportSize(dm,childA,&size)); 85 PetscCall(DMPlexGetSupport(dm,childA,&supp)); 86 87 /* find a point sA in supp(childA) that has the same parent */ 88 for (i = 0; i < size; i++) { 89 PetscInt sParent; 90 91 sA = supp[i]; 92 if (sA == parent) continue; 93 PetscCall(DMPlexGetTreeParent(dm,sA,&sParent,NULL)); 94 if (sParent == parent) { 95 break; 96 } 97 } 98 PetscCheck(i != size,PETSC_COMM_SELF,PETSC_ERR_PLIB,"could not find support in children"); 99 /* find out which point sB is in an equivalent position to sA under 100 * parentOrientB */ 101 PetscCall(DMPlexReferenceTreeGetChildSymmetry_Default(dm,parent,parentOrientA,0,sA,parentOrientB,&sOrientB,&sB)); 102 PetscCall(DMPlexGetConeSize(dm,sA,&sConeSize)); 103 PetscCall(DMPlexGetCone(dm,sA,&coneA)); 104 PetscCall(DMPlexGetCone(dm,sB,&coneB)); 105 PetscCall(DMPlexGetConeOrientation(dm,sA,&oA)); 106 PetscCall(DMPlexGetConeOrientation(dm,sB,&oB)); 107 /* step through the cone of sA in natural order */ 108 for (i = 0; i < sConeSize; i++) { 109 if (coneA[i] == childA) { 110 /* if childA is at position i in coneA, 111 * then we want the point that is at sOrientB*i in coneB */ 112 PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize -(sOrientB+1) - i) % sConeSize); 113 if (childB) *childB = coneB[j]; 114 if (childOrientB) { 115 DMPolytopeType ct; 116 PetscInt oBtrue; 117 118 PetscCall(DMPlexGetConeSize(dm,childA,&coneSize)); 119 /* compose sOrientB and oB[j] */ 120 PetscCheck(coneSize == 0 || coneSize == 2,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected a vertex or an edge"); 121 ct = coneSize ? DM_POLYTOPE_SEGMENT : DM_POLYTOPE_POINT; 122 /* we may have to flip an edge */ 123 oBtrue = (sOrientB >= 0) ? oB[j] : DMPolytopeTypeComposeOrientation(ct, -1, oB[j]); 124 oBtrue = DMPolytopeConvertNewOrientation_Internal(ct, oBtrue); 125 ABswap = DihedralSwap(coneSize,DMPolytopeConvertNewOrientation_Internal(ct, oA[i]),oBtrue); 126 *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap); 127 } 128 break; 129 } 130 } 131 PetscCheck(i != sConeSize,PETSC_COMM_SELF,PETSC_ERR_PLIB,"support cone mismatch"); 132 PetscFunctionReturn(0); 133 } 134 /* get the cone size and symmetry swap */ 135 PetscCall(DMPlexGetConeSize(dm,parent,&coneSize)); 136 ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB); 137 if (dim == 2) { 138 /* orientations refer to cones: we want them to refer to vertices: 139 * if it's a rotation, they are the same, but if the order is reversed, a 140 * permutation that puts side i first does *not* put vertex i first */ 141 oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1); 142 oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1); 143 ABswapVert = DihedralSwap(coneSize, oAvert, oBvert); 144 } else { 145 ABswapVert = ABswap; 146 } 147 if (childB) { 148 /* assume that each child corresponds to a vertex, in the same order */ 149 PetscInt p, posA = -1, numChildren, i; 150 const PetscInt *children; 151 152 /* count which position the child is in */ 153 PetscCall(DMPlexGetTreeChildren(dm,parent,&numChildren,&children)); 154 for (i = 0; i < numChildren; i++) { 155 p = children[i]; 156 if (p == childA) { 157 posA = i; 158 break; 159 } 160 } 161 if (posA >= coneSize) { 162 /* this is the triangle in the middle of a uniformly refined triangle: it 163 * is invariant */ 164 PetscCheck(dim == 2 && posA == 3,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Expected a middle triangle, got something else"); 165 *childB = childA; 166 } 167 else { 168 /* figure out position B by applying ABswapVert */ 169 PetscInt posB; 170 171 posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize -(ABswapVert + 1) - posA) % coneSize); 172 if (childB) *childB = children[posB]; 173 } 174 } 175 if (childOrientB) *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap); 176 PetscFunctionReturn(0); 177 } 178 179 /*@ 180 DMPlexReferenceTreeGetChildSymmetry - Given a reference tree, transform a childid and orientation from one parent frame to another 181 182 Input Parameters: 183 + dm - the reference tree DMPlex object 184 . parent - the parent point 185 . parentOrientA - the reference orientation for describing the parent 186 . childOrientA - the reference orientation for describing the child 187 . childA - the reference childID for describing the child 188 - parentOrientB - the new orientation for describing the parent 189 190 Output Parameters: 191 + childOrientB - if not NULL, set to the new oreintation for describing the child 192 - childB - if not NULL, the new childID for describing the child 193 194 Level: developer 195 196 .seealso: `DMPlexGetReferenceTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetTree()` 197 @*/ 198 PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB) 199 { 200 DM_Plex *mesh = (DM_Plex *)dm->data; 201 202 PetscFunctionBegin; 203 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 204 PetscCheck(mesh->getchildsymmetry,PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlexReferenceTreeGetChildSymmetry not implemented"); 205 PetscCall(mesh->getchildsymmetry(dm,parent,parentOrientA,childOrientA,childA,parentOrientB,childOrientB,childB)); 206 PetscFunctionReturn(0); 207 } 208 209 static PetscErrorCode DMPlexSetTree_Internal(DM,PetscSection,PetscInt*,PetscInt*,PetscBool,PetscBool); 210 211 PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[]) 212 { 213 PetscFunctionBegin; 214 PetscCall(DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_TRUE,PETSC_FALSE)); 215 PetscFunctionReturn(0); 216 } 217 218 PetscErrorCode DMPlexCreateReferenceTree_Union(DM K, DM Kref, const char *labelName, DM *ref) 219 { 220 MPI_Comm comm; 221 PetscInt dim, p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs; 222 PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts; 223 DMLabel identity, identityRef; 224 PetscSection unionSection, unionConeSection, parentSection; 225 PetscScalar *unionCoords; 226 IS perm; 227 228 PetscFunctionBegin; 229 comm = PetscObjectComm((PetscObject)K); 230 PetscCall(DMGetDimension(K, &dim)); 231 PetscCall(DMPlexGetChart(K, &pStart, &pEnd)); 232 PetscCall(DMGetLabel(K, labelName, &identity)); 233 PetscCall(DMGetLabel(Kref, labelName, &identityRef)); 234 PetscCall(DMPlexGetChart(Kref, &pRefStart, &pRefEnd)); 235 PetscCall(PetscSectionCreate(comm, &unionSection)); 236 PetscCall(PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart))); 237 /* count points that will go in the union */ 238 for (p = pStart; p < pEnd; p++) { 239 PetscCall(PetscSectionSetDof(unionSection, p - pStart, 1)); 240 } 241 for (p = pRefStart; p < pRefEnd; p++) { 242 PetscInt q, qSize; 243 PetscCall(DMLabelGetValue(identityRef, p, &q)); 244 PetscCall(DMLabelGetStratumSize(identityRef, q, &qSize)); 245 if (qSize > 1) { 246 PetscCall(PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1)); 247 } 248 } 249 PetscCall(PetscMalloc1(pEnd - pStart + pRefEnd - pRefStart,&permvals)); 250 offset = 0; 251 /* stratify points in the union by topological dimension */ 252 for (d = 0; d <= dim; d++) { 253 PetscInt cStart, cEnd, c; 254 255 PetscCall(DMPlexGetHeightStratum(K, d, &cStart, &cEnd)); 256 for (c = cStart; c < cEnd; c++) { 257 permvals[offset++] = c; 258 } 259 260 PetscCall(DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd)); 261 for (c = cStart; c < cEnd; c++) { 262 permvals[offset++] = c + (pEnd - pStart); 263 } 264 } 265 PetscCall(ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm)); 266 PetscCall(PetscSectionSetPermutation(unionSection,perm)); 267 PetscCall(PetscSectionSetUp(unionSection)); 268 PetscCall(PetscSectionGetStorageSize(unionSection,&numUnionPoints)); 269 PetscCall(PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints)); 270 /* count dimension points */ 271 for (d = 0; d <= dim; d++) { 272 PetscInt cStart, cOff, cOff2; 273 PetscCall(DMPlexGetHeightStratum(K,d,&cStart,NULL)); 274 PetscCall(PetscSectionGetOffset(unionSection,cStart-pStart,&cOff)); 275 if (d < dim) { 276 PetscCall(DMPlexGetHeightStratum(K,d+1,&cStart,NULL)); 277 PetscCall(PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2)); 278 } 279 else { 280 cOff2 = numUnionPoints; 281 } 282 numDimPoints[dim - d] = cOff2 - cOff; 283 } 284 PetscCall(PetscSectionCreate(comm, &unionConeSection)); 285 PetscCall(PetscSectionSetChart(unionConeSection, 0, numUnionPoints)); 286 /* count the cones in the union */ 287 for (p = pStart; p < pEnd; p++) { 288 PetscInt dof, uOff; 289 290 PetscCall(DMPlexGetConeSize(K, p, &dof)); 291 PetscCall(PetscSectionGetOffset(unionSection, p - pStart,&uOff)); 292 PetscCall(PetscSectionSetDof(unionConeSection, uOff, dof)); 293 coneSizes[uOff] = dof; 294 } 295 for (p = pRefStart; p < pRefEnd; p++) { 296 PetscInt dof, uDof, uOff; 297 298 PetscCall(DMPlexGetConeSize(Kref, p, &dof)); 299 PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof)); 300 PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff)); 301 if (uDof) { 302 PetscCall(PetscSectionSetDof(unionConeSection, uOff, dof)); 303 coneSizes[uOff] = dof; 304 } 305 } 306 PetscCall(PetscSectionSetUp(unionConeSection)); 307 PetscCall(PetscSectionGetStorageSize(unionConeSection,&numCones)); 308 PetscCall(PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations)); 309 /* write the cones in the union */ 310 for (p = pStart; p < pEnd; p++) { 311 PetscInt dof, uOff, c, cOff; 312 const PetscInt *cone, *orientation; 313 314 PetscCall(DMPlexGetConeSize(K, p, &dof)); 315 PetscCall(DMPlexGetCone(K, p, &cone)); 316 PetscCall(DMPlexGetConeOrientation(K, p, &orientation)); 317 PetscCall(PetscSectionGetOffset(unionSection, p - pStart,&uOff)); 318 PetscCall(PetscSectionGetOffset(unionConeSection,uOff,&cOff)); 319 for (c = 0; c < dof; c++) { 320 PetscInt e, eOff; 321 e = cone[c]; 322 PetscCall(PetscSectionGetOffset(unionSection, e - pStart, &eOff)); 323 unionCones[cOff + c] = eOff; 324 unionOrientations[cOff + c] = orientation[c]; 325 } 326 } 327 for (p = pRefStart; p < pRefEnd; p++) { 328 PetscInt dof, uDof, uOff, c, cOff; 329 const PetscInt *cone, *orientation; 330 331 PetscCall(DMPlexGetConeSize(Kref, p, &dof)); 332 PetscCall(DMPlexGetCone(Kref, p, &cone)); 333 PetscCall(DMPlexGetConeOrientation(Kref, p, &orientation)); 334 PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof)); 335 PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff)); 336 if (uDof) { 337 PetscCall(PetscSectionGetOffset(unionConeSection,uOff,&cOff)); 338 for (c = 0; c < dof; c++) { 339 PetscInt e, eOff, eDof; 340 341 e = cone[c]; 342 PetscCall(PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof)); 343 if (eDof) { 344 PetscCall(PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff)); 345 } 346 else { 347 PetscCall(DMLabelGetValue(identityRef, e, &e)); 348 PetscCall(PetscSectionGetOffset(unionSection, e - pStart, &eOff)); 349 } 350 unionCones[cOff + c] = eOff; 351 unionOrientations[cOff + c] = orientation[c]; 352 } 353 } 354 } 355 /* get the coordinates */ 356 { 357 PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff; 358 PetscSection KcoordsSec, KrefCoordsSec; 359 Vec KcoordsVec, KrefCoordsVec; 360 PetscScalar *Kcoords; 361 362 PetscCall(DMGetCoordinateSection(K, &KcoordsSec)); 363 PetscCall(DMGetCoordinatesLocal(K, &KcoordsVec)); 364 PetscCall(DMGetCoordinateSection(Kref, &KrefCoordsSec)); 365 PetscCall(DMGetCoordinatesLocal(Kref, &KrefCoordsVec)); 366 367 numVerts = numDimPoints[0]; 368 PetscCall(PetscMalloc1(numVerts * dim,&unionCoords)); 369 PetscCall(DMPlexGetDepthStratum(K,0,&vStart,&vEnd)); 370 371 offset = 0; 372 for (v = vStart; v < vEnd; v++) { 373 PetscCall(PetscSectionGetOffset(unionSection,v - pStart,&vOff)); 374 PetscCall(VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords)); 375 for (d = 0; d < dim; d++) { 376 unionCoords[offset * dim + d] = Kcoords[d]; 377 } 378 offset++; 379 } 380 PetscCall(DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd)); 381 for (v = vRefStart; v < vRefEnd; v++) { 382 PetscCall(PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof)); 383 PetscCall(PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff)); 384 PetscCall(VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords)); 385 if (vDof) { 386 for (d = 0; d < dim; d++) { 387 unionCoords[offset * dim + d] = Kcoords[d]; 388 } 389 offset++; 390 } 391 } 392 } 393 PetscCall(DMCreate(comm,ref)); 394 PetscCall(DMSetType(*ref,DMPLEX)); 395 PetscCall(DMSetDimension(*ref,dim)); 396 PetscCall(DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords)); 397 /* set the tree */ 398 PetscCall(PetscSectionCreate(comm,&parentSection)); 399 PetscCall(PetscSectionSetChart(parentSection,0,numUnionPoints)); 400 for (p = pRefStart; p < pRefEnd; p++) { 401 PetscInt uDof, uOff; 402 403 PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof)); 404 PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff)); 405 if (uDof) PetscCall(PetscSectionSetDof(parentSection,uOff,1)); 406 } 407 PetscCall(PetscSectionSetUp(parentSection)); 408 PetscCall(PetscSectionGetStorageSize(parentSection,&parentSize)); 409 PetscCall(PetscMalloc2(parentSize,&parents,parentSize,&childIDs)); 410 for (p = pRefStart; p < pRefEnd; p++) { 411 PetscInt uDof, uOff; 412 413 PetscCall(PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof)); 414 PetscCall(PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff)); 415 if (uDof) { 416 PetscInt pOff, parent, parentU; 417 PetscCall(PetscSectionGetOffset(parentSection,uOff,&pOff)); 418 PetscCall(DMLabelGetValue(identityRef,p,&parent)); 419 PetscCall(PetscSectionGetOffset(unionSection, parent - pStart,&parentU)); 420 parents[pOff] = parentU; 421 childIDs[pOff] = uOff; 422 } 423 } 424 PetscCall(DMPlexCreateReferenceTree_SetTree(*ref,parentSection,parents,childIDs)); 425 PetscCall(PetscSectionDestroy(&parentSection)); 426 PetscCall(PetscFree2(parents,childIDs)); 427 428 /* clean up */ 429 PetscCall(PetscSectionDestroy(&unionSection)); 430 PetscCall(PetscSectionDestroy(&unionConeSection)); 431 PetscCall(ISDestroy(&perm)); 432 PetscCall(PetscFree(unionCoords)); 433 PetscCall(PetscFree2(unionCones,unionOrientations)); 434 PetscCall(PetscFree2(coneSizes,numDimPoints)); 435 PetscFunctionReturn(0); 436 } 437 438 /*@ 439 DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement. 440 441 Collective 442 443 Input Parameters: 444 + comm - the MPI communicator 445 . dim - the spatial dimension 446 - simplex - Flag for simplex, otherwise use a tensor-product cell 447 448 Output Parameters: 449 . ref - the reference tree DMPlex object 450 451 Level: intermediate 452 453 .seealso: `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()` 454 @*/ 455 PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref) 456 { 457 DM_Plex *mesh; 458 DM K, Kref; 459 PetscInt p, pStart, pEnd; 460 DMLabel identity; 461 462 PetscFunctionBegin; 463 #if 1 464 comm = PETSC_COMM_SELF; 465 #endif 466 /* create a reference element */ 467 PetscCall(DMPlexCreateReferenceCell(comm, DMPolytopeTypeSimpleShape(dim, simplex), &K)); 468 PetscCall(DMCreateLabel(K, "identity")); 469 PetscCall(DMGetLabel(K, "identity", &identity)); 470 PetscCall(DMPlexGetChart(K, &pStart, &pEnd)); 471 for (p = pStart; p < pEnd; p++) { 472 PetscCall(DMLabelSetValue(identity, p, p)); 473 } 474 /* refine it */ 475 PetscCall(DMRefine(K,comm,&Kref)); 476 477 /* the reference tree is the union of these two, without duplicating 478 * points that appear in both */ 479 PetscCall(DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref)); 480 mesh = (DM_Plex *) (*ref)->data; 481 mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default; 482 PetscCall(DMDestroy(&K)); 483 PetscCall(DMDestroy(&Kref)); 484 PetscFunctionReturn(0); 485 } 486 487 static PetscErrorCode DMPlexTreeSymmetrize(DM dm) 488 { 489 DM_Plex *mesh = (DM_Plex *)dm->data; 490 PetscSection childSec, pSec; 491 PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT; 492 PetscInt *offsets, *children, pStart, pEnd; 493 494 PetscFunctionBegin; 495 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 496 PetscCall(PetscSectionDestroy(&mesh->childSection)); 497 PetscCall(PetscFree(mesh->children)); 498 pSec = mesh->parentSection; 499 if (!pSec) PetscFunctionReturn(0); 500 PetscCall(PetscSectionGetStorageSize(pSec,&pSize)); 501 for (p = 0; p < pSize; p++) { 502 PetscInt par = mesh->parents[p]; 503 504 parMax = PetscMax(parMax,par+1); 505 parMin = PetscMin(parMin,par); 506 } 507 if (parMin > parMax) { 508 parMin = -1; 509 parMax = -1; 510 } 511 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec)); 512 PetscCall(PetscSectionSetChart(childSec,parMin,parMax)); 513 for (p = 0; p < pSize; p++) { 514 PetscInt par = mesh->parents[p]; 515 516 PetscCall(PetscSectionAddDof(childSec,par,1)); 517 } 518 PetscCall(PetscSectionSetUp(childSec)); 519 PetscCall(PetscSectionGetStorageSize(childSec,&cSize)); 520 PetscCall(PetscMalloc1(cSize,&children)); 521 PetscCall(PetscCalloc1(parMax-parMin,&offsets)); 522 PetscCall(PetscSectionGetChart(pSec,&pStart,&pEnd)); 523 for (p = pStart; p < pEnd; p++) { 524 PetscInt dof, off, i; 525 526 PetscCall(PetscSectionGetDof(pSec,p,&dof)); 527 PetscCall(PetscSectionGetOffset(pSec,p,&off)); 528 for (i = 0; i < dof; i++) { 529 PetscInt par = mesh->parents[off + i], cOff; 530 531 PetscCall(PetscSectionGetOffset(childSec,par,&cOff)); 532 children[cOff + offsets[par-parMin]++] = p; 533 } 534 } 535 mesh->childSection = childSec; 536 mesh->children = children; 537 PetscCall(PetscFree(offsets)); 538 PetscFunctionReturn(0); 539 } 540 541 static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew) 542 { 543 PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL; 544 const PetscInt *vals; 545 PetscSection secNew; 546 PetscBool anyNew, globalAnyNew; 547 PetscBool compress; 548 549 PetscFunctionBegin; 550 PetscCall(PetscSectionGetChart(section,&pStart,&pEnd)); 551 PetscCall(ISGetLocalSize(is,&size)); 552 PetscCall(ISGetIndices(is,&vals)); 553 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew)); 554 PetscCall(PetscSectionSetChart(secNew,pStart,pEnd)); 555 for (i = 0; i < size; i++) { 556 PetscInt dof; 557 558 p = vals[i]; 559 if (p < pStart || p >= pEnd) continue; 560 PetscCall(PetscSectionGetDof(section, p, &dof)); 561 if (dof) break; 562 } 563 if (i == size) { 564 PetscCall(PetscSectionSetUp(secNew)); 565 anyNew = PETSC_FALSE; 566 compress = PETSC_FALSE; 567 sizeNew = 0; 568 } 569 else { 570 anyNew = PETSC_TRUE; 571 for (p = pStart; p < pEnd; p++) { 572 PetscInt dof, off; 573 574 PetscCall(PetscSectionGetDof(section, p, &dof)); 575 PetscCall(PetscSectionGetOffset(section, p, &off)); 576 for (i = 0; i < dof; i++) { 577 PetscInt q = vals[off + i], qDof = 0; 578 579 if (q >= pStart && q < pEnd) { 580 PetscCall(PetscSectionGetDof(section, q, &qDof)); 581 } 582 if (qDof) PetscCall(PetscSectionAddDof(secNew, p, qDof)); 583 else { 584 PetscCall(PetscSectionAddDof(secNew, p, 1)); 585 } 586 } 587 } 588 PetscCall(PetscSectionSetUp(secNew)); 589 PetscCall(PetscSectionGetStorageSize(secNew,&sizeNew)); 590 PetscCall(PetscMalloc1(sizeNew,&valsNew)); 591 compress = PETSC_FALSE; 592 for (p = pStart; p < pEnd; p++) { 593 PetscInt dof, off, count, offNew, dofNew; 594 595 PetscCall(PetscSectionGetDof(section, p, &dof)); 596 PetscCall(PetscSectionGetOffset(section, p, &off)); 597 PetscCall(PetscSectionGetDof(secNew, p, &dofNew)); 598 PetscCall(PetscSectionGetOffset(secNew, p, &offNew)); 599 count = 0; 600 for (i = 0; i < dof; i++) { 601 PetscInt q = vals[off + i], qDof = 0, qOff = 0, j; 602 603 if (q >= pStart && q < pEnd) { 604 PetscCall(PetscSectionGetDof(section, q, &qDof)); 605 PetscCall(PetscSectionGetOffset(section, q, &qOff)); 606 } 607 if (qDof) { 608 PetscInt oldCount = count; 609 610 for (j = 0; j < qDof; j++) { 611 PetscInt k, r = vals[qOff + j]; 612 613 for (k = 0; k < oldCount; k++) { 614 if (valsNew[offNew + k] == r) { 615 break; 616 } 617 } 618 if (k == oldCount) { 619 valsNew[offNew + count++] = r; 620 } 621 } 622 } 623 else { 624 PetscInt k, oldCount = count; 625 626 for (k = 0; k < oldCount; k++) { 627 if (valsNew[offNew + k] == q) { 628 break; 629 } 630 } 631 if (k == oldCount) { 632 valsNew[offNew + count++] = q; 633 } 634 } 635 } 636 if (count < dofNew) { 637 PetscCall(PetscSectionSetDof(secNew, p, count)); 638 compress = PETSC_TRUE; 639 } 640 } 641 } 642 PetscCall(ISRestoreIndices(is,&vals)); 643 PetscCall(MPIU_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew))); 644 if (!globalAnyNew) { 645 PetscCall(PetscSectionDestroy(&secNew)); 646 *sectionNew = NULL; 647 *isNew = NULL; 648 } 649 else { 650 PetscBool globalCompress; 651 652 PetscCall(MPIU_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew))); 653 if (compress) { 654 PetscSection secComp; 655 PetscInt *valsComp = NULL; 656 657 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp)); 658 PetscCall(PetscSectionSetChart(secComp,pStart,pEnd)); 659 for (p = pStart; p < pEnd; p++) { 660 PetscInt dof; 661 662 PetscCall(PetscSectionGetDof(secNew, p, &dof)); 663 PetscCall(PetscSectionSetDof(secComp, p, dof)); 664 } 665 PetscCall(PetscSectionSetUp(secComp)); 666 PetscCall(PetscSectionGetStorageSize(secComp,&sizeNew)); 667 PetscCall(PetscMalloc1(sizeNew,&valsComp)); 668 for (p = pStart; p < pEnd; p++) { 669 PetscInt dof, off, offNew, j; 670 671 PetscCall(PetscSectionGetDof(secNew, p, &dof)); 672 PetscCall(PetscSectionGetOffset(secNew, p, &off)); 673 PetscCall(PetscSectionGetOffset(secComp, p, &offNew)); 674 for (j = 0; j < dof; j++) { 675 valsComp[offNew + j] = valsNew[off + j]; 676 } 677 } 678 PetscCall(PetscSectionDestroy(&secNew)); 679 secNew = secComp; 680 PetscCall(PetscFree(valsNew)); 681 valsNew = valsComp; 682 } 683 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew)); 684 } 685 PetscFunctionReturn(0); 686 } 687 688 static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm) 689 { 690 PetscInt p, pStart, pEnd, *anchors, size; 691 PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT; 692 PetscSection aSec; 693 DMLabel canonLabel; 694 IS aIS; 695 696 PetscFunctionBegin; 697 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 698 PetscCall(DMPlexGetChart(dm,&pStart,&pEnd)); 699 PetscCall(DMGetLabel(dm,"canonical",&canonLabel)); 700 for (p = pStart; p < pEnd; p++) { 701 PetscInt parent; 702 703 if (canonLabel) { 704 PetscInt canon; 705 706 PetscCall(DMLabelGetValue(canonLabel,p,&canon)); 707 if (p != canon) continue; 708 } 709 PetscCall(DMPlexGetTreeParent(dm,p,&parent,NULL)); 710 if (parent != p) { 711 aMin = PetscMin(aMin,p); 712 aMax = PetscMax(aMax,p+1); 713 } 714 } 715 if (aMin > aMax) { 716 aMin = -1; 717 aMax = -1; 718 } 719 PetscCall(PetscSectionCreate(PETSC_COMM_SELF,&aSec)); 720 PetscCall(PetscSectionSetChart(aSec,aMin,aMax)); 721 for (p = aMin; p < aMax; p++) { 722 PetscInt parent, ancestor = p; 723 724 if (canonLabel) { 725 PetscInt canon; 726 727 PetscCall(DMLabelGetValue(canonLabel,p,&canon)); 728 if (p != canon) continue; 729 } 730 PetscCall(DMPlexGetTreeParent(dm,p,&parent,NULL)); 731 while (parent != ancestor) { 732 ancestor = parent; 733 PetscCall(DMPlexGetTreeParent(dm,ancestor,&parent,NULL)); 734 } 735 if (ancestor != p) { 736 PetscInt closureSize, *closure = NULL; 737 738 PetscCall(DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure)); 739 PetscCall(PetscSectionSetDof(aSec,p,closureSize)); 740 PetscCall(DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure)); 741 } 742 } 743 PetscCall(PetscSectionSetUp(aSec)); 744 PetscCall(PetscSectionGetStorageSize(aSec,&size)); 745 PetscCall(PetscMalloc1(size,&anchors)); 746 for (p = aMin; p < aMax; p++) { 747 PetscInt parent, ancestor = p; 748 749 if (canonLabel) { 750 PetscInt canon; 751 752 PetscCall(DMLabelGetValue(canonLabel,p,&canon)); 753 if (p != canon) continue; 754 } 755 PetscCall(DMPlexGetTreeParent(dm,p,&parent,NULL)); 756 while (parent != ancestor) { 757 ancestor = parent; 758 PetscCall(DMPlexGetTreeParent(dm,ancestor,&parent,NULL)); 759 } 760 if (ancestor != p) { 761 PetscInt j, closureSize, *closure = NULL, aOff; 762 763 PetscCall(PetscSectionGetOffset(aSec,p,&aOff)); 764 765 PetscCall(DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure)); 766 for (j = 0; j < closureSize; j++) { 767 anchors[aOff + j] = closure[2*j]; 768 } 769 PetscCall(DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure)); 770 } 771 } 772 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,size,anchors,PETSC_OWN_POINTER,&aIS)); 773 { 774 PetscSection aSecNew = aSec; 775 IS aISNew = aIS; 776 777 PetscCall(PetscObjectReference((PetscObject)aSec)); 778 PetscCall(PetscObjectReference((PetscObject)aIS)); 779 while (aSecNew) { 780 PetscCall(PetscSectionDestroy(&aSec)); 781 PetscCall(ISDestroy(&aIS)); 782 aSec = aSecNew; 783 aIS = aISNew; 784 aSecNew = NULL; 785 aISNew = NULL; 786 PetscCall(AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew)); 787 } 788 } 789 PetscCall(DMPlexSetAnchors(dm,aSec,aIS)); 790 PetscCall(PetscSectionDestroy(&aSec)); 791 PetscCall(ISDestroy(&aIS)); 792 PetscFunctionReturn(0); 793 } 794 795 static PetscErrorCode DMPlexGetTrueSupportSize(DM dm,PetscInt p,PetscInt *dof,PetscInt *numTrueSupp) 796 { 797 PetscFunctionBegin; 798 if (numTrueSupp[p] == -1) { 799 PetscInt i, alldof; 800 const PetscInt *supp; 801 PetscInt count = 0; 802 803 PetscCall(DMPlexGetSupportSize(dm,p,&alldof)); 804 PetscCall(DMPlexGetSupport(dm,p,&supp)); 805 for (i = 0; i < alldof; i++) { 806 PetscInt q = supp[i], numCones, j; 807 const PetscInt *cone; 808 809 PetscCall(DMPlexGetConeSize(dm,q,&numCones)); 810 PetscCall(DMPlexGetCone(dm,q,&cone)); 811 for (j = 0; j < numCones; j++) { 812 if (cone[j] == p) break; 813 } 814 if (j < numCones) count++; 815 } 816 numTrueSupp[p] = count; 817 } 818 *dof = numTrueSupp[p]; 819 PetscFunctionReturn(0); 820 } 821 822 static PetscErrorCode DMPlexTreeExchangeSupports(DM dm) 823 { 824 DM_Plex *mesh = (DM_Plex *)dm->data; 825 PetscSection newSupportSection; 826 PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth; 827 PetscInt *numTrueSupp; 828 PetscInt *offsets; 829 830 PetscFunctionBegin; 831 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 832 /* symmetrize the hierarchy */ 833 PetscCall(DMPlexGetDepth(dm,&depth)); 834 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)),&newSupportSection)); 835 PetscCall(DMPlexGetChart(dm,&pStart,&pEnd)); 836 PetscCall(PetscSectionSetChart(newSupportSection,pStart,pEnd)); 837 PetscCall(PetscCalloc1(pEnd,&offsets)); 838 PetscCall(PetscMalloc1(pEnd,&numTrueSupp)); 839 for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1; 840 /* if a point is in the (true) support of q, it should be in the support of 841 * parent(q) */ 842 for (d = 0; d <= depth; d++) { 843 PetscCall(DMPlexGetHeightStratum(dm,d,&pStart,&pEnd)); 844 for (p = pStart; p < pEnd; ++p) { 845 PetscInt dof, q, qdof, parent; 846 847 PetscCall(DMPlexGetTrueSupportSize(dm,p,&dof,numTrueSupp)); 848 PetscCall(PetscSectionAddDof(newSupportSection, p, dof)); 849 q = p; 850 PetscCall(DMPlexGetTreeParent(dm,q,&parent,NULL)); 851 while (parent != q && parent >= pStart && parent < pEnd) { 852 q = parent; 853 854 PetscCall(DMPlexGetTrueSupportSize(dm,q,&qdof,numTrueSupp)); 855 PetscCall(PetscSectionAddDof(newSupportSection,p,qdof)); 856 PetscCall(PetscSectionAddDof(newSupportSection,q,dof)); 857 PetscCall(DMPlexGetTreeParent(dm,q,&parent,NULL)); 858 } 859 } 860 } 861 PetscCall(PetscSectionSetUp(newSupportSection)); 862 PetscCall(PetscSectionGetStorageSize(newSupportSection,&newSize)); 863 PetscCall(PetscMalloc1(newSize,&newSupports)); 864 for (d = 0; d <= depth; d++) { 865 PetscCall(DMPlexGetHeightStratum(dm,d,&pStart,&pEnd)); 866 for (p = pStart; p < pEnd; p++) { 867 PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent; 868 869 PetscCall(PetscSectionGetDof(mesh->supportSection, p, &dof)); 870 PetscCall(PetscSectionGetOffset(mesh->supportSection, p, &off)); 871 PetscCall(PetscSectionGetDof(newSupportSection, p, &newDof)); 872 PetscCall(PetscSectionGetOffset(newSupportSection, p, &newOff)); 873 for (i = 0; i < dof; i++) { 874 PetscInt numCones, j; 875 const PetscInt *cone; 876 PetscInt q = mesh->supports[off + i]; 877 878 PetscCall(DMPlexGetConeSize(dm,q,&numCones)); 879 PetscCall(DMPlexGetCone(dm,q,&cone)); 880 for (j = 0; j < numCones; j++) { 881 if (cone[j] == p) break; 882 } 883 if (j < numCones) newSupports[newOff+offsets[p]++] = q; 884 } 885 886 q = p; 887 PetscCall(DMPlexGetTreeParent(dm,q,&parent,NULL)); 888 while (parent != q && parent >= pStart && parent < pEnd) { 889 q = parent; 890 PetscCall(PetscSectionGetDof(mesh->supportSection, q, &qdof)); 891 PetscCall(PetscSectionGetOffset(mesh->supportSection, q, &qoff)); 892 PetscCall(PetscSectionGetOffset(newSupportSection, q, &newqOff)); 893 for (i = 0; i < qdof; i++) { 894 PetscInt numCones, j; 895 const PetscInt *cone; 896 PetscInt r = mesh->supports[qoff + i]; 897 898 PetscCall(DMPlexGetConeSize(dm,r,&numCones)); 899 PetscCall(DMPlexGetCone(dm,r,&cone)); 900 for (j = 0; j < numCones; j++) { 901 if (cone[j] == q) break; 902 } 903 if (j < numCones) newSupports[newOff+offsets[p]++] = r; 904 } 905 for (i = 0; i < dof; i++) { 906 PetscInt numCones, j; 907 const PetscInt *cone; 908 PetscInt r = mesh->supports[off + i]; 909 910 PetscCall(DMPlexGetConeSize(dm,r,&numCones)); 911 PetscCall(DMPlexGetCone(dm,r,&cone)); 912 for (j = 0; j < numCones; j++) { 913 if (cone[j] == p) break; 914 } 915 if (j < numCones) newSupports[newqOff+offsets[q]++] = r; 916 } 917 PetscCall(DMPlexGetTreeParent(dm,q,&parent,NULL)); 918 } 919 } 920 } 921 PetscCall(PetscSectionDestroy(&mesh->supportSection)); 922 mesh->supportSection = newSupportSection; 923 PetscCall(PetscFree(mesh->supports)); 924 mesh->supports = newSupports; 925 PetscCall(PetscFree(offsets)); 926 PetscCall(PetscFree(numTrueSupp)); 927 928 PetscFunctionReturn(0); 929 } 930 931 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM,PetscSection,PetscSection,Mat); 932 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM,PetscSection,PetscSection,Mat); 933 934 static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports) 935 { 936 DM_Plex *mesh = (DM_Plex *)dm->data; 937 DM refTree; 938 PetscInt size; 939 940 PetscFunctionBegin; 941 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 942 PetscValidHeaderSpecific(parentSection, PETSC_SECTION_CLASSID, 2); 943 PetscCall(PetscObjectReference((PetscObject)parentSection)); 944 PetscCall(PetscSectionDestroy(&mesh->parentSection)); 945 mesh->parentSection = parentSection; 946 PetscCall(PetscSectionGetStorageSize(parentSection,&size)); 947 if (parents != mesh->parents) { 948 PetscCall(PetscFree(mesh->parents)); 949 PetscCall(PetscMalloc1(size,&mesh->parents)); 950 PetscCall(PetscArraycpy(mesh->parents, parents, size)); 951 } 952 if (childIDs != mesh->childIDs) { 953 PetscCall(PetscFree(mesh->childIDs)); 954 PetscCall(PetscMalloc1(size,&mesh->childIDs)); 955 PetscCall(PetscArraycpy(mesh->childIDs, childIDs, size)); 956 } 957 PetscCall(DMPlexGetReferenceTree(dm,&refTree)); 958 if (refTree) { 959 DMLabel canonLabel; 960 961 PetscCall(DMGetLabel(refTree,"canonical",&canonLabel)); 962 if (canonLabel) { 963 PetscInt i; 964 965 for (i = 0; i < size; i++) { 966 PetscInt canon; 967 PetscCall(DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon)); 968 if (canon >= 0) { 969 mesh->childIDs[i] = canon; 970 } 971 } 972 } 973 mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference; 974 } else { 975 mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct; 976 } 977 PetscCall(DMPlexTreeSymmetrize(dm)); 978 if (computeCanonical) { 979 PetscInt d, dim; 980 981 /* add the canonical label */ 982 PetscCall(DMGetDimension(dm,&dim)); 983 PetscCall(DMCreateLabel(dm,"canonical")); 984 for (d = 0; d <= dim; d++) { 985 PetscInt p, dStart, dEnd, canon = -1, cNumChildren; 986 const PetscInt *cChildren; 987 988 PetscCall(DMPlexGetDepthStratum(dm,d,&dStart,&dEnd)); 989 for (p = dStart; p < dEnd; p++) { 990 PetscCall(DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren)); 991 if (cNumChildren) { 992 canon = p; 993 break; 994 } 995 } 996 if (canon == -1) continue; 997 for (p = dStart; p < dEnd; p++) { 998 PetscInt numChildren, i; 999 const PetscInt *children; 1000 1001 PetscCall(DMPlexGetTreeChildren(dm,p,&numChildren,&children)); 1002 if (numChildren) { 1003 PetscCheck(numChildren == cNumChildren,PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"All parent points in a stratum should have the same number of children: %" PetscInt_FMT " != %" PetscInt_FMT, numChildren, cNumChildren); 1004 PetscCall(DMSetLabelValue(dm,"canonical",p,canon)); 1005 for (i = 0; i < numChildren; i++) { 1006 PetscCall(DMSetLabelValue(dm,"canonical",children[i],cChildren[i])); 1007 } 1008 } 1009 } 1010 } 1011 } 1012 if (exchangeSupports) PetscCall(DMPlexTreeExchangeSupports(dm)); 1013 mesh->createanchors = DMPlexCreateAnchors_Tree; 1014 /* reset anchors */ 1015 PetscCall(DMPlexSetAnchors(dm,NULL,NULL)); 1016 PetscFunctionReturn(0); 1017 } 1018 1019 /*@ 1020 DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates 1021 the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its 1022 tree root. 1023 1024 Collective on dm 1025 1026 Input Parameters: 1027 + dm - the DMPlex object 1028 . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 1029 offset indexes the parent and childID list; the reference count of parentSection is incremented 1030 . parents - a list of the point parents; copied, can be destroyed 1031 - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 1032 the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed 1033 1034 Level: intermediate 1035 1036 .seealso: `DMPlexGetTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetAnchors()`, `DMPlexGetTreeParent()`, `DMPlexGetTreeChildren()` 1037 @*/ 1038 PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[]) 1039 { 1040 PetscFunctionBegin; 1041 PetscCall(DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE)); 1042 PetscFunctionReturn(0); 1043 } 1044 1045 /*@ 1046 DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points. 1047 Collective on dm 1048 1049 Input Parameter: 1050 . dm - the DMPlex object 1051 1052 Output Parameters: 1053 + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 1054 offset indexes the parent and childID list 1055 . parents - a list of the point parents 1056 . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 1057 the child corresponds to the point in the reference tree with index childID 1058 . childSection - the inverse of the parent section 1059 - children - a list of the point children 1060 1061 Level: intermediate 1062 1063 .seealso: `DMPlexSetTree()`, `DMPlexSetReferenceTree()`, `DMPlexSetAnchors()`, `DMPlexGetTreeParent()`, `DMPlexGetTreeChildren()` 1064 @*/ 1065 PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[]) 1066 { 1067 DM_Plex *mesh = (DM_Plex *)dm->data; 1068 1069 PetscFunctionBegin; 1070 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1071 if (parentSection) *parentSection = mesh->parentSection; 1072 if (parents) *parents = mesh->parents; 1073 if (childIDs) *childIDs = mesh->childIDs; 1074 if (childSection) *childSection = mesh->childSection; 1075 if (children) *children = mesh->children; 1076 PetscFunctionReturn(0); 1077 } 1078 1079 /*@ 1080 DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the DAG) 1081 1082 Input Parameters: 1083 + dm - the DMPlex object 1084 - point - the query point 1085 1086 Output Parameters: 1087 + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent 1088 - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point 1089 does not have a parent 1090 1091 Level: intermediate 1092 1093 .seealso: `DMPlexSetTree()`, `DMPlexGetTree()`, `DMPlexGetTreeChildren()` 1094 @*/ 1095 PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID) 1096 { 1097 DM_Plex *mesh = (DM_Plex *)dm->data; 1098 PetscSection pSec; 1099 1100 PetscFunctionBegin; 1101 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1102 pSec = mesh->parentSection; 1103 if (pSec && point >= pSec->pStart && point < pSec->pEnd) { 1104 PetscInt dof; 1105 1106 PetscCall(PetscSectionGetDof (pSec, point, &dof)); 1107 if (dof) { 1108 PetscInt off; 1109 1110 PetscCall(PetscSectionGetOffset (pSec, point, &off)); 1111 if (parent) *parent = mesh->parents[off]; 1112 if (childID) *childID = mesh->childIDs[off]; 1113 PetscFunctionReturn(0); 1114 } 1115 } 1116 if (parent) { 1117 *parent = point; 1118 } 1119 if (childID) { 1120 *childID = 0; 1121 } 1122 PetscFunctionReturn(0); 1123 } 1124 1125 /*@C 1126 DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the DAG) 1127 1128 Input Parameters: 1129 + dm - the DMPlex object 1130 - point - the query point 1131 1132 Output Parameters: 1133 + numChildren - if not NULL, set to the number of children 1134 - children - if not NULL, set to a list children, or set to NULL if the point has no children 1135 1136 Level: intermediate 1137 1138 Fortran Notes: 1139 Since it returns an array, this routine is only available in Fortran 90, and you must 1140 include petsc.h90 in your code. 1141 1142 .seealso: `DMPlexSetTree()`, `DMPlexGetTree()`, `DMPlexGetTreeParent()` 1143 @*/ 1144 PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[]) 1145 { 1146 DM_Plex *mesh = (DM_Plex *)dm->data; 1147 PetscSection childSec; 1148 PetscInt dof = 0; 1149 1150 PetscFunctionBegin; 1151 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1152 childSec = mesh->childSection; 1153 if (childSec && point >= childSec->pStart && point < childSec->pEnd) { 1154 PetscCall(PetscSectionGetDof (childSec, point, &dof)); 1155 } 1156 if (numChildren) *numChildren = dof; 1157 if (children) { 1158 if (dof) { 1159 PetscInt off; 1160 1161 PetscCall(PetscSectionGetOffset (childSec, point, &off)); 1162 *children = &mesh->children[off]; 1163 } 1164 else { 1165 *children = NULL; 1166 } 1167 } 1168 PetscFunctionReturn(0); 1169 } 1170 1171 static PetscErrorCode EvaluateBasis(PetscSpace space, PetscInt nBasis, PetscInt nFunctionals, PetscInt nComps, PetscInt nPoints, const PetscInt *pointsPerFn, const PetscReal *points, const PetscReal *weights, PetscReal *work, Mat basisAtPoints) 1172 { 1173 PetscInt f, b, p, c, offset, qPoints; 1174 1175 PetscFunctionBegin; 1176 PetscCall(PetscSpaceEvaluate(space,nPoints,points,work,NULL,NULL)); 1177 for (f = 0, offset = 0; f < nFunctionals; f++) { 1178 qPoints = pointsPerFn[f]; 1179 for (b = 0; b < nBasis; b++) { 1180 PetscScalar val = 0.; 1181 1182 for (p = 0; p < qPoints; p++) { 1183 for (c = 0; c < nComps; c++) { 1184 val += work[((offset + p) * nBasis + b) * nComps + c] * weights[(offset + p) * nComps + c]; 1185 } 1186 } 1187 PetscCall(MatSetValue(basisAtPoints,b,f,val,INSERT_VALUES)); 1188 } 1189 offset += qPoints; 1190 } 1191 PetscCall(MatAssemblyBegin(basisAtPoints,MAT_FINAL_ASSEMBLY)); 1192 PetscCall(MatAssemblyEnd(basisAtPoints,MAT_FINAL_ASSEMBLY)); 1193 PetscFunctionReturn(0); 1194 } 1195 1196 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat) 1197 { 1198 PetscDS ds; 1199 PetscInt spdim; 1200 PetscInt numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd; 1201 const PetscInt *anchors; 1202 PetscSection aSec; 1203 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent; 1204 IS aIS; 1205 1206 PetscFunctionBegin; 1207 PetscCall(DMPlexGetChart(dm,&pStart,&pEnd)); 1208 PetscCall(DMGetDS(dm,&ds)); 1209 PetscCall(PetscDSGetNumFields(ds,&numFields)); 1210 PetscCall(DMPlexGetHeightStratum(dm,0,&cStart,&cEnd)); 1211 PetscCall(DMPlexGetAnchors(dm,&aSec,&aIS)); 1212 PetscCall(ISGetIndices(aIS,&anchors)); 1213 PetscCall(PetscSectionGetChart(cSec,&conStart,&conEnd)); 1214 PetscCall(DMGetDimension(dm,&spdim)); 1215 PetscCall(PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent)); 1216 1217 for (f = 0; f < numFields; f++) { 1218 PetscObject disc; 1219 PetscClassId id; 1220 PetscSpace bspace; 1221 PetscDualSpace dspace; 1222 PetscInt i, j, k, nPoints, Nc, offset; 1223 PetscInt fSize, maxDof; 1224 PetscReal *weights, *pointsRef, *pointsReal, *work; 1225 PetscScalar *scwork; 1226 const PetscScalar *X; 1227 PetscInt *sizes, *workIndRow, *workIndCol; 1228 Mat Amat, Bmat, Xmat; 1229 const PetscInt *numDof = NULL; 1230 const PetscInt ***perms = NULL; 1231 const PetscScalar ***flips = NULL; 1232 1233 PetscCall(PetscDSGetDiscretization(ds,f,&disc)); 1234 PetscCall(PetscObjectGetClassId(disc,&id)); 1235 if (id == PETSCFE_CLASSID) { 1236 PetscFE fe = (PetscFE) disc; 1237 1238 PetscCall(PetscFEGetBasisSpace(fe,&bspace)); 1239 PetscCall(PetscFEGetDualSpace(fe,&dspace)); 1240 PetscCall(PetscDualSpaceGetDimension(dspace,&fSize)); 1241 PetscCall(PetscFEGetNumComponents(fe,&Nc)); 1242 } 1243 else if (id == PETSCFV_CLASSID) { 1244 PetscFV fv = (PetscFV) disc; 1245 1246 PetscCall(PetscFVGetNumComponents(fv,&Nc)); 1247 PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)fv),&bspace)); 1248 PetscCall(PetscSpaceSetType(bspace,PETSCSPACEPOLYNOMIAL)); 1249 PetscCall(PetscSpaceSetDegree(bspace,0,PETSC_DETERMINE)); 1250 PetscCall(PetscSpaceSetNumComponents(bspace,Nc)); 1251 PetscCall(PetscSpaceSetNumVariables(bspace,spdim)); 1252 PetscCall(PetscSpaceSetUp(bspace)); 1253 PetscCall(PetscFVGetDualSpace(fv,&dspace)); 1254 PetscCall(PetscDualSpaceGetDimension(dspace,&fSize)); 1255 } 1256 else SETERRQ(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id); 1257 PetscCall(PetscDualSpaceGetNumDof(dspace,&numDof)); 1258 for (i = 0, maxDof = 0; i <= spdim; i++) {maxDof = PetscMax(maxDof,numDof[i]);} 1259 PetscCall(PetscDualSpaceGetSymmetries(dspace,&perms,&flips)); 1260 1261 PetscCall(MatCreate(PETSC_COMM_SELF,&Amat)); 1262 PetscCall(MatSetSizes(Amat,fSize,fSize,fSize,fSize)); 1263 PetscCall(MatSetType(Amat,MATSEQDENSE)); 1264 PetscCall(MatSetUp(Amat)); 1265 PetscCall(MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat)); 1266 PetscCall(MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat)); 1267 nPoints = 0; 1268 for (i = 0; i < fSize; i++) { 1269 PetscInt qPoints, thisNc; 1270 PetscQuadrature quad; 1271 1272 PetscCall(PetscDualSpaceGetFunctional(dspace,i,&quad)); 1273 PetscCall(PetscQuadratureGetData(quad,NULL,&thisNc,&qPoints,NULL,NULL)); 1274 PetscCheck(thisNc == Nc,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %" PetscInt_FMT " does not much basis dim %" PetscInt_FMT,thisNc,Nc); 1275 nPoints += qPoints; 1276 } 1277 PetscCall(PetscMalloc7(fSize,&sizes,nPoints*Nc,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal,nPoints*fSize*Nc,&work,maxDof,&workIndRow,maxDof,&workIndCol)); 1278 PetscCall(PetscMalloc1(maxDof * maxDof,&scwork)); 1279 offset = 0; 1280 for (i = 0; i < fSize; i++) { 1281 PetscInt qPoints; 1282 const PetscReal *p, *w; 1283 PetscQuadrature quad; 1284 1285 PetscCall(PetscDualSpaceGetFunctional(dspace,i,&quad)); 1286 PetscCall(PetscQuadratureGetData(quad,NULL,NULL,&qPoints,&p,&w)); 1287 PetscCall(PetscArraycpy(weights+Nc*offset,w,Nc*qPoints)); 1288 PetscCall(PetscArraycpy(pointsRef+spdim*offset,p,spdim*qPoints)); 1289 sizes[i] = qPoints; 1290 offset += qPoints; 1291 } 1292 PetscCall(EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsRef,weights,work,Amat)); 1293 PetscCall(MatLUFactor(Amat,NULL,NULL,NULL)); 1294 for (c = cStart; c < cEnd; c++) { 1295 PetscInt parent; 1296 PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL; 1297 PetscInt *childOffsets, *parentOffsets; 1298 1299 PetscCall(DMPlexGetTreeParent(dm,c,&parent,NULL)); 1300 if (parent == c) continue; 1301 PetscCall(DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure)); 1302 for (i = 0; i < closureSize; i++) { 1303 PetscInt p = closure[2*i]; 1304 PetscInt conDof; 1305 1306 if (p < conStart || p >= conEnd) continue; 1307 if (numFields) { 1308 PetscCall(PetscSectionGetFieldDof(cSec,p,f,&conDof)); 1309 } 1310 else { 1311 PetscCall(PetscSectionGetDof(cSec,p,&conDof)); 1312 } 1313 if (conDof) break; 1314 } 1315 if (i == closureSize) { 1316 PetscCall(DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure)); 1317 continue; 1318 } 1319 1320 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ)); 1321 PetscCall(DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent)); 1322 for (i = 0; i < nPoints; i++) { 1323 const PetscReal xi0[3] = {-1.,-1.,-1.}; 1324 1325 CoordinatesRefToReal(spdim, spdim, xi0, v0, J, &pointsRef[i*spdim],vtmp); 1326 CoordinatesRealToRef(spdim, spdim, xi0, v0parent, invJparent, vtmp, &pointsReal[i*spdim]); 1327 } 1328 PetscCall(EvaluateBasis(bspace,fSize,fSize,Nc,nPoints,sizes,pointsReal,weights,work,Bmat)); 1329 PetscCall(MatMatSolve(Amat,Bmat,Xmat)); 1330 PetscCall(MatDenseGetArrayRead(Xmat,&X)); 1331 PetscCall(DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP)); 1332 PetscCall(PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets)); 1333 childOffsets[0] = 0; 1334 for (i = 0; i < closureSize; i++) { 1335 PetscInt p = closure[2*i]; 1336 PetscInt dof; 1337 1338 if (numFields) { 1339 PetscCall(PetscSectionGetFieldDof(section,p,f,&dof)); 1340 } 1341 else { 1342 PetscCall(PetscSectionGetDof(section,p,&dof)); 1343 } 1344 childOffsets[i+1]=childOffsets[i]+dof; 1345 } 1346 parentOffsets[0] = 0; 1347 for (i = 0; i < closureSizeP; i++) { 1348 PetscInt p = closureP[2*i]; 1349 PetscInt dof; 1350 1351 if (numFields) { 1352 PetscCall(PetscSectionGetFieldDof(section,p,f,&dof)); 1353 } 1354 else { 1355 PetscCall(PetscSectionGetDof(section,p,&dof)); 1356 } 1357 parentOffsets[i+1]=parentOffsets[i]+dof; 1358 } 1359 for (i = 0; i < closureSize; i++) { 1360 PetscInt conDof, conOff, aDof, aOff, nWork; 1361 PetscInt p = closure[2*i]; 1362 PetscInt o = closure[2*i+1]; 1363 const PetscInt *perm; 1364 const PetscScalar *flip; 1365 1366 if (p < conStart || p >= conEnd) continue; 1367 if (numFields) { 1368 PetscCall(PetscSectionGetFieldDof(cSec,p,f,&conDof)); 1369 PetscCall(PetscSectionGetFieldOffset(cSec,p,f,&conOff)); 1370 } 1371 else { 1372 PetscCall(PetscSectionGetDof(cSec,p,&conDof)); 1373 PetscCall(PetscSectionGetOffset(cSec,p,&conOff)); 1374 } 1375 if (!conDof) continue; 1376 perm = (perms && perms[i]) ? perms[i][o] : NULL; 1377 flip = (flips && flips[i]) ? flips[i][o] : NULL; 1378 PetscCall(PetscSectionGetDof(aSec,p,&aDof)); 1379 PetscCall(PetscSectionGetOffset(aSec,p,&aOff)); 1380 nWork = childOffsets[i+1]-childOffsets[i]; 1381 for (k = 0; k < aDof; k++) { 1382 PetscInt a = anchors[aOff + k]; 1383 PetscInt aSecDof, aSecOff; 1384 1385 if (numFields) { 1386 PetscCall(PetscSectionGetFieldDof(section,a,f,&aSecDof)); 1387 PetscCall(PetscSectionGetFieldOffset(section,a,f,&aSecOff)); 1388 } 1389 else { 1390 PetscCall(PetscSectionGetDof(section,a,&aSecDof)); 1391 PetscCall(PetscSectionGetOffset(section,a,&aSecOff)); 1392 } 1393 if (!aSecDof) continue; 1394 1395 for (j = 0; j < closureSizeP; j++) { 1396 PetscInt q = closureP[2*j]; 1397 PetscInt oq = closureP[2*j+1]; 1398 1399 if (q == a) { 1400 PetscInt r, s, nWorkP; 1401 const PetscInt *permP; 1402 const PetscScalar *flipP; 1403 1404 permP = (perms && perms[j]) ? perms[j][oq] : NULL; 1405 flipP = (flips && flips[j]) ? flips[j][oq] : NULL; 1406 nWorkP = parentOffsets[j+1]-parentOffsets[j]; 1407 /* get a copy of the child-to-anchor portion of the matrix, and transpose so that rows correspond to the 1408 * child and columns correspond to the anchor: BUT the maxrix returned by MatDenseGetArrayRead() is 1409 * column-major, so transpose-transpose = do nothing */ 1410 for (r = 0; r < nWork; r++) { 1411 for (s = 0; s < nWorkP; s++) { 1412 scwork[r * nWorkP + s] = X[fSize * (r + childOffsets[i]) + (s + parentOffsets[j])]; 1413 } 1414 } 1415 for (r = 0; r < nWork; r++) {workIndRow[perm ? perm[r] : r] = conOff + r;} 1416 for (s = 0; s < nWorkP; s++) {workIndCol[permP ? permP[s] : s] = aSecOff + s;} 1417 if (flip) { 1418 for (r = 0; r < nWork; r++) { 1419 for (s = 0; s < nWorkP; s++) { 1420 scwork[r * nWorkP + s] *= flip[r]; 1421 } 1422 } 1423 } 1424 if (flipP) { 1425 for (r = 0; r < nWork; r++) { 1426 for (s = 0; s < nWorkP; s++) { 1427 scwork[r * nWorkP + s] *= flipP[s]; 1428 } 1429 } 1430 } 1431 PetscCall(MatSetValues(cMat,nWork,workIndRow,nWorkP,workIndCol,scwork,INSERT_VALUES)); 1432 break; 1433 } 1434 } 1435 } 1436 } 1437 PetscCall(MatDenseRestoreArrayRead(Xmat,&X)); 1438 PetscCall(PetscFree2(childOffsets,parentOffsets)); 1439 PetscCall(DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure)); 1440 PetscCall(DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP)); 1441 } 1442 PetscCall(MatDestroy(&Amat)); 1443 PetscCall(MatDestroy(&Bmat)); 1444 PetscCall(MatDestroy(&Xmat)); 1445 PetscCall(PetscFree(scwork)); 1446 PetscCall(PetscFree7(sizes,weights,pointsRef,pointsReal,work,workIndRow,workIndCol)); 1447 if (id == PETSCFV_CLASSID) { 1448 PetscCall(PetscSpaceDestroy(&bspace)); 1449 } 1450 } 1451 PetscCall(MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY)); 1452 PetscCall(MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY)); 1453 PetscCall(PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent)); 1454 PetscCall(ISRestoreIndices(aIS,&anchors)); 1455 1456 PetscFunctionReturn(0); 1457 } 1458 1459 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN) 1460 { 1461 Mat refCmat; 1462 PetscDS ds; 1463 PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN; 1464 PetscScalar ***refPointFieldMats; 1465 PetscSection refConSec, refAnSec, refSection; 1466 IS refAnIS; 1467 const PetscInt *refAnchors; 1468 const PetscInt **perms; 1469 const PetscScalar **flips; 1470 1471 PetscFunctionBegin; 1472 PetscCall(DMGetDS(refTree,&ds)); 1473 PetscCall(PetscDSGetNumFields(ds,&numFields)); 1474 maxFields = PetscMax(1,numFields); 1475 PetscCall(DMGetDefaultConstraints(refTree,&refConSec,&refCmat,NULL)); 1476 PetscCall(DMPlexGetAnchors(refTree,&refAnSec,&refAnIS)); 1477 PetscCall(ISGetIndices(refAnIS,&refAnchors)); 1478 PetscCall(DMGetLocalSection(refTree,&refSection)); 1479 PetscCall(PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd)); 1480 PetscCall(PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats)); 1481 PetscCall(PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN)); 1482 PetscCall(PetscSectionGetMaxDof(refConSec,&maxDof)); 1483 PetscCall(PetscSectionGetMaxDof(refAnSec,&maxAnDof)); 1484 PetscCall(PetscMalloc1(maxDof,&rows)); 1485 PetscCall(PetscMalloc1(maxDof*maxAnDof,&cols)); 1486 for (p = pRefStart; p < pRefEnd; p++) { 1487 PetscInt parent, closureSize, *closure = NULL, pDof; 1488 1489 PetscCall(DMPlexGetTreeParent(refTree,p,&parent,NULL)); 1490 PetscCall(PetscSectionGetDof(refConSec,p,&pDof)); 1491 if (!pDof || parent == p) continue; 1492 1493 PetscCall(PetscMalloc1(maxFields,&refPointFieldMats[p-pRefStart])); 1494 PetscCall(PetscCalloc1(maxFields,&refPointFieldN[p-pRefStart])); 1495 PetscCall(DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure)); 1496 for (f = 0; f < maxFields; f++) { 1497 PetscInt cDof, cOff, numCols, r, i; 1498 1499 if (f < numFields) { 1500 PetscCall(PetscSectionGetFieldDof(refConSec,p,f,&cDof)); 1501 PetscCall(PetscSectionGetFieldOffset(refConSec,p,f,&cOff)); 1502 PetscCall(PetscSectionGetFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips)); 1503 } else { 1504 PetscCall(PetscSectionGetDof(refConSec,p,&cDof)); 1505 PetscCall(PetscSectionGetOffset(refConSec,p,&cOff)); 1506 PetscCall(PetscSectionGetPointSyms(refSection,closureSize,closure,&perms,&flips)); 1507 } 1508 1509 for (r = 0; r < cDof; r++) { 1510 rows[r] = cOff + r; 1511 } 1512 numCols = 0; 1513 for (i = 0; i < closureSize; i++) { 1514 PetscInt q = closure[2*i]; 1515 PetscInt aDof, aOff, j; 1516 const PetscInt *perm = perms ? perms[i] : NULL; 1517 1518 if (numFields) { 1519 PetscCall(PetscSectionGetFieldDof(refSection,q,f,&aDof)); 1520 PetscCall(PetscSectionGetFieldOffset(refSection,q,f,&aOff)); 1521 } 1522 else { 1523 PetscCall(PetscSectionGetDof(refSection,q,&aDof)); 1524 PetscCall(PetscSectionGetOffset(refSection,q,&aOff)); 1525 } 1526 1527 for (j = 0; j < aDof; j++) { 1528 cols[numCols++] = aOff + (perm ? perm[j] : j); 1529 } 1530 } 1531 refPointFieldN[p-pRefStart][f] = numCols; 1532 PetscCall(PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f])); 1533 PetscCall(MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f])); 1534 if (flips) { 1535 PetscInt colOff = 0; 1536 1537 for (i = 0; i < closureSize; i++) { 1538 PetscInt q = closure[2*i]; 1539 PetscInt aDof, aOff, j; 1540 const PetscScalar *flip = flips ? flips[i] : NULL; 1541 1542 if (numFields) { 1543 PetscCall(PetscSectionGetFieldDof(refSection,q,f,&aDof)); 1544 PetscCall(PetscSectionGetFieldOffset(refSection,q,f,&aOff)); 1545 } 1546 else { 1547 PetscCall(PetscSectionGetDof(refSection,q,&aDof)); 1548 PetscCall(PetscSectionGetOffset(refSection,q,&aOff)); 1549 } 1550 if (flip) { 1551 PetscInt k; 1552 for (k = 0; k < cDof; k++) { 1553 for (j = 0; j < aDof; j++) { 1554 refPointFieldMats[p-pRefStart][f][k * numCols + colOff + j] *= flip[j]; 1555 } 1556 } 1557 } 1558 colOff += aDof; 1559 } 1560 } 1561 if (numFields) { 1562 PetscCall(PetscSectionRestoreFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips)); 1563 } else { 1564 PetscCall(PetscSectionRestorePointSyms(refSection,closureSize,closure,&perms,&flips)); 1565 } 1566 } 1567 PetscCall(DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure)); 1568 } 1569 *childrenMats = refPointFieldMats; 1570 *childrenN = refPointFieldN; 1571 PetscCall(ISRestoreIndices(refAnIS,&refAnchors)); 1572 PetscCall(PetscFree(rows)); 1573 PetscCall(PetscFree(cols)); 1574 PetscFunctionReturn(0); 1575 } 1576 1577 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN) 1578 { 1579 PetscDS ds; 1580 PetscInt **refPointFieldN; 1581 PetscScalar ***refPointFieldMats; 1582 PetscInt numFields, maxFields, pRefStart, pRefEnd, p, f; 1583 PetscSection refConSec; 1584 1585 PetscFunctionBegin; 1586 refPointFieldN = *childrenN; 1587 *childrenN = NULL; 1588 refPointFieldMats = *childrenMats; 1589 *childrenMats = NULL; 1590 PetscCall(DMGetDS(refTree,&ds)); 1591 PetscCall(PetscDSGetNumFields(ds,&numFields)); 1592 maxFields = PetscMax(1,numFields); 1593 PetscCall(DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL)); 1594 PetscCall(PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd)); 1595 for (p = pRefStart; p < pRefEnd; p++) { 1596 PetscInt parent, pDof; 1597 1598 PetscCall(DMPlexGetTreeParent(refTree,p,&parent,NULL)); 1599 PetscCall(PetscSectionGetDof(refConSec,p,&pDof)); 1600 if (!pDof || parent == p) continue; 1601 1602 for (f = 0; f < maxFields; f++) { 1603 PetscInt cDof; 1604 1605 if (numFields) { 1606 PetscCall(PetscSectionGetFieldDof(refConSec,p,f,&cDof)); 1607 } 1608 else { 1609 PetscCall(PetscSectionGetDof(refConSec,p,&cDof)); 1610 } 1611 1612 PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f])); 1613 } 1614 PetscCall(PetscFree(refPointFieldMats[p - pRefStart])); 1615 PetscCall(PetscFree(refPointFieldN[p - pRefStart])); 1616 } 1617 PetscCall(PetscFree(refPointFieldMats)); 1618 PetscCall(PetscFree(refPointFieldN)); 1619 PetscFunctionReturn(0); 1620 } 1621 1622 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat) 1623 { 1624 DM refTree; 1625 PetscDS ds; 1626 Mat refCmat; 1627 PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN; 1628 PetscScalar ***refPointFieldMats, *pointWork; 1629 PetscSection refConSec, refAnSec, anSec; 1630 IS refAnIS, anIS; 1631 const PetscInt *anchors; 1632 1633 PetscFunctionBegin; 1634 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1635 PetscCall(DMGetDS(dm,&ds)); 1636 PetscCall(PetscDSGetNumFields(ds,&numFields)); 1637 maxFields = PetscMax(1,numFields); 1638 PetscCall(DMPlexGetReferenceTree(dm,&refTree)); 1639 PetscCall(DMCopyDisc(dm,refTree)); 1640 PetscCall(DMGetDefaultConstraints(refTree,&refConSec,&refCmat,NULL)); 1641 PetscCall(DMPlexGetAnchors(refTree,&refAnSec,&refAnIS)); 1642 PetscCall(DMPlexGetAnchors(dm,&anSec,&anIS)); 1643 PetscCall(ISGetIndices(anIS,&anchors)); 1644 PetscCall(PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd)); 1645 PetscCall(PetscSectionGetChart(conSec,&conStart,&conEnd)); 1646 PetscCall(PetscSectionGetMaxDof(refConSec,&maxDof)); 1647 PetscCall(PetscSectionGetMaxDof(refAnSec,&maxAnDof)); 1648 PetscCall(PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork)); 1649 1650 /* step 1: get submats for every constrained point in the reference tree */ 1651 PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN)); 1652 1653 /* step 2: compute the preorder */ 1654 PetscCall(DMPlexGetChart(dm,&pStart,&pEnd)); 1655 PetscCall(PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm)); 1656 for (p = pStart; p < pEnd; p++) { 1657 perm[p - pStart] = p; 1658 iperm[p - pStart] = p-pStart; 1659 } 1660 for (p = 0; p < pEnd - pStart;) { 1661 PetscInt point = perm[p]; 1662 PetscInt parent; 1663 1664 PetscCall(DMPlexGetTreeParent(dm,point,&parent,NULL)); 1665 if (parent == point) { 1666 p++; 1667 } 1668 else { 1669 PetscInt size, closureSize, *closure = NULL, i; 1670 1671 PetscCall(DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure)); 1672 for (i = 0; i < closureSize; i++) { 1673 PetscInt q = closure[2*i]; 1674 if (iperm[q-pStart] > iperm[point-pStart]) { 1675 /* swap */ 1676 perm[p] = q; 1677 perm[iperm[q-pStart]] = point; 1678 iperm[point-pStart] = iperm[q-pStart]; 1679 iperm[q-pStart] = p; 1680 break; 1681 } 1682 } 1683 size = closureSize; 1684 PetscCall(DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure)); 1685 if (i == size) { 1686 p++; 1687 } 1688 } 1689 } 1690 1691 /* step 3: fill the constraint matrix */ 1692 /* we are going to use a preorder progressive fill strategy. Mat doesn't 1693 * allow progressive fill without assembly, so we are going to set up the 1694 * values outside of the Mat first. 1695 */ 1696 { 1697 PetscInt nRows, row, nnz; 1698 PetscBool done; 1699 const PetscInt *ia, *ja; 1700 PetscScalar *vals; 1701 1702 PetscCall(MatGetRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done)); 1703 PetscCheck(done,PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not get RowIJ of constraint matrix"); 1704 nnz = ia[nRows]; 1705 /* malloc and then zero rows right before we fill them: this way valgrind 1706 * can tell if we are doing progressive fill in the wrong order */ 1707 PetscCall(PetscMalloc1(nnz,&vals)); 1708 for (p = 0; p < pEnd - pStart; p++) { 1709 PetscInt parent, childid, closureSize, *closure = NULL; 1710 PetscInt point = perm[p], pointDof; 1711 1712 PetscCall(DMPlexGetTreeParent(dm,point,&parent,&childid)); 1713 if ((point < conStart) || (point >= conEnd) || (parent == point)) continue; 1714 PetscCall(PetscSectionGetDof(conSec,point,&pointDof)); 1715 if (!pointDof) continue; 1716 PetscCall(DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure)); 1717 for (f = 0; f < maxFields; f++) { 1718 PetscInt cDof, cOff, numCols, numFillCols, i, r, matOffset, offset; 1719 PetscScalar *pointMat; 1720 const PetscInt **perms; 1721 const PetscScalar **flips; 1722 1723 if (numFields) { 1724 PetscCall(PetscSectionGetFieldDof(conSec,point,f,&cDof)); 1725 PetscCall(PetscSectionGetFieldOffset(conSec,point,f,&cOff)); 1726 } 1727 else { 1728 PetscCall(PetscSectionGetDof(conSec,point,&cDof)); 1729 PetscCall(PetscSectionGetOffset(conSec,point,&cOff)); 1730 } 1731 if (!cDof) continue; 1732 if (numFields) PetscCall(PetscSectionGetFieldPointSyms(section,f,closureSize,closure,&perms,&flips)); 1733 else PetscCall(PetscSectionGetPointSyms(section,closureSize,closure,&perms,&flips)); 1734 1735 /* make sure that every row for this point is the same size */ 1736 if (PetscDefined(USE_DEBUG)) { 1737 for (r = 0; r < cDof; r++) { 1738 if (cDof > 1 && r) { 1739 PetscCheck((ia[cOff+r+1]-ia[cOff+r]) == (ia[cOff+r]-ia[cOff+r-1]),PETSC_COMM_SELF,PETSC_ERR_PLIB,"Two point rows have different nnz: %" PetscInt_FMT " vs. %" PetscInt_FMT, (ia[cOff+r+1]-ia[cOff+r]), (ia[cOff+r]-ia[cOff+r-1])); 1740 } 1741 } 1742 } 1743 /* zero rows */ 1744 for (i = ia[cOff] ; i< ia[cOff+cDof];i++) { 1745 vals[i] = 0.; 1746 } 1747 matOffset = ia[cOff]; 1748 numFillCols = ia[cOff+1] - matOffset; 1749 pointMat = refPointFieldMats[childid-pRefStart][f]; 1750 numCols = refPointFieldN[childid-pRefStart][f]; 1751 offset = 0; 1752 for (i = 0; i < closureSize; i++) { 1753 PetscInt q = closure[2*i]; 1754 PetscInt aDof, aOff, j, k, qConDof, qConOff; 1755 const PetscInt *perm = perms ? perms[i] : NULL; 1756 const PetscScalar *flip = flips ? flips[i] : NULL; 1757 1758 qConDof = qConOff = 0; 1759 if (numFields) { 1760 PetscCall(PetscSectionGetFieldDof(section,q,f,&aDof)); 1761 PetscCall(PetscSectionGetFieldOffset(section,q,f,&aOff)); 1762 if (q >= conStart && q < conEnd) { 1763 PetscCall(PetscSectionGetFieldDof(conSec,q,f,&qConDof)); 1764 PetscCall(PetscSectionGetFieldOffset(conSec,q,f,&qConOff)); 1765 } 1766 } 1767 else { 1768 PetscCall(PetscSectionGetDof(section,q,&aDof)); 1769 PetscCall(PetscSectionGetOffset(section,q,&aOff)); 1770 if (q >= conStart && q < conEnd) { 1771 PetscCall(PetscSectionGetDof(conSec,q,&qConDof)); 1772 PetscCall(PetscSectionGetOffset(conSec,q,&qConOff)); 1773 } 1774 } 1775 if (!aDof) continue; 1776 if (qConDof) { 1777 /* this point has anchors: its rows of the matrix should already 1778 * be filled, thanks to preordering */ 1779 /* first multiply into pointWork, then set in matrix */ 1780 PetscInt aMatOffset = ia[qConOff]; 1781 PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset; 1782 for (r = 0; r < cDof; r++) { 1783 for (j = 0; j < aNumFillCols; j++) { 1784 PetscScalar inVal = 0; 1785 for (k = 0; k < aDof; k++) { 1786 PetscInt col = perm ? perm[k] : k; 1787 1788 inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j] * (flip ? flip[col] : 1.); 1789 } 1790 pointWork[r * aNumFillCols + j] = inVal; 1791 } 1792 } 1793 /* assume that the columns are sorted, spend less time searching */ 1794 for (j = 0, k = 0; j < aNumFillCols; j++) { 1795 PetscInt col = ja[aMatOffset + j]; 1796 for (;k < numFillCols; k++) { 1797 if (ja[matOffset + k] == col) { 1798 break; 1799 } 1800 } 1801 PetscCheck(k != numFillCols,PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%" PetscInt_FMT ", %" PetscInt_FMT ")", cOff, col); 1802 for (r = 0; r < cDof; r++) { 1803 vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j]; 1804 } 1805 } 1806 } 1807 else { 1808 /* find where to put this portion of pointMat into the matrix */ 1809 for (k = 0; k < numFillCols; k++) { 1810 if (ja[matOffset + k] == aOff) { 1811 break; 1812 } 1813 } 1814 PetscCheck(k != numFillCols,PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%" PetscInt_FMT ", %" PetscInt_FMT ")", cOff, aOff); 1815 for (r = 0; r < cDof; r++) { 1816 for (j = 0; j < aDof; j++) { 1817 PetscInt col = perm ? perm[j] : j; 1818 1819 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.); 1820 } 1821 } 1822 } 1823 offset += aDof; 1824 } 1825 if (numFields) { 1826 PetscCall(PetscSectionRestoreFieldPointSyms(section,f,closureSize,closure,&perms,&flips)); 1827 } else { 1828 PetscCall(PetscSectionRestorePointSyms(section,closureSize,closure,&perms,&flips)); 1829 } 1830 } 1831 PetscCall(DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure)); 1832 } 1833 for (row = 0; row < nRows; row++) { 1834 PetscCall(MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES)); 1835 } 1836 PetscCall(MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done)); 1837 PetscCheck(done,PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix"); 1838 PetscCall(MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY)); 1839 PetscCall(MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY)); 1840 PetscCall(PetscFree(vals)); 1841 } 1842 1843 /* clean up */ 1844 PetscCall(ISRestoreIndices(anIS,&anchors)); 1845 PetscCall(PetscFree2(perm,iperm)); 1846 PetscCall(PetscFree(pointWork)); 1847 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN)); 1848 PetscFunctionReturn(0); 1849 } 1850 1851 /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of 1852 * a non-conforming mesh. Local refinement comes later */ 1853 PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm) 1854 { 1855 DM K; 1856 PetscMPIInt rank; 1857 PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd; 1858 PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations; 1859 PetscInt *Kembedding; 1860 PetscInt *cellClosure=NULL, nc; 1861 PetscScalar *newVertexCoords; 1862 PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset; 1863 PetscSection parentSection; 1864 1865 PetscFunctionBegin; 1866 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank)); 1867 PetscCall(DMGetDimension(dm,&dim)); 1868 PetscCall(DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm)); 1869 PetscCall(DMSetDimension(*ncdm,dim)); 1870 1871 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 1872 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection)); 1873 PetscCall(DMPlexGetReferenceTree(dm,&K)); 1874 if (rank == 0) { 1875 /* compute the new charts */ 1876 PetscCall(PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd)); 1877 offset = 0; 1878 for (d = 0; d <= dim; d++) { 1879 PetscInt pOldCount, kStart, kEnd, k; 1880 1881 pNewStart[d] = offset; 1882 PetscCall(DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d])); 1883 PetscCall(DMPlexGetHeightStratum(K,d,&kStart,&kEnd)); 1884 pOldCount = pOldEnd[d] - pOldStart[d]; 1885 /* adding the new points */ 1886 pNewCount[d] = pOldCount + kEnd - kStart; 1887 if (!d) { 1888 /* removing the cell */ 1889 pNewCount[d]--; 1890 } 1891 for (k = kStart; k < kEnd; k++) { 1892 PetscInt parent; 1893 PetscCall(DMPlexGetTreeParent(K,k,&parent,NULL)); 1894 if (parent == k) { 1895 /* avoid double counting points that won't actually be new */ 1896 pNewCount[d]--; 1897 } 1898 } 1899 pNewEnd[d] = pNewStart[d] + pNewCount[d]; 1900 offset = pNewEnd[d]; 1901 1902 } 1903 PetscCheck(cell >= pOldStart[0] && cell < pOldEnd[0],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"%" PetscInt_FMT " not in cell range [%" PetscInt_FMT ", %" PetscInt_FMT ")", cell, pOldStart[0], pOldEnd[0]); 1904 /* get the current closure of the cell that we are removing */ 1905 PetscCall(DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure)); 1906 1907 PetscCall(PetscMalloc1(pNewEnd[dim],&newConeSizes)); 1908 { 1909 DMPolytopeType pct, qct; 1910 PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j; 1911 1912 PetscCall(DMPlexGetChart(K,&kStart,&kEnd)); 1913 PetscCall(PetscMalloc4(kEnd-kStart,&Kembedding,kEnd-kStart,&perm,kEnd-kStart,&iperm,kEnd-kStart,&preOrient)); 1914 1915 for (k = kStart; k < kEnd; k++) { 1916 perm[k - kStart] = k; 1917 iperm [k - kStart] = k - kStart; 1918 preOrient[k - kStart] = 0; 1919 } 1920 1921 PetscCall(DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK)); 1922 for (j = 1; j < closureSizeK; j++) { 1923 PetscInt parentOrientA = closureK[2*j+1]; 1924 PetscInt parentOrientB = cellClosure[2*j+1]; 1925 PetscInt p, q; 1926 1927 p = closureK[2*j]; 1928 q = cellClosure[2*j]; 1929 PetscCall(DMPlexGetCellType(K, p, &pct)); 1930 PetscCall(DMPlexGetCellType(dm, q, &qct)); 1931 for (d = 0; d <= dim; d++) { 1932 if (q >= pOldStart[d] && q < pOldEnd[d]) { 1933 Kembedding[p] = (q - pOldStart[d]) + pNewStart[d]; 1934 } 1935 } 1936 parentOrientA = DMPolytopeConvertNewOrientation_Internal(pct, parentOrientA); 1937 parentOrientB = DMPolytopeConvertNewOrientation_Internal(qct, parentOrientB); 1938 if (parentOrientA != parentOrientB) { 1939 PetscInt numChildren, i; 1940 const PetscInt *children; 1941 1942 PetscCall(DMPlexGetTreeChildren(K,p,&numChildren,&children)); 1943 for (i = 0; i < numChildren; i++) { 1944 PetscInt kPerm, oPerm; 1945 1946 k = children[i]; 1947 PetscCall(DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm)); 1948 /* perm = what refTree position I'm in */ 1949 perm[kPerm-kStart] = k; 1950 /* iperm = who is at this position */ 1951 iperm[k-kStart] = kPerm-kStart; 1952 preOrient[kPerm-kStart] = oPerm; 1953 } 1954 } 1955 } 1956 PetscCall(DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK)); 1957 } 1958 PetscCall(PetscSectionSetChart(parentSection,0,pNewEnd[dim])); 1959 offset = 0; 1960 numNewCones = 0; 1961 for (d = 0; d <= dim; d++) { 1962 PetscInt kStart, kEnd, k; 1963 PetscInt p; 1964 PetscInt size; 1965 1966 for (p = pOldStart[d]; p < pOldEnd[d]; p++) { 1967 /* skip cell 0 */ 1968 if (p == cell) continue; 1969 /* old cones to new cones */ 1970 PetscCall(DMPlexGetConeSize(dm,p,&size)); 1971 newConeSizes[offset++] = size; 1972 numNewCones += size; 1973 } 1974 1975 PetscCall(DMPlexGetHeightStratum(K,d,&kStart,&kEnd)); 1976 for (k = kStart; k < kEnd; k++) { 1977 PetscInt kParent; 1978 1979 PetscCall(DMPlexGetTreeParent(K,k,&kParent,NULL)); 1980 if (kParent != k) { 1981 Kembedding[k] = offset; 1982 PetscCall(DMPlexGetConeSize(K,k,&size)); 1983 newConeSizes[offset++] = size; 1984 numNewCones += size; 1985 if (kParent != 0) { 1986 PetscCall(PetscSectionSetDof(parentSection,Kembedding[k],1)); 1987 } 1988 } 1989 } 1990 } 1991 1992 PetscCall(PetscSectionSetUp(parentSection)); 1993 PetscCall(PetscSectionGetStorageSize(parentSection,&numPointsWithParents)); 1994 PetscCall(PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations)); 1995 PetscCall(PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs)); 1996 1997 /* fill new cones */ 1998 offset = 0; 1999 for (d = 0; d <= dim; d++) { 2000 PetscInt kStart, kEnd, k, l; 2001 PetscInt p; 2002 PetscInt size; 2003 const PetscInt *cone, *orientation; 2004 2005 for (p = pOldStart[d]; p < pOldEnd[d]; p++) { 2006 /* skip cell 0 */ 2007 if (p == cell) continue; 2008 /* old cones to new cones */ 2009 PetscCall(DMPlexGetConeSize(dm,p,&size)); 2010 PetscCall(DMPlexGetCone(dm,p,&cone)); 2011 PetscCall(DMPlexGetConeOrientation(dm,p,&orientation)); 2012 for (l = 0; l < size; l++) { 2013 newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1]; 2014 newOrientations[offset++] = orientation[l]; 2015 } 2016 } 2017 2018 PetscCall(DMPlexGetHeightStratum(K,d,&kStart,&kEnd)); 2019 for (k = kStart; k < kEnd; k++) { 2020 PetscInt kPerm = perm[k], kParent; 2021 PetscInt preO = preOrient[k]; 2022 2023 PetscCall(DMPlexGetTreeParent(K,k,&kParent,NULL)); 2024 if (kParent != k) { 2025 /* embed new cones */ 2026 PetscCall(DMPlexGetConeSize(K,k,&size)); 2027 PetscCall(DMPlexGetCone(K,kPerm,&cone)); 2028 PetscCall(DMPlexGetConeOrientation(K,kPerm,&orientation)); 2029 for (l = 0; l < size; l++) { 2030 PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size); 2031 PetscInt newO, lSize, oTrue; 2032 DMPolytopeType ct = DM_NUM_POLYTOPES; 2033 2034 q = iperm[cone[m]]; 2035 newCones[offset] = Kembedding[q]; 2036 PetscCall(DMPlexGetConeSize(K,q,&lSize)); 2037 if (lSize == 2) ct = DM_POLYTOPE_SEGMENT; 2038 else if (lSize == 4) ct = DM_POLYTOPE_QUADRILATERAL; 2039 oTrue = DMPolytopeConvertNewOrientation_Internal(ct, orientation[m]); 2040 oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2); 2041 newO = DihedralCompose(lSize,oTrue,preOrient[q]); 2042 newOrientations[offset++] = DMPolytopeConvertOldOrientation_Internal(ct, newO); 2043 } 2044 if (kParent != 0) { 2045 PetscInt newPoint = Kembedding[kParent]; 2046 PetscCall(PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset)); 2047 parents[pOffset] = newPoint; 2048 childIDs[pOffset] = k; 2049 } 2050 } 2051 } 2052 } 2053 2054 PetscCall(PetscMalloc1(dim*(pNewEnd[dim]-pNewStart[dim]),&newVertexCoords)); 2055 2056 /* fill coordinates */ 2057 offset = 0; 2058 { 2059 PetscInt kStart, kEnd, l; 2060 PetscSection vSection; 2061 PetscInt v; 2062 Vec coords; 2063 PetscScalar *coordvals; 2064 PetscInt dof, off; 2065 PetscReal v0[3], J[9], detJ; 2066 2067 if (PetscDefined(USE_DEBUG)) { 2068 PetscInt k; 2069 PetscCall(DMPlexGetHeightStratum(K,0,&kStart,&kEnd)); 2070 for (k = kStart; k < kEnd; k++) { 2071 PetscCall(DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ)); 2072 PetscCheck(detJ > 0.,PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %" PetscInt_FMT " has bad determinant",k); 2073 } 2074 } 2075 PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ)); 2076 PetscCall(DMGetCoordinateSection(dm,&vSection)); 2077 PetscCall(DMGetCoordinatesLocal(dm,&coords)); 2078 PetscCall(VecGetArray(coords,&coordvals)); 2079 for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) { 2080 2081 PetscCall(PetscSectionGetDof(vSection,v,&dof)); 2082 PetscCall(PetscSectionGetOffset(vSection,v,&off)); 2083 for (l = 0; l < dof; l++) { 2084 newVertexCoords[offset++] = coordvals[off + l]; 2085 } 2086 } 2087 PetscCall(VecRestoreArray(coords,&coordvals)); 2088 2089 PetscCall(DMGetCoordinateSection(K,&vSection)); 2090 PetscCall(DMGetCoordinatesLocal(K,&coords)); 2091 PetscCall(VecGetArray(coords,&coordvals)); 2092 PetscCall(DMPlexGetDepthStratum(K,0,&kStart,&kEnd)); 2093 for (v = kStart; v < kEnd; v++) { 2094 PetscReal coord[3], newCoord[3]; 2095 PetscInt vPerm = perm[v]; 2096 PetscInt kParent; 2097 const PetscReal xi0[3] = {-1.,-1.,-1.}; 2098 2099 PetscCall(DMPlexGetTreeParent(K,v,&kParent,NULL)); 2100 if (kParent != v) { 2101 /* this is a new vertex */ 2102 PetscCall(PetscSectionGetOffset(vSection,vPerm,&off)); 2103 for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]); 2104 CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord); 2105 for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l]; 2106 offset += dim; 2107 } 2108 } 2109 PetscCall(VecRestoreArray(coords,&coordvals)); 2110 } 2111 2112 /* need to reverse the order of pNewCount: vertices first, cells last */ 2113 for (d = 0; d < (dim + 1) / 2; d++) { 2114 PetscInt tmp; 2115 2116 tmp = pNewCount[d]; 2117 pNewCount[d] = pNewCount[dim - d]; 2118 pNewCount[dim - d] = tmp; 2119 } 2120 2121 PetscCall(DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords)); 2122 PetscCall(DMPlexSetReferenceTree(*ncdm,K)); 2123 PetscCall(DMPlexSetTree(*ncdm,parentSection,parents,childIDs)); 2124 2125 /* clean up */ 2126 PetscCall(DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure)); 2127 PetscCall(PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd)); 2128 PetscCall(PetscFree(newConeSizes)); 2129 PetscCall(PetscFree2(newCones,newOrientations)); 2130 PetscCall(PetscFree(newVertexCoords)); 2131 PetscCall(PetscFree2(parents,childIDs)); 2132 PetscCall(PetscFree4(Kembedding,perm,iperm,preOrient)); 2133 } 2134 else { 2135 PetscInt p, counts[4]; 2136 PetscInt *coneSizes, *cones, *orientations; 2137 Vec coordVec; 2138 PetscScalar *coords; 2139 2140 for (d = 0; d <= dim; d++) { 2141 PetscInt dStart, dEnd; 2142 2143 PetscCall(DMPlexGetDepthStratum(dm,d,&dStart,&dEnd)); 2144 counts[d] = dEnd - dStart; 2145 } 2146 PetscCall(PetscMalloc1(pEnd-pStart,&coneSizes)); 2147 for (p = pStart; p < pEnd; p++) { 2148 PetscCall(DMPlexGetConeSize(dm,p,&coneSizes[p-pStart])); 2149 } 2150 PetscCall(DMPlexGetCones(dm, &cones)); 2151 PetscCall(DMPlexGetConeOrientations(dm, &orientations)); 2152 PetscCall(DMGetCoordinatesLocal(dm,&coordVec)); 2153 PetscCall(VecGetArray(coordVec,&coords)); 2154 2155 PetscCall(PetscSectionSetChart(parentSection,pStart,pEnd)); 2156 PetscCall(PetscSectionSetUp(parentSection)); 2157 PetscCall(DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL)); 2158 PetscCall(DMPlexSetReferenceTree(*ncdm,K)); 2159 PetscCall(DMPlexSetTree(*ncdm,parentSection,NULL,NULL)); 2160 PetscCall(VecRestoreArray(coordVec,&coords)); 2161 } 2162 PetscCall(PetscSectionDestroy(&parentSection)); 2163 2164 PetscFunctionReturn(0); 2165 } 2166 2167 PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 2168 { 2169 PetscSF coarseToFineEmbedded; 2170 PetscSection globalCoarse, globalFine; 2171 PetscSection localCoarse, localFine; 2172 PetscSection aSec, cSec; 2173 PetscSection rootIndicesSec, rootMatricesSec; 2174 PetscSection leafIndicesSec, leafMatricesSec; 2175 PetscInt *rootIndices, *leafIndices; 2176 PetscScalar *rootMatrices, *leafMatrices; 2177 IS aIS; 2178 const PetscInt *anchors; 2179 Mat cMat; 2180 PetscInt numFields, maxFields; 2181 PetscInt pStartC, pEndC, pStartF, pEndF, p; 2182 PetscInt aStart, aEnd, cStart, cEnd; 2183 PetscInt *maxChildIds; 2184 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 2185 const PetscInt ***perms; 2186 const PetscScalar ***flips; 2187 2188 PetscFunctionBegin; 2189 PetscCall(DMPlexGetChart(coarse,&pStartC,&pEndC)); 2190 PetscCall(DMPlexGetChart(fine,&pStartF,&pEndF)); 2191 PetscCall(DMGetGlobalSection(fine,&globalFine)); 2192 { /* winnow fine points that don't have global dofs out of the sf */ 2193 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l; 2194 const PetscInt *leaves; 2195 2196 PetscCall(PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL)); 2197 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 2198 p = leaves ? leaves[l] : l; 2199 PetscCall(PetscSectionGetDof(globalFine,p,&dof)); 2200 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&cdof)); 2201 if ((dof - cdof) > 0) { 2202 numPointsWithDofs++; 2203 } 2204 } 2205 PetscCall(PetscMalloc1(numPointsWithDofs,&pointsWithDofs)); 2206 for (l = 0, offset = 0; l < nleaves; l++) { 2207 p = leaves ? leaves[l] : l; 2208 PetscCall(PetscSectionGetDof(globalFine,p,&dof)); 2209 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&cdof)); 2210 if ((dof - cdof) > 0) { 2211 pointsWithDofs[offset++] = l; 2212 } 2213 } 2214 PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded)); 2215 PetscCall(PetscFree(pointsWithDofs)); 2216 } 2217 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 2218 PetscCall(PetscMalloc1(pEndC-pStartC,&maxChildIds)); 2219 for (p = pStartC; p < pEndC; p++) { 2220 maxChildIds[p - pStartC] = -2; 2221 } 2222 PetscCall(PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX)); 2223 PetscCall(PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX)); 2224 2225 PetscCall(DMGetLocalSection(coarse,&localCoarse)); 2226 PetscCall(DMGetGlobalSection(coarse,&globalCoarse)); 2227 2228 PetscCall(DMPlexGetAnchors(coarse,&aSec,&aIS)); 2229 PetscCall(ISGetIndices(aIS,&anchors)); 2230 PetscCall(PetscSectionGetChart(aSec,&aStart,&aEnd)); 2231 2232 PetscCall(DMGetDefaultConstraints(coarse,&cSec,&cMat,NULL)); 2233 PetscCall(PetscSectionGetChart(cSec,&cStart,&cEnd)); 2234 2235 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 2236 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec)); 2237 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec)); 2238 PetscCall(PetscSectionSetChart(rootIndicesSec,pStartC,pEndC)); 2239 PetscCall(PetscSectionSetChart(rootMatricesSec,pStartC,pEndC)); 2240 PetscCall(PetscSectionGetNumFields(localCoarse,&numFields)); 2241 maxFields = PetscMax(1,numFields); 2242 PetscCall(PetscMalloc7(maxFields+1,&offsets,maxFields+1,&offsetsCopy,maxFields+1,&newOffsets,maxFields+1,&newOffsetsCopy,maxFields+1,&rowOffsets,maxFields+1,&numD,maxFields+1,&numO)); 2243 PetscCall(PetscMalloc2(maxFields+1,(PetscInt****)&perms,maxFields+1,(PetscScalar****)&flips)); 2244 PetscCall(PetscMemzero((void *) perms, (maxFields+1) * sizeof(const PetscInt **))); 2245 PetscCall(PetscMemzero((void *) flips, (maxFields+1) * sizeof(const PetscScalar **))); 2246 2247 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 2248 PetscInt dof, matSize = 0; 2249 PetscInt aDof = 0; 2250 PetscInt cDof = 0; 2251 PetscInt maxChildId = maxChildIds[p - pStartC]; 2252 PetscInt numRowIndices = 0; 2253 PetscInt numColIndices = 0; 2254 PetscInt f; 2255 2256 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 2257 if (dof < 0) { 2258 dof = -(dof + 1); 2259 } 2260 if (p >= aStart && p < aEnd) { 2261 PetscCall(PetscSectionGetDof(aSec,p,&aDof)); 2262 } 2263 if (p >= cStart && p < cEnd) { 2264 PetscCall(PetscSectionGetDof(cSec,p,&cDof)); 2265 } 2266 for (f = 0; f <= numFields; f++) offsets[f] = 0; 2267 for (f = 0; f <= numFields; f++) newOffsets[f] = 0; 2268 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 2269 PetscInt *closure = NULL, closureSize, cl; 2270 2271 PetscCall(DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure)); 2272 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 2273 PetscInt c = closure[2 * cl], clDof; 2274 2275 PetscCall(PetscSectionGetDof(localCoarse,c,&clDof)); 2276 numRowIndices += clDof; 2277 for (f = 0; f < numFields; f++) { 2278 PetscCall(PetscSectionGetFieldDof(localCoarse,c,f,&clDof)); 2279 offsets[f + 1] += clDof; 2280 } 2281 } 2282 for (f = 0; f < numFields; f++) { 2283 offsets[f + 1] += offsets[f]; 2284 newOffsets[f + 1] = offsets[f + 1]; 2285 } 2286 /* get the number of indices needed and their field offsets */ 2287 PetscCall(DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE)); 2288 PetscCall(DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure)); 2289 if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */ 2290 numColIndices = numRowIndices; 2291 matSize = 0; 2292 } 2293 else if (numFields) { /* we send one submat for each field: sum their sizes */ 2294 matSize = 0; 2295 for (f = 0; f < numFields; f++) { 2296 PetscInt numRow, numCol; 2297 2298 numRow = offsets[f + 1] - offsets[f]; 2299 numCol = newOffsets[f + 1] - newOffsets[f]; 2300 matSize += numRow * numCol; 2301 } 2302 } 2303 else { 2304 matSize = numRowIndices * numColIndices; 2305 } 2306 } else if (maxChildId == -1) { 2307 if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */ 2308 PetscInt aOff, a; 2309 2310 PetscCall(PetscSectionGetOffset(aSec,p,&aOff)); 2311 for (f = 0; f < numFields; f++) { 2312 PetscInt fDof; 2313 2314 PetscCall(PetscSectionGetFieldDof(localCoarse,p,f,&fDof)); 2315 offsets[f+1] = fDof; 2316 } 2317 for (a = 0; a < aDof; a++) { 2318 PetscInt anchor = anchors[a + aOff], aLocalDof; 2319 2320 PetscCall(PetscSectionGetDof(localCoarse,anchor,&aLocalDof)); 2321 numColIndices += aLocalDof; 2322 for (f = 0; f < numFields; f++) { 2323 PetscInt fDof; 2324 2325 PetscCall(PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof)); 2326 newOffsets[f+1] += fDof; 2327 } 2328 } 2329 if (numFields) { 2330 matSize = 0; 2331 for (f = 0; f < numFields; f++) { 2332 matSize += offsets[f+1] * newOffsets[f+1]; 2333 } 2334 } 2335 else { 2336 matSize = numColIndices * dof; 2337 } 2338 } 2339 else { /* no children, and no constraints on dofs: just get the global indices */ 2340 numColIndices = dof; 2341 matSize = 0; 2342 } 2343 } 2344 /* we will pack the column indices with the field offsets */ 2345 PetscCall(PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0)); 2346 PetscCall(PetscSectionSetDof(rootMatricesSec,p,matSize)); 2347 } 2348 PetscCall(PetscSectionSetUp(rootIndicesSec)); 2349 PetscCall(PetscSectionSetUp(rootMatricesSec)); 2350 { 2351 PetscInt numRootIndices, numRootMatrices; 2352 2353 PetscCall(PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices)); 2354 PetscCall(PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices)); 2355 PetscCall(PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices)); 2356 for (p = pStartC; p < pEndC; p++) { 2357 PetscInt numRowIndices, numColIndices, matSize, dof; 2358 PetscInt pIndOff, pMatOff, f; 2359 PetscInt *pInd; 2360 PetscInt maxChildId = maxChildIds[p - pStartC]; 2361 PetscScalar *pMat = NULL; 2362 2363 PetscCall(PetscSectionGetDof(rootIndicesSec,p,&numColIndices)); 2364 if (!numColIndices) { 2365 continue; 2366 } 2367 for (f = 0; f <= numFields; f++) { 2368 offsets[f] = 0; 2369 newOffsets[f] = 0; 2370 offsetsCopy[f] = 0; 2371 newOffsetsCopy[f] = 0; 2372 } 2373 numColIndices -= 2 * numFields; 2374 PetscCall(PetscSectionGetOffset(rootIndicesSec,p,&pIndOff)); 2375 pInd = &(rootIndices[pIndOff]); 2376 PetscCall(PetscSectionGetDof(rootMatricesSec,p,&matSize)); 2377 if (matSize) { 2378 PetscCall(PetscSectionGetOffset(rootMatricesSec,p,&pMatOff)); 2379 pMat = &rootMatrices[pMatOff]; 2380 } 2381 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 2382 if (dof < 0) { 2383 dof = -(dof + 1); 2384 } 2385 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 2386 PetscInt i, j; 2387 PetscInt numRowIndices = matSize / numColIndices; 2388 2389 if (!numRowIndices) { /* don't need to calculate the mat, just the indices */ 2390 PetscInt numIndices, *indices; 2391 PetscCall(DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,PETSC_TRUE,&numIndices,&indices,offsets,NULL)); 2392 PetscCheck(numIndices == numColIndices,PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations"); 2393 for (i = 0; i < numColIndices; i++) { 2394 pInd[i] = indices[i]; 2395 } 2396 for (i = 0; i < numFields; i++) { 2397 pInd[numColIndices + i] = offsets[i+1]; 2398 pInd[numColIndices + numFields + i] = offsets[i+1]; 2399 } 2400 PetscCall(DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,PETSC_TRUE,&numIndices,&indices,offsets,NULL)); 2401 } 2402 else { 2403 PetscInt closureSize, *closure = NULL, cl; 2404 PetscScalar *pMatIn, *pMatModified; 2405 PetscInt numPoints,*points; 2406 2407 PetscCall(DMGetWorkArray(coarse,numRowIndices * numRowIndices,MPIU_SCALAR,&pMatIn)); 2408 for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */ 2409 for (j = 0; j < numRowIndices; j++) { 2410 pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.; 2411 } 2412 } 2413 PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure)); 2414 for (f = 0; f < maxFields; f++) { 2415 if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f])); 2416 else PetscCall(PetscSectionGetPointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f])); 2417 } 2418 if (numFields) { 2419 for (cl = 0; cl < closureSize; cl++) { 2420 PetscInt c = closure[2 * cl]; 2421 2422 for (f = 0; f < numFields; f++) { 2423 PetscInt fDof; 2424 2425 PetscCall(PetscSectionGetFieldDof(localCoarse,c,f,&fDof)); 2426 offsets[f + 1] += fDof; 2427 } 2428 } 2429 for (f = 0; f < numFields; f++) { 2430 offsets[f + 1] += offsets[f]; 2431 newOffsets[f + 1] = offsets[f + 1]; 2432 } 2433 } 2434 /* TODO : flips here ? */ 2435 /* apply hanging node constraints on the right, get the new points and the new offsets */ 2436 PetscCall(DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,perms,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE)); 2437 for (f = 0; f < maxFields; f++) { 2438 if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f])); 2439 else PetscCall(PetscSectionRestorePointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f])); 2440 } 2441 for (f = 0; f < maxFields; f++) { 2442 if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f])); 2443 else PetscCall(PetscSectionGetPointSyms(localCoarse,numPoints,points,&perms[f],&flips[f])); 2444 } 2445 if (!numFields) { 2446 for (i = 0; i < numRowIndices * numColIndices; i++) { 2447 pMat[i] = pMatModified[i]; 2448 } 2449 } 2450 else { 2451 PetscInt i, j, count; 2452 for (f = 0, count = 0; f < numFields; f++) { 2453 for (i = offsets[f]; i < offsets[f+1]; i++) { 2454 for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) { 2455 pMat[count] = pMatModified[i * numColIndices + j]; 2456 } 2457 } 2458 } 2459 } 2460 PetscCall(DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatModified)); 2461 PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure)); 2462 PetscCall(DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatIn)); 2463 if (numFields) { 2464 for (f = 0; f < numFields; f++) { 2465 pInd[numColIndices + f] = offsets[f+1]; 2466 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2467 } 2468 for (cl = 0; cl < numPoints; cl++) { 2469 PetscInt globalOff, c = points[2*cl]; 2470 PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff)); 2471 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, NULL, pInd)); 2472 } 2473 } else { 2474 for (cl = 0; cl < numPoints; cl++) { 2475 PetscInt c = points[2*cl], globalOff; 2476 const PetscInt *perm = perms[0] ? perms[0][cl] : NULL; 2477 2478 PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff)); 2479 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perm, NULL, pInd)); 2480 } 2481 } 2482 for (f = 0; f < maxFields; f++) { 2483 if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f])); 2484 else PetscCall(PetscSectionRestorePointSyms(localCoarse,numPoints,points,&perms[f],&flips[f])); 2485 } 2486 PetscCall(DMRestoreWorkArray(coarse,numPoints,MPIU_SCALAR,&points)); 2487 } 2488 } 2489 else if (matSize) { 2490 PetscInt cOff; 2491 PetscInt *rowIndices, *colIndices, a, aDof, aOff; 2492 2493 numRowIndices = matSize / numColIndices; 2494 PetscCheck(numRowIndices == dof,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Miscounted dofs"); 2495 PetscCall(DMGetWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices)); 2496 PetscCall(DMGetWorkArray(coarse,numColIndices,MPIU_INT,&colIndices)); 2497 PetscCall(PetscSectionGetOffset(cSec,p,&cOff)); 2498 PetscCall(PetscSectionGetDof(aSec,p,&aDof)); 2499 PetscCall(PetscSectionGetOffset(aSec,p,&aOff)); 2500 if (numFields) { 2501 for (f = 0; f < numFields; f++) { 2502 PetscInt fDof; 2503 2504 PetscCall(PetscSectionGetFieldDof(cSec,p,f,&fDof)); 2505 offsets[f + 1] = fDof; 2506 for (a = 0; a < aDof; a++) { 2507 PetscInt anchor = anchors[a + aOff]; 2508 PetscCall(PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof)); 2509 newOffsets[f + 1] += fDof; 2510 } 2511 } 2512 for (f = 0; f < numFields; f++) { 2513 offsets[f + 1] += offsets[f]; 2514 offsetsCopy[f + 1] = offsets[f + 1]; 2515 newOffsets[f + 1] += newOffsets[f]; 2516 newOffsetsCopy[f + 1] = newOffsets[f + 1]; 2517 } 2518 PetscCall(DMPlexGetIndicesPointFields_Internal(cSec,PETSC_TRUE,p,cOff,offsetsCopy,PETSC_TRUE,NULL,-1, NULL,rowIndices)); 2519 for (a = 0; a < aDof; a++) { 2520 PetscInt anchor = anchors[a + aOff], lOff; 2521 PetscCall(PetscSectionGetOffset(localCoarse,anchor,&lOff)); 2522 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_TRUE,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,-1, NULL,colIndices)); 2523 } 2524 } 2525 else { 2526 PetscCall(DMPlexGetIndicesPoint_Internal(cSec,PETSC_TRUE,p,cOff,offsetsCopy,PETSC_TRUE,NULL, NULL,rowIndices)); 2527 for (a = 0; a < aDof; a++) { 2528 PetscInt anchor = anchors[a + aOff], lOff; 2529 PetscCall(PetscSectionGetOffset(localCoarse,anchor,&lOff)); 2530 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_TRUE,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL, NULL,colIndices)); 2531 } 2532 } 2533 if (numFields) { 2534 PetscInt count, a; 2535 2536 for (f = 0, count = 0; f < numFields; f++) { 2537 PetscInt iSize = offsets[f + 1] - offsets[f]; 2538 PetscInt jSize = newOffsets[f + 1] - newOffsets[f]; 2539 PetscCall(MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count])); 2540 count += iSize * jSize; 2541 pInd[numColIndices + f] = offsets[f+1]; 2542 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2543 } 2544 for (a = 0; a < aDof; a++) { 2545 PetscInt anchor = anchors[a + aOff]; 2546 PetscInt gOff; 2547 PetscCall(PetscSectionGetOffset(globalCoarse,anchor,&gOff)); 2548 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,-1, NULL,pInd)); 2549 } 2550 } 2551 else { 2552 PetscInt a; 2553 PetscCall(MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat)); 2554 for (a = 0; a < aDof; a++) { 2555 PetscInt anchor = anchors[a + aOff]; 2556 PetscInt gOff; 2557 PetscCall(PetscSectionGetOffset(globalCoarse,anchor,&gOff)); 2558 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL, NULL,pInd)); 2559 } 2560 } 2561 PetscCall(DMRestoreWorkArray(coarse,numColIndices,MPIU_INT,&colIndices)); 2562 PetscCall(DMRestoreWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices)); 2563 } 2564 else { 2565 PetscInt gOff; 2566 2567 PetscCall(PetscSectionGetOffset(globalCoarse,p,&gOff)); 2568 if (numFields) { 2569 for (f = 0; f < numFields; f++) { 2570 PetscInt fDof; 2571 PetscCall(PetscSectionGetFieldDof(localCoarse,p,f,&fDof)); 2572 offsets[f + 1] = fDof + offsets[f]; 2573 } 2574 for (f = 0; f < numFields; f++) { 2575 pInd[numColIndices + f] = offsets[f+1]; 2576 pInd[numColIndices + numFields + f] = offsets[f+1]; 2577 } 2578 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1, NULL,pInd)); 2579 } else { 2580 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL, NULL,pInd)); 2581 } 2582 } 2583 } 2584 PetscCall(PetscFree(maxChildIds)); 2585 } 2586 { 2587 PetscSF indicesSF, matricesSF; 2588 PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices; 2589 2590 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec)); 2591 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafMatricesSec)); 2592 PetscCall(PetscSFDistributeSection(coarseToFineEmbedded,rootIndicesSec,&remoteOffsetsIndices,leafIndicesSec)); 2593 PetscCall(PetscSFDistributeSection(coarseToFineEmbedded,rootMatricesSec,&remoteOffsetsMatrices,leafMatricesSec)); 2594 PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded,rootIndicesSec,remoteOffsetsIndices,leafIndicesSec,&indicesSF)); 2595 PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded,rootMatricesSec,remoteOffsetsMatrices,leafMatricesSec,&matricesSF)); 2596 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 2597 PetscCall(PetscFree(remoteOffsetsIndices)); 2598 PetscCall(PetscFree(remoteOffsetsMatrices)); 2599 PetscCall(PetscSectionGetStorageSize(leafIndicesSec,&numLeafIndices)); 2600 PetscCall(PetscSectionGetStorageSize(leafMatricesSec,&numLeafMatrices)); 2601 PetscCall(PetscMalloc2(numLeafIndices,&leafIndices,numLeafMatrices,&leafMatrices)); 2602 PetscCall(PetscSFBcastBegin(indicesSF,MPIU_INT,rootIndices,leafIndices,MPI_REPLACE)); 2603 PetscCall(PetscSFBcastBegin(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices,MPI_REPLACE)); 2604 PetscCall(PetscSFBcastEnd(indicesSF,MPIU_INT,rootIndices,leafIndices,MPI_REPLACE)); 2605 PetscCall(PetscSFBcastEnd(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices,MPI_REPLACE)); 2606 PetscCall(PetscSFDestroy(&matricesSF)); 2607 PetscCall(PetscSFDestroy(&indicesSF)); 2608 PetscCall(PetscFree2(rootIndices,rootMatrices)); 2609 PetscCall(PetscSectionDestroy(&rootIndicesSec)); 2610 PetscCall(PetscSectionDestroy(&rootMatricesSec)); 2611 } 2612 /* count to preallocate */ 2613 PetscCall(DMGetLocalSection(fine,&localFine)); 2614 { 2615 PetscInt nGlobal; 2616 PetscInt *dnnz, *onnz; 2617 PetscLayout rowMap, colMap; 2618 PetscInt rowStart, rowEnd, colStart, colEnd; 2619 PetscInt maxDof; 2620 PetscInt *rowIndices; 2621 DM refTree; 2622 PetscInt **refPointFieldN; 2623 PetscScalar ***refPointFieldMats; 2624 PetscSection refConSec, refAnSec; 2625 PetscInt pRefStart,pRefEnd,maxConDof,maxColumns,leafStart,leafEnd; 2626 PetscScalar *pointWork; 2627 2628 PetscCall(PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal)); 2629 PetscCall(PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz)); 2630 PetscCall(MatGetLayouts(mat,&rowMap,&colMap)); 2631 PetscCall(PetscLayoutSetUp(rowMap)); 2632 PetscCall(PetscLayoutSetUp(colMap)); 2633 PetscCall(PetscLayoutGetRange(rowMap,&rowStart,&rowEnd)); 2634 PetscCall(PetscLayoutGetRange(colMap,&colStart,&colEnd)); 2635 PetscCall(PetscSectionGetMaxDof(localFine,&maxDof)); 2636 PetscCall(PetscSectionGetChart(leafIndicesSec,&leafStart,&leafEnd)); 2637 PetscCall(DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices)); 2638 for (p = leafStart; p < leafEnd; p++) { 2639 PetscInt gDof, gcDof, gOff; 2640 PetscInt numColIndices, pIndOff, *pInd; 2641 PetscInt matSize; 2642 PetscInt i; 2643 2644 PetscCall(PetscSectionGetDof(globalFine,p,&gDof)); 2645 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&gcDof)); 2646 if ((gDof - gcDof) <= 0) { 2647 continue; 2648 } 2649 PetscCall(PetscSectionGetOffset(globalFine,p,&gOff)); 2650 PetscCheck(gOff >= 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"I though having global dofs meant a non-negative offset"); 2651 PetscCheck(gOff >= rowStart && (gOff + gDof - gcDof) <= rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"I thought the row map would constrain the global dofs"); 2652 PetscCall(PetscSectionGetDof(leafIndicesSec,p,&numColIndices)); 2653 PetscCall(PetscSectionGetOffset(leafIndicesSec,p,&pIndOff)); 2654 numColIndices -= 2 * numFields; 2655 PetscCheck(numColIndices > 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"global fine dof with no dofs to interpolate from"); 2656 pInd = &leafIndices[pIndOff]; 2657 offsets[0] = 0; 2658 offsetsCopy[0] = 0; 2659 newOffsets[0] = 0; 2660 newOffsetsCopy[0] = 0; 2661 if (numFields) { 2662 PetscInt f; 2663 for (f = 0; f < numFields; f++) { 2664 PetscInt rowDof; 2665 2666 PetscCall(PetscSectionGetFieldDof(localFine,p,f,&rowDof)); 2667 offsets[f + 1] = offsets[f] + rowDof; 2668 offsetsCopy[f + 1] = offsets[f + 1]; 2669 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2670 numD[f] = 0; 2671 numO[f] = 0; 2672 } 2673 PetscCall(DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,rowIndices)); 2674 for (f = 0; f < numFields; f++) { 2675 PetscInt colOffset = newOffsets[f]; 2676 PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f]; 2677 2678 for (i = 0; i < numFieldCols; i++) { 2679 PetscInt gInd = pInd[i + colOffset]; 2680 2681 if (gInd >= colStart && gInd < colEnd) { 2682 numD[f]++; 2683 } 2684 else if (gInd >= 0) { /* negative means non-entry */ 2685 numO[f]++; 2686 } 2687 } 2688 } 2689 } 2690 else { 2691 PetscCall(DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,rowIndices)); 2692 numD[0] = 0; 2693 numO[0] = 0; 2694 for (i = 0; i < numColIndices; i++) { 2695 PetscInt gInd = pInd[i]; 2696 2697 if (gInd >= colStart && gInd < colEnd) { 2698 numD[0]++; 2699 } 2700 else if (gInd >= 0) { /* negative means non-entry */ 2701 numO[0]++; 2702 } 2703 } 2704 } 2705 PetscCall(PetscSectionGetDof(leafMatricesSec,p,&matSize)); 2706 if (!matSize) { /* incoming matrix is identity */ 2707 PetscInt childId; 2708 2709 childId = childIds[p-pStartF]; 2710 if (childId < 0) { /* no child interpolation: one nnz per */ 2711 if (numFields) { 2712 PetscInt f; 2713 for (f = 0; f < numFields; f++) { 2714 PetscInt numRows = offsets[f+1] - offsets[f], row; 2715 for (row = 0; row < numRows; row++) { 2716 PetscInt gIndCoarse = pInd[newOffsets[f] + row]; 2717 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2718 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2719 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2720 dnnz[gIndFine - rowStart] = 1; 2721 } 2722 else if (gIndCoarse >= 0) { /* remote */ 2723 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2724 onnz[gIndFine - rowStart] = 1; 2725 } 2726 else { /* constrained */ 2727 PetscCheck(gIndFine < 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2728 } 2729 } 2730 } 2731 } 2732 else { 2733 PetscInt i; 2734 for (i = 0; i < gDof; i++) { 2735 PetscInt gIndCoarse = pInd[i]; 2736 PetscInt gIndFine = rowIndices[i]; 2737 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2738 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2739 dnnz[gIndFine - rowStart] = 1; 2740 } 2741 else if (gIndCoarse >= 0) { /* remote */ 2742 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2743 onnz[gIndFine - rowStart] = 1; 2744 } 2745 else { /* constrained */ 2746 PetscCheck(gIndFine < 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2747 } 2748 } 2749 } 2750 } 2751 else { /* interpolate from all */ 2752 if (numFields) { 2753 PetscInt f; 2754 for (f = 0; f < numFields; f++) { 2755 PetscInt numRows = offsets[f+1] - offsets[f], row; 2756 for (row = 0; row < numRows; row++) { 2757 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2758 if (gIndFine >= 0) { 2759 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2760 dnnz[gIndFine - rowStart] = numD[f]; 2761 onnz[gIndFine - rowStart] = numO[f]; 2762 } 2763 } 2764 } 2765 } 2766 else { 2767 PetscInt i; 2768 for (i = 0; i < gDof; i++) { 2769 PetscInt gIndFine = rowIndices[i]; 2770 if (gIndFine >= 0) { 2771 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2772 dnnz[gIndFine - rowStart] = numD[0]; 2773 onnz[gIndFine - rowStart] = numO[0]; 2774 } 2775 } 2776 } 2777 } 2778 } 2779 else { /* interpolate from all */ 2780 if (numFields) { 2781 PetscInt f; 2782 for (f = 0; f < numFields; f++) { 2783 PetscInt numRows = offsets[f+1] - offsets[f], row; 2784 for (row = 0; row < numRows; row++) { 2785 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2786 if (gIndFine >= 0) { 2787 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2788 dnnz[gIndFine - rowStart] = numD[f]; 2789 onnz[gIndFine - rowStart] = numO[f]; 2790 } 2791 } 2792 } 2793 } 2794 else { /* every dof get a full row */ 2795 PetscInt i; 2796 for (i = 0; i < gDof; i++) { 2797 PetscInt gIndFine = rowIndices[i]; 2798 if (gIndFine >= 0) { 2799 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2800 dnnz[gIndFine - rowStart] = numD[0]; 2801 onnz[gIndFine - rowStart] = numO[0]; 2802 } 2803 } 2804 } 2805 } 2806 } 2807 PetscCall(MatXAIJSetPreallocation(mat,1,dnnz,onnz,NULL,NULL)); 2808 PetscCall(PetscFree2(dnnz,onnz)); 2809 2810 PetscCall(DMPlexGetReferenceTree(fine,&refTree)); 2811 PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN)); 2812 PetscCall(DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL)); 2813 PetscCall(DMPlexGetAnchors(refTree,&refAnSec,NULL)); 2814 PetscCall(PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd)); 2815 PetscCall(PetscSectionGetMaxDof(refConSec,&maxConDof)); 2816 PetscCall(PetscSectionGetMaxDof(leafIndicesSec,&maxColumns)); 2817 PetscCall(PetscMalloc1(maxConDof*maxColumns,&pointWork)); 2818 for (p = leafStart; p < leafEnd; p++) { 2819 PetscInt gDof, gcDof, gOff; 2820 PetscInt numColIndices, pIndOff, *pInd; 2821 PetscInt matSize; 2822 PetscInt childId; 2823 2824 PetscCall(PetscSectionGetDof(globalFine,p,&gDof)); 2825 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&gcDof)); 2826 if ((gDof - gcDof) <= 0) { 2827 continue; 2828 } 2829 childId = childIds[p-pStartF]; 2830 PetscCall(PetscSectionGetOffset(globalFine,p,&gOff)); 2831 PetscCall(PetscSectionGetDof(leafIndicesSec,p,&numColIndices)); 2832 PetscCall(PetscSectionGetOffset(leafIndicesSec,p,&pIndOff)); 2833 numColIndices -= 2 * numFields; 2834 pInd = &leafIndices[pIndOff]; 2835 offsets[0] = 0; 2836 offsetsCopy[0] = 0; 2837 newOffsets[0] = 0; 2838 newOffsetsCopy[0] = 0; 2839 rowOffsets[0] = 0; 2840 if (numFields) { 2841 PetscInt f; 2842 for (f = 0; f < numFields; f++) { 2843 PetscInt rowDof; 2844 2845 PetscCall(PetscSectionGetFieldDof(localFine,p,f,&rowDof)); 2846 offsets[f + 1] = offsets[f] + rowDof; 2847 offsetsCopy[f + 1] = offsets[f + 1]; 2848 rowOffsets[f + 1] = pInd[numColIndices + f]; 2849 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2850 } 2851 PetscCall(DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,rowIndices)); 2852 } 2853 else { 2854 PetscCall(DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,rowIndices)); 2855 } 2856 PetscCall(PetscSectionGetDof(leafMatricesSec,p,&matSize)); 2857 if (!matSize) { /* incoming matrix is identity */ 2858 if (childId < 0) { /* no child interpolation: scatter */ 2859 if (numFields) { 2860 PetscInt f; 2861 for (f = 0; f < numFields; f++) { 2862 PetscInt numRows = offsets[f+1] - offsets[f], row; 2863 for (row = 0; row < numRows; row++) { 2864 PetscCall(MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES)); 2865 } 2866 } 2867 } 2868 else { 2869 PetscInt numRows = gDof, row; 2870 for (row = 0; row < numRows; row++) { 2871 PetscCall(MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES)); 2872 } 2873 } 2874 } 2875 else { /* interpolate from all */ 2876 if (numFields) { 2877 PetscInt f; 2878 for (f = 0; f < numFields; f++) { 2879 PetscInt numRows = offsets[f+1] - offsets[f]; 2880 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 2881 PetscCall(MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES)); 2882 } 2883 } 2884 else { 2885 PetscCall(MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES)); 2886 } 2887 } 2888 } 2889 else { /* interpolate from all */ 2890 PetscInt pMatOff; 2891 PetscScalar *pMat; 2892 2893 PetscCall(PetscSectionGetOffset(leafMatricesSec,p,&pMatOff)); 2894 pMat = &leafMatrices[pMatOff]; 2895 if (childId < 0) { /* copy the incoming matrix */ 2896 if (numFields) { 2897 PetscInt f, count; 2898 for (f = 0, count = 0; f < numFields; f++) { 2899 PetscInt numRows = offsets[f+1]-offsets[f]; 2900 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2901 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2902 PetscScalar *inMat = &pMat[count]; 2903 2904 PetscCall(MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES)); 2905 count += numCols * numInRows; 2906 } 2907 } 2908 else { 2909 PetscCall(MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES)); 2910 } 2911 } 2912 else { /* multiply the incoming matrix by the child interpolation */ 2913 if (numFields) { 2914 PetscInt f, count; 2915 for (f = 0, count = 0; f < numFields; f++) { 2916 PetscInt numRows = offsets[f+1]-offsets[f]; 2917 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2918 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2919 PetscScalar *inMat = &pMat[count]; 2920 PetscInt i, j, k; 2921 PetscCheck(refPointFieldN[childId - pRefStart][f] == numInRows,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2922 for (i = 0; i < numRows; i++) { 2923 for (j = 0; j < numCols; j++) { 2924 PetscScalar val = 0.; 2925 for (k = 0; k < numInRows; k++) { 2926 val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j]; 2927 } 2928 pointWork[i * numCols + j] = val; 2929 } 2930 } 2931 PetscCall(MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES)); 2932 count += numCols * numInRows; 2933 } 2934 } 2935 else { /* every dof gets a full row */ 2936 PetscInt numRows = gDof; 2937 PetscInt numCols = numColIndices; 2938 PetscInt numInRows = matSize / numColIndices; 2939 PetscInt i, j, k; 2940 PetscCheck(refPointFieldN[childId - pRefStart][0] == numInRows,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2941 for (i = 0; i < numRows; i++) { 2942 for (j = 0; j < numCols; j++) { 2943 PetscScalar val = 0.; 2944 for (k = 0; k < numInRows; k++) { 2945 val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j]; 2946 } 2947 pointWork[i * numCols + j] = val; 2948 } 2949 } 2950 PetscCall(MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES)); 2951 } 2952 } 2953 } 2954 } 2955 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN)); 2956 PetscCall(DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices)); 2957 PetscCall(PetscFree(pointWork)); 2958 } 2959 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 2960 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 2961 PetscCall(PetscSectionDestroy(&leafIndicesSec)); 2962 PetscCall(PetscSectionDestroy(&leafMatricesSec)); 2963 PetscCall(PetscFree2(leafIndices,leafMatrices)); 2964 PetscCall(PetscFree2(*(PetscInt****)&perms,*(PetscScalar****)&flips)); 2965 PetscCall(PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO)); 2966 PetscCall(ISRestoreIndices(aIS,&anchors)); 2967 PetscFunctionReturn(0); 2968 } 2969 2970 /* 2971 * Assuming a nodal basis (w.r.t. the dual basis) basis: 2972 * 2973 * for each coarse dof \phi^c_i: 2974 * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i: 2975 * for each fine dof \phi^f_j; 2976 * a_{i,j} = 0; 2977 * for each fine dof \phi^f_k: 2978 * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l 2979 * [^^^ this is = \phi^c_i ^^^] 2980 */ 2981 PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj) 2982 { 2983 PetscDS ds; 2984 PetscSection section, cSection; 2985 DMLabel canonical, depth; 2986 Mat cMat, mat; 2987 PetscInt *nnz; 2988 PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd; 2989 PetscInt m, n; 2990 PetscScalar *pointScalar; 2991 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent; 2992 2993 PetscFunctionBegin; 2994 PetscCall(DMGetLocalSection(refTree,§ion)); 2995 PetscCall(DMGetDimension(refTree, &dim)); 2996 PetscCall(PetscMalloc6(dim,&v0,dim,&v0parent,dim,&vtmp,dim*dim,&J,dim*dim,&Jparent,dim*dim,&invJ)); 2997 PetscCall(PetscMalloc2(dim,&pointScalar,dim,&pointRef)); 2998 PetscCall(DMGetDS(refTree,&ds)); 2999 PetscCall(PetscDSGetNumFields(ds,&numFields)); 3000 PetscCall(PetscSectionGetNumFields(section,&numSecFields)); 3001 PetscCall(DMGetLabel(refTree,"canonical",&canonical)); 3002 PetscCall(DMGetLabel(refTree,"depth",&depth)); 3003 PetscCall(DMGetDefaultConstraints(refTree,&cSection,&cMat,NULL)); 3004 PetscCall(DMPlexGetChart(refTree, &pStart, &pEnd)); 3005 PetscCall(DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd)); 3006 PetscCall(MatGetSize(cMat,&n,&m)); /* the injector has transpose sizes from the constraint matrix */ 3007 /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */ 3008 PetscCall(PetscCalloc1(m,&nnz)); 3009 for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */ 3010 const PetscInt *children; 3011 PetscInt numChildren; 3012 PetscInt i, numChildDof, numSelfDof; 3013 3014 if (canonical) { 3015 PetscInt pCanonical; 3016 PetscCall(DMLabelGetValue(canonical,p,&pCanonical)); 3017 if (p != pCanonical) continue; 3018 } 3019 PetscCall(DMPlexGetTreeChildren(refTree,p,&numChildren,&children)); 3020 if (!numChildren) continue; 3021 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3022 PetscInt child = children[i]; 3023 PetscInt dof; 3024 3025 PetscCall(PetscSectionGetDof(section,child,&dof)); 3026 numChildDof += dof; 3027 } 3028 PetscCall(PetscSectionGetDof(section,p,&numSelfDof)); 3029 if (!numChildDof || !numSelfDof) continue; 3030 for (f = 0; f < numFields; f++) { 3031 PetscInt selfOff; 3032 3033 if (numSecFields) { /* count the dofs for just this field */ 3034 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3035 PetscInt child = children[i]; 3036 PetscInt dof; 3037 3038 PetscCall(PetscSectionGetFieldDof(section,child,f,&dof)); 3039 numChildDof += dof; 3040 } 3041 PetscCall(PetscSectionGetFieldDof(section,p,f,&numSelfDof)); 3042 PetscCall(PetscSectionGetFieldOffset(section,p,f,&selfOff)); 3043 } 3044 else { 3045 PetscCall(PetscSectionGetOffset(section,p,&selfOff)); 3046 } 3047 for (i = 0; i < numSelfDof; i++) { 3048 nnz[selfOff + i] = numChildDof; 3049 } 3050 } 3051 } 3052 PetscCall(MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat)); 3053 PetscCall(PetscFree(nnz)); 3054 /* Setp 2: compute entries */ 3055 for (p = pStart; p < pEnd; p++) { 3056 const PetscInt *children; 3057 PetscInt numChildren; 3058 PetscInt i, numChildDof, numSelfDof; 3059 3060 /* same conditions about when entries occur */ 3061 if (canonical) { 3062 PetscInt pCanonical; 3063 PetscCall(DMLabelGetValue(canonical,p,&pCanonical)); 3064 if (p != pCanonical) continue; 3065 } 3066 PetscCall(DMPlexGetTreeChildren(refTree,p,&numChildren,&children)); 3067 if (!numChildren) continue; 3068 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3069 PetscInt child = children[i]; 3070 PetscInt dof; 3071 3072 PetscCall(PetscSectionGetDof(section,child,&dof)); 3073 numChildDof += dof; 3074 } 3075 PetscCall(PetscSectionGetDof(section,p,&numSelfDof)); 3076 if (!numChildDof || !numSelfDof) continue; 3077 3078 for (f = 0; f < numFields; f++) { 3079 PetscInt pI = -1, cI = -1; 3080 PetscInt selfOff, Nc, parentCell; 3081 PetscInt cellShapeOff; 3082 PetscObject disc; 3083 PetscDualSpace dsp; 3084 PetscClassId classId; 3085 PetscScalar *pointMat; 3086 PetscInt *matRows, *matCols; 3087 PetscInt pO = PETSC_MIN_INT; 3088 const PetscInt *depthNumDof; 3089 3090 if (numSecFields) { 3091 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3092 PetscInt child = children[i]; 3093 PetscInt dof; 3094 3095 PetscCall(PetscSectionGetFieldDof(section,child,f,&dof)); 3096 numChildDof += dof; 3097 } 3098 PetscCall(PetscSectionGetFieldDof(section,p,f,&numSelfDof)); 3099 PetscCall(PetscSectionGetFieldOffset(section,p,f,&selfOff)); 3100 } 3101 else { 3102 PetscCall(PetscSectionGetOffset(section,p,&selfOff)); 3103 } 3104 3105 /* find a cell whose closure contains p */ 3106 if (p >= cStart && p < cEnd) { 3107 parentCell = p; 3108 } 3109 else { 3110 PetscInt *star = NULL; 3111 PetscInt numStar; 3112 3113 parentCell = -1; 3114 PetscCall(DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star)); 3115 for (i = numStar - 1; i >= 0; i--) { 3116 PetscInt c = star[2 * i]; 3117 3118 if (c >= cStart && c < cEnd) { 3119 parentCell = c; 3120 break; 3121 } 3122 } 3123 PetscCall(DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star)); 3124 } 3125 /* determine the offset of p's shape functions within parentCell's shape functions */ 3126 PetscCall(PetscDSGetDiscretization(ds,f,&disc)); 3127 PetscCall(PetscObjectGetClassId(disc,&classId)); 3128 if (classId == PETSCFE_CLASSID) { 3129 PetscCall(PetscFEGetDualSpace((PetscFE)disc,&dsp)); 3130 } 3131 else if (classId == PETSCFV_CLASSID) { 3132 PetscCall(PetscFVGetDualSpace((PetscFV)disc,&dsp)); 3133 } 3134 else { 3135 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object"); 3136 } 3137 PetscCall(PetscDualSpaceGetNumDof(dsp,&depthNumDof)); 3138 PetscCall(PetscDualSpaceGetNumComponents(dsp,&Nc)); 3139 { 3140 PetscInt *closure = NULL; 3141 PetscInt numClosure; 3142 3143 PetscCall(DMPlexGetTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure)); 3144 for (i = 0, pI = -1, cellShapeOff = 0; i < numClosure; i++) { 3145 PetscInt point = closure[2 * i], pointDepth; 3146 3147 pO = closure[2 * i + 1]; 3148 if (point == p) { 3149 pI = i; 3150 break; 3151 } 3152 PetscCall(DMLabelGetValue(depth,point,&pointDepth)); 3153 cellShapeOff += depthNumDof[pointDepth]; 3154 } 3155 PetscCall(DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure)); 3156 } 3157 3158 PetscCall(DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat)); 3159 PetscCall(DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows)); 3160 matCols = matRows + numSelfDof; 3161 for (i = 0; i < numSelfDof; i++) { 3162 matRows[i] = selfOff + i; 3163 } 3164 for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.; 3165 { 3166 PetscInt colOff = 0; 3167 3168 for (i = 0; i < numChildren; i++) { 3169 PetscInt child = children[i]; 3170 PetscInt dof, off, j; 3171 3172 if (numSecFields) { 3173 PetscCall(PetscSectionGetFieldDof(cSection,child,f,&dof)); 3174 PetscCall(PetscSectionGetFieldOffset(cSection,child,f,&off)); 3175 } 3176 else { 3177 PetscCall(PetscSectionGetDof(cSection,child,&dof)); 3178 PetscCall(PetscSectionGetOffset(cSection,child,&off)); 3179 } 3180 3181 for (j = 0; j < dof; j++) { 3182 matCols[colOff++] = off + j; 3183 } 3184 } 3185 } 3186 if (classId == PETSCFE_CLASSID) { 3187 PetscFE fe = (PetscFE) disc; 3188 PetscInt fSize; 3189 const PetscInt ***perms; 3190 const PetscScalar ***flips; 3191 const PetscInt *pperms; 3192 3193 PetscCall(PetscFEGetDualSpace(fe,&dsp)); 3194 PetscCall(PetscDualSpaceGetDimension(dsp,&fSize)); 3195 PetscCall(PetscDualSpaceGetSymmetries(dsp, &perms, &flips)); 3196 pperms = perms ? perms[pI] ? perms[pI][pO] : NULL : NULL; 3197 for (i = 0; i < numSelfDof; i++) { /* for every shape function */ 3198 PetscQuadrature q; 3199 PetscInt dim, thisNc, numPoints, j, k; 3200 const PetscReal *points; 3201 const PetscReal *weights; 3202 PetscInt *closure = NULL; 3203 PetscInt numClosure; 3204 PetscInt iCell = pperms ? pperms[i] : i; 3205 PetscInt parentCellShapeDof = cellShapeOff + iCell; 3206 PetscTabulation Tparent; 3207 3208 PetscCall(PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q)); 3209 PetscCall(PetscQuadratureGetData(q,&dim,&thisNc,&numPoints,&points,&weights)); 3210 PetscCheck(thisNc == Nc,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %" PetscInt_FMT " does not much basis dim %" PetscInt_FMT,thisNc,Nc); 3211 PetscCall(PetscFECreateTabulation(fe,1,numPoints,points,0,&Tparent)); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */ 3212 for (j = 0; j < numPoints; j++) { 3213 PetscInt childCell = -1; 3214 PetscReal *parentValAtPoint; 3215 const PetscReal xi0[3] = {-1.,-1.,-1.}; 3216 const PetscReal *pointReal = &points[dim * j]; 3217 const PetscScalar *point; 3218 PetscTabulation Tchild; 3219 PetscInt childCellShapeOff, pointMatOff; 3220 #if defined(PETSC_USE_COMPLEX) 3221 PetscInt d; 3222 3223 for (d = 0; d < dim; d++) { 3224 pointScalar[d] = points[dim * j + d]; 3225 } 3226 point = pointScalar; 3227 #else 3228 point = pointReal; 3229 #endif 3230 3231 parentValAtPoint = &Tparent->T[0][(fSize * j + parentCellShapeDof) * Nc]; 3232 3233 for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/ 3234 PetscInt child = children[k]; 3235 PetscInt *star = NULL; 3236 PetscInt numStar, s; 3237 3238 PetscCall(DMPlexGetTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star)); 3239 for (s = numStar - 1; s >= 0; s--) { 3240 PetscInt c = star[2 * s]; 3241 3242 if (c < cStart || c >= cEnd) continue; 3243 PetscCall(DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell)); 3244 if (childCell >= 0) break; 3245 } 3246 PetscCall(DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star)); 3247 if (childCell >= 0) break; 3248 } 3249 PetscCheck(childCell >= 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not locate quadrature point"); 3250 PetscCall(DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ)); 3251 PetscCall(DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent)); 3252 CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp); 3253 CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef); 3254 3255 PetscCall(PetscFECreateTabulation(fe,1,1,pointRef,0,&Tchild)); 3256 PetscCall(DMPlexGetTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure)); 3257 for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */ 3258 PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT; 3259 PetscInt l; 3260 const PetscInt *cperms; 3261 3262 PetscCall(DMLabelGetValue(depth,child,&childDepth)); 3263 childDof = depthNumDof[childDepth]; 3264 for (l = 0, cI = -1, childCellShapeOff = 0; l < numClosure; l++) { 3265 PetscInt point = closure[2 * l]; 3266 PetscInt pointDepth; 3267 3268 childO = closure[2 * l + 1]; 3269 if (point == child) { 3270 cI = l; 3271 break; 3272 } 3273 PetscCall(DMLabelGetValue(depth,point,&pointDepth)); 3274 childCellShapeOff += depthNumDof[pointDepth]; 3275 } 3276 if (l == numClosure) { 3277 pointMatOff += childDof; 3278 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */ 3279 } 3280 cperms = perms ? perms[cI] ? perms[cI][childO] : NULL : NULL; 3281 for (l = 0; l < childDof; l++) { 3282 PetscInt lCell = cperms ? cperms[l] : l; 3283 PetscInt childCellDof = childCellShapeOff + lCell; 3284 PetscReal *childValAtPoint; 3285 PetscReal val = 0.; 3286 3287 childValAtPoint = &Tchild->T[0][childCellDof * Nc]; 3288 for (m = 0; m < Nc; m++) { 3289 val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m]; 3290 } 3291 3292 pointMat[i * numChildDof + pointMatOff + l] += val; 3293 } 3294 pointMatOff += childDof; 3295 } 3296 PetscCall(DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure)); 3297 PetscCall(PetscTabulationDestroy(&Tchild)); 3298 } 3299 PetscCall(PetscTabulationDestroy(&Tparent)); 3300 } 3301 } 3302 else { /* just the volume-weighted averages of the children */ 3303 PetscReal parentVol; 3304 PetscInt childCell; 3305 3306 PetscCall(DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL)); 3307 for (i = 0, childCell = 0; i < numChildren; i++) { 3308 PetscInt child = children[i], j; 3309 PetscReal childVol; 3310 3311 if (child < cStart || child >= cEnd) continue; 3312 PetscCall(DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL)); 3313 for (j = 0; j < Nc; j++) { 3314 pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol; 3315 } 3316 childCell++; 3317 } 3318 } 3319 /* Insert pointMat into mat */ 3320 PetscCall(MatSetValues(mat,numSelfDof,matRows,numChildDof,matCols,pointMat,INSERT_VALUES)); 3321 PetscCall(DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows)); 3322 PetscCall(DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat)); 3323 } 3324 } 3325 PetscCall(PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ)); 3326 PetscCall(PetscFree2(pointScalar,pointRef)); 3327 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 3328 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 3329 *inj = mat; 3330 PetscFunctionReturn(0); 3331 } 3332 3333 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3334 { 3335 PetscDS ds; 3336 PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof; 3337 PetscScalar ***refPointFieldMats; 3338 PetscSection refConSec, refSection; 3339 3340 PetscFunctionBegin; 3341 PetscCall(DMGetDS(refTree,&ds)); 3342 PetscCall(PetscDSGetNumFields(ds,&numFields)); 3343 PetscCall(DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL)); 3344 PetscCall(DMGetLocalSection(refTree,&refSection)); 3345 PetscCall(PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd)); 3346 PetscCall(PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats)); 3347 PetscCall(PetscSectionGetMaxDof(refConSec,&maxDof)); 3348 PetscCall(PetscMalloc1(maxDof,&rows)); 3349 PetscCall(PetscMalloc1(maxDof*maxDof,&cols)); 3350 for (p = pRefStart; p < pRefEnd; p++) { 3351 PetscInt parent, pDof, parentDof; 3352 3353 PetscCall(DMPlexGetTreeParent(refTree,p,&parent,NULL)); 3354 PetscCall(PetscSectionGetDof(refConSec,p,&pDof)); 3355 PetscCall(PetscSectionGetDof(refSection,parent,&parentDof)); 3356 if (!pDof || !parentDof || parent == p) continue; 3357 3358 PetscCall(PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart])); 3359 for (f = 0; f < numFields; f++) { 3360 PetscInt cDof, cOff, numCols, r; 3361 3362 if (numFields > 1) { 3363 PetscCall(PetscSectionGetFieldDof(refConSec,p,f,&cDof)); 3364 PetscCall(PetscSectionGetFieldOffset(refConSec,p,f,&cOff)); 3365 } 3366 else { 3367 PetscCall(PetscSectionGetDof(refConSec,p,&cDof)); 3368 PetscCall(PetscSectionGetOffset(refConSec,p,&cOff)); 3369 } 3370 3371 for (r = 0; r < cDof; r++) { 3372 rows[r] = cOff + r; 3373 } 3374 numCols = 0; 3375 { 3376 PetscInt aDof, aOff, j; 3377 3378 if (numFields > 1) { 3379 PetscCall(PetscSectionGetFieldDof(refSection,parent,f,&aDof)); 3380 PetscCall(PetscSectionGetFieldOffset(refSection,parent,f,&aOff)); 3381 } 3382 else { 3383 PetscCall(PetscSectionGetDof(refSection,parent,&aDof)); 3384 PetscCall(PetscSectionGetOffset(refSection,parent,&aOff)); 3385 } 3386 3387 for (j = 0; j < aDof; j++) { 3388 cols[numCols++] = aOff + j; 3389 } 3390 } 3391 PetscCall(PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f])); 3392 /* transpose of constraint matrix */ 3393 PetscCall(MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f])); 3394 } 3395 } 3396 *childrenMats = refPointFieldMats; 3397 PetscCall(PetscFree(rows)); 3398 PetscCall(PetscFree(cols)); 3399 PetscFunctionReturn(0); 3400 } 3401 3402 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3403 { 3404 PetscDS ds; 3405 PetscScalar ***refPointFieldMats; 3406 PetscInt numFields, pRefStart, pRefEnd, p, f; 3407 PetscSection refConSec, refSection; 3408 3409 PetscFunctionBegin; 3410 refPointFieldMats = *childrenMats; 3411 *childrenMats = NULL; 3412 PetscCall(DMGetDS(refTree,&ds)); 3413 PetscCall(DMGetLocalSection(refTree,&refSection)); 3414 PetscCall(PetscDSGetNumFields(ds,&numFields)); 3415 PetscCall(DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL)); 3416 PetscCall(PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd)); 3417 for (p = pRefStart; p < pRefEnd; p++) { 3418 PetscInt parent, pDof, parentDof; 3419 3420 PetscCall(DMPlexGetTreeParent(refTree,p,&parent,NULL)); 3421 PetscCall(PetscSectionGetDof(refConSec,p,&pDof)); 3422 PetscCall(PetscSectionGetDof(refSection,parent,&parentDof)); 3423 if (!pDof || !parentDof || parent == p) continue; 3424 3425 for (f = 0; f < numFields; f++) { 3426 PetscInt cDof; 3427 3428 if (numFields > 1) { 3429 PetscCall(PetscSectionGetFieldDof(refConSec,p,f,&cDof)); 3430 } 3431 else { 3432 PetscCall(PetscSectionGetDof(refConSec,p,&cDof)); 3433 } 3434 3435 PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f])); 3436 } 3437 PetscCall(PetscFree(refPointFieldMats[p - pRefStart])); 3438 } 3439 PetscCall(PetscFree(refPointFieldMats)); 3440 PetscFunctionReturn(0); 3441 } 3442 3443 static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree,Mat *injRef) 3444 { 3445 Mat cMatRef; 3446 PetscObject injRefObj; 3447 3448 PetscFunctionBegin; 3449 PetscCall(DMGetDefaultConstraints(refTree,NULL,&cMatRef,NULL)); 3450 PetscCall(PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj)); 3451 *injRef = (Mat) injRefObj; 3452 if (!*injRef) { 3453 PetscCall(DMPlexComputeInjectorReferenceTree(refTree,injRef)); 3454 PetscCall(PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)*injRef)); 3455 /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */ 3456 PetscCall(PetscObjectDereference((PetscObject)*injRef)); 3457 } 3458 PetscFunctionReturn(0); 3459 } 3460 3461 static PetscErrorCode DMPlexTransferInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, const PetscInt *childIds, Vec fineVec, PetscInt numFields, PetscInt *offsets, PetscSection *rootMultiSec, PetscSection *multiLeafSec, PetscInt **gatheredIndices, PetscScalar **gatheredValues) 3462 { 3463 PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti; 3464 PetscSection globalCoarse, globalFine; 3465 PetscSection localCoarse, localFine, leafIndicesSec; 3466 PetscSection multiRootSec, rootIndicesSec; 3467 PetscInt *leafInds, *rootInds = NULL; 3468 const PetscInt *rootDegrees; 3469 PetscScalar *leafVals = NULL, *rootVals = NULL; 3470 PetscSF coarseToFineEmbedded; 3471 3472 PetscFunctionBegin; 3473 PetscCall(DMPlexGetChart(coarse,&pStartC,&pEndC)); 3474 PetscCall(DMPlexGetChart(fine,&pStartF,&pEndF)); 3475 PetscCall(DMGetLocalSection(fine,&localFine)); 3476 PetscCall(DMGetGlobalSection(fine,&globalFine)); 3477 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec)); 3478 PetscCall(PetscSectionSetChart(leafIndicesSec,pStartF, pEndF)); 3479 PetscCall(PetscSectionGetMaxDof(localFine,&maxDof)); 3480 { /* winnow fine points that don't have global dofs out of the sf */ 3481 PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices; 3482 const PetscInt *leaves; 3483 3484 PetscCall(PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL)); 3485 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3486 p = leaves ? leaves[l] : l; 3487 PetscCall(PetscSectionGetDof(globalFine,p,&dof)); 3488 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&cdof)); 3489 if ((dof - cdof) > 0) { 3490 numPointsWithDofs++; 3491 3492 PetscCall(PetscSectionGetDof(localFine,p,&dof)); 3493 PetscCall(PetscSectionSetDof(leafIndicesSec,p,dof + 1)); 3494 } 3495 } 3496 PetscCall(PetscMalloc1(numPointsWithDofs,&pointsWithDofs)); 3497 PetscCall(PetscSectionSetUp(leafIndicesSec)); 3498 PetscCall(PetscSectionGetStorageSize(leafIndicesSec,&numIndices)); 3499 PetscCall(PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1),&leafInds)); 3500 if (gatheredValues) PetscCall(PetscMalloc1(numIndices,&leafVals)); 3501 for (l = 0, offset = 0; l < nleaves; l++) { 3502 p = leaves ? leaves[l] : l; 3503 PetscCall(PetscSectionGetDof(globalFine,p,&dof)); 3504 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&cdof)); 3505 if ((dof - cdof) > 0) { 3506 PetscInt off, gOff; 3507 PetscInt *pInd; 3508 PetscScalar *pVal = NULL; 3509 3510 pointsWithDofs[offset++] = l; 3511 3512 PetscCall(PetscSectionGetOffset(leafIndicesSec,p,&off)); 3513 3514 pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds; 3515 if (gatheredValues) { 3516 PetscInt i; 3517 3518 pVal = &leafVals[off + 1]; 3519 for (i = 0; i < dof; i++) pVal[i] = 0.; 3520 } 3521 PetscCall(PetscSectionGetOffset(globalFine,p,&gOff)); 3522 3523 offsets[0] = 0; 3524 if (numFields) { 3525 PetscInt f; 3526 3527 for (f = 0; f < numFields; f++) { 3528 PetscInt fDof; 3529 PetscCall(PetscSectionGetFieldDof(localFine,p,f,&fDof)); 3530 offsets[f + 1] = fDof + offsets[f]; 3531 } 3532 PetscCall(DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1, NULL,pInd)); 3533 } else { 3534 PetscCall(DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL, NULL,pInd)); 3535 } 3536 if (gatheredValues) PetscCall(VecGetValues(fineVec,dof,pInd,pVal)); 3537 } 3538 } 3539 PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded)); 3540 PetscCall(PetscFree(pointsWithDofs)); 3541 } 3542 3543 PetscCall(DMPlexGetChart(coarse,&pStartC,&pEndC)); 3544 PetscCall(DMGetLocalSection(coarse,&localCoarse)); 3545 PetscCall(DMGetGlobalSection(coarse,&globalCoarse)); 3546 3547 { /* there may be the case where an sf root has a parent: broadcast parents back to children */ 3548 MPI_Datatype threeInt; 3549 PetscMPIInt rank; 3550 PetscInt (*parentNodeAndIdCoarse)[3]; 3551 PetscInt (*parentNodeAndIdFine)[3]; 3552 PetscInt p, nleaves, nleavesToParents; 3553 PetscSF pointSF, sfToParents; 3554 const PetscInt *ilocal; 3555 const PetscSFNode *iremote; 3556 PetscSFNode *iremoteToParents; 3557 PetscInt *ilocalToParents; 3558 3559 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank)); 3560 PetscCallMPI(MPI_Type_contiguous(3,MPIU_INT,&threeInt)); 3561 PetscCallMPI(MPI_Type_commit(&threeInt)); 3562 PetscCall(PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine)); 3563 PetscCall(DMGetPointSF(coarse,&pointSF)); 3564 PetscCall(PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote)); 3565 for (p = pStartC; p < pEndC; p++) { 3566 PetscInt parent, childId; 3567 PetscCall(DMPlexGetTreeParent(coarse,p,&parent,&childId)); 3568 parentNodeAndIdCoarse[p - pStartC][0] = rank; 3569 parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC; 3570 parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId; 3571 if (nleaves > 0) { 3572 PetscInt leaf = -1; 3573 3574 if (ilocal) { 3575 PetscCall(PetscFindInt(parent,nleaves,ilocal,&leaf)); 3576 } 3577 else { 3578 leaf = p - pStartC; 3579 } 3580 if (leaf >= 0) { 3581 parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank; 3582 parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index; 3583 } 3584 } 3585 } 3586 for (p = pStartF; p < pEndF; p++) { 3587 parentNodeAndIdFine[p - pStartF][0] = -1; 3588 parentNodeAndIdFine[p - pStartF][1] = -1; 3589 parentNodeAndIdFine[p - pStartF][2] = -1; 3590 } 3591 PetscCall(PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine,MPI_REPLACE)); 3592 PetscCall(PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine,MPI_REPLACE)); 3593 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3594 PetscInt dof; 3595 3596 PetscCall(PetscSectionGetDof(leafIndicesSec,p,&dof)); 3597 if (dof) { 3598 PetscInt off; 3599 3600 PetscCall(PetscSectionGetOffset(leafIndicesSec,p,&off)); 3601 if (gatheredIndices) { 3602 leafInds[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3603 } else if (gatheredValues) { 3604 leafVals[off] = (PetscScalar) PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3605 } 3606 } 3607 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3608 nleavesToParents++; 3609 } 3610 } 3611 PetscCall(PetscMalloc1(nleavesToParents,&ilocalToParents)); 3612 PetscCall(PetscMalloc1(nleavesToParents,&iremoteToParents)); 3613 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3614 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3615 ilocalToParents[nleavesToParents] = p - pStartF; 3616 iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p-pStartF][0]; 3617 iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1]; 3618 nleavesToParents++; 3619 } 3620 } 3621 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents)); 3622 PetscCall(PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER)); 3623 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 3624 3625 coarseToFineEmbedded = sfToParents; 3626 3627 PetscCall(PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine)); 3628 PetscCallMPI(MPI_Type_free(&threeInt)); 3629 } 3630 3631 { /* winnow out coarse points that don't have dofs */ 3632 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3633 PetscSF sfDofsOnly; 3634 3635 for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) { 3636 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 3637 PetscCall(PetscSectionGetConstraintDof(globalCoarse,p,&cdof)); 3638 if ((dof - cdof) > 0) { 3639 numPointsWithDofs++; 3640 } 3641 } 3642 PetscCall(PetscMalloc1(numPointsWithDofs,&pointsWithDofs)); 3643 for (p = pStartC, offset = 0; p < pEndC; p++) { 3644 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 3645 PetscCall(PetscSectionGetConstraintDof(globalCoarse,p,&cdof)); 3646 if ((dof - cdof) > 0) { 3647 pointsWithDofs[offset++] = p - pStartC; 3648 } 3649 } 3650 PetscCall(PetscSFCreateEmbeddedRootSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly)); 3651 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 3652 PetscCall(PetscFree(pointsWithDofs)); 3653 coarseToFineEmbedded = sfDofsOnly; 3654 } 3655 3656 /* communicate back to the coarse mesh which coarse points have children (that may require injection) */ 3657 PetscCall(PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees)); 3658 PetscCall(PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees)); 3659 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec)); 3660 PetscCall(PetscSectionSetChart(multiRootSec,pStartC,pEndC)); 3661 for (p = pStartC; p < pEndC; p++) { 3662 PetscCall(PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC])); 3663 } 3664 PetscCall(PetscSectionSetUp(multiRootSec)); 3665 PetscCall(PetscSectionGetStorageSize(multiRootSec,&numMulti)); 3666 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec)); 3667 { /* distribute the leaf section */ 3668 PetscSF multi, multiInv, indicesSF; 3669 PetscInt *remoteOffsets, numRootIndices; 3670 3671 PetscCall(PetscSFGetMultiSF(coarseToFineEmbedded,&multi)); 3672 PetscCall(PetscSFCreateInverseSF(multi,&multiInv)); 3673 PetscCall(PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec)); 3674 PetscCall(PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF)); 3675 PetscCall(PetscFree(remoteOffsets)); 3676 PetscCall(PetscSFDestroy(&multiInv)); 3677 PetscCall(PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices)); 3678 if (gatheredIndices) { 3679 PetscCall(PetscMalloc1(numRootIndices,&rootInds)); 3680 PetscCall(PetscSFBcastBegin(indicesSF,MPIU_INT,leafInds,rootInds,MPI_REPLACE)); 3681 PetscCall(PetscSFBcastEnd(indicesSF,MPIU_INT,leafInds,rootInds,MPI_REPLACE)); 3682 } 3683 if (gatheredValues) { 3684 PetscCall(PetscMalloc1(numRootIndices,&rootVals)); 3685 PetscCall(PetscSFBcastBegin(indicesSF,MPIU_SCALAR,leafVals,rootVals,MPI_REPLACE)); 3686 PetscCall(PetscSFBcastEnd(indicesSF,MPIU_SCALAR,leafVals,rootVals,MPI_REPLACE)); 3687 } 3688 PetscCall(PetscSFDestroy(&indicesSF)); 3689 } 3690 PetscCall(PetscSectionDestroy(&leafIndicesSec)); 3691 PetscCall(PetscFree(leafInds)); 3692 PetscCall(PetscFree(leafVals)); 3693 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 3694 *rootMultiSec = multiRootSec; 3695 *multiLeafSec = rootIndicesSec; 3696 if (gatheredIndices) *gatheredIndices = rootInds; 3697 if (gatheredValues) *gatheredValues = rootVals; 3698 PetscFunctionReturn(0); 3699 } 3700 3701 PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 3702 { 3703 DM refTree; 3704 PetscSection multiRootSec, rootIndicesSec; 3705 PetscSection globalCoarse, globalFine; 3706 PetscSection localCoarse, localFine; 3707 PetscSection cSecRef; 3708 PetscInt *rootIndices = NULL, *parentIndices, pRefStart, pRefEnd; 3709 Mat injRef; 3710 PetscInt numFields, maxDof; 3711 PetscInt pStartC, pEndC, pStartF, pEndF, p; 3712 PetscInt *offsets, *offsetsCopy, *rowOffsets; 3713 PetscLayout rowMap, colMap; 3714 PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO; 3715 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 3716 3717 PetscFunctionBegin; 3718 3719 /* get the templates for the fine-to-coarse injection from the reference tree */ 3720 PetscCall(DMPlexGetReferenceTree(coarse,&refTree)); 3721 PetscCall(DMGetDefaultConstraints(refTree,&cSecRef,NULL,NULL)); 3722 PetscCall(PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd)); 3723 PetscCall(DMPlexReferenceTreeGetInjector(refTree,&injRef)); 3724 3725 PetscCall(DMPlexGetChart(fine,&pStartF,&pEndF)); 3726 PetscCall(DMGetLocalSection(fine,&localFine)); 3727 PetscCall(DMGetGlobalSection(fine,&globalFine)); 3728 PetscCall(PetscSectionGetNumFields(localFine,&numFields)); 3729 PetscCall(DMPlexGetChart(coarse,&pStartC,&pEndC)); 3730 PetscCall(DMGetLocalSection(coarse,&localCoarse)); 3731 PetscCall(DMGetGlobalSection(coarse,&globalCoarse)); 3732 PetscCall(PetscSectionGetMaxDof(localCoarse,&maxDof)); 3733 { 3734 PetscInt maxFields = PetscMax(1,numFields) + 1; 3735 PetscCall(PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets)); 3736 } 3737 3738 PetscCall(DMPlexTransferInjectorTree(coarse,fine,coarseToFine,childIds,NULL,numFields,offsets,&multiRootSec,&rootIndicesSec,&rootIndices,NULL)); 3739 3740 PetscCall(PetscMalloc1(maxDof,&parentIndices)); 3741 3742 /* count indices */ 3743 PetscCall(MatGetLayouts(mat,&rowMap,&colMap)); 3744 PetscCall(PetscLayoutSetUp(rowMap)); 3745 PetscCall(PetscLayoutSetUp(colMap)); 3746 PetscCall(PetscLayoutGetRange(rowMap,&rowStart,&rowEnd)); 3747 PetscCall(PetscLayoutGetRange(colMap,&colStart,&colEnd)); 3748 PetscCall(PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO)); 3749 for (p = pStartC; p < pEndC; p++) { 3750 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3751 3752 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 3753 PetscCall(PetscSectionGetConstraintDof(globalCoarse,p,&cdof)); 3754 if ((dof - cdof) <= 0) continue; 3755 PetscCall(PetscSectionGetOffset(globalCoarse,p,&gOff)); 3756 3757 rowOffsets[0] = 0; 3758 offsetsCopy[0] = 0; 3759 if (numFields) { 3760 PetscInt f; 3761 3762 for (f = 0; f < numFields; f++) { 3763 PetscInt fDof; 3764 PetscCall(PetscSectionGetFieldDof(localCoarse,p,f,&fDof)); 3765 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3766 } 3767 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,parentIndices)); 3768 } else { 3769 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,parentIndices)); 3770 rowOffsets[1] = offsetsCopy[0]; 3771 } 3772 3773 PetscCall(PetscSectionGetDof(multiRootSec,p,&numLeaves)); 3774 PetscCall(PetscSectionGetOffset(multiRootSec,p,&leafStart)); 3775 leafEnd = leafStart + numLeaves; 3776 for (l = leafStart; l < leafEnd; l++) { 3777 PetscInt numIndices, childId, offset; 3778 const PetscInt *childIndices; 3779 3780 PetscCall(PetscSectionGetDof(rootIndicesSec,l,&numIndices)); 3781 PetscCall(PetscSectionGetOffset(rootIndicesSec,l,&offset)); 3782 childId = rootIndices[offset++]; 3783 childIndices = &rootIndices[offset]; 3784 numIndices--; 3785 3786 if (childId == -1) { /* equivalent points: scatter */ 3787 PetscInt i; 3788 3789 for (i = 0; i < numIndices; i++) { 3790 PetscInt colIndex = childIndices[i]; 3791 PetscInt rowIndex = parentIndices[i]; 3792 if (rowIndex < 0) continue; 3793 PetscCheck(colIndex >= 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unconstrained fine and constrained coarse"); 3794 if (colIndex >= colStart && colIndex < colEnd) { 3795 nnzD[rowIndex - rowStart] = 1; 3796 } 3797 else { 3798 nnzO[rowIndex - rowStart] = 1; 3799 } 3800 } 3801 } 3802 else { 3803 PetscInt parentId, f, lim; 3804 3805 PetscCall(DMPlexGetTreeParent(refTree,childId,&parentId,NULL)); 3806 3807 lim = PetscMax(1,numFields); 3808 offsets[0] = 0; 3809 if (numFields) { 3810 PetscInt f; 3811 3812 for (f = 0; f < numFields; f++) { 3813 PetscInt fDof; 3814 PetscCall(PetscSectionGetFieldDof(cSecRef,childId,f,&fDof)); 3815 3816 offsets[f + 1] = fDof + offsets[f]; 3817 } 3818 } 3819 else { 3820 PetscInt cDof; 3821 3822 PetscCall(PetscSectionGetDof(cSecRef,childId,&cDof)); 3823 offsets[1] = cDof; 3824 } 3825 for (f = 0; f < lim; f++) { 3826 PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1]; 3827 PetscInt childStart = offsets[f], childEnd = offsets[f + 1]; 3828 PetscInt i, numD = 0, numO = 0; 3829 3830 for (i = childStart; i < childEnd; i++) { 3831 PetscInt colIndex = childIndices[i]; 3832 3833 if (colIndex < 0) continue; 3834 if (colIndex >= colStart && colIndex < colEnd) { 3835 numD++; 3836 } 3837 else { 3838 numO++; 3839 } 3840 } 3841 for (i = parentStart; i < parentEnd; i++) { 3842 PetscInt rowIndex = parentIndices[i]; 3843 3844 if (rowIndex < 0) continue; 3845 nnzD[rowIndex - rowStart] += numD; 3846 nnzO[rowIndex - rowStart] += numO; 3847 } 3848 } 3849 } 3850 } 3851 } 3852 /* preallocate */ 3853 PetscCall(MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL)); 3854 PetscCall(PetscFree2(nnzD,nnzO)); 3855 /* insert values */ 3856 PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats)); 3857 for (p = pStartC; p < pEndC; p++) { 3858 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3859 3860 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 3861 PetscCall(PetscSectionGetConstraintDof(globalCoarse,p,&cdof)); 3862 if ((dof - cdof) <= 0) continue; 3863 PetscCall(PetscSectionGetOffset(globalCoarse,p,&gOff)); 3864 3865 rowOffsets[0] = 0; 3866 offsetsCopy[0] = 0; 3867 if (numFields) { 3868 PetscInt f; 3869 3870 for (f = 0; f < numFields; f++) { 3871 PetscInt fDof; 3872 PetscCall(PetscSectionGetFieldDof(localCoarse,p,f,&fDof)); 3873 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3874 } 3875 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,parentIndices)); 3876 } else { 3877 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,parentIndices)); 3878 rowOffsets[1] = offsetsCopy[0]; 3879 } 3880 3881 PetscCall(PetscSectionGetDof(multiRootSec,p,&numLeaves)); 3882 PetscCall(PetscSectionGetOffset(multiRootSec,p,&leafStart)); 3883 leafEnd = leafStart + numLeaves; 3884 for (l = leafStart; l < leafEnd; l++) { 3885 PetscInt numIndices, childId, offset; 3886 const PetscInt *childIndices; 3887 3888 PetscCall(PetscSectionGetDof(rootIndicesSec,l,&numIndices)); 3889 PetscCall(PetscSectionGetOffset(rootIndicesSec,l,&offset)); 3890 childId = rootIndices[offset++]; 3891 childIndices = &rootIndices[offset]; 3892 numIndices--; 3893 3894 if (childId == -1) { /* equivalent points: scatter */ 3895 PetscInt i; 3896 3897 for (i = 0; i < numIndices; i++) { 3898 PetscCall(MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES)); 3899 } 3900 } 3901 else { 3902 PetscInt parentId, f, lim; 3903 3904 PetscCall(DMPlexGetTreeParent(refTree,childId,&parentId,NULL)); 3905 3906 lim = PetscMax(1,numFields); 3907 offsets[0] = 0; 3908 if (numFields) { 3909 PetscInt f; 3910 3911 for (f = 0; f < numFields; f++) { 3912 PetscInt fDof; 3913 PetscCall(PetscSectionGetFieldDof(cSecRef,childId,f,&fDof)); 3914 3915 offsets[f + 1] = fDof + offsets[f]; 3916 } 3917 } 3918 else { 3919 PetscInt cDof; 3920 3921 PetscCall(PetscSectionGetDof(cSecRef,childId,&cDof)); 3922 offsets[1] = cDof; 3923 } 3924 for (f = 0; f < lim; f++) { 3925 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 3926 PetscInt *rowIndices = &parentIndices[rowOffsets[f]]; 3927 const PetscInt *colIndices = &childIndices[offsets[f]]; 3928 3929 PetscCall(MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES)); 3930 } 3931 } 3932 } 3933 } 3934 PetscCall(PetscSectionDestroy(&multiRootSec)); 3935 PetscCall(PetscSectionDestroy(&rootIndicesSec)); 3936 PetscCall(PetscFree(parentIndices)); 3937 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats)); 3938 PetscCall(PetscFree(rootIndices)); 3939 PetscCall(PetscFree3(offsets,offsetsCopy,rowOffsets)); 3940 3941 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 3942 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 3943 PetscFunctionReturn(0); 3944 } 3945 3946 static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom) 3947 { 3948 PetscSF coarseToFineEmbedded; 3949 PetscSection globalCoarse, globalFine; 3950 PetscSection localCoarse, localFine; 3951 PetscSection aSec, cSec; 3952 PetscSection rootValuesSec; 3953 PetscSection leafValuesSec; 3954 PetscScalar *rootValues, *leafValues; 3955 IS aIS; 3956 const PetscInt *anchors; 3957 Mat cMat; 3958 PetscInt numFields; 3959 PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd; 3960 PetscInt aStart, aEnd, cStart, cEnd; 3961 PetscInt *maxChildIds; 3962 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 3963 PetscFV fv = NULL; 3964 PetscInt dim, numFVcomps = -1, fvField = -1; 3965 DM cellDM = NULL, gradDM = NULL; 3966 const PetscScalar *cellGeomArray = NULL; 3967 const PetscScalar *gradArray = NULL; 3968 3969 PetscFunctionBegin; 3970 PetscCall(VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE)); 3971 PetscCall(DMPlexGetChart(coarse,&pStartC,&pEndC)); 3972 PetscCall(DMPlexGetSimplexOrBoxCells(coarse,0,&cellStart,&cellEnd)); 3973 PetscCall(DMPlexGetChart(fine,&pStartF,&pEndF)); 3974 PetscCall(DMGetGlobalSection(fine,&globalFine)); 3975 PetscCall(DMGetCoordinateDim(coarse,&dim)); 3976 { /* winnow fine points that don't have global dofs out of the sf */ 3977 PetscInt nleaves, l; 3978 const PetscInt *leaves; 3979 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3980 3981 PetscCall(PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL)); 3982 3983 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3984 PetscInt p = leaves ? leaves[l] : l; 3985 3986 PetscCall(PetscSectionGetDof(globalFine,p,&dof)); 3987 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&cdof)); 3988 if ((dof - cdof) > 0) { 3989 numPointsWithDofs++; 3990 } 3991 } 3992 PetscCall(PetscMalloc1(numPointsWithDofs,&pointsWithDofs)); 3993 for (l = 0, offset = 0; l < nleaves; l++) { 3994 PetscInt p = leaves ? leaves[l] : l; 3995 3996 PetscCall(PetscSectionGetDof(globalFine,p,&dof)); 3997 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&cdof)); 3998 if ((dof - cdof) > 0) { 3999 pointsWithDofs[offset++] = l; 4000 } 4001 } 4002 PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded)); 4003 PetscCall(PetscFree(pointsWithDofs)); 4004 } 4005 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 4006 PetscCall(PetscMalloc1(pEndC-pStartC,&maxChildIds)); 4007 for (p = pStartC; p < pEndC; p++) { 4008 maxChildIds[p - pStartC] = -2; 4009 } 4010 PetscCall(PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX)); 4011 PetscCall(PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX)); 4012 4013 PetscCall(DMGetLocalSection(coarse,&localCoarse)); 4014 PetscCall(DMGetGlobalSection(coarse,&globalCoarse)); 4015 4016 PetscCall(DMPlexGetAnchors(coarse,&aSec,&aIS)); 4017 PetscCall(ISGetIndices(aIS,&anchors)); 4018 PetscCall(PetscSectionGetChart(aSec,&aStart,&aEnd)); 4019 4020 PetscCall(DMGetDefaultConstraints(coarse,&cSec,&cMat,NULL)); 4021 PetscCall(PetscSectionGetChart(cSec,&cStart,&cEnd)); 4022 4023 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 4024 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootValuesSec)); 4025 PetscCall(PetscSectionSetChart(rootValuesSec,pStartC,pEndC)); 4026 PetscCall(PetscSectionGetNumFields(localCoarse,&numFields)); 4027 { 4028 PetscInt maxFields = PetscMax(1,numFields) + 1; 4029 PetscCall(PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO)); 4030 } 4031 if (grad) { 4032 PetscInt i; 4033 4034 PetscCall(VecGetDM(cellGeom,&cellDM)); 4035 PetscCall(VecGetArrayRead(cellGeom,&cellGeomArray)); 4036 PetscCall(VecGetDM(grad,&gradDM)); 4037 PetscCall(VecGetArrayRead(grad,&gradArray)); 4038 for (i = 0; i < PetscMax(1,numFields); i++) { 4039 PetscObject obj; 4040 PetscClassId id; 4041 4042 PetscCall(DMGetField(coarse, i, NULL, &obj)); 4043 PetscCall(PetscObjectGetClassId(obj,&id)); 4044 if (id == PETSCFV_CLASSID) { 4045 fv = (PetscFV) obj; 4046 PetscCall(PetscFVGetNumComponents(fv,&numFVcomps)); 4047 fvField = i; 4048 break; 4049 } 4050 } 4051 } 4052 4053 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 4054 PetscInt dof; 4055 PetscInt maxChildId = maxChildIds[p - pStartC]; 4056 PetscInt numValues = 0; 4057 4058 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 4059 if (dof < 0) { 4060 dof = -(dof + 1); 4061 } 4062 offsets[0] = 0; 4063 newOffsets[0] = 0; 4064 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 4065 PetscInt *closure = NULL, closureSize, cl; 4066 4067 PetscCall(DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure)); 4068 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 4069 PetscInt c = closure[2 * cl], clDof; 4070 4071 PetscCall(PetscSectionGetDof(localCoarse,c,&clDof)); 4072 numValues += clDof; 4073 } 4074 PetscCall(DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure)); 4075 } 4076 else if (maxChildId == -1) { 4077 PetscCall(PetscSectionGetDof(localCoarse,p,&numValues)); 4078 } 4079 /* we will pack the column indices with the field offsets */ 4080 if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) { 4081 /* also send the centroid, and the gradient */ 4082 numValues += dim * (1 + numFVcomps); 4083 } 4084 PetscCall(PetscSectionSetDof(rootValuesSec,p,numValues)); 4085 } 4086 PetscCall(PetscSectionSetUp(rootValuesSec)); 4087 { 4088 PetscInt numRootValues; 4089 const PetscScalar *coarseArray; 4090 4091 PetscCall(PetscSectionGetStorageSize(rootValuesSec,&numRootValues)); 4092 PetscCall(PetscMalloc1(numRootValues,&rootValues)); 4093 PetscCall(VecGetArrayRead(vecCoarseLocal,&coarseArray)); 4094 for (p = pStartC; p < pEndC; p++) { 4095 PetscInt numValues; 4096 PetscInt pValOff; 4097 PetscScalar *pVal; 4098 PetscInt maxChildId = maxChildIds[p - pStartC]; 4099 4100 PetscCall(PetscSectionGetDof(rootValuesSec,p,&numValues)); 4101 if (!numValues) { 4102 continue; 4103 } 4104 PetscCall(PetscSectionGetOffset(rootValuesSec,p,&pValOff)); 4105 pVal = &(rootValues[pValOff]); 4106 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 4107 PetscInt closureSize = numValues; 4108 PetscCall(DMPlexVecGetClosure(coarse,NULL,vecCoarseLocal,p,&closureSize,&pVal)); 4109 if (grad && p >= cellStart && p < cellEnd) { 4110 PetscFVCellGeom *cg; 4111 PetscScalar *gradVals = NULL; 4112 PetscInt i; 4113 4114 pVal += (numValues - dim * (1 + numFVcomps)); 4115 4116 PetscCall(DMPlexPointLocalRead(cellDM,p,cellGeomArray,(void *) &cg)); 4117 for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i]; 4118 pVal += dim; 4119 PetscCall(DMPlexPointGlobalRead(gradDM,p,gradArray,(void *) &gradVals)); 4120 for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i]; 4121 } 4122 } 4123 else if (maxChildId == -1) { 4124 PetscInt lDof, lOff, i; 4125 4126 PetscCall(PetscSectionGetDof(localCoarse,p,&lDof)); 4127 PetscCall(PetscSectionGetOffset(localCoarse,p,&lOff)); 4128 for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i]; 4129 } 4130 } 4131 PetscCall(VecRestoreArrayRead(vecCoarseLocal,&coarseArray)); 4132 PetscCall(PetscFree(maxChildIds)); 4133 } 4134 { 4135 PetscSF valuesSF; 4136 PetscInt *remoteOffsetsValues, numLeafValues; 4137 4138 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafValuesSec)); 4139 PetscCall(PetscSFDistributeSection(coarseToFineEmbedded,rootValuesSec,&remoteOffsetsValues,leafValuesSec)); 4140 PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded,rootValuesSec,remoteOffsetsValues,leafValuesSec,&valuesSF)); 4141 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 4142 PetscCall(PetscFree(remoteOffsetsValues)); 4143 PetscCall(PetscSectionGetStorageSize(leafValuesSec,&numLeafValues)); 4144 PetscCall(PetscMalloc1(numLeafValues,&leafValues)); 4145 PetscCall(PetscSFBcastBegin(valuesSF,MPIU_SCALAR,rootValues,leafValues,MPI_REPLACE)); 4146 PetscCall(PetscSFBcastEnd(valuesSF,MPIU_SCALAR,rootValues,leafValues,MPI_REPLACE)); 4147 PetscCall(PetscSFDestroy(&valuesSF)); 4148 PetscCall(PetscFree(rootValues)); 4149 PetscCall(PetscSectionDestroy(&rootValuesSec)); 4150 } 4151 PetscCall(DMGetLocalSection(fine,&localFine)); 4152 { 4153 PetscInt maxDof; 4154 PetscInt *rowIndices; 4155 DM refTree; 4156 PetscInt **refPointFieldN; 4157 PetscScalar ***refPointFieldMats; 4158 PetscSection refConSec, refAnSec; 4159 PetscInt pRefStart,pRefEnd,leafStart,leafEnd; 4160 PetscScalar *pointWork; 4161 4162 PetscCall(PetscSectionGetMaxDof(localFine,&maxDof)); 4163 PetscCall(DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices)); 4164 PetscCall(DMGetWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork)); 4165 PetscCall(DMPlexGetReferenceTree(fine,&refTree)); 4166 PetscCall(DMCopyDisc(fine,refTree)); 4167 PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN)); 4168 PetscCall(DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL)); 4169 PetscCall(DMPlexGetAnchors(refTree,&refAnSec,NULL)); 4170 PetscCall(PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd)); 4171 PetscCall(PetscSectionGetChart(leafValuesSec,&leafStart,&leafEnd)); 4172 PetscCall(DMPlexGetSimplexOrBoxCells(fine,0,&cellStart,&cellEnd)); 4173 for (p = leafStart; p < leafEnd; p++) { 4174 PetscInt gDof, gcDof, gOff, lDof; 4175 PetscInt numValues, pValOff; 4176 PetscInt childId; 4177 const PetscScalar *pVal; 4178 const PetscScalar *fvGradData = NULL; 4179 4180 PetscCall(PetscSectionGetDof(globalFine,p,&gDof)); 4181 PetscCall(PetscSectionGetDof(localFine,p,&lDof)); 4182 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&gcDof)); 4183 if ((gDof - gcDof) <= 0) { 4184 continue; 4185 } 4186 PetscCall(PetscSectionGetOffset(globalFine,p,&gOff)); 4187 PetscCall(PetscSectionGetDof(leafValuesSec,p,&numValues)); 4188 if (!numValues) continue; 4189 PetscCall(PetscSectionGetOffset(leafValuesSec,p,&pValOff)); 4190 pVal = &leafValues[pValOff]; 4191 offsets[0] = 0; 4192 offsetsCopy[0] = 0; 4193 newOffsets[0] = 0; 4194 newOffsetsCopy[0] = 0; 4195 childId = cids[p - pStartF]; 4196 if (numFields) { 4197 PetscInt f; 4198 for (f = 0; f < numFields; f++) { 4199 PetscInt rowDof; 4200 4201 PetscCall(PetscSectionGetFieldDof(localFine,p,f,&rowDof)); 4202 offsets[f + 1] = offsets[f] + rowDof; 4203 offsetsCopy[f + 1] = offsets[f + 1]; 4204 /* TODO: closure indices */ 4205 newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]); 4206 } 4207 PetscCall(DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,NULL,rowIndices)); 4208 } 4209 else { 4210 offsets[0] = 0; 4211 offsets[1] = lDof; 4212 newOffsets[0] = 0; 4213 newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0]; 4214 PetscCall(DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,NULL,rowIndices)); 4215 } 4216 if (childId == -1) { /* no child interpolation: one nnz per */ 4217 PetscCall(VecSetValues(vecFine,numValues,rowIndices,pVal,INSERT_VALUES)); 4218 } else { 4219 PetscInt f; 4220 4221 if (grad && p >= cellStart && p < cellEnd) { 4222 numValues -= (dim * (1 + numFVcomps)); 4223 fvGradData = &pVal[numValues]; 4224 } 4225 for (f = 0; f < PetscMax(1,numFields); f++) { 4226 const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f]; 4227 PetscInt numRows = offsets[f+1] - offsets[f]; 4228 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 4229 const PetscScalar *cVal = &pVal[newOffsets[f]]; 4230 PetscScalar *rVal = &pointWork[offsets[f]]; 4231 PetscInt i, j; 4232 4233 #if 0 4234 PetscCall(PetscInfo(coarse,"childId %" PetscInt_FMT ", numRows %" PetscInt_FMT ", numCols %" PetscInt_FMT ", refPointFieldN %" PetscInt_FMT " maxDof %" PetscInt_FMT "\n",childId,numRows,numCols,refPointFieldN[childId - pRefStart][f], maxDof)); 4235 #endif 4236 for (i = 0; i < numRows; i++) { 4237 PetscScalar val = 0.; 4238 for (j = 0; j < numCols; j++) { 4239 val += childMat[i * numCols + j] * cVal[j]; 4240 } 4241 rVal[i] = val; 4242 } 4243 if (f == fvField && p >= cellStart && p < cellEnd) { 4244 PetscReal centroid[3]; 4245 PetscScalar diff[3]; 4246 const PetscScalar *parentCentroid = &fvGradData[0]; 4247 const PetscScalar *gradient = &fvGradData[dim]; 4248 4249 PetscCall(DMPlexComputeCellGeometryFVM(fine,p,NULL,centroid,NULL)); 4250 for (i = 0; i < dim; i++) { 4251 diff[i] = centroid[i] - parentCentroid[i]; 4252 } 4253 for (i = 0; i < numFVcomps; i++) { 4254 PetscScalar val = 0.; 4255 4256 for (j = 0; j < dim; j++) { 4257 val += gradient[dim * i + j] * diff[j]; 4258 } 4259 rVal[i] += val; 4260 } 4261 } 4262 PetscCall(VecSetValues(vecFine,numRows,&rowIndices[offsets[f]],rVal,INSERT_VALUES)); 4263 } 4264 } 4265 } 4266 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN)); 4267 PetscCall(DMRestoreWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork)); 4268 PetscCall(DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices)); 4269 } 4270 PetscCall(PetscFree(leafValues)); 4271 PetscCall(PetscSectionDestroy(&leafValuesSec)); 4272 PetscCall(PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO)); 4273 PetscCall(ISRestoreIndices(aIS,&anchors)); 4274 PetscFunctionReturn(0); 4275 } 4276 4277 static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids) 4278 { 4279 DM refTree; 4280 PetscSection multiRootSec, rootIndicesSec; 4281 PetscSection globalCoarse, globalFine; 4282 PetscSection localCoarse, localFine; 4283 PetscSection cSecRef; 4284 PetscInt *parentIndices, pRefStart, pRefEnd; 4285 PetscScalar *rootValues, *parentValues; 4286 Mat injRef; 4287 PetscInt numFields, maxDof; 4288 PetscInt pStartC, pEndC, pStartF, pEndF, p; 4289 PetscInt *offsets, *offsetsCopy, *rowOffsets; 4290 PetscLayout rowMap, colMap; 4291 PetscInt rowStart, rowEnd, colStart, colEnd; 4292 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 4293 4294 PetscFunctionBegin; 4295 4296 /* get the templates for the fine-to-coarse injection from the reference tree */ 4297 PetscCall(VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE)); 4298 PetscCall(VecSetOption(vecCoarse,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE)); 4299 PetscCall(DMPlexGetReferenceTree(coarse,&refTree)); 4300 PetscCall(DMCopyDisc(coarse,refTree)); 4301 PetscCall(DMGetDefaultConstraints(refTree,&cSecRef,NULL,NULL)); 4302 PetscCall(PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd)); 4303 PetscCall(DMPlexReferenceTreeGetInjector(refTree,&injRef)); 4304 4305 PetscCall(DMPlexGetChart(fine,&pStartF,&pEndF)); 4306 PetscCall(DMGetLocalSection(fine,&localFine)); 4307 PetscCall(DMGetGlobalSection(fine,&globalFine)); 4308 PetscCall(PetscSectionGetNumFields(localFine,&numFields)); 4309 PetscCall(DMPlexGetChart(coarse,&pStartC,&pEndC)); 4310 PetscCall(DMGetLocalSection(coarse,&localCoarse)); 4311 PetscCall(DMGetGlobalSection(coarse,&globalCoarse)); 4312 PetscCall(PetscSectionGetMaxDof(localCoarse,&maxDof)); 4313 { 4314 PetscInt maxFields = PetscMax(1,numFields) + 1; 4315 PetscCall(PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets)); 4316 } 4317 4318 PetscCall(DMPlexTransferInjectorTree(coarse,fine,coarseToFine,cids,vecFine,numFields,offsets,&multiRootSec,&rootIndicesSec,NULL,&rootValues)); 4319 4320 PetscCall(PetscMalloc2(maxDof,&parentIndices,maxDof,&parentValues)); 4321 4322 /* count indices */ 4323 PetscCall(VecGetLayout(vecFine,&colMap)); 4324 PetscCall(VecGetLayout(vecCoarse,&rowMap)); 4325 PetscCall(PetscLayoutSetUp(rowMap)); 4326 PetscCall(PetscLayoutSetUp(colMap)); 4327 PetscCall(PetscLayoutGetRange(rowMap,&rowStart,&rowEnd)); 4328 PetscCall(PetscLayoutGetRange(colMap,&colStart,&colEnd)); 4329 /* insert values */ 4330 PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats)); 4331 for (p = pStartC; p < pEndC; p++) { 4332 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 4333 PetscBool contribute = PETSC_FALSE; 4334 4335 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 4336 PetscCall(PetscSectionGetConstraintDof(globalCoarse,p,&cdof)); 4337 if ((dof - cdof) <= 0) continue; 4338 PetscCall(PetscSectionGetDof(localCoarse,p,&dof)); 4339 PetscCall(PetscSectionGetOffset(globalCoarse,p,&gOff)); 4340 4341 rowOffsets[0] = 0; 4342 offsetsCopy[0] = 0; 4343 if (numFields) { 4344 PetscInt f; 4345 4346 for (f = 0; f < numFields; f++) { 4347 PetscInt fDof; 4348 PetscCall(PetscSectionGetFieldDof(localCoarse,p,f,&fDof)); 4349 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 4350 } 4351 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,NULL,parentIndices)); 4352 } else { 4353 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,NULL,parentIndices)); 4354 rowOffsets[1] = offsetsCopy[0]; 4355 } 4356 4357 PetscCall(PetscSectionGetDof(multiRootSec,p,&numLeaves)); 4358 PetscCall(PetscSectionGetOffset(multiRootSec,p,&leafStart)); 4359 leafEnd = leafStart + numLeaves; 4360 for (l = 0; l < dof; l++) parentValues[l] = 0.; 4361 for (l = leafStart; l < leafEnd; l++) { 4362 PetscInt numIndices, childId, offset; 4363 const PetscScalar *childValues; 4364 4365 PetscCall(PetscSectionGetDof(rootIndicesSec,l,&numIndices)); 4366 PetscCall(PetscSectionGetOffset(rootIndicesSec,l,&offset)); 4367 childId = (PetscInt) PetscRealPart(rootValues[offset++]); 4368 childValues = &rootValues[offset]; 4369 numIndices--; 4370 4371 if (childId == -2) { /* skip */ 4372 continue; 4373 } else if (childId == -1) { /* equivalent points: scatter */ 4374 PetscInt m; 4375 4376 contribute = PETSC_TRUE; 4377 for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m]; 4378 } else { /* contributions from children: sum with injectors from reference tree */ 4379 PetscInt parentId, f, lim; 4380 4381 contribute = PETSC_TRUE; 4382 PetscCall(DMPlexGetTreeParent(refTree,childId,&parentId,NULL)); 4383 4384 lim = PetscMax(1,numFields); 4385 offsets[0] = 0; 4386 if (numFields) { 4387 PetscInt f; 4388 4389 for (f = 0; f < numFields; f++) { 4390 PetscInt fDof; 4391 PetscCall(PetscSectionGetFieldDof(cSecRef,childId,f,&fDof)); 4392 4393 offsets[f + 1] = fDof + offsets[f]; 4394 } 4395 } 4396 else { 4397 PetscInt cDof; 4398 4399 PetscCall(PetscSectionGetDof(cSecRef,childId,&cDof)); 4400 offsets[1] = cDof; 4401 } 4402 for (f = 0; f < lim; f++) { 4403 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 4404 PetscInt n = offsets[f+1]-offsets[f]; 4405 PetscInt m = rowOffsets[f+1]-rowOffsets[f]; 4406 PetscInt i, j; 4407 const PetscScalar *colValues = &childValues[offsets[f]]; 4408 4409 for (i = 0; i < m; i++) { 4410 PetscScalar val = 0.; 4411 for (j = 0; j < n; j++) { 4412 val += childMat[n * i + j] * colValues[j]; 4413 } 4414 parentValues[rowOffsets[f] + i] += val; 4415 } 4416 } 4417 } 4418 } 4419 if (contribute) PetscCall(VecSetValues(vecCoarse,dof,parentIndices,parentValues,INSERT_VALUES)); 4420 } 4421 PetscCall(PetscSectionDestroy(&multiRootSec)); 4422 PetscCall(PetscSectionDestroy(&rootIndicesSec)); 4423 PetscCall(PetscFree2(parentIndices,parentValues)); 4424 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats)); 4425 PetscCall(PetscFree(rootValues)); 4426 PetscCall(PetscFree3(offsets,offsetsCopy,rowOffsets)); 4427 PetscFunctionReturn(0); 4428 } 4429 4430 /*@ 4431 DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening 4432 that can be represented by a common reference tree used by both. This routine can be used for a combination of 4433 coarsening and refinement at the same time. 4434 4435 collective 4436 4437 Input Parameters: 4438 + dmIn - The DMPlex mesh for the input vector 4439 . vecIn - The input vector 4440 . sfRefine - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in 4441 the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn 4442 . sfCoarsen - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in 4443 the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn 4444 . cidsRefine - The childIds of the points in dmOut. These childIds relate back to the reference tree: childid[j] = k implies 4445 that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference 4446 tree was refined from its parent. childid[j] = -1 indicates that the point j in dmOut is exactly 4447 equivalent to its root in dmIn, so no interpolation is necessary. childid[j] = -2 indicates that this 4448 point j in dmOut is not a leaf of sfRefine. 4449 . cidsCoarsen - The childIds of the points in dmIn. These childIds relate back to the reference tree: childid[j] = k implies 4450 that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference 4451 tree coarsens to its parent. childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen. 4452 . useBCs - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer. 4453 - time - Used if boundary values are time dependent. 4454 4455 Output Parameters: 4456 . vecOut - Using interpolation and injection operators calculated on the reference tree, the transferred 4457 projection of vecIn from dmIn to dmOut. Note that any field discretized with a PetscFV finite volume 4458 method that uses gradient reconstruction will use reconstructed gradients when interpolating from 4459 coarse points to fine points. 4460 4461 Level: developer 4462 4463 .seealso: `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`, `PetscFVGetComputeGradients()` 4464 @*/ 4465 PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time) 4466 { 4467 PetscFunctionBegin; 4468 PetscCall(VecSet(vecOut,0.0)); 4469 if (sfRefine) { 4470 Vec vecInLocal; 4471 DM dmGrad = NULL; 4472 Vec faceGeom = NULL, cellGeom = NULL, grad = NULL; 4473 4474 PetscCall(DMGetLocalVector(dmIn,&vecInLocal)); 4475 PetscCall(VecSet(vecInLocal,0.0)); 4476 { 4477 PetscInt numFields, i; 4478 4479 PetscCall(DMGetNumFields(dmIn, &numFields)); 4480 for (i = 0; i < numFields; i++) { 4481 PetscObject obj; 4482 PetscClassId classid; 4483 4484 PetscCall(DMGetField(dmIn, i, NULL, &obj)); 4485 PetscCall(PetscObjectGetClassId(obj, &classid)); 4486 if (classid == PETSCFV_CLASSID) { 4487 PetscCall(DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad)); 4488 break; 4489 } 4490 } 4491 } 4492 if (useBCs) PetscCall(DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL)); 4493 PetscCall(DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal)); 4494 PetscCall(DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal)); 4495 if (dmGrad) { 4496 PetscCall(DMGetGlobalVector(dmGrad,&grad)); 4497 PetscCall(DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad)); 4498 } 4499 PetscCall(DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfRefine,cidsRefine,grad,cellGeom)); 4500 PetscCall(DMRestoreLocalVector(dmIn,&vecInLocal)); 4501 if (dmGrad) { 4502 PetscCall(DMRestoreGlobalVector(dmGrad,&grad)); 4503 } 4504 } 4505 if (sfCoarsen) PetscCall(DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfCoarsen,cidsCoarsen)); 4506 PetscCall(VecAssemblyBegin(vecOut)); 4507 PetscCall(VecAssemblyEnd(vecOut)); 4508 PetscFunctionReturn(0); 4509 } 4510