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