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 PetscCall(DMGetCoordinatesLocalSetUp(dm)); 1875 if (rank == 0) { 1876 /* compute the new charts */ 1877 PetscCall(PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd)); 1878 offset = 0; 1879 for (d = 0; d <= dim; d++) { 1880 PetscInt pOldCount, kStart, kEnd, k; 1881 1882 pNewStart[d] = offset; 1883 PetscCall(DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d])); 1884 PetscCall(DMPlexGetHeightStratum(K,d,&kStart,&kEnd)); 1885 pOldCount = pOldEnd[d] - pOldStart[d]; 1886 /* adding the new points */ 1887 pNewCount[d] = pOldCount + kEnd - kStart; 1888 if (!d) { 1889 /* removing the cell */ 1890 pNewCount[d]--; 1891 } 1892 for (k = kStart; k < kEnd; k++) { 1893 PetscInt parent; 1894 PetscCall(DMPlexGetTreeParent(K,k,&parent,NULL)); 1895 if (parent == k) { 1896 /* avoid double counting points that won't actually be new */ 1897 pNewCount[d]--; 1898 } 1899 } 1900 pNewEnd[d] = pNewStart[d] + pNewCount[d]; 1901 offset = pNewEnd[d]; 1902 1903 } 1904 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]); 1905 /* get the current closure of the cell that we are removing */ 1906 PetscCall(DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure)); 1907 1908 PetscCall(PetscMalloc1(pNewEnd[dim],&newConeSizes)); 1909 { 1910 DMPolytopeType pct, qct; 1911 PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j; 1912 1913 PetscCall(DMPlexGetChart(K,&kStart,&kEnd)); 1914 PetscCall(PetscMalloc4(kEnd-kStart,&Kembedding,kEnd-kStart,&perm,kEnd-kStart,&iperm,kEnd-kStart,&preOrient)); 1915 1916 for (k = kStart; k < kEnd; k++) { 1917 perm[k - kStart] = k; 1918 iperm [k - kStart] = k - kStart; 1919 preOrient[k - kStart] = 0; 1920 } 1921 1922 PetscCall(DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK)); 1923 for (j = 1; j < closureSizeK; j++) { 1924 PetscInt parentOrientA = closureK[2*j+1]; 1925 PetscInt parentOrientB = cellClosure[2*j+1]; 1926 PetscInt p, q; 1927 1928 p = closureK[2*j]; 1929 q = cellClosure[2*j]; 1930 PetscCall(DMPlexGetCellType(K, p, &pct)); 1931 PetscCall(DMPlexGetCellType(dm, q, &qct)); 1932 for (d = 0; d <= dim; d++) { 1933 if (q >= pOldStart[d] && q < pOldEnd[d]) { 1934 Kembedding[p] = (q - pOldStart[d]) + pNewStart[d]; 1935 } 1936 } 1937 parentOrientA = DMPolytopeConvertNewOrientation_Internal(pct, parentOrientA); 1938 parentOrientB = DMPolytopeConvertNewOrientation_Internal(qct, parentOrientB); 1939 if (parentOrientA != parentOrientB) { 1940 PetscInt numChildren, i; 1941 const PetscInt *children; 1942 1943 PetscCall(DMPlexGetTreeChildren(K,p,&numChildren,&children)); 1944 for (i = 0; i < numChildren; i++) { 1945 PetscInt kPerm, oPerm; 1946 1947 k = children[i]; 1948 PetscCall(DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm)); 1949 /* perm = what refTree position I'm in */ 1950 perm[kPerm-kStart] = k; 1951 /* iperm = who is at this position */ 1952 iperm[k-kStart] = kPerm-kStart; 1953 preOrient[kPerm-kStart] = oPerm; 1954 } 1955 } 1956 } 1957 PetscCall(DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK)); 1958 } 1959 PetscCall(PetscSectionSetChart(parentSection,0,pNewEnd[dim])); 1960 offset = 0; 1961 numNewCones = 0; 1962 for (d = 0; d <= dim; d++) { 1963 PetscInt kStart, kEnd, k; 1964 PetscInt p; 1965 PetscInt size; 1966 1967 for (p = pOldStart[d]; p < pOldEnd[d]; p++) { 1968 /* skip cell 0 */ 1969 if (p == cell) continue; 1970 /* old cones to new cones */ 1971 PetscCall(DMPlexGetConeSize(dm,p,&size)); 1972 newConeSizes[offset++] = size; 1973 numNewCones += size; 1974 } 1975 1976 PetscCall(DMPlexGetHeightStratum(K,d,&kStart,&kEnd)); 1977 for (k = kStart; k < kEnd; k++) { 1978 PetscInt kParent; 1979 1980 PetscCall(DMPlexGetTreeParent(K,k,&kParent,NULL)); 1981 if (kParent != k) { 1982 Kembedding[k] = offset; 1983 PetscCall(DMPlexGetConeSize(K,k,&size)); 1984 newConeSizes[offset++] = size; 1985 numNewCones += size; 1986 if (kParent != 0) { 1987 PetscCall(PetscSectionSetDof(parentSection,Kembedding[k],1)); 1988 } 1989 } 1990 } 1991 } 1992 1993 PetscCall(PetscSectionSetUp(parentSection)); 1994 PetscCall(PetscSectionGetStorageSize(parentSection,&numPointsWithParents)); 1995 PetscCall(PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations)); 1996 PetscCall(PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs)); 1997 1998 /* fill new cones */ 1999 offset = 0; 2000 for (d = 0; d <= dim; d++) { 2001 PetscInt kStart, kEnd, k, l; 2002 PetscInt p; 2003 PetscInt size; 2004 const PetscInt *cone, *orientation; 2005 2006 for (p = pOldStart[d]; p < pOldEnd[d]; p++) { 2007 /* skip cell 0 */ 2008 if (p == cell) continue; 2009 /* old cones to new cones */ 2010 PetscCall(DMPlexGetConeSize(dm,p,&size)); 2011 PetscCall(DMPlexGetCone(dm,p,&cone)); 2012 PetscCall(DMPlexGetConeOrientation(dm,p,&orientation)); 2013 for (l = 0; l < size; l++) { 2014 newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1]; 2015 newOrientations[offset++] = orientation[l]; 2016 } 2017 } 2018 2019 PetscCall(DMPlexGetHeightStratum(K,d,&kStart,&kEnd)); 2020 for (k = kStart; k < kEnd; k++) { 2021 PetscInt kPerm = perm[k], kParent; 2022 PetscInt preO = preOrient[k]; 2023 2024 PetscCall(DMPlexGetTreeParent(K,k,&kParent,NULL)); 2025 if (kParent != k) { 2026 /* embed new cones */ 2027 PetscCall(DMPlexGetConeSize(K,k,&size)); 2028 PetscCall(DMPlexGetCone(K,kPerm,&cone)); 2029 PetscCall(DMPlexGetConeOrientation(K,kPerm,&orientation)); 2030 for (l = 0; l < size; l++) { 2031 PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size); 2032 PetscInt newO, lSize, oTrue; 2033 DMPolytopeType ct = DM_NUM_POLYTOPES; 2034 2035 q = iperm[cone[m]]; 2036 newCones[offset] = Kembedding[q]; 2037 PetscCall(DMPlexGetConeSize(K,q,&lSize)); 2038 if (lSize == 2) ct = DM_POLYTOPE_SEGMENT; 2039 else if (lSize == 4) ct = DM_POLYTOPE_QUADRILATERAL; 2040 oTrue = DMPolytopeConvertNewOrientation_Internal(ct, orientation[m]); 2041 oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2); 2042 newO = DihedralCompose(lSize,oTrue,preOrient[q]); 2043 newOrientations[offset++] = DMPolytopeConvertOldOrientation_Internal(ct, newO); 2044 } 2045 if (kParent != 0) { 2046 PetscInt newPoint = Kembedding[kParent]; 2047 PetscCall(PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset)); 2048 parents[pOffset] = newPoint; 2049 childIDs[pOffset] = k; 2050 } 2051 } 2052 } 2053 } 2054 2055 PetscCall(PetscMalloc1(dim*(pNewEnd[dim]-pNewStart[dim]),&newVertexCoords)); 2056 2057 /* fill coordinates */ 2058 offset = 0; 2059 { 2060 PetscInt kStart, kEnd, l; 2061 PetscSection vSection; 2062 PetscInt v; 2063 Vec coords; 2064 PetscScalar *coordvals; 2065 PetscInt dof, off; 2066 PetscReal v0[3], J[9], detJ; 2067 2068 if (PetscDefined(USE_DEBUG)) { 2069 PetscInt k; 2070 PetscCall(DMPlexGetHeightStratum(K,0,&kStart,&kEnd)); 2071 for (k = kStart; k < kEnd; k++) { 2072 PetscCall(DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ)); 2073 PetscCheck(detJ > 0.,PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %" PetscInt_FMT " has bad determinant",k); 2074 } 2075 } 2076 PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ)); 2077 PetscCall(DMGetCoordinateSection(dm,&vSection)); 2078 PetscCall(DMGetCoordinatesLocal(dm,&coords)); 2079 PetscCall(VecGetArray(coords,&coordvals)); 2080 for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) { 2081 2082 PetscCall(PetscSectionGetDof(vSection,v,&dof)); 2083 PetscCall(PetscSectionGetOffset(vSection,v,&off)); 2084 for (l = 0; l < dof; l++) { 2085 newVertexCoords[offset++] = coordvals[off + l]; 2086 } 2087 } 2088 PetscCall(VecRestoreArray(coords,&coordvals)); 2089 2090 PetscCall(DMGetCoordinateSection(K,&vSection)); 2091 PetscCall(DMGetCoordinatesLocal(K,&coords)); 2092 PetscCall(VecGetArray(coords,&coordvals)); 2093 PetscCall(DMPlexGetDepthStratum(K,0,&kStart,&kEnd)); 2094 for (v = kStart; v < kEnd; v++) { 2095 PetscReal coord[3], newCoord[3]; 2096 PetscInt vPerm = perm[v]; 2097 PetscInt kParent; 2098 const PetscReal xi0[3] = {-1.,-1.,-1.}; 2099 2100 PetscCall(DMPlexGetTreeParent(K,v,&kParent,NULL)); 2101 if (kParent != v) { 2102 /* this is a new vertex */ 2103 PetscCall(PetscSectionGetOffset(vSection,vPerm,&off)); 2104 for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]); 2105 CoordinatesRefToReal(dim, dim, xi0, v0, J, coord, newCoord); 2106 for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l]; 2107 offset += dim; 2108 } 2109 } 2110 PetscCall(VecRestoreArray(coords,&coordvals)); 2111 } 2112 2113 /* need to reverse the order of pNewCount: vertices first, cells last */ 2114 for (d = 0; d < (dim + 1) / 2; d++) { 2115 PetscInt tmp; 2116 2117 tmp = pNewCount[d]; 2118 pNewCount[d] = pNewCount[dim - d]; 2119 pNewCount[dim - d] = tmp; 2120 } 2121 2122 PetscCall(DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords)); 2123 PetscCall(DMPlexSetReferenceTree(*ncdm,K)); 2124 PetscCall(DMPlexSetTree(*ncdm,parentSection,parents,childIDs)); 2125 2126 /* clean up */ 2127 PetscCall(DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure)); 2128 PetscCall(PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd)); 2129 PetscCall(PetscFree(newConeSizes)); 2130 PetscCall(PetscFree2(newCones,newOrientations)); 2131 PetscCall(PetscFree(newVertexCoords)); 2132 PetscCall(PetscFree2(parents,childIDs)); 2133 PetscCall(PetscFree4(Kembedding,perm,iperm,preOrient)); 2134 } 2135 else { 2136 PetscInt p, counts[4]; 2137 PetscInt *coneSizes, *cones, *orientations; 2138 Vec coordVec; 2139 PetscScalar *coords; 2140 2141 for (d = 0; d <= dim; d++) { 2142 PetscInt dStart, dEnd; 2143 2144 PetscCall(DMPlexGetDepthStratum(dm,d,&dStart,&dEnd)); 2145 counts[d] = dEnd - dStart; 2146 } 2147 PetscCall(PetscMalloc1(pEnd-pStart,&coneSizes)); 2148 for (p = pStart; p < pEnd; p++) { 2149 PetscCall(DMPlexGetConeSize(dm,p,&coneSizes[p-pStart])); 2150 } 2151 PetscCall(DMPlexGetCones(dm, &cones)); 2152 PetscCall(DMPlexGetConeOrientations(dm, &orientations)); 2153 PetscCall(DMGetCoordinatesLocal(dm,&coordVec)); 2154 PetscCall(VecGetArray(coordVec,&coords)); 2155 2156 PetscCall(PetscSectionSetChart(parentSection,pStart,pEnd)); 2157 PetscCall(PetscSectionSetUp(parentSection)); 2158 PetscCall(DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL)); 2159 PetscCall(DMPlexSetReferenceTree(*ncdm,K)); 2160 PetscCall(DMPlexSetTree(*ncdm,parentSection,NULL,NULL)); 2161 PetscCall(VecRestoreArray(coordVec,&coords)); 2162 } 2163 PetscCall(PetscSectionDestroy(&parentSection)); 2164 2165 PetscFunctionReturn(0); 2166 } 2167 2168 PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 2169 { 2170 PetscSF coarseToFineEmbedded; 2171 PetscSection globalCoarse, globalFine; 2172 PetscSection localCoarse, localFine; 2173 PetscSection aSec, cSec; 2174 PetscSection rootIndicesSec, rootMatricesSec; 2175 PetscSection leafIndicesSec, leafMatricesSec; 2176 PetscInt *rootIndices, *leafIndices; 2177 PetscScalar *rootMatrices, *leafMatrices; 2178 IS aIS; 2179 const PetscInt *anchors; 2180 Mat cMat; 2181 PetscInt numFields, maxFields; 2182 PetscInt pStartC, pEndC, pStartF, pEndF, p; 2183 PetscInt aStart, aEnd, cStart, cEnd; 2184 PetscInt *maxChildIds; 2185 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 2186 const PetscInt ***perms; 2187 const PetscScalar ***flips; 2188 2189 PetscFunctionBegin; 2190 PetscCall(DMPlexGetChart(coarse,&pStartC,&pEndC)); 2191 PetscCall(DMPlexGetChart(fine,&pStartF,&pEndF)); 2192 PetscCall(DMGetGlobalSection(fine,&globalFine)); 2193 { /* winnow fine points that don't have global dofs out of the sf */ 2194 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l; 2195 const PetscInt *leaves; 2196 2197 PetscCall(PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL)); 2198 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 2199 p = leaves ? leaves[l] : l; 2200 PetscCall(PetscSectionGetDof(globalFine,p,&dof)); 2201 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&cdof)); 2202 if ((dof - cdof) > 0) { 2203 numPointsWithDofs++; 2204 } 2205 } 2206 PetscCall(PetscMalloc1(numPointsWithDofs,&pointsWithDofs)); 2207 for (l = 0, offset = 0; l < nleaves; l++) { 2208 p = leaves ? leaves[l] : l; 2209 PetscCall(PetscSectionGetDof(globalFine,p,&dof)); 2210 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&cdof)); 2211 if ((dof - cdof) > 0) { 2212 pointsWithDofs[offset++] = l; 2213 } 2214 } 2215 PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded)); 2216 PetscCall(PetscFree(pointsWithDofs)); 2217 } 2218 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 2219 PetscCall(PetscMalloc1(pEndC-pStartC,&maxChildIds)); 2220 for (p = pStartC; p < pEndC; p++) { 2221 maxChildIds[p - pStartC] = -2; 2222 } 2223 PetscCall(PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX)); 2224 PetscCall(PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX)); 2225 2226 PetscCall(DMGetLocalSection(coarse,&localCoarse)); 2227 PetscCall(DMGetGlobalSection(coarse,&globalCoarse)); 2228 2229 PetscCall(DMPlexGetAnchors(coarse,&aSec,&aIS)); 2230 PetscCall(ISGetIndices(aIS,&anchors)); 2231 PetscCall(PetscSectionGetChart(aSec,&aStart,&aEnd)); 2232 2233 PetscCall(DMGetDefaultConstraints(coarse,&cSec,&cMat,NULL)); 2234 PetscCall(PetscSectionGetChart(cSec,&cStart,&cEnd)); 2235 2236 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 2237 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec)); 2238 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec)); 2239 PetscCall(PetscSectionSetChart(rootIndicesSec,pStartC,pEndC)); 2240 PetscCall(PetscSectionSetChart(rootMatricesSec,pStartC,pEndC)); 2241 PetscCall(PetscSectionGetNumFields(localCoarse,&numFields)); 2242 maxFields = PetscMax(1,numFields); 2243 PetscCall(PetscMalloc7(maxFields+1,&offsets,maxFields+1,&offsetsCopy,maxFields+1,&newOffsets,maxFields+1,&newOffsetsCopy,maxFields+1,&rowOffsets,maxFields+1,&numD,maxFields+1,&numO)); 2244 PetscCall(PetscMalloc2(maxFields+1,(PetscInt****)&perms,maxFields+1,(PetscScalar****)&flips)); 2245 PetscCall(PetscMemzero((void *) perms, (maxFields+1) * sizeof(const PetscInt **))); 2246 PetscCall(PetscMemzero((void *) flips, (maxFields+1) * sizeof(const PetscScalar **))); 2247 2248 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 2249 PetscInt dof, matSize = 0; 2250 PetscInt aDof = 0; 2251 PetscInt cDof = 0; 2252 PetscInt maxChildId = maxChildIds[p - pStartC]; 2253 PetscInt numRowIndices = 0; 2254 PetscInt numColIndices = 0; 2255 PetscInt f; 2256 2257 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 2258 if (dof < 0) { 2259 dof = -(dof + 1); 2260 } 2261 if (p >= aStart && p < aEnd) { 2262 PetscCall(PetscSectionGetDof(aSec,p,&aDof)); 2263 } 2264 if (p >= cStart && p < cEnd) { 2265 PetscCall(PetscSectionGetDof(cSec,p,&cDof)); 2266 } 2267 for (f = 0; f <= numFields; f++) offsets[f] = 0; 2268 for (f = 0; f <= numFields; f++) newOffsets[f] = 0; 2269 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 2270 PetscInt *closure = NULL, closureSize, cl; 2271 2272 PetscCall(DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure)); 2273 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 2274 PetscInt c = closure[2 * cl], clDof; 2275 2276 PetscCall(PetscSectionGetDof(localCoarse,c,&clDof)); 2277 numRowIndices += clDof; 2278 for (f = 0; f < numFields; f++) { 2279 PetscCall(PetscSectionGetFieldDof(localCoarse,c,f,&clDof)); 2280 offsets[f + 1] += clDof; 2281 } 2282 } 2283 for (f = 0; f < numFields; f++) { 2284 offsets[f + 1] += offsets[f]; 2285 newOffsets[f + 1] = offsets[f + 1]; 2286 } 2287 /* get the number of indices needed and their field offsets */ 2288 PetscCall(DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE)); 2289 PetscCall(DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure)); 2290 if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */ 2291 numColIndices = numRowIndices; 2292 matSize = 0; 2293 } 2294 else if (numFields) { /* we send one submat for each field: sum their sizes */ 2295 matSize = 0; 2296 for (f = 0; f < numFields; f++) { 2297 PetscInt numRow, numCol; 2298 2299 numRow = offsets[f + 1] - offsets[f]; 2300 numCol = newOffsets[f + 1] - newOffsets[f]; 2301 matSize += numRow * numCol; 2302 } 2303 } 2304 else { 2305 matSize = numRowIndices * numColIndices; 2306 } 2307 } else if (maxChildId == -1) { 2308 if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */ 2309 PetscInt aOff, a; 2310 2311 PetscCall(PetscSectionGetOffset(aSec,p,&aOff)); 2312 for (f = 0; f < numFields; f++) { 2313 PetscInt fDof; 2314 2315 PetscCall(PetscSectionGetFieldDof(localCoarse,p,f,&fDof)); 2316 offsets[f+1] = fDof; 2317 } 2318 for (a = 0; a < aDof; a++) { 2319 PetscInt anchor = anchors[a + aOff], aLocalDof; 2320 2321 PetscCall(PetscSectionGetDof(localCoarse,anchor,&aLocalDof)); 2322 numColIndices += aLocalDof; 2323 for (f = 0; f < numFields; f++) { 2324 PetscInt fDof; 2325 2326 PetscCall(PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof)); 2327 newOffsets[f+1] += fDof; 2328 } 2329 } 2330 if (numFields) { 2331 matSize = 0; 2332 for (f = 0; f < numFields; f++) { 2333 matSize += offsets[f+1] * newOffsets[f+1]; 2334 } 2335 } 2336 else { 2337 matSize = numColIndices * dof; 2338 } 2339 } 2340 else { /* no children, and no constraints on dofs: just get the global indices */ 2341 numColIndices = dof; 2342 matSize = 0; 2343 } 2344 } 2345 /* we will pack the column indices with the field offsets */ 2346 PetscCall(PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0)); 2347 PetscCall(PetscSectionSetDof(rootMatricesSec,p,matSize)); 2348 } 2349 PetscCall(PetscSectionSetUp(rootIndicesSec)); 2350 PetscCall(PetscSectionSetUp(rootMatricesSec)); 2351 { 2352 PetscInt numRootIndices, numRootMatrices; 2353 2354 PetscCall(PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices)); 2355 PetscCall(PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices)); 2356 PetscCall(PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices)); 2357 for (p = pStartC; p < pEndC; p++) { 2358 PetscInt numRowIndices, numColIndices, matSize, dof; 2359 PetscInt pIndOff, pMatOff, f; 2360 PetscInt *pInd; 2361 PetscInt maxChildId = maxChildIds[p - pStartC]; 2362 PetscScalar *pMat = NULL; 2363 2364 PetscCall(PetscSectionGetDof(rootIndicesSec,p,&numColIndices)); 2365 if (!numColIndices) { 2366 continue; 2367 } 2368 for (f = 0; f <= numFields; f++) { 2369 offsets[f] = 0; 2370 newOffsets[f] = 0; 2371 offsetsCopy[f] = 0; 2372 newOffsetsCopy[f] = 0; 2373 } 2374 numColIndices -= 2 * numFields; 2375 PetscCall(PetscSectionGetOffset(rootIndicesSec,p,&pIndOff)); 2376 pInd = &(rootIndices[pIndOff]); 2377 PetscCall(PetscSectionGetDof(rootMatricesSec,p,&matSize)); 2378 if (matSize) { 2379 PetscCall(PetscSectionGetOffset(rootMatricesSec,p,&pMatOff)); 2380 pMat = &rootMatrices[pMatOff]; 2381 } 2382 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 2383 if (dof < 0) { 2384 dof = -(dof + 1); 2385 } 2386 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 2387 PetscInt i, j; 2388 PetscInt numRowIndices = matSize / numColIndices; 2389 2390 if (!numRowIndices) { /* don't need to calculate the mat, just the indices */ 2391 PetscInt numIndices, *indices; 2392 PetscCall(DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,PETSC_TRUE,&numIndices,&indices,offsets,NULL)); 2393 PetscCheck(numIndices == numColIndices,PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations"); 2394 for (i = 0; i < numColIndices; i++) { 2395 pInd[i] = indices[i]; 2396 } 2397 for (i = 0; i < numFields; i++) { 2398 pInd[numColIndices + i] = offsets[i+1]; 2399 pInd[numColIndices + numFields + i] = offsets[i+1]; 2400 } 2401 PetscCall(DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,PETSC_TRUE,&numIndices,&indices,offsets,NULL)); 2402 } 2403 else { 2404 PetscInt closureSize, *closure = NULL, cl; 2405 PetscScalar *pMatIn, *pMatModified; 2406 PetscInt numPoints,*points; 2407 2408 PetscCall(DMGetWorkArray(coarse,numRowIndices * numRowIndices,MPIU_SCALAR,&pMatIn)); 2409 for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */ 2410 for (j = 0; j < numRowIndices; j++) { 2411 pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.; 2412 } 2413 } 2414 PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure)); 2415 for (f = 0; f < maxFields; f++) { 2416 if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f])); 2417 else PetscCall(PetscSectionGetPointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f])); 2418 } 2419 if (numFields) { 2420 for (cl = 0; cl < closureSize; cl++) { 2421 PetscInt c = closure[2 * cl]; 2422 2423 for (f = 0; f < numFields; f++) { 2424 PetscInt fDof; 2425 2426 PetscCall(PetscSectionGetFieldDof(localCoarse,c,f,&fDof)); 2427 offsets[f + 1] += fDof; 2428 } 2429 } 2430 for (f = 0; f < numFields; f++) { 2431 offsets[f + 1] += offsets[f]; 2432 newOffsets[f + 1] = offsets[f + 1]; 2433 } 2434 } 2435 /* TODO : flips here ? */ 2436 /* apply hanging node constraints on the right, get the new points and the new offsets */ 2437 PetscCall(DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,perms,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE)); 2438 for (f = 0; f < maxFields; f++) { 2439 if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f])); 2440 else PetscCall(PetscSectionRestorePointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f])); 2441 } 2442 for (f = 0; f < maxFields; f++) { 2443 if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f])); 2444 else PetscCall(PetscSectionGetPointSyms(localCoarse,numPoints,points,&perms[f],&flips[f])); 2445 } 2446 if (!numFields) { 2447 for (i = 0; i < numRowIndices * numColIndices; i++) { 2448 pMat[i] = pMatModified[i]; 2449 } 2450 } 2451 else { 2452 PetscInt i, j, count; 2453 for (f = 0, count = 0; f < numFields; f++) { 2454 for (i = offsets[f]; i < offsets[f+1]; i++) { 2455 for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) { 2456 pMat[count] = pMatModified[i * numColIndices + j]; 2457 } 2458 } 2459 } 2460 } 2461 PetscCall(DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatModified)); 2462 PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure)); 2463 PetscCall(DMRestoreWorkArray(coarse,numRowIndices * numColIndices,MPIU_SCALAR,&pMatIn)); 2464 if (numFields) { 2465 for (f = 0; f < numFields; f++) { 2466 pInd[numColIndices + f] = offsets[f+1]; 2467 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2468 } 2469 for (cl = 0; cl < numPoints; cl++) { 2470 PetscInt globalOff, c = points[2*cl]; 2471 PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff)); 2472 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, NULL, pInd)); 2473 } 2474 } else { 2475 for (cl = 0; cl < numPoints; cl++) { 2476 PetscInt c = points[2*cl], globalOff; 2477 const PetscInt *perm = perms[0] ? perms[0][cl] : NULL; 2478 2479 PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff)); 2480 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perm, NULL, pInd)); 2481 } 2482 } 2483 for (f = 0; f < maxFields; f++) { 2484 if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f])); 2485 else PetscCall(PetscSectionRestorePointSyms(localCoarse,numPoints,points,&perms[f],&flips[f])); 2486 } 2487 PetscCall(DMRestoreWorkArray(coarse,numPoints,MPIU_SCALAR,&points)); 2488 } 2489 } 2490 else if (matSize) { 2491 PetscInt cOff; 2492 PetscInt *rowIndices, *colIndices, a, aDof, aOff; 2493 2494 numRowIndices = matSize / numColIndices; 2495 PetscCheck(numRowIndices == dof,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Miscounted dofs"); 2496 PetscCall(DMGetWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices)); 2497 PetscCall(DMGetWorkArray(coarse,numColIndices,MPIU_INT,&colIndices)); 2498 PetscCall(PetscSectionGetOffset(cSec,p,&cOff)); 2499 PetscCall(PetscSectionGetDof(aSec,p,&aDof)); 2500 PetscCall(PetscSectionGetOffset(aSec,p,&aOff)); 2501 if (numFields) { 2502 for (f = 0; f < numFields; f++) { 2503 PetscInt fDof; 2504 2505 PetscCall(PetscSectionGetFieldDof(cSec,p,f,&fDof)); 2506 offsets[f + 1] = fDof; 2507 for (a = 0; a < aDof; a++) { 2508 PetscInt anchor = anchors[a + aOff]; 2509 PetscCall(PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof)); 2510 newOffsets[f + 1] += fDof; 2511 } 2512 } 2513 for (f = 0; f < numFields; f++) { 2514 offsets[f + 1] += offsets[f]; 2515 offsetsCopy[f + 1] = offsets[f + 1]; 2516 newOffsets[f + 1] += newOffsets[f]; 2517 newOffsetsCopy[f + 1] = newOffsets[f + 1]; 2518 } 2519 PetscCall(DMPlexGetIndicesPointFields_Internal(cSec,PETSC_TRUE,p,cOff,offsetsCopy,PETSC_TRUE,NULL,-1, NULL,rowIndices)); 2520 for (a = 0; a < aDof; a++) { 2521 PetscInt anchor = anchors[a + aOff], lOff; 2522 PetscCall(PetscSectionGetOffset(localCoarse,anchor,&lOff)); 2523 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_TRUE,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,-1, NULL,colIndices)); 2524 } 2525 } 2526 else { 2527 PetscCall(DMPlexGetIndicesPoint_Internal(cSec,PETSC_TRUE,p,cOff,offsetsCopy,PETSC_TRUE,NULL, NULL,rowIndices)); 2528 for (a = 0; a < aDof; a++) { 2529 PetscInt anchor = anchors[a + aOff], lOff; 2530 PetscCall(PetscSectionGetOffset(localCoarse,anchor,&lOff)); 2531 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_TRUE,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL, NULL,colIndices)); 2532 } 2533 } 2534 if (numFields) { 2535 PetscInt count, a; 2536 2537 for (f = 0, count = 0; f < numFields; f++) { 2538 PetscInt iSize = offsets[f + 1] - offsets[f]; 2539 PetscInt jSize = newOffsets[f + 1] - newOffsets[f]; 2540 PetscCall(MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count])); 2541 count += iSize * jSize; 2542 pInd[numColIndices + f] = offsets[f+1]; 2543 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2544 } 2545 for (a = 0; a < aDof; a++) { 2546 PetscInt anchor = anchors[a + aOff]; 2547 PetscInt gOff; 2548 PetscCall(PetscSectionGetOffset(globalCoarse,anchor,&gOff)); 2549 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,-1, NULL,pInd)); 2550 } 2551 } 2552 else { 2553 PetscInt a; 2554 PetscCall(MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat)); 2555 for (a = 0; a < aDof; a++) { 2556 PetscInt anchor = anchors[a + aOff]; 2557 PetscInt gOff; 2558 PetscCall(PetscSectionGetOffset(globalCoarse,anchor,&gOff)); 2559 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL, NULL,pInd)); 2560 } 2561 } 2562 PetscCall(DMRestoreWorkArray(coarse,numColIndices,MPIU_INT,&colIndices)); 2563 PetscCall(DMRestoreWorkArray(coarse,numRowIndices,MPIU_INT,&rowIndices)); 2564 } 2565 else { 2566 PetscInt gOff; 2567 2568 PetscCall(PetscSectionGetOffset(globalCoarse,p,&gOff)); 2569 if (numFields) { 2570 for (f = 0; f < numFields; f++) { 2571 PetscInt fDof; 2572 PetscCall(PetscSectionGetFieldDof(localCoarse,p,f,&fDof)); 2573 offsets[f + 1] = fDof + offsets[f]; 2574 } 2575 for (f = 0; f < numFields; f++) { 2576 pInd[numColIndices + f] = offsets[f+1]; 2577 pInd[numColIndices + numFields + f] = offsets[f+1]; 2578 } 2579 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1, NULL,pInd)); 2580 } else { 2581 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL, NULL,pInd)); 2582 } 2583 } 2584 } 2585 PetscCall(PetscFree(maxChildIds)); 2586 } 2587 { 2588 PetscSF indicesSF, matricesSF; 2589 PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices; 2590 2591 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec)); 2592 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafMatricesSec)); 2593 PetscCall(PetscSFDistributeSection(coarseToFineEmbedded,rootIndicesSec,&remoteOffsetsIndices,leafIndicesSec)); 2594 PetscCall(PetscSFDistributeSection(coarseToFineEmbedded,rootMatricesSec,&remoteOffsetsMatrices,leafMatricesSec)); 2595 PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded,rootIndicesSec,remoteOffsetsIndices,leafIndicesSec,&indicesSF)); 2596 PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded,rootMatricesSec,remoteOffsetsMatrices,leafMatricesSec,&matricesSF)); 2597 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 2598 PetscCall(PetscFree(remoteOffsetsIndices)); 2599 PetscCall(PetscFree(remoteOffsetsMatrices)); 2600 PetscCall(PetscSectionGetStorageSize(leafIndicesSec,&numLeafIndices)); 2601 PetscCall(PetscSectionGetStorageSize(leafMatricesSec,&numLeafMatrices)); 2602 PetscCall(PetscMalloc2(numLeafIndices,&leafIndices,numLeafMatrices,&leafMatrices)); 2603 PetscCall(PetscSFBcastBegin(indicesSF,MPIU_INT,rootIndices,leafIndices,MPI_REPLACE)); 2604 PetscCall(PetscSFBcastBegin(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices,MPI_REPLACE)); 2605 PetscCall(PetscSFBcastEnd(indicesSF,MPIU_INT,rootIndices,leafIndices,MPI_REPLACE)); 2606 PetscCall(PetscSFBcastEnd(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices,MPI_REPLACE)); 2607 PetscCall(PetscSFDestroy(&matricesSF)); 2608 PetscCall(PetscSFDestroy(&indicesSF)); 2609 PetscCall(PetscFree2(rootIndices,rootMatrices)); 2610 PetscCall(PetscSectionDestroy(&rootIndicesSec)); 2611 PetscCall(PetscSectionDestroy(&rootMatricesSec)); 2612 } 2613 /* count to preallocate */ 2614 PetscCall(DMGetLocalSection(fine,&localFine)); 2615 { 2616 PetscInt nGlobal; 2617 PetscInt *dnnz, *onnz; 2618 PetscLayout rowMap, colMap; 2619 PetscInt rowStart, rowEnd, colStart, colEnd; 2620 PetscInt maxDof; 2621 PetscInt *rowIndices; 2622 DM refTree; 2623 PetscInt **refPointFieldN; 2624 PetscScalar ***refPointFieldMats; 2625 PetscSection refConSec, refAnSec; 2626 PetscInt pRefStart,pRefEnd,maxConDof,maxColumns,leafStart,leafEnd; 2627 PetscScalar *pointWork; 2628 2629 PetscCall(PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal)); 2630 PetscCall(PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz)); 2631 PetscCall(MatGetLayouts(mat,&rowMap,&colMap)); 2632 PetscCall(PetscLayoutSetUp(rowMap)); 2633 PetscCall(PetscLayoutSetUp(colMap)); 2634 PetscCall(PetscLayoutGetRange(rowMap,&rowStart,&rowEnd)); 2635 PetscCall(PetscLayoutGetRange(colMap,&colStart,&colEnd)); 2636 PetscCall(PetscSectionGetMaxDof(localFine,&maxDof)); 2637 PetscCall(PetscSectionGetChart(leafIndicesSec,&leafStart,&leafEnd)); 2638 PetscCall(DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices)); 2639 for (p = leafStart; p < leafEnd; p++) { 2640 PetscInt gDof, gcDof, gOff; 2641 PetscInt numColIndices, pIndOff, *pInd; 2642 PetscInt matSize; 2643 PetscInt i; 2644 2645 PetscCall(PetscSectionGetDof(globalFine,p,&gDof)); 2646 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&gcDof)); 2647 if ((gDof - gcDof) <= 0) { 2648 continue; 2649 } 2650 PetscCall(PetscSectionGetOffset(globalFine,p,&gOff)); 2651 PetscCheck(gOff >= 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"I though having global dofs meant a non-negative offset"); 2652 PetscCheck(gOff >= rowStart && (gOff + gDof - gcDof) <= rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"I thought the row map would constrain the global dofs"); 2653 PetscCall(PetscSectionGetDof(leafIndicesSec,p,&numColIndices)); 2654 PetscCall(PetscSectionGetOffset(leafIndicesSec,p,&pIndOff)); 2655 numColIndices -= 2 * numFields; 2656 PetscCheck(numColIndices > 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"global fine dof with no dofs to interpolate from"); 2657 pInd = &leafIndices[pIndOff]; 2658 offsets[0] = 0; 2659 offsetsCopy[0] = 0; 2660 newOffsets[0] = 0; 2661 newOffsetsCopy[0] = 0; 2662 if (numFields) { 2663 PetscInt f; 2664 for (f = 0; f < numFields; f++) { 2665 PetscInt rowDof; 2666 2667 PetscCall(PetscSectionGetFieldDof(localFine,p,f,&rowDof)); 2668 offsets[f + 1] = offsets[f] + rowDof; 2669 offsetsCopy[f + 1] = offsets[f + 1]; 2670 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2671 numD[f] = 0; 2672 numO[f] = 0; 2673 } 2674 PetscCall(DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,rowIndices)); 2675 for (f = 0; f < numFields; f++) { 2676 PetscInt colOffset = newOffsets[f]; 2677 PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f]; 2678 2679 for (i = 0; i < numFieldCols; i++) { 2680 PetscInt gInd = pInd[i + colOffset]; 2681 2682 if (gInd >= colStart && gInd < colEnd) { 2683 numD[f]++; 2684 } 2685 else if (gInd >= 0) { /* negative means non-entry */ 2686 numO[f]++; 2687 } 2688 } 2689 } 2690 } 2691 else { 2692 PetscCall(DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,rowIndices)); 2693 numD[0] = 0; 2694 numO[0] = 0; 2695 for (i = 0; i < numColIndices; i++) { 2696 PetscInt gInd = pInd[i]; 2697 2698 if (gInd >= colStart && gInd < colEnd) { 2699 numD[0]++; 2700 } 2701 else if (gInd >= 0) { /* negative means non-entry */ 2702 numO[0]++; 2703 } 2704 } 2705 } 2706 PetscCall(PetscSectionGetDof(leafMatricesSec,p,&matSize)); 2707 if (!matSize) { /* incoming matrix is identity */ 2708 PetscInt childId; 2709 2710 childId = childIds[p-pStartF]; 2711 if (childId < 0) { /* no child interpolation: one nnz per */ 2712 if (numFields) { 2713 PetscInt f; 2714 for (f = 0; f < numFields; f++) { 2715 PetscInt numRows = offsets[f+1] - offsets[f], row; 2716 for (row = 0; row < numRows; row++) { 2717 PetscInt gIndCoarse = pInd[newOffsets[f] + row]; 2718 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2719 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2720 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2721 dnnz[gIndFine - rowStart] = 1; 2722 } 2723 else if (gIndCoarse >= 0) { /* remote */ 2724 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2725 onnz[gIndFine - rowStart] = 1; 2726 } 2727 else { /* constrained */ 2728 PetscCheck(gIndFine < 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2729 } 2730 } 2731 } 2732 } 2733 else { 2734 PetscInt i; 2735 for (i = 0; i < gDof; i++) { 2736 PetscInt gIndCoarse = pInd[i]; 2737 PetscInt gIndFine = rowIndices[i]; 2738 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2739 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2740 dnnz[gIndFine - rowStart] = 1; 2741 } 2742 else if (gIndCoarse >= 0) { /* remote */ 2743 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2744 onnz[gIndFine - rowStart] = 1; 2745 } 2746 else { /* constrained */ 2747 PetscCheck(gIndFine < 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2748 } 2749 } 2750 } 2751 } 2752 else { /* interpolate from all */ 2753 if (numFields) { 2754 PetscInt f; 2755 for (f = 0; f < numFields; f++) { 2756 PetscInt numRows = offsets[f+1] - offsets[f], row; 2757 for (row = 0; row < numRows; row++) { 2758 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2759 if (gIndFine >= 0) { 2760 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2761 dnnz[gIndFine - rowStart] = numD[f]; 2762 onnz[gIndFine - rowStart] = numO[f]; 2763 } 2764 } 2765 } 2766 } 2767 else { 2768 PetscInt i; 2769 for (i = 0; i < gDof; i++) { 2770 PetscInt gIndFine = rowIndices[i]; 2771 if (gIndFine >= 0) { 2772 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2773 dnnz[gIndFine - rowStart] = numD[0]; 2774 onnz[gIndFine - rowStart] = numO[0]; 2775 } 2776 } 2777 } 2778 } 2779 } 2780 else { /* interpolate from all */ 2781 if (numFields) { 2782 PetscInt f; 2783 for (f = 0; f < numFields; f++) { 2784 PetscInt numRows = offsets[f+1] - offsets[f], row; 2785 for (row = 0; row < numRows; row++) { 2786 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2787 if (gIndFine >= 0) { 2788 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2789 dnnz[gIndFine - rowStart] = numD[f]; 2790 onnz[gIndFine - rowStart] = numO[f]; 2791 } 2792 } 2793 } 2794 } 2795 else { /* every dof get a full row */ 2796 PetscInt i; 2797 for (i = 0; i < gDof; i++) { 2798 PetscInt gIndFine = rowIndices[i]; 2799 if (gIndFine >= 0) { 2800 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2801 dnnz[gIndFine - rowStart] = numD[0]; 2802 onnz[gIndFine - rowStart] = numO[0]; 2803 } 2804 } 2805 } 2806 } 2807 } 2808 PetscCall(MatXAIJSetPreallocation(mat,1,dnnz,onnz,NULL,NULL)); 2809 PetscCall(PetscFree2(dnnz,onnz)); 2810 2811 PetscCall(DMPlexGetReferenceTree(fine,&refTree)); 2812 PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN)); 2813 PetscCall(DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL)); 2814 PetscCall(DMPlexGetAnchors(refTree,&refAnSec,NULL)); 2815 PetscCall(PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd)); 2816 PetscCall(PetscSectionGetMaxDof(refConSec,&maxConDof)); 2817 PetscCall(PetscSectionGetMaxDof(leafIndicesSec,&maxColumns)); 2818 PetscCall(PetscMalloc1(maxConDof*maxColumns,&pointWork)); 2819 for (p = leafStart; p < leafEnd; p++) { 2820 PetscInt gDof, gcDof, gOff; 2821 PetscInt numColIndices, pIndOff, *pInd; 2822 PetscInt matSize; 2823 PetscInt childId; 2824 2825 PetscCall(PetscSectionGetDof(globalFine,p,&gDof)); 2826 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&gcDof)); 2827 if ((gDof - gcDof) <= 0) { 2828 continue; 2829 } 2830 childId = childIds[p-pStartF]; 2831 PetscCall(PetscSectionGetOffset(globalFine,p,&gOff)); 2832 PetscCall(PetscSectionGetDof(leafIndicesSec,p,&numColIndices)); 2833 PetscCall(PetscSectionGetOffset(leafIndicesSec,p,&pIndOff)); 2834 numColIndices -= 2 * numFields; 2835 pInd = &leafIndices[pIndOff]; 2836 offsets[0] = 0; 2837 offsetsCopy[0] = 0; 2838 newOffsets[0] = 0; 2839 newOffsetsCopy[0] = 0; 2840 rowOffsets[0] = 0; 2841 if (numFields) { 2842 PetscInt f; 2843 for (f = 0; f < numFields; f++) { 2844 PetscInt rowDof; 2845 2846 PetscCall(PetscSectionGetFieldDof(localFine,p,f,&rowDof)); 2847 offsets[f + 1] = offsets[f] + rowDof; 2848 offsetsCopy[f + 1] = offsets[f + 1]; 2849 rowOffsets[f + 1] = pInd[numColIndices + f]; 2850 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2851 } 2852 PetscCall(DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,rowIndices)); 2853 } 2854 else { 2855 PetscCall(DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,rowIndices)); 2856 } 2857 PetscCall(PetscSectionGetDof(leafMatricesSec,p,&matSize)); 2858 if (!matSize) { /* incoming matrix is identity */ 2859 if (childId < 0) { /* no child interpolation: scatter */ 2860 if (numFields) { 2861 PetscInt f; 2862 for (f = 0; f < numFields; f++) { 2863 PetscInt numRows = offsets[f+1] - offsets[f], row; 2864 for (row = 0; row < numRows; row++) { 2865 PetscCall(MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES)); 2866 } 2867 } 2868 } 2869 else { 2870 PetscInt numRows = gDof, row; 2871 for (row = 0; row < numRows; row++) { 2872 PetscCall(MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES)); 2873 } 2874 } 2875 } 2876 else { /* interpolate from all */ 2877 if (numFields) { 2878 PetscInt f; 2879 for (f = 0; f < numFields; f++) { 2880 PetscInt numRows = offsets[f+1] - offsets[f]; 2881 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 2882 PetscCall(MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES)); 2883 } 2884 } 2885 else { 2886 PetscCall(MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES)); 2887 } 2888 } 2889 } 2890 else { /* interpolate from all */ 2891 PetscInt pMatOff; 2892 PetscScalar *pMat; 2893 2894 PetscCall(PetscSectionGetOffset(leafMatricesSec,p,&pMatOff)); 2895 pMat = &leafMatrices[pMatOff]; 2896 if (childId < 0) { /* copy the incoming matrix */ 2897 if (numFields) { 2898 PetscInt f, count; 2899 for (f = 0, count = 0; f < numFields; f++) { 2900 PetscInt numRows = offsets[f+1]-offsets[f]; 2901 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2902 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2903 PetscScalar *inMat = &pMat[count]; 2904 2905 PetscCall(MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES)); 2906 count += numCols * numInRows; 2907 } 2908 } 2909 else { 2910 PetscCall(MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES)); 2911 } 2912 } 2913 else { /* multiply the incoming matrix by the child interpolation */ 2914 if (numFields) { 2915 PetscInt f, count; 2916 for (f = 0, count = 0; f < numFields; f++) { 2917 PetscInt numRows = offsets[f+1]-offsets[f]; 2918 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2919 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2920 PetscScalar *inMat = &pMat[count]; 2921 PetscInt i, j, k; 2922 PetscCheck(refPointFieldN[childId - pRefStart][f] == numInRows,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2923 for (i = 0; i < numRows; i++) { 2924 for (j = 0; j < numCols; j++) { 2925 PetscScalar val = 0.; 2926 for (k = 0; k < numInRows; k++) { 2927 val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j]; 2928 } 2929 pointWork[i * numCols + j] = val; 2930 } 2931 } 2932 PetscCall(MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES)); 2933 count += numCols * numInRows; 2934 } 2935 } 2936 else { /* every dof gets a full row */ 2937 PetscInt numRows = gDof; 2938 PetscInt numCols = numColIndices; 2939 PetscInt numInRows = matSize / numColIndices; 2940 PetscInt i, j, k; 2941 PetscCheck(refPointFieldN[childId - pRefStart][0] == numInRows,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2942 for (i = 0; i < numRows; i++) { 2943 for (j = 0; j < numCols; j++) { 2944 PetscScalar val = 0.; 2945 for (k = 0; k < numInRows; k++) { 2946 val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j]; 2947 } 2948 pointWork[i * numCols + j] = val; 2949 } 2950 } 2951 PetscCall(MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES)); 2952 } 2953 } 2954 } 2955 } 2956 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN)); 2957 PetscCall(DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices)); 2958 PetscCall(PetscFree(pointWork)); 2959 } 2960 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 2961 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 2962 PetscCall(PetscSectionDestroy(&leafIndicesSec)); 2963 PetscCall(PetscSectionDestroy(&leafMatricesSec)); 2964 PetscCall(PetscFree2(leafIndices,leafMatrices)); 2965 PetscCall(PetscFree2(*(PetscInt****)&perms,*(PetscScalar****)&flips)); 2966 PetscCall(PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO)); 2967 PetscCall(ISRestoreIndices(aIS,&anchors)); 2968 PetscFunctionReturn(0); 2969 } 2970 2971 /* 2972 * Assuming a nodal basis (w.r.t. the dual basis) basis: 2973 * 2974 * for each coarse dof \phi^c_i: 2975 * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i: 2976 * for each fine dof \phi^f_j; 2977 * a_{i,j} = 0; 2978 * for each fine dof \phi^f_k: 2979 * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l 2980 * [^^^ this is = \phi^c_i ^^^] 2981 */ 2982 PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj) 2983 { 2984 PetscDS ds; 2985 PetscSection section, cSection; 2986 DMLabel canonical, depth; 2987 Mat cMat, mat; 2988 PetscInt *nnz; 2989 PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd; 2990 PetscInt m, n; 2991 PetscScalar *pointScalar; 2992 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent; 2993 2994 PetscFunctionBegin; 2995 PetscCall(DMGetLocalSection(refTree,§ion)); 2996 PetscCall(DMGetDimension(refTree, &dim)); 2997 PetscCall(PetscMalloc6(dim,&v0,dim,&v0parent,dim,&vtmp,dim*dim,&J,dim*dim,&Jparent,dim*dim,&invJ)); 2998 PetscCall(PetscMalloc2(dim,&pointScalar,dim,&pointRef)); 2999 PetscCall(DMGetDS(refTree,&ds)); 3000 PetscCall(PetscDSGetNumFields(ds,&numFields)); 3001 PetscCall(PetscSectionGetNumFields(section,&numSecFields)); 3002 PetscCall(DMGetLabel(refTree,"canonical",&canonical)); 3003 PetscCall(DMGetLabel(refTree,"depth",&depth)); 3004 PetscCall(DMGetDefaultConstraints(refTree,&cSection,&cMat,NULL)); 3005 PetscCall(DMPlexGetChart(refTree, &pStart, &pEnd)); 3006 PetscCall(DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd)); 3007 PetscCall(MatGetSize(cMat,&n,&m)); /* the injector has transpose sizes from the constraint matrix */ 3008 /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */ 3009 PetscCall(PetscCalloc1(m,&nnz)); 3010 for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */ 3011 const PetscInt *children; 3012 PetscInt numChildren; 3013 PetscInt i, numChildDof, numSelfDof; 3014 3015 if (canonical) { 3016 PetscInt pCanonical; 3017 PetscCall(DMLabelGetValue(canonical,p,&pCanonical)); 3018 if (p != pCanonical) continue; 3019 } 3020 PetscCall(DMPlexGetTreeChildren(refTree,p,&numChildren,&children)); 3021 if (!numChildren) continue; 3022 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3023 PetscInt child = children[i]; 3024 PetscInt dof; 3025 3026 PetscCall(PetscSectionGetDof(section,child,&dof)); 3027 numChildDof += dof; 3028 } 3029 PetscCall(PetscSectionGetDof(section,p,&numSelfDof)); 3030 if (!numChildDof || !numSelfDof) continue; 3031 for (f = 0; f < numFields; f++) { 3032 PetscInt selfOff; 3033 3034 if (numSecFields) { /* count the dofs for just this field */ 3035 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3036 PetscInt child = children[i]; 3037 PetscInt dof; 3038 3039 PetscCall(PetscSectionGetFieldDof(section,child,f,&dof)); 3040 numChildDof += dof; 3041 } 3042 PetscCall(PetscSectionGetFieldDof(section,p,f,&numSelfDof)); 3043 PetscCall(PetscSectionGetFieldOffset(section,p,f,&selfOff)); 3044 } 3045 else { 3046 PetscCall(PetscSectionGetOffset(section,p,&selfOff)); 3047 } 3048 for (i = 0; i < numSelfDof; i++) { 3049 nnz[selfOff + i] = numChildDof; 3050 } 3051 } 3052 } 3053 PetscCall(MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat)); 3054 PetscCall(PetscFree(nnz)); 3055 /* Setp 2: compute entries */ 3056 for (p = pStart; p < pEnd; p++) { 3057 const PetscInt *children; 3058 PetscInt numChildren; 3059 PetscInt i, numChildDof, numSelfDof; 3060 3061 /* same conditions about when entries occur */ 3062 if (canonical) { 3063 PetscInt pCanonical; 3064 PetscCall(DMLabelGetValue(canonical,p,&pCanonical)); 3065 if (p != pCanonical) continue; 3066 } 3067 PetscCall(DMPlexGetTreeChildren(refTree,p,&numChildren,&children)); 3068 if (!numChildren) continue; 3069 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3070 PetscInt child = children[i]; 3071 PetscInt dof; 3072 3073 PetscCall(PetscSectionGetDof(section,child,&dof)); 3074 numChildDof += dof; 3075 } 3076 PetscCall(PetscSectionGetDof(section,p,&numSelfDof)); 3077 if (!numChildDof || !numSelfDof) continue; 3078 3079 for (f = 0; f < numFields; f++) { 3080 PetscInt pI = -1, cI = -1; 3081 PetscInt selfOff, Nc, parentCell; 3082 PetscInt cellShapeOff; 3083 PetscObject disc; 3084 PetscDualSpace dsp; 3085 PetscClassId classId; 3086 PetscScalar *pointMat; 3087 PetscInt *matRows, *matCols; 3088 PetscInt pO = PETSC_MIN_INT; 3089 const PetscInt *depthNumDof; 3090 3091 if (numSecFields) { 3092 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3093 PetscInt child = children[i]; 3094 PetscInt dof; 3095 3096 PetscCall(PetscSectionGetFieldDof(section,child,f,&dof)); 3097 numChildDof += dof; 3098 } 3099 PetscCall(PetscSectionGetFieldDof(section,p,f,&numSelfDof)); 3100 PetscCall(PetscSectionGetFieldOffset(section,p,f,&selfOff)); 3101 } 3102 else { 3103 PetscCall(PetscSectionGetOffset(section,p,&selfOff)); 3104 } 3105 3106 /* find a cell whose closure contains p */ 3107 if (p >= cStart && p < cEnd) { 3108 parentCell = p; 3109 } 3110 else { 3111 PetscInt *star = NULL; 3112 PetscInt numStar; 3113 3114 parentCell = -1; 3115 PetscCall(DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star)); 3116 for (i = numStar - 1; i >= 0; i--) { 3117 PetscInt c = star[2 * i]; 3118 3119 if (c >= cStart && c < cEnd) { 3120 parentCell = c; 3121 break; 3122 } 3123 } 3124 PetscCall(DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star)); 3125 } 3126 /* determine the offset of p's shape functions within parentCell's shape functions */ 3127 PetscCall(PetscDSGetDiscretization(ds,f,&disc)); 3128 PetscCall(PetscObjectGetClassId(disc,&classId)); 3129 if (classId == PETSCFE_CLASSID) { 3130 PetscCall(PetscFEGetDualSpace((PetscFE)disc,&dsp)); 3131 } 3132 else if (classId == PETSCFV_CLASSID) { 3133 PetscCall(PetscFVGetDualSpace((PetscFV)disc,&dsp)); 3134 } 3135 else { 3136 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object"); 3137 } 3138 PetscCall(PetscDualSpaceGetNumDof(dsp,&depthNumDof)); 3139 PetscCall(PetscDualSpaceGetNumComponents(dsp,&Nc)); 3140 { 3141 PetscInt *closure = NULL; 3142 PetscInt numClosure; 3143 3144 PetscCall(DMPlexGetTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure)); 3145 for (i = 0, pI = -1, cellShapeOff = 0; i < numClosure; i++) { 3146 PetscInt point = closure[2 * i], pointDepth; 3147 3148 pO = closure[2 * i + 1]; 3149 if (point == p) { 3150 pI = i; 3151 break; 3152 } 3153 PetscCall(DMLabelGetValue(depth,point,&pointDepth)); 3154 cellShapeOff += depthNumDof[pointDepth]; 3155 } 3156 PetscCall(DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure)); 3157 } 3158 3159 PetscCall(DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat)); 3160 PetscCall(DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows)); 3161 matCols = matRows + numSelfDof; 3162 for (i = 0; i < numSelfDof; i++) { 3163 matRows[i] = selfOff + i; 3164 } 3165 for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.; 3166 { 3167 PetscInt colOff = 0; 3168 3169 for (i = 0; i < numChildren; i++) { 3170 PetscInt child = children[i]; 3171 PetscInt dof, off, j; 3172 3173 if (numSecFields) { 3174 PetscCall(PetscSectionGetFieldDof(cSection,child,f,&dof)); 3175 PetscCall(PetscSectionGetFieldOffset(cSection,child,f,&off)); 3176 } 3177 else { 3178 PetscCall(PetscSectionGetDof(cSection,child,&dof)); 3179 PetscCall(PetscSectionGetOffset(cSection,child,&off)); 3180 } 3181 3182 for (j = 0; j < dof; j++) { 3183 matCols[colOff++] = off + j; 3184 } 3185 } 3186 } 3187 if (classId == PETSCFE_CLASSID) { 3188 PetscFE fe = (PetscFE) disc; 3189 PetscInt fSize; 3190 const PetscInt ***perms; 3191 const PetscScalar ***flips; 3192 const PetscInt *pperms; 3193 3194 PetscCall(PetscFEGetDualSpace(fe,&dsp)); 3195 PetscCall(PetscDualSpaceGetDimension(dsp,&fSize)); 3196 PetscCall(PetscDualSpaceGetSymmetries(dsp, &perms, &flips)); 3197 pperms = perms ? perms[pI] ? perms[pI][pO] : NULL : NULL; 3198 for (i = 0; i < numSelfDof; i++) { /* for every shape function */ 3199 PetscQuadrature q; 3200 PetscInt dim, thisNc, numPoints, j, k; 3201 const PetscReal *points; 3202 const PetscReal *weights; 3203 PetscInt *closure = NULL; 3204 PetscInt numClosure; 3205 PetscInt iCell = pperms ? pperms[i] : i; 3206 PetscInt parentCellShapeDof = cellShapeOff + iCell; 3207 PetscTabulation Tparent; 3208 3209 PetscCall(PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q)); 3210 PetscCall(PetscQuadratureGetData(q,&dim,&thisNc,&numPoints,&points,&weights)); 3211 PetscCheck(thisNc == Nc,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %" PetscInt_FMT " does not much basis dim %" PetscInt_FMT,thisNc,Nc); 3212 PetscCall(PetscFECreateTabulation(fe,1,numPoints,points,0,&Tparent)); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */ 3213 for (j = 0; j < numPoints; j++) { 3214 PetscInt childCell = -1; 3215 PetscReal *parentValAtPoint; 3216 const PetscReal xi0[3] = {-1.,-1.,-1.}; 3217 const PetscReal *pointReal = &points[dim * j]; 3218 const PetscScalar *point; 3219 PetscTabulation Tchild; 3220 PetscInt childCellShapeOff, pointMatOff; 3221 #if defined(PETSC_USE_COMPLEX) 3222 PetscInt d; 3223 3224 for (d = 0; d < dim; d++) { 3225 pointScalar[d] = points[dim * j + d]; 3226 } 3227 point = pointScalar; 3228 #else 3229 point = pointReal; 3230 #endif 3231 3232 parentValAtPoint = &Tparent->T[0][(fSize * j + parentCellShapeDof) * Nc]; 3233 3234 for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/ 3235 PetscInt child = children[k]; 3236 PetscInt *star = NULL; 3237 PetscInt numStar, s; 3238 3239 PetscCall(DMPlexGetTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star)); 3240 for (s = numStar - 1; s >= 0; s--) { 3241 PetscInt c = star[2 * s]; 3242 3243 if (c < cStart || c >= cEnd) continue; 3244 PetscCall(DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell)); 3245 if (childCell >= 0) break; 3246 } 3247 PetscCall(DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star)); 3248 if (childCell >= 0) break; 3249 } 3250 PetscCheck(childCell >= 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not locate quadrature point"); 3251 PetscCall(DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ)); 3252 PetscCall(DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent)); 3253 CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp); 3254 CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef); 3255 3256 PetscCall(PetscFECreateTabulation(fe,1,1,pointRef,0,&Tchild)); 3257 PetscCall(DMPlexGetTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure)); 3258 for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */ 3259 PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT; 3260 PetscInt l; 3261 const PetscInt *cperms; 3262 3263 PetscCall(DMLabelGetValue(depth,child,&childDepth)); 3264 childDof = depthNumDof[childDepth]; 3265 for (l = 0, cI = -1, childCellShapeOff = 0; l < numClosure; l++) { 3266 PetscInt point = closure[2 * l]; 3267 PetscInt pointDepth; 3268 3269 childO = closure[2 * l + 1]; 3270 if (point == child) { 3271 cI = l; 3272 break; 3273 } 3274 PetscCall(DMLabelGetValue(depth,point,&pointDepth)); 3275 childCellShapeOff += depthNumDof[pointDepth]; 3276 } 3277 if (l == numClosure) { 3278 pointMatOff += childDof; 3279 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */ 3280 } 3281 cperms = perms ? perms[cI] ? perms[cI][childO] : NULL : NULL; 3282 for (l = 0; l < childDof; l++) { 3283 PetscInt lCell = cperms ? cperms[l] : l; 3284 PetscInt childCellDof = childCellShapeOff + lCell; 3285 PetscReal *childValAtPoint; 3286 PetscReal val = 0.; 3287 3288 childValAtPoint = &Tchild->T[0][childCellDof * Nc]; 3289 for (m = 0; m < Nc; m++) { 3290 val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m]; 3291 } 3292 3293 pointMat[i * numChildDof + pointMatOff + l] += val; 3294 } 3295 pointMatOff += childDof; 3296 } 3297 PetscCall(DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure)); 3298 PetscCall(PetscTabulationDestroy(&Tchild)); 3299 } 3300 PetscCall(PetscTabulationDestroy(&Tparent)); 3301 } 3302 } 3303 else { /* just the volume-weighted averages of the children */ 3304 PetscReal parentVol; 3305 PetscInt childCell; 3306 3307 PetscCall(DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL)); 3308 for (i = 0, childCell = 0; i < numChildren; i++) { 3309 PetscInt child = children[i], j; 3310 PetscReal childVol; 3311 3312 if (child < cStart || child >= cEnd) continue; 3313 PetscCall(DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL)); 3314 for (j = 0; j < Nc; j++) { 3315 pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol; 3316 } 3317 childCell++; 3318 } 3319 } 3320 /* Insert pointMat into mat */ 3321 PetscCall(MatSetValues(mat,numSelfDof,matRows,numChildDof,matCols,pointMat,INSERT_VALUES)); 3322 PetscCall(DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows)); 3323 PetscCall(DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat)); 3324 } 3325 } 3326 PetscCall(PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ)); 3327 PetscCall(PetscFree2(pointScalar,pointRef)); 3328 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 3329 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 3330 *inj = mat; 3331 PetscFunctionReturn(0); 3332 } 3333 3334 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3335 { 3336 PetscDS ds; 3337 PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof; 3338 PetscScalar ***refPointFieldMats; 3339 PetscSection refConSec, refSection; 3340 3341 PetscFunctionBegin; 3342 PetscCall(DMGetDS(refTree,&ds)); 3343 PetscCall(PetscDSGetNumFields(ds,&numFields)); 3344 PetscCall(DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL)); 3345 PetscCall(DMGetLocalSection(refTree,&refSection)); 3346 PetscCall(PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd)); 3347 PetscCall(PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats)); 3348 PetscCall(PetscSectionGetMaxDof(refConSec,&maxDof)); 3349 PetscCall(PetscMalloc1(maxDof,&rows)); 3350 PetscCall(PetscMalloc1(maxDof*maxDof,&cols)); 3351 for (p = pRefStart; p < pRefEnd; p++) { 3352 PetscInt parent, pDof, parentDof; 3353 3354 PetscCall(DMPlexGetTreeParent(refTree,p,&parent,NULL)); 3355 PetscCall(PetscSectionGetDof(refConSec,p,&pDof)); 3356 PetscCall(PetscSectionGetDof(refSection,parent,&parentDof)); 3357 if (!pDof || !parentDof || parent == p) continue; 3358 3359 PetscCall(PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart])); 3360 for (f = 0; f < numFields; f++) { 3361 PetscInt cDof, cOff, numCols, r; 3362 3363 if (numFields > 1) { 3364 PetscCall(PetscSectionGetFieldDof(refConSec,p,f,&cDof)); 3365 PetscCall(PetscSectionGetFieldOffset(refConSec,p,f,&cOff)); 3366 } 3367 else { 3368 PetscCall(PetscSectionGetDof(refConSec,p,&cDof)); 3369 PetscCall(PetscSectionGetOffset(refConSec,p,&cOff)); 3370 } 3371 3372 for (r = 0; r < cDof; r++) { 3373 rows[r] = cOff + r; 3374 } 3375 numCols = 0; 3376 { 3377 PetscInt aDof, aOff, j; 3378 3379 if (numFields > 1) { 3380 PetscCall(PetscSectionGetFieldDof(refSection,parent,f,&aDof)); 3381 PetscCall(PetscSectionGetFieldOffset(refSection,parent,f,&aOff)); 3382 } 3383 else { 3384 PetscCall(PetscSectionGetDof(refSection,parent,&aDof)); 3385 PetscCall(PetscSectionGetOffset(refSection,parent,&aOff)); 3386 } 3387 3388 for (j = 0; j < aDof; j++) { 3389 cols[numCols++] = aOff + j; 3390 } 3391 } 3392 PetscCall(PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f])); 3393 /* transpose of constraint matrix */ 3394 PetscCall(MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f])); 3395 } 3396 } 3397 *childrenMats = refPointFieldMats; 3398 PetscCall(PetscFree(rows)); 3399 PetscCall(PetscFree(cols)); 3400 PetscFunctionReturn(0); 3401 } 3402 3403 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3404 { 3405 PetscDS ds; 3406 PetscScalar ***refPointFieldMats; 3407 PetscInt numFields, pRefStart, pRefEnd, p, f; 3408 PetscSection refConSec, refSection; 3409 3410 PetscFunctionBegin; 3411 refPointFieldMats = *childrenMats; 3412 *childrenMats = NULL; 3413 PetscCall(DMGetDS(refTree,&ds)); 3414 PetscCall(DMGetLocalSection(refTree,&refSection)); 3415 PetscCall(PetscDSGetNumFields(ds,&numFields)); 3416 PetscCall(DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL)); 3417 PetscCall(PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd)); 3418 for (p = pRefStart; p < pRefEnd; p++) { 3419 PetscInt parent, pDof, parentDof; 3420 3421 PetscCall(DMPlexGetTreeParent(refTree,p,&parent,NULL)); 3422 PetscCall(PetscSectionGetDof(refConSec,p,&pDof)); 3423 PetscCall(PetscSectionGetDof(refSection,parent,&parentDof)); 3424 if (!pDof || !parentDof || parent == p) continue; 3425 3426 for (f = 0; f < numFields; f++) { 3427 PetscInt cDof; 3428 3429 if (numFields > 1) { 3430 PetscCall(PetscSectionGetFieldDof(refConSec,p,f,&cDof)); 3431 } 3432 else { 3433 PetscCall(PetscSectionGetDof(refConSec,p,&cDof)); 3434 } 3435 3436 PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f])); 3437 } 3438 PetscCall(PetscFree(refPointFieldMats[p - pRefStart])); 3439 } 3440 PetscCall(PetscFree(refPointFieldMats)); 3441 PetscFunctionReturn(0); 3442 } 3443 3444 static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree,Mat *injRef) 3445 { 3446 Mat cMatRef; 3447 PetscObject injRefObj; 3448 3449 PetscFunctionBegin; 3450 PetscCall(DMGetDefaultConstraints(refTree,NULL,&cMatRef,NULL)); 3451 PetscCall(PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj)); 3452 *injRef = (Mat) injRefObj; 3453 if (!*injRef) { 3454 PetscCall(DMPlexComputeInjectorReferenceTree(refTree,injRef)); 3455 PetscCall(PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)*injRef)); 3456 /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */ 3457 PetscCall(PetscObjectDereference((PetscObject)*injRef)); 3458 } 3459 PetscFunctionReturn(0); 3460 } 3461 3462 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) 3463 { 3464 PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti; 3465 PetscSection globalCoarse, globalFine; 3466 PetscSection localCoarse, localFine, leafIndicesSec; 3467 PetscSection multiRootSec, rootIndicesSec; 3468 PetscInt *leafInds, *rootInds = NULL; 3469 const PetscInt *rootDegrees; 3470 PetscScalar *leafVals = NULL, *rootVals = NULL; 3471 PetscSF coarseToFineEmbedded; 3472 3473 PetscFunctionBegin; 3474 PetscCall(DMPlexGetChart(coarse,&pStartC,&pEndC)); 3475 PetscCall(DMPlexGetChart(fine,&pStartF,&pEndF)); 3476 PetscCall(DMGetLocalSection(fine,&localFine)); 3477 PetscCall(DMGetGlobalSection(fine,&globalFine)); 3478 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec)); 3479 PetscCall(PetscSectionSetChart(leafIndicesSec,pStartF, pEndF)); 3480 PetscCall(PetscSectionGetMaxDof(localFine,&maxDof)); 3481 { /* winnow fine points that don't have global dofs out of the sf */ 3482 PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices; 3483 const PetscInt *leaves; 3484 3485 PetscCall(PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL)); 3486 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3487 p = leaves ? leaves[l] : l; 3488 PetscCall(PetscSectionGetDof(globalFine,p,&dof)); 3489 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&cdof)); 3490 if ((dof - cdof) > 0) { 3491 numPointsWithDofs++; 3492 3493 PetscCall(PetscSectionGetDof(localFine,p,&dof)); 3494 PetscCall(PetscSectionSetDof(leafIndicesSec,p,dof + 1)); 3495 } 3496 } 3497 PetscCall(PetscMalloc1(numPointsWithDofs,&pointsWithDofs)); 3498 PetscCall(PetscSectionSetUp(leafIndicesSec)); 3499 PetscCall(PetscSectionGetStorageSize(leafIndicesSec,&numIndices)); 3500 PetscCall(PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1),&leafInds)); 3501 if (gatheredValues) PetscCall(PetscMalloc1(numIndices,&leafVals)); 3502 for (l = 0, offset = 0; l < nleaves; l++) { 3503 p = leaves ? leaves[l] : l; 3504 PetscCall(PetscSectionGetDof(globalFine,p,&dof)); 3505 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&cdof)); 3506 if ((dof - cdof) > 0) { 3507 PetscInt off, gOff; 3508 PetscInt *pInd; 3509 PetscScalar *pVal = NULL; 3510 3511 pointsWithDofs[offset++] = l; 3512 3513 PetscCall(PetscSectionGetOffset(leafIndicesSec,p,&off)); 3514 3515 pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds; 3516 if (gatheredValues) { 3517 PetscInt i; 3518 3519 pVal = &leafVals[off + 1]; 3520 for (i = 0; i < dof; i++) pVal[i] = 0.; 3521 } 3522 PetscCall(PetscSectionGetOffset(globalFine,p,&gOff)); 3523 3524 offsets[0] = 0; 3525 if (numFields) { 3526 PetscInt f; 3527 3528 for (f = 0; f < numFields; f++) { 3529 PetscInt fDof; 3530 PetscCall(PetscSectionGetFieldDof(localFine,p,f,&fDof)); 3531 offsets[f + 1] = fDof + offsets[f]; 3532 } 3533 PetscCall(DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1, NULL,pInd)); 3534 } else { 3535 PetscCall(DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL, NULL,pInd)); 3536 } 3537 if (gatheredValues) PetscCall(VecGetValues(fineVec,dof,pInd,pVal)); 3538 } 3539 } 3540 PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded)); 3541 PetscCall(PetscFree(pointsWithDofs)); 3542 } 3543 3544 PetscCall(DMPlexGetChart(coarse,&pStartC,&pEndC)); 3545 PetscCall(DMGetLocalSection(coarse,&localCoarse)); 3546 PetscCall(DMGetGlobalSection(coarse,&globalCoarse)); 3547 3548 { /* there may be the case where an sf root has a parent: broadcast parents back to children */ 3549 MPI_Datatype threeInt; 3550 PetscMPIInt rank; 3551 PetscInt (*parentNodeAndIdCoarse)[3]; 3552 PetscInt (*parentNodeAndIdFine)[3]; 3553 PetscInt p, nleaves, nleavesToParents; 3554 PetscSF pointSF, sfToParents; 3555 const PetscInt *ilocal; 3556 const PetscSFNode *iremote; 3557 PetscSFNode *iremoteToParents; 3558 PetscInt *ilocalToParents; 3559 3560 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank)); 3561 PetscCallMPI(MPI_Type_contiguous(3,MPIU_INT,&threeInt)); 3562 PetscCallMPI(MPI_Type_commit(&threeInt)); 3563 PetscCall(PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine)); 3564 PetscCall(DMGetPointSF(coarse,&pointSF)); 3565 PetscCall(PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote)); 3566 for (p = pStartC; p < pEndC; p++) { 3567 PetscInt parent, childId; 3568 PetscCall(DMPlexGetTreeParent(coarse,p,&parent,&childId)); 3569 parentNodeAndIdCoarse[p - pStartC][0] = rank; 3570 parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC; 3571 parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId; 3572 if (nleaves > 0) { 3573 PetscInt leaf = -1; 3574 3575 if (ilocal) { 3576 PetscCall(PetscFindInt(parent,nleaves,ilocal,&leaf)); 3577 } 3578 else { 3579 leaf = p - pStartC; 3580 } 3581 if (leaf >= 0) { 3582 parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank; 3583 parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index; 3584 } 3585 } 3586 } 3587 for (p = pStartF; p < pEndF; p++) { 3588 parentNodeAndIdFine[p - pStartF][0] = -1; 3589 parentNodeAndIdFine[p - pStartF][1] = -1; 3590 parentNodeAndIdFine[p - pStartF][2] = -1; 3591 } 3592 PetscCall(PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine,MPI_REPLACE)); 3593 PetscCall(PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine,MPI_REPLACE)); 3594 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3595 PetscInt dof; 3596 3597 PetscCall(PetscSectionGetDof(leafIndicesSec,p,&dof)); 3598 if (dof) { 3599 PetscInt off; 3600 3601 PetscCall(PetscSectionGetOffset(leafIndicesSec,p,&off)); 3602 if (gatheredIndices) { 3603 leafInds[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3604 } else if (gatheredValues) { 3605 leafVals[off] = (PetscScalar) PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3606 } 3607 } 3608 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3609 nleavesToParents++; 3610 } 3611 } 3612 PetscCall(PetscMalloc1(nleavesToParents,&ilocalToParents)); 3613 PetscCall(PetscMalloc1(nleavesToParents,&iremoteToParents)); 3614 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3615 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3616 ilocalToParents[nleavesToParents] = p - pStartF; 3617 iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p-pStartF][0]; 3618 iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1]; 3619 nleavesToParents++; 3620 } 3621 } 3622 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents)); 3623 PetscCall(PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER)); 3624 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 3625 3626 coarseToFineEmbedded = sfToParents; 3627 3628 PetscCall(PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine)); 3629 PetscCallMPI(MPI_Type_free(&threeInt)); 3630 } 3631 3632 { /* winnow out coarse points that don't have dofs */ 3633 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3634 PetscSF sfDofsOnly; 3635 3636 for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) { 3637 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 3638 PetscCall(PetscSectionGetConstraintDof(globalCoarse,p,&cdof)); 3639 if ((dof - cdof) > 0) { 3640 numPointsWithDofs++; 3641 } 3642 } 3643 PetscCall(PetscMalloc1(numPointsWithDofs,&pointsWithDofs)); 3644 for (p = pStartC, offset = 0; p < pEndC; p++) { 3645 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 3646 PetscCall(PetscSectionGetConstraintDof(globalCoarse,p,&cdof)); 3647 if ((dof - cdof) > 0) { 3648 pointsWithDofs[offset++] = p - pStartC; 3649 } 3650 } 3651 PetscCall(PetscSFCreateEmbeddedRootSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly)); 3652 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 3653 PetscCall(PetscFree(pointsWithDofs)); 3654 coarseToFineEmbedded = sfDofsOnly; 3655 } 3656 3657 /* communicate back to the coarse mesh which coarse points have children (that may require injection) */ 3658 PetscCall(PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees)); 3659 PetscCall(PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees)); 3660 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec)); 3661 PetscCall(PetscSectionSetChart(multiRootSec,pStartC,pEndC)); 3662 for (p = pStartC; p < pEndC; p++) { 3663 PetscCall(PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC])); 3664 } 3665 PetscCall(PetscSectionSetUp(multiRootSec)); 3666 PetscCall(PetscSectionGetStorageSize(multiRootSec,&numMulti)); 3667 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec)); 3668 { /* distribute the leaf section */ 3669 PetscSF multi, multiInv, indicesSF; 3670 PetscInt *remoteOffsets, numRootIndices; 3671 3672 PetscCall(PetscSFGetMultiSF(coarseToFineEmbedded,&multi)); 3673 PetscCall(PetscSFCreateInverseSF(multi,&multiInv)); 3674 PetscCall(PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec)); 3675 PetscCall(PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF)); 3676 PetscCall(PetscFree(remoteOffsets)); 3677 PetscCall(PetscSFDestroy(&multiInv)); 3678 PetscCall(PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices)); 3679 if (gatheredIndices) { 3680 PetscCall(PetscMalloc1(numRootIndices,&rootInds)); 3681 PetscCall(PetscSFBcastBegin(indicesSF,MPIU_INT,leafInds,rootInds,MPI_REPLACE)); 3682 PetscCall(PetscSFBcastEnd(indicesSF,MPIU_INT,leafInds,rootInds,MPI_REPLACE)); 3683 } 3684 if (gatheredValues) { 3685 PetscCall(PetscMalloc1(numRootIndices,&rootVals)); 3686 PetscCall(PetscSFBcastBegin(indicesSF,MPIU_SCALAR,leafVals,rootVals,MPI_REPLACE)); 3687 PetscCall(PetscSFBcastEnd(indicesSF,MPIU_SCALAR,leafVals,rootVals,MPI_REPLACE)); 3688 } 3689 PetscCall(PetscSFDestroy(&indicesSF)); 3690 } 3691 PetscCall(PetscSectionDestroy(&leafIndicesSec)); 3692 PetscCall(PetscFree(leafInds)); 3693 PetscCall(PetscFree(leafVals)); 3694 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 3695 *rootMultiSec = multiRootSec; 3696 *multiLeafSec = rootIndicesSec; 3697 if (gatheredIndices) *gatheredIndices = rootInds; 3698 if (gatheredValues) *gatheredValues = rootVals; 3699 PetscFunctionReturn(0); 3700 } 3701 3702 PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 3703 { 3704 DM refTree; 3705 PetscSection multiRootSec, rootIndicesSec; 3706 PetscSection globalCoarse, globalFine; 3707 PetscSection localCoarse, localFine; 3708 PetscSection cSecRef; 3709 PetscInt *rootIndices = NULL, *parentIndices, pRefStart, pRefEnd; 3710 Mat injRef; 3711 PetscInt numFields, maxDof; 3712 PetscInt pStartC, pEndC, pStartF, pEndF, p; 3713 PetscInt *offsets, *offsetsCopy, *rowOffsets; 3714 PetscLayout rowMap, colMap; 3715 PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO; 3716 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 3717 3718 PetscFunctionBegin; 3719 3720 /* get the templates for the fine-to-coarse injection from the reference tree */ 3721 PetscCall(DMPlexGetReferenceTree(coarse,&refTree)); 3722 PetscCall(DMGetDefaultConstraints(refTree,&cSecRef,NULL,NULL)); 3723 PetscCall(PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd)); 3724 PetscCall(DMPlexReferenceTreeGetInjector(refTree,&injRef)); 3725 3726 PetscCall(DMPlexGetChart(fine,&pStartF,&pEndF)); 3727 PetscCall(DMGetLocalSection(fine,&localFine)); 3728 PetscCall(DMGetGlobalSection(fine,&globalFine)); 3729 PetscCall(PetscSectionGetNumFields(localFine,&numFields)); 3730 PetscCall(DMPlexGetChart(coarse,&pStartC,&pEndC)); 3731 PetscCall(DMGetLocalSection(coarse,&localCoarse)); 3732 PetscCall(DMGetGlobalSection(coarse,&globalCoarse)); 3733 PetscCall(PetscSectionGetMaxDof(localCoarse,&maxDof)); 3734 { 3735 PetscInt maxFields = PetscMax(1,numFields) + 1; 3736 PetscCall(PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets)); 3737 } 3738 3739 PetscCall(DMPlexTransferInjectorTree(coarse,fine,coarseToFine,childIds,NULL,numFields,offsets,&multiRootSec,&rootIndicesSec,&rootIndices,NULL)); 3740 3741 PetscCall(PetscMalloc1(maxDof,&parentIndices)); 3742 3743 /* count indices */ 3744 PetscCall(MatGetLayouts(mat,&rowMap,&colMap)); 3745 PetscCall(PetscLayoutSetUp(rowMap)); 3746 PetscCall(PetscLayoutSetUp(colMap)); 3747 PetscCall(PetscLayoutGetRange(rowMap,&rowStart,&rowEnd)); 3748 PetscCall(PetscLayoutGetRange(colMap,&colStart,&colEnd)); 3749 PetscCall(PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO)); 3750 for (p = pStartC; p < pEndC; p++) { 3751 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3752 3753 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 3754 PetscCall(PetscSectionGetConstraintDof(globalCoarse,p,&cdof)); 3755 if ((dof - cdof) <= 0) continue; 3756 PetscCall(PetscSectionGetOffset(globalCoarse,p,&gOff)); 3757 3758 rowOffsets[0] = 0; 3759 offsetsCopy[0] = 0; 3760 if (numFields) { 3761 PetscInt f; 3762 3763 for (f = 0; f < numFields; f++) { 3764 PetscInt fDof; 3765 PetscCall(PetscSectionGetFieldDof(localCoarse,p,f,&fDof)); 3766 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3767 } 3768 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,parentIndices)); 3769 } else { 3770 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,parentIndices)); 3771 rowOffsets[1] = offsetsCopy[0]; 3772 } 3773 3774 PetscCall(PetscSectionGetDof(multiRootSec,p,&numLeaves)); 3775 PetscCall(PetscSectionGetOffset(multiRootSec,p,&leafStart)); 3776 leafEnd = leafStart + numLeaves; 3777 for (l = leafStart; l < leafEnd; l++) { 3778 PetscInt numIndices, childId, offset; 3779 const PetscInt *childIndices; 3780 3781 PetscCall(PetscSectionGetDof(rootIndicesSec,l,&numIndices)); 3782 PetscCall(PetscSectionGetOffset(rootIndicesSec,l,&offset)); 3783 childId = rootIndices[offset++]; 3784 childIndices = &rootIndices[offset]; 3785 numIndices--; 3786 3787 if (childId == -1) { /* equivalent points: scatter */ 3788 PetscInt i; 3789 3790 for (i = 0; i < numIndices; i++) { 3791 PetscInt colIndex = childIndices[i]; 3792 PetscInt rowIndex = parentIndices[i]; 3793 if (rowIndex < 0) continue; 3794 PetscCheck(colIndex >= 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unconstrained fine and constrained coarse"); 3795 if (colIndex >= colStart && colIndex < colEnd) { 3796 nnzD[rowIndex - rowStart] = 1; 3797 } 3798 else { 3799 nnzO[rowIndex - rowStart] = 1; 3800 } 3801 } 3802 } 3803 else { 3804 PetscInt parentId, f, lim; 3805 3806 PetscCall(DMPlexGetTreeParent(refTree,childId,&parentId,NULL)); 3807 3808 lim = PetscMax(1,numFields); 3809 offsets[0] = 0; 3810 if (numFields) { 3811 PetscInt f; 3812 3813 for (f = 0; f < numFields; f++) { 3814 PetscInt fDof; 3815 PetscCall(PetscSectionGetFieldDof(cSecRef,childId,f,&fDof)); 3816 3817 offsets[f + 1] = fDof + offsets[f]; 3818 } 3819 } 3820 else { 3821 PetscInt cDof; 3822 3823 PetscCall(PetscSectionGetDof(cSecRef,childId,&cDof)); 3824 offsets[1] = cDof; 3825 } 3826 for (f = 0; f < lim; f++) { 3827 PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1]; 3828 PetscInt childStart = offsets[f], childEnd = offsets[f + 1]; 3829 PetscInt i, numD = 0, numO = 0; 3830 3831 for (i = childStart; i < childEnd; i++) { 3832 PetscInt colIndex = childIndices[i]; 3833 3834 if (colIndex < 0) continue; 3835 if (colIndex >= colStart && colIndex < colEnd) { 3836 numD++; 3837 } 3838 else { 3839 numO++; 3840 } 3841 } 3842 for (i = parentStart; i < parentEnd; i++) { 3843 PetscInt rowIndex = parentIndices[i]; 3844 3845 if (rowIndex < 0) continue; 3846 nnzD[rowIndex - rowStart] += numD; 3847 nnzO[rowIndex - rowStart] += numO; 3848 } 3849 } 3850 } 3851 } 3852 } 3853 /* preallocate */ 3854 PetscCall(MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL)); 3855 PetscCall(PetscFree2(nnzD,nnzO)); 3856 /* insert values */ 3857 PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats)); 3858 for (p = pStartC; p < pEndC; p++) { 3859 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3860 3861 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 3862 PetscCall(PetscSectionGetConstraintDof(globalCoarse,p,&cdof)); 3863 if ((dof - cdof) <= 0) continue; 3864 PetscCall(PetscSectionGetOffset(globalCoarse,p,&gOff)); 3865 3866 rowOffsets[0] = 0; 3867 offsetsCopy[0] = 0; 3868 if (numFields) { 3869 PetscInt f; 3870 3871 for (f = 0; f < numFields; f++) { 3872 PetscInt fDof; 3873 PetscCall(PetscSectionGetFieldDof(localCoarse,p,f,&fDof)); 3874 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3875 } 3876 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1, NULL,parentIndices)); 3877 } else { 3878 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL, NULL,parentIndices)); 3879 rowOffsets[1] = offsetsCopy[0]; 3880 } 3881 3882 PetscCall(PetscSectionGetDof(multiRootSec,p,&numLeaves)); 3883 PetscCall(PetscSectionGetOffset(multiRootSec,p,&leafStart)); 3884 leafEnd = leafStart + numLeaves; 3885 for (l = leafStart; l < leafEnd; l++) { 3886 PetscInt numIndices, childId, offset; 3887 const PetscInt *childIndices; 3888 3889 PetscCall(PetscSectionGetDof(rootIndicesSec,l,&numIndices)); 3890 PetscCall(PetscSectionGetOffset(rootIndicesSec,l,&offset)); 3891 childId = rootIndices[offset++]; 3892 childIndices = &rootIndices[offset]; 3893 numIndices--; 3894 3895 if (childId == -1) { /* equivalent points: scatter */ 3896 PetscInt i; 3897 3898 for (i = 0; i < numIndices; i++) { 3899 PetscCall(MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES)); 3900 } 3901 } 3902 else { 3903 PetscInt parentId, f, lim; 3904 3905 PetscCall(DMPlexGetTreeParent(refTree,childId,&parentId,NULL)); 3906 3907 lim = PetscMax(1,numFields); 3908 offsets[0] = 0; 3909 if (numFields) { 3910 PetscInt f; 3911 3912 for (f = 0; f < numFields; f++) { 3913 PetscInt fDof; 3914 PetscCall(PetscSectionGetFieldDof(cSecRef,childId,f,&fDof)); 3915 3916 offsets[f + 1] = fDof + offsets[f]; 3917 } 3918 } 3919 else { 3920 PetscInt cDof; 3921 3922 PetscCall(PetscSectionGetDof(cSecRef,childId,&cDof)); 3923 offsets[1] = cDof; 3924 } 3925 for (f = 0; f < lim; f++) { 3926 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 3927 PetscInt *rowIndices = &parentIndices[rowOffsets[f]]; 3928 const PetscInt *colIndices = &childIndices[offsets[f]]; 3929 3930 PetscCall(MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES)); 3931 } 3932 } 3933 } 3934 } 3935 PetscCall(PetscSectionDestroy(&multiRootSec)); 3936 PetscCall(PetscSectionDestroy(&rootIndicesSec)); 3937 PetscCall(PetscFree(parentIndices)); 3938 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats)); 3939 PetscCall(PetscFree(rootIndices)); 3940 PetscCall(PetscFree3(offsets,offsetsCopy,rowOffsets)); 3941 3942 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 3943 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 3944 PetscFunctionReturn(0); 3945 } 3946 3947 static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom) 3948 { 3949 PetscSF coarseToFineEmbedded; 3950 PetscSection globalCoarse, globalFine; 3951 PetscSection localCoarse, localFine; 3952 PetscSection aSec, cSec; 3953 PetscSection rootValuesSec; 3954 PetscSection leafValuesSec; 3955 PetscScalar *rootValues, *leafValues; 3956 IS aIS; 3957 const PetscInt *anchors; 3958 Mat cMat; 3959 PetscInt numFields; 3960 PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd; 3961 PetscInt aStart, aEnd, cStart, cEnd; 3962 PetscInt *maxChildIds; 3963 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 3964 PetscFV fv = NULL; 3965 PetscInt dim, numFVcomps = -1, fvField = -1; 3966 DM cellDM = NULL, gradDM = NULL; 3967 const PetscScalar *cellGeomArray = NULL; 3968 const PetscScalar *gradArray = NULL; 3969 3970 PetscFunctionBegin; 3971 PetscCall(VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE)); 3972 PetscCall(DMPlexGetChart(coarse,&pStartC,&pEndC)); 3973 PetscCall(DMPlexGetSimplexOrBoxCells(coarse,0,&cellStart,&cellEnd)); 3974 PetscCall(DMPlexGetChart(fine,&pStartF,&pEndF)); 3975 PetscCall(DMGetGlobalSection(fine,&globalFine)); 3976 PetscCall(DMGetCoordinateDim(coarse,&dim)); 3977 { /* winnow fine points that don't have global dofs out of the sf */ 3978 PetscInt nleaves, l; 3979 const PetscInt *leaves; 3980 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3981 3982 PetscCall(PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL)); 3983 3984 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3985 PetscInt p = leaves ? leaves[l] : l; 3986 3987 PetscCall(PetscSectionGetDof(globalFine,p,&dof)); 3988 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&cdof)); 3989 if ((dof - cdof) > 0) { 3990 numPointsWithDofs++; 3991 } 3992 } 3993 PetscCall(PetscMalloc1(numPointsWithDofs,&pointsWithDofs)); 3994 for (l = 0, offset = 0; l < nleaves; l++) { 3995 PetscInt p = leaves ? leaves[l] : l; 3996 3997 PetscCall(PetscSectionGetDof(globalFine,p,&dof)); 3998 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&cdof)); 3999 if ((dof - cdof) > 0) { 4000 pointsWithDofs[offset++] = l; 4001 } 4002 } 4003 PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded)); 4004 PetscCall(PetscFree(pointsWithDofs)); 4005 } 4006 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 4007 PetscCall(PetscMalloc1(pEndC-pStartC,&maxChildIds)); 4008 for (p = pStartC; p < pEndC; p++) { 4009 maxChildIds[p - pStartC] = -2; 4010 } 4011 PetscCall(PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX)); 4012 PetscCall(PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX)); 4013 4014 PetscCall(DMGetLocalSection(coarse,&localCoarse)); 4015 PetscCall(DMGetGlobalSection(coarse,&globalCoarse)); 4016 4017 PetscCall(DMPlexGetAnchors(coarse,&aSec,&aIS)); 4018 PetscCall(ISGetIndices(aIS,&anchors)); 4019 PetscCall(PetscSectionGetChart(aSec,&aStart,&aEnd)); 4020 4021 PetscCall(DMGetDefaultConstraints(coarse,&cSec,&cMat,NULL)); 4022 PetscCall(PetscSectionGetChart(cSec,&cStart,&cEnd)); 4023 4024 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 4025 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootValuesSec)); 4026 PetscCall(PetscSectionSetChart(rootValuesSec,pStartC,pEndC)); 4027 PetscCall(PetscSectionGetNumFields(localCoarse,&numFields)); 4028 { 4029 PetscInt maxFields = PetscMax(1,numFields) + 1; 4030 PetscCall(PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO)); 4031 } 4032 if (grad) { 4033 PetscInt i; 4034 4035 PetscCall(VecGetDM(cellGeom,&cellDM)); 4036 PetscCall(VecGetArrayRead(cellGeom,&cellGeomArray)); 4037 PetscCall(VecGetDM(grad,&gradDM)); 4038 PetscCall(VecGetArrayRead(grad,&gradArray)); 4039 for (i = 0; i < PetscMax(1,numFields); i++) { 4040 PetscObject obj; 4041 PetscClassId id; 4042 4043 PetscCall(DMGetField(coarse, i, NULL, &obj)); 4044 PetscCall(PetscObjectGetClassId(obj,&id)); 4045 if (id == PETSCFV_CLASSID) { 4046 fv = (PetscFV) obj; 4047 PetscCall(PetscFVGetNumComponents(fv,&numFVcomps)); 4048 fvField = i; 4049 break; 4050 } 4051 } 4052 } 4053 4054 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 4055 PetscInt dof; 4056 PetscInt maxChildId = maxChildIds[p - pStartC]; 4057 PetscInt numValues = 0; 4058 4059 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 4060 if (dof < 0) { 4061 dof = -(dof + 1); 4062 } 4063 offsets[0] = 0; 4064 newOffsets[0] = 0; 4065 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 4066 PetscInt *closure = NULL, closureSize, cl; 4067 4068 PetscCall(DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure)); 4069 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 4070 PetscInt c = closure[2 * cl], clDof; 4071 4072 PetscCall(PetscSectionGetDof(localCoarse,c,&clDof)); 4073 numValues += clDof; 4074 } 4075 PetscCall(DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure)); 4076 } 4077 else if (maxChildId == -1) { 4078 PetscCall(PetscSectionGetDof(localCoarse,p,&numValues)); 4079 } 4080 /* we will pack the column indices with the field offsets */ 4081 if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) { 4082 /* also send the centroid, and the gradient */ 4083 numValues += dim * (1 + numFVcomps); 4084 } 4085 PetscCall(PetscSectionSetDof(rootValuesSec,p,numValues)); 4086 } 4087 PetscCall(PetscSectionSetUp(rootValuesSec)); 4088 { 4089 PetscInt numRootValues; 4090 const PetscScalar *coarseArray; 4091 4092 PetscCall(PetscSectionGetStorageSize(rootValuesSec,&numRootValues)); 4093 PetscCall(PetscMalloc1(numRootValues,&rootValues)); 4094 PetscCall(VecGetArrayRead(vecCoarseLocal,&coarseArray)); 4095 for (p = pStartC; p < pEndC; p++) { 4096 PetscInt numValues; 4097 PetscInt pValOff; 4098 PetscScalar *pVal; 4099 PetscInt maxChildId = maxChildIds[p - pStartC]; 4100 4101 PetscCall(PetscSectionGetDof(rootValuesSec,p,&numValues)); 4102 if (!numValues) { 4103 continue; 4104 } 4105 PetscCall(PetscSectionGetOffset(rootValuesSec,p,&pValOff)); 4106 pVal = &(rootValues[pValOff]); 4107 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 4108 PetscInt closureSize = numValues; 4109 PetscCall(DMPlexVecGetClosure(coarse,NULL,vecCoarseLocal,p,&closureSize,&pVal)); 4110 if (grad && p >= cellStart && p < cellEnd) { 4111 PetscFVCellGeom *cg; 4112 PetscScalar *gradVals = NULL; 4113 PetscInt i; 4114 4115 pVal += (numValues - dim * (1 + numFVcomps)); 4116 4117 PetscCall(DMPlexPointLocalRead(cellDM,p,cellGeomArray,(void *) &cg)); 4118 for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i]; 4119 pVal += dim; 4120 PetscCall(DMPlexPointGlobalRead(gradDM,p,gradArray,(void *) &gradVals)); 4121 for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i]; 4122 } 4123 } 4124 else if (maxChildId == -1) { 4125 PetscInt lDof, lOff, i; 4126 4127 PetscCall(PetscSectionGetDof(localCoarse,p,&lDof)); 4128 PetscCall(PetscSectionGetOffset(localCoarse,p,&lOff)); 4129 for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i]; 4130 } 4131 } 4132 PetscCall(VecRestoreArrayRead(vecCoarseLocal,&coarseArray)); 4133 PetscCall(PetscFree(maxChildIds)); 4134 } 4135 { 4136 PetscSF valuesSF; 4137 PetscInt *remoteOffsetsValues, numLeafValues; 4138 4139 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafValuesSec)); 4140 PetscCall(PetscSFDistributeSection(coarseToFineEmbedded,rootValuesSec,&remoteOffsetsValues,leafValuesSec)); 4141 PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded,rootValuesSec,remoteOffsetsValues,leafValuesSec,&valuesSF)); 4142 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 4143 PetscCall(PetscFree(remoteOffsetsValues)); 4144 PetscCall(PetscSectionGetStorageSize(leafValuesSec,&numLeafValues)); 4145 PetscCall(PetscMalloc1(numLeafValues,&leafValues)); 4146 PetscCall(PetscSFBcastBegin(valuesSF,MPIU_SCALAR,rootValues,leafValues,MPI_REPLACE)); 4147 PetscCall(PetscSFBcastEnd(valuesSF,MPIU_SCALAR,rootValues,leafValues,MPI_REPLACE)); 4148 PetscCall(PetscSFDestroy(&valuesSF)); 4149 PetscCall(PetscFree(rootValues)); 4150 PetscCall(PetscSectionDestroy(&rootValuesSec)); 4151 } 4152 PetscCall(DMGetLocalSection(fine,&localFine)); 4153 { 4154 PetscInt maxDof; 4155 PetscInt *rowIndices; 4156 DM refTree; 4157 PetscInt **refPointFieldN; 4158 PetscScalar ***refPointFieldMats; 4159 PetscSection refConSec, refAnSec; 4160 PetscInt pRefStart,pRefEnd,leafStart,leafEnd; 4161 PetscScalar *pointWork; 4162 4163 PetscCall(PetscSectionGetMaxDof(localFine,&maxDof)); 4164 PetscCall(DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices)); 4165 PetscCall(DMGetWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork)); 4166 PetscCall(DMPlexGetReferenceTree(fine,&refTree)); 4167 PetscCall(DMCopyDisc(fine,refTree)); 4168 PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN)); 4169 PetscCall(DMGetDefaultConstraints(refTree,&refConSec,NULL,NULL)); 4170 PetscCall(DMPlexGetAnchors(refTree,&refAnSec,NULL)); 4171 PetscCall(PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd)); 4172 PetscCall(PetscSectionGetChart(leafValuesSec,&leafStart,&leafEnd)); 4173 PetscCall(DMPlexGetSimplexOrBoxCells(fine,0,&cellStart,&cellEnd)); 4174 for (p = leafStart; p < leafEnd; p++) { 4175 PetscInt gDof, gcDof, gOff, lDof; 4176 PetscInt numValues, pValOff; 4177 PetscInt childId; 4178 const PetscScalar *pVal; 4179 const PetscScalar *fvGradData = NULL; 4180 4181 PetscCall(PetscSectionGetDof(globalFine,p,&gDof)); 4182 PetscCall(PetscSectionGetDof(localFine,p,&lDof)); 4183 PetscCall(PetscSectionGetConstraintDof(globalFine,p,&gcDof)); 4184 if ((gDof - gcDof) <= 0) { 4185 continue; 4186 } 4187 PetscCall(PetscSectionGetOffset(globalFine,p,&gOff)); 4188 PetscCall(PetscSectionGetDof(leafValuesSec,p,&numValues)); 4189 if (!numValues) continue; 4190 PetscCall(PetscSectionGetOffset(leafValuesSec,p,&pValOff)); 4191 pVal = &leafValues[pValOff]; 4192 offsets[0] = 0; 4193 offsetsCopy[0] = 0; 4194 newOffsets[0] = 0; 4195 newOffsetsCopy[0] = 0; 4196 childId = cids[p - pStartF]; 4197 if (numFields) { 4198 PetscInt f; 4199 for (f = 0; f < numFields; f++) { 4200 PetscInt rowDof; 4201 4202 PetscCall(PetscSectionGetFieldDof(localFine,p,f,&rowDof)); 4203 offsets[f + 1] = offsets[f] + rowDof; 4204 offsetsCopy[f + 1] = offsets[f + 1]; 4205 /* TODO: closure indices */ 4206 newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]); 4207 } 4208 PetscCall(DMPlexGetIndicesPointFields_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,NULL,rowIndices)); 4209 } 4210 else { 4211 offsets[0] = 0; 4212 offsets[1] = lDof; 4213 newOffsets[0] = 0; 4214 newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0]; 4215 PetscCall(DMPlexGetIndicesPoint_Internal(localFine,PETSC_FALSE,p,gOff,offsetsCopy,PETSC_FALSE,NULL,NULL,rowIndices)); 4216 } 4217 if (childId == -1) { /* no child interpolation: one nnz per */ 4218 PetscCall(VecSetValues(vecFine,numValues,rowIndices,pVal,INSERT_VALUES)); 4219 } else { 4220 PetscInt f; 4221 4222 if (grad && p >= cellStart && p < cellEnd) { 4223 numValues -= (dim * (1 + numFVcomps)); 4224 fvGradData = &pVal[numValues]; 4225 } 4226 for (f = 0; f < PetscMax(1,numFields); f++) { 4227 const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f]; 4228 PetscInt numRows = offsets[f+1] - offsets[f]; 4229 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 4230 const PetscScalar *cVal = &pVal[newOffsets[f]]; 4231 PetscScalar *rVal = &pointWork[offsets[f]]; 4232 PetscInt i, j; 4233 4234 #if 0 4235 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)); 4236 #endif 4237 for (i = 0; i < numRows; i++) { 4238 PetscScalar val = 0.; 4239 for (j = 0; j < numCols; j++) { 4240 val += childMat[i * numCols + j] * cVal[j]; 4241 } 4242 rVal[i] = val; 4243 } 4244 if (f == fvField && p >= cellStart && p < cellEnd) { 4245 PetscReal centroid[3]; 4246 PetscScalar diff[3]; 4247 const PetscScalar *parentCentroid = &fvGradData[0]; 4248 const PetscScalar *gradient = &fvGradData[dim]; 4249 4250 PetscCall(DMPlexComputeCellGeometryFVM(fine,p,NULL,centroid,NULL)); 4251 for (i = 0; i < dim; i++) { 4252 diff[i] = centroid[i] - parentCentroid[i]; 4253 } 4254 for (i = 0; i < numFVcomps; i++) { 4255 PetscScalar val = 0.; 4256 4257 for (j = 0; j < dim; j++) { 4258 val += gradient[dim * i + j] * diff[j]; 4259 } 4260 rVal[i] += val; 4261 } 4262 } 4263 PetscCall(VecSetValues(vecFine,numRows,&rowIndices[offsets[f]],rVal,INSERT_VALUES)); 4264 } 4265 } 4266 } 4267 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN)); 4268 PetscCall(DMRestoreWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork)); 4269 PetscCall(DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices)); 4270 } 4271 PetscCall(PetscFree(leafValues)); 4272 PetscCall(PetscSectionDestroy(&leafValuesSec)); 4273 PetscCall(PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO)); 4274 PetscCall(ISRestoreIndices(aIS,&anchors)); 4275 PetscFunctionReturn(0); 4276 } 4277 4278 static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids) 4279 { 4280 DM refTree; 4281 PetscSection multiRootSec, rootIndicesSec; 4282 PetscSection globalCoarse, globalFine; 4283 PetscSection localCoarse, localFine; 4284 PetscSection cSecRef; 4285 PetscInt *parentIndices, pRefStart, pRefEnd; 4286 PetscScalar *rootValues, *parentValues; 4287 Mat injRef; 4288 PetscInt numFields, maxDof; 4289 PetscInt pStartC, pEndC, pStartF, pEndF, p; 4290 PetscInt *offsets, *offsetsCopy, *rowOffsets; 4291 PetscLayout rowMap, colMap; 4292 PetscInt rowStart, rowEnd, colStart, colEnd; 4293 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 4294 4295 PetscFunctionBegin; 4296 4297 /* get the templates for the fine-to-coarse injection from the reference tree */ 4298 PetscCall(VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE)); 4299 PetscCall(VecSetOption(vecCoarse,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE)); 4300 PetscCall(DMPlexGetReferenceTree(coarse,&refTree)); 4301 PetscCall(DMCopyDisc(coarse,refTree)); 4302 PetscCall(DMGetDefaultConstraints(refTree,&cSecRef,NULL,NULL)); 4303 PetscCall(PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd)); 4304 PetscCall(DMPlexReferenceTreeGetInjector(refTree,&injRef)); 4305 4306 PetscCall(DMPlexGetChart(fine,&pStartF,&pEndF)); 4307 PetscCall(DMGetLocalSection(fine,&localFine)); 4308 PetscCall(DMGetGlobalSection(fine,&globalFine)); 4309 PetscCall(PetscSectionGetNumFields(localFine,&numFields)); 4310 PetscCall(DMPlexGetChart(coarse,&pStartC,&pEndC)); 4311 PetscCall(DMGetLocalSection(coarse,&localCoarse)); 4312 PetscCall(DMGetGlobalSection(coarse,&globalCoarse)); 4313 PetscCall(PetscSectionGetMaxDof(localCoarse,&maxDof)); 4314 { 4315 PetscInt maxFields = PetscMax(1,numFields) + 1; 4316 PetscCall(PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets)); 4317 } 4318 4319 PetscCall(DMPlexTransferInjectorTree(coarse,fine,coarseToFine,cids,vecFine,numFields,offsets,&multiRootSec,&rootIndicesSec,NULL,&rootValues)); 4320 4321 PetscCall(PetscMalloc2(maxDof,&parentIndices,maxDof,&parentValues)); 4322 4323 /* count indices */ 4324 PetscCall(VecGetLayout(vecFine,&colMap)); 4325 PetscCall(VecGetLayout(vecCoarse,&rowMap)); 4326 PetscCall(PetscLayoutSetUp(rowMap)); 4327 PetscCall(PetscLayoutSetUp(colMap)); 4328 PetscCall(PetscLayoutGetRange(rowMap,&rowStart,&rowEnd)); 4329 PetscCall(PetscLayoutGetRange(colMap,&colStart,&colEnd)); 4330 /* insert values */ 4331 PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats)); 4332 for (p = pStartC; p < pEndC; p++) { 4333 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 4334 PetscBool contribute = PETSC_FALSE; 4335 4336 PetscCall(PetscSectionGetDof(globalCoarse,p,&dof)); 4337 PetscCall(PetscSectionGetConstraintDof(globalCoarse,p,&cdof)); 4338 if ((dof - cdof) <= 0) continue; 4339 PetscCall(PetscSectionGetDof(localCoarse,p,&dof)); 4340 PetscCall(PetscSectionGetOffset(globalCoarse,p,&gOff)); 4341 4342 rowOffsets[0] = 0; 4343 offsetsCopy[0] = 0; 4344 if (numFields) { 4345 PetscInt f; 4346 4347 for (f = 0; f < numFields; f++) { 4348 PetscInt fDof; 4349 PetscCall(PetscSectionGetFieldDof(localCoarse,p,f,&fDof)); 4350 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 4351 } 4352 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,NULL,parentIndices)); 4353 } else { 4354 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse,PETSC_FALSE,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,NULL,parentIndices)); 4355 rowOffsets[1] = offsetsCopy[0]; 4356 } 4357 4358 PetscCall(PetscSectionGetDof(multiRootSec,p,&numLeaves)); 4359 PetscCall(PetscSectionGetOffset(multiRootSec,p,&leafStart)); 4360 leafEnd = leafStart + numLeaves; 4361 for (l = 0; l < dof; l++) parentValues[l] = 0.; 4362 for (l = leafStart; l < leafEnd; l++) { 4363 PetscInt numIndices, childId, offset; 4364 const PetscScalar *childValues; 4365 4366 PetscCall(PetscSectionGetDof(rootIndicesSec,l,&numIndices)); 4367 PetscCall(PetscSectionGetOffset(rootIndicesSec,l,&offset)); 4368 childId = (PetscInt) PetscRealPart(rootValues[offset++]); 4369 childValues = &rootValues[offset]; 4370 numIndices--; 4371 4372 if (childId == -2) { /* skip */ 4373 continue; 4374 } else if (childId == -1) { /* equivalent points: scatter */ 4375 PetscInt m; 4376 4377 contribute = PETSC_TRUE; 4378 for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m]; 4379 } else { /* contributions from children: sum with injectors from reference tree */ 4380 PetscInt parentId, f, lim; 4381 4382 contribute = PETSC_TRUE; 4383 PetscCall(DMPlexGetTreeParent(refTree,childId,&parentId,NULL)); 4384 4385 lim = PetscMax(1,numFields); 4386 offsets[0] = 0; 4387 if (numFields) { 4388 PetscInt f; 4389 4390 for (f = 0; f < numFields; f++) { 4391 PetscInt fDof; 4392 PetscCall(PetscSectionGetFieldDof(cSecRef,childId,f,&fDof)); 4393 4394 offsets[f + 1] = fDof + offsets[f]; 4395 } 4396 } 4397 else { 4398 PetscInt cDof; 4399 4400 PetscCall(PetscSectionGetDof(cSecRef,childId,&cDof)); 4401 offsets[1] = cDof; 4402 } 4403 for (f = 0; f < lim; f++) { 4404 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 4405 PetscInt n = offsets[f+1]-offsets[f]; 4406 PetscInt m = rowOffsets[f+1]-rowOffsets[f]; 4407 PetscInt i, j; 4408 const PetscScalar *colValues = &childValues[offsets[f]]; 4409 4410 for (i = 0; i < m; i++) { 4411 PetscScalar val = 0.; 4412 for (j = 0; j < n; j++) { 4413 val += childMat[n * i + j] * colValues[j]; 4414 } 4415 parentValues[rowOffsets[f] + i] += val; 4416 } 4417 } 4418 } 4419 } 4420 if (contribute) PetscCall(VecSetValues(vecCoarse,dof,parentIndices,parentValues,INSERT_VALUES)); 4421 } 4422 PetscCall(PetscSectionDestroy(&multiRootSec)); 4423 PetscCall(PetscSectionDestroy(&rootIndicesSec)); 4424 PetscCall(PetscFree2(parentIndices,parentValues)); 4425 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats)); 4426 PetscCall(PetscFree(rootValues)); 4427 PetscCall(PetscFree3(offsets,offsetsCopy,rowOffsets)); 4428 PetscFunctionReturn(0); 4429 } 4430 4431 /*@ 4432 DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening 4433 that can be represented by a common reference tree used by both. This routine can be used for a combination of 4434 coarsening and refinement at the same time. 4435 4436 collective 4437 4438 Input Parameters: 4439 + dmIn - The DMPlex mesh for the input vector 4440 . vecIn - The input vector 4441 . sfRefine - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in 4442 the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn 4443 . sfCoarsen - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in 4444 the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn 4445 . cidsRefine - The childIds of the points in dmOut. These childIds relate back to the reference tree: childid[j] = k implies 4446 that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference 4447 tree was refined from its parent. childid[j] = -1 indicates that the point j in dmOut is exactly 4448 equivalent to its root in dmIn, so no interpolation is necessary. childid[j] = -2 indicates that this 4449 point j in dmOut is not a leaf of sfRefine. 4450 . cidsCoarsen - The childIds of the points in dmIn. These childIds relate back to the reference tree: childid[j] = k implies 4451 that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference 4452 tree coarsens to its parent. childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen. 4453 . useBCs - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer. 4454 - time - Used if boundary values are time dependent. 4455 4456 Output Parameters: 4457 . vecOut - Using interpolation and injection operators calculated on the reference tree, the transferred 4458 projection of vecIn from dmIn to dmOut. Note that any field discretized with a PetscFV finite volume 4459 method that uses gradient reconstruction will use reconstructed gradients when interpolating from 4460 coarse points to fine points. 4461 4462 Level: developer 4463 4464 .seealso: `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`, `PetscFVGetComputeGradients()` 4465 @*/ 4466 PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time) 4467 { 4468 PetscFunctionBegin; 4469 PetscCall(VecSet(vecOut,0.0)); 4470 if (sfRefine) { 4471 Vec vecInLocal; 4472 DM dmGrad = NULL; 4473 Vec faceGeom = NULL, cellGeom = NULL, grad = NULL; 4474 4475 PetscCall(DMGetLocalVector(dmIn,&vecInLocal)); 4476 PetscCall(VecSet(vecInLocal,0.0)); 4477 { 4478 PetscInt numFields, i; 4479 4480 PetscCall(DMGetNumFields(dmIn, &numFields)); 4481 for (i = 0; i < numFields; i++) { 4482 PetscObject obj; 4483 PetscClassId classid; 4484 4485 PetscCall(DMGetField(dmIn, i, NULL, &obj)); 4486 PetscCall(PetscObjectGetClassId(obj, &classid)); 4487 if (classid == PETSCFV_CLASSID) { 4488 PetscCall(DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad)); 4489 break; 4490 } 4491 } 4492 } 4493 if (useBCs) PetscCall(DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL)); 4494 PetscCall(DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal)); 4495 PetscCall(DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal)); 4496 if (dmGrad) { 4497 PetscCall(DMGetGlobalVector(dmGrad,&grad)); 4498 PetscCall(DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad)); 4499 } 4500 PetscCall(DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfRefine,cidsRefine,grad,cellGeom)); 4501 PetscCall(DMRestoreLocalVector(dmIn,&vecInLocal)); 4502 if (dmGrad) { 4503 PetscCall(DMRestoreGlobalVector(dmGrad,&grad)); 4504 } 4505 } 4506 if (sfCoarsen) PetscCall(DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfCoarsen,cidsCoarsen)); 4507 PetscCall(VecAssemblyBegin(vecOut)); 4508 PetscCall(VecAssemblyEnd(vecOut)); 4509 PetscFunctionReturn(0); 4510 } 4511