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 = 0, 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 2266 if (matSize == 0) { /* don't need to calculate the mat, just the indices */ 2267 PetscInt numIndices, *indices; 2268 PetscCall(DMPlexGetClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL)); 2269 PetscCheck(numIndices == numColIndices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "mismatching constraint indices calculations"); 2270 for (i = 0; i < numColIndices; i++) pInd[i] = indices[i]; 2271 for (i = 0; i < numFields; i++) { 2272 pInd[numColIndices + i] = offsets[i + 1]; 2273 pInd[numColIndices + numFields + i] = offsets[i + 1]; 2274 } 2275 PetscCall(DMPlexRestoreClosureIndices(coarse, localCoarse, globalCoarse, p, PETSC_TRUE, &numIndices, &indices, offsets, NULL)); 2276 } else { 2277 PetscInt closureSize, *closure = NULL, cl; 2278 PetscScalar *pMatIn, *pMatModified; 2279 PetscInt numPoints, *points; 2280 2281 { 2282 PetscInt *closure = NULL, closureSize, cl; 2283 2284 PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure)); 2285 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 2286 PetscInt c = closure[2 * cl], clDof; 2287 2288 PetscCall(PetscSectionGetDof(localCoarse, c, &clDof)); 2289 numRowIndices += clDof; 2290 } 2291 PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure)); 2292 } 2293 2294 PetscCall(DMGetWorkArray(coarse, numRowIndices * numRowIndices, MPIU_SCALAR, &pMatIn)); 2295 for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */ 2296 for (j = 0; j < numRowIndices; j++) pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.; 2297 } 2298 PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure)); 2299 for (f = 0; f < maxFields; f++) { 2300 if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f])); 2301 else PetscCall(PetscSectionGetPointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f])); 2302 } 2303 if (numFields) { 2304 for (cl = 0; cl < closureSize; cl++) { 2305 PetscInt c = closure[2 * cl]; 2306 2307 for (f = 0; f < numFields; f++) { 2308 PetscInt fDof; 2309 2310 PetscCall(PetscSectionGetFieldDof(localCoarse, c, f, &fDof)); 2311 offsets[f + 1] += fDof; 2312 } 2313 } 2314 for (f = 0; f < numFields; f++) { 2315 offsets[f + 1] += offsets[f]; 2316 newOffsets[f + 1] = offsets[f + 1]; 2317 } 2318 } 2319 /* TODO : flips here ? */ 2320 /* apply hanging node constraints on the right, get the new points and the new offsets */ 2321 PetscCall(DMPlexAnchorsModifyMat(coarse, localCoarse, closureSize, numRowIndices, closure, perms, pMatIn, &numPoints, NULL, &points, &pMatModified, newOffsets, PETSC_FALSE)); 2322 for (f = 0; f < maxFields; f++) { 2323 if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse, f, closureSize, closure, &perms[f], &flips[f])); 2324 else PetscCall(PetscSectionRestorePointSyms(localCoarse, closureSize, closure, &perms[f], &flips[f])); 2325 } 2326 for (f = 0; f < maxFields; f++) { 2327 if (numFields) PetscCall(PetscSectionGetFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f])); 2328 else PetscCall(PetscSectionGetPointSyms(localCoarse, numPoints, points, &perms[f], &flips[f])); 2329 } 2330 if (!numFields) { 2331 for (i = 0; i < numRowIndices * numColIndices; i++) pMat[i] = pMatModified[i]; 2332 } else { 2333 PetscInt i, j, count; 2334 for (f = 0, count = 0; f < numFields; f++) { 2335 for (i = offsets[f]; i < offsets[f + 1]; i++) { 2336 for (j = newOffsets[f]; j < newOffsets[f + 1]; j++, count++) pMat[count] = pMatModified[i * numColIndices + j]; 2337 } 2338 } 2339 } 2340 PetscCall(DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatModified)); 2341 PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure)); 2342 PetscCall(DMRestoreWorkArray(coarse, numRowIndices * numColIndices, MPIU_SCALAR, &pMatIn)); 2343 if (numFields) { 2344 for (f = 0; f < numFields; f++) { 2345 pInd[numColIndices + f] = offsets[f + 1]; 2346 pInd[numColIndices + numFields + f] = newOffsets[f + 1]; 2347 } 2348 for (cl = 0; cl < numPoints; cl++) { 2349 PetscInt globalOff, c = points[2 * cl]; 2350 PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff)); 2351 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, NULL, pInd)); 2352 } 2353 } else { 2354 for (cl = 0; cl < numPoints; cl++) { 2355 PetscInt c = points[2 * cl], globalOff; 2356 const PetscInt *perm = perms[0] ? perms[0][cl] : NULL; 2357 2358 PetscCall(PetscSectionGetOffset(globalCoarse, c, &globalOff)); 2359 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, c, globalOff < 0 ? -(globalOff + 1) : globalOff, newOffsets, PETSC_FALSE, perm, NULL, pInd)); 2360 } 2361 } 2362 for (f = 0; f < maxFields; f++) { 2363 if (numFields) PetscCall(PetscSectionRestoreFieldPointSyms(localCoarse, f, numPoints, points, &perms[f], &flips[f])); 2364 else PetscCall(PetscSectionRestorePointSyms(localCoarse, numPoints, points, &perms[f], &flips[f])); 2365 } 2366 PetscCall(DMRestoreWorkArray(coarse, numPoints, MPIU_SCALAR, &points)); 2367 } 2368 } else if (matSize) { 2369 PetscInt cOff; 2370 PetscInt *rowIndices, *colIndices, a, aDof = 0, aOff; 2371 2372 numRowIndices = dof; 2373 PetscCall(DMGetWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices)); 2374 PetscCall(DMGetWorkArray(coarse, numColIndices, MPIU_INT, &colIndices)); 2375 PetscCall(PetscSectionGetOffset(cSec, p, &cOff)); 2376 PetscCall(PetscSectionGetDof(aSec, p, &aDof)); 2377 PetscCall(PetscSectionGetOffset(aSec, p, &aOff)); 2378 if (numFields) { 2379 for (f = 0; f < numFields; f++) { 2380 PetscInt fDof; 2381 2382 PetscCall(PetscSectionGetFieldDof(cSec, p, f, &fDof)); 2383 offsets[f + 1] = fDof; 2384 for (a = 0; a < aDof; a++) { 2385 PetscInt anchor = anchors[a + aOff]; 2386 PetscCall(PetscSectionGetFieldDof(localCoarse, anchor, f, &fDof)); 2387 newOffsets[f + 1] += fDof; 2388 } 2389 } 2390 for (f = 0; f < numFields; f++) { 2391 offsets[f + 1] += offsets[f]; 2392 offsetsCopy[f + 1] = offsets[f + 1]; 2393 newOffsets[f + 1] += newOffsets[f]; 2394 newOffsetsCopy[f + 1] = newOffsets[f + 1]; 2395 } 2396 PetscCall(DMPlexGetIndicesPointFields_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, -1, NULL, rowIndices)); 2397 for (a = 0; a < aDof; a++) { 2398 PetscInt anchor = anchors[a + aOff], lOff; 2399 PetscCall(PetscSectionGetOffset(localCoarse, anchor, &lOff)); 2400 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, -1, NULL, colIndices)); 2401 } 2402 } else { 2403 PetscCall(DMPlexGetIndicesPoint_Internal(cSec, PETSC_TRUE, p, cOff, offsetsCopy, PETSC_TRUE, NULL, NULL, rowIndices)); 2404 for (a = 0; a < aDof; a++) { 2405 PetscInt anchor = anchors[a + aOff], lOff; 2406 PetscCall(PetscSectionGetOffset(localCoarse, anchor, &lOff)); 2407 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_TRUE, anchor, lOff, newOffsetsCopy, PETSC_TRUE, NULL, NULL, colIndices)); 2408 } 2409 } 2410 if (numFields) { 2411 PetscInt count, a; 2412 2413 for (f = 0, count = 0; f < numFields; f++) { 2414 PetscInt iSize = offsets[f + 1] - offsets[f]; 2415 PetscInt jSize = newOffsets[f + 1] - newOffsets[f]; 2416 PetscCall(MatGetValues(cMat, iSize, &rowIndices[offsets[f]], jSize, &colIndices[newOffsets[f]], &pMat[count])); 2417 count += iSize * jSize; 2418 pInd[numColIndices + f] = offsets[f + 1]; 2419 pInd[numColIndices + numFields + f] = newOffsets[f + 1]; 2420 } 2421 for (a = 0; a < aDof; a++) { 2422 PetscInt anchor = anchors[a + aOff]; 2423 PetscInt gOff; 2424 PetscCall(PetscSectionGetOffset(globalCoarse, anchor, &gOff)); 2425 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, -1, NULL, pInd)); 2426 } 2427 } else { 2428 PetscInt a; 2429 PetscCall(MatGetValues(cMat, numRowIndices, rowIndices, numColIndices, colIndices, pMat)); 2430 for (a = 0; a < aDof; a++) { 2431 PetscInt anchor = anchors[a + aOff]; 2432 PetscInt gOff; 2433 PetscCall(PetscSectionGetOffset(globalCoarse, anchor, &gOff)); 2434 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, anchor, gOff < 0 ? -(gOff + 1) : gOff, newOffsets, PETSC_FALSE, NULL, NULL, pInd)); 2435 } 2436 } 2437 PetscCall(DMRestoreWorkArray(coarse, numColIndices, MPIU_INT, &colIndices)); 2438 PetscCall(DMRestoreWorkArray(coarse, numRowIndices, MPIU_INT, &rowIndices)); 2439 } else { 2440 PetscInt gOff; 2441 2442 PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff)); 2443 if (numFields) { 2444 for (f = 0; f < numFields; f++) { 2445 PetscInt fDof; 2446 PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof)); 2447 offsets[f + 1] = fDof + offsets[f]; 2448 } 2449 for (f = 0; f < numFields; f++) { 2450 pInd[numColIndices + f] = offsets[f + 1]; 2451 pInd[numColIndices + numFields + f] = offsets[f + 1]; 2452 } 2453 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd)); 2454 } else { 2455 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd)); 2456 } 2457 } 2458 } 2459 PetscCall(PetscFree(maxChildIds)); 2460 } 2461 { 2462 PetscSF indicesSF, matricesSF; 2463 PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices; 2464 2465 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec)); 2466 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafMatricesSec)); 2467 PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootIndicesSec, &remoteOffsetsIndices, leafIndicesSec)); 2468 PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootMatricesSec, &remoteOffsetsMatrices, leafMatricesSec)); 2469 PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootIndicesSec, remoteOffsetsIndices, leafIndicesSec, &indicesSF)); 2470 PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootMatricesSec, remoteOffsetsMatrices, leafMatricesSec, &matricesSF)); 2471 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 2472 PetscCall(PetscFree(remoteOffsetsIndices)); 2473 PetscCall(PetscFree(remoteOffsetsMatrices)); 2474 PetscCall(PetscSectionGetStorageSize(leafIndicesSec, &numLeafIndices)); 2475 PetscCall(PetscSectionGetStorageSize(leafMatricesSec, &numLeafMatrices)); 2476 PetscCall(PetscMalloc2(numLeafIndices, &leafIndices, numLeafMatrices, &leafMatrices)); 2477 PetscCall(PetscSFBcastBegin(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE)); 2478 PetscCall(PetscSFBcastBegin(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE)); 2479 PetscCall(PetscSFBcastEnd(indicesSF, MPIU_INT, rootIndices, leafIndices, MPI_REPLACE)); 2480 PetscCall(PetscSFBcastEnd(matricesSF, MPIU_SCALAR, rootMatrices, leafMatrices, MPI_REPLACE)); 2481 PetscCall(PetscSFDestroy(&matricesSF)); 2482 PetscCall(PetscSFDestroy(&indicesSF)); 2483 PetscCall(PetscFree2(rootIndices, rootMatrices)); 2484 PetscCall(PetscSectionDestroy(&rootIndicesSec)); 2485 PetscCall(PetscSectionDestroy(&rootMatricesSec)); 2486 } 2487 /* count to preallocate */ 2488 PetscCall(DMGetLocalSection(fine, &localFine)); 2489 { 2490 PetscInt nGlobal; 2491 PetscInt *dnnz, *onnz; 2492 PetscLayout rowMap, colMap; 2493 PetscInt rowStart, rowEnd, colStart, colEnd; 2494 PetscInt maxDof; 2495 PetscInt *rowIndices; 2496 DM refTree; 2497 PetscInt **refPointFieldN; 2498 PetscScalar ***refPointFieldMats; 2499 PetscSection refConSec, refAnSec; 2500 PetscInt pRefStart, pRefEnd, maxConDof, maxColumns, leafStart, leafEnd; 2501 PetscScalar *pointWork; 2502 2503 PetscCall(PetscSectionGetConstrainedStorageSize(globalFine, &nGlobal)); 2504 PetscCall(PetscCalloc2(nGlobal, &dnnz, nGlobal, &onnz)); 2505 PetscCall(MatGetLayouts(mat, &rowMap, &colMap)); 2506 PetscCall(PetscLayoutSetUp(rowMap)); 2507 PetscCall(PetscLayoutSetUp(colMap)); 2508 PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd)); 2509 PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd)); 2510 PetscCall(PetscSectionGetMaxDof(localFine, &maxDof)); 2511 PetscCall(PetscSectionGetChart(leafIndicesSec, &leafStart, &leafEnd)); 2512 PetscCall(DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices)); 2513 for (p = leafStart; p < leafEnd; p++) { 2514 PetscInt gDof, gcDof, gOff; 2515 PetscInt numColIndices, pIndOff, *pInd; 2516 PetscInt matSize; 2517 PetscInt i; 2518 2519 PetscCall(PetscSectionGetDof(globalFine, p, &gDof)); 2520 PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof)); 2521 if ((gDof - gcDof) <= 0) continue; 2522 PetscCall(PetscSectionGetOffset(globalFine, p, &gOff)); 2523 PetscCheck(gOff >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "I though having global dofs meant a non-negative offset"); 2524 PetscCheck(gOff >= rowStart && (gOff + gDof - gcDof) <= rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "I thought the row map would constrain the global dofs"); 2525 PetscCall(PetscSectionGetDof(leafIndicesSec, p, &numColIndices)); 2526 PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &pIndOff)); 2527 numColIndices -= 2 * numFields; 2528 PetscCheck(numColIndices > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "global fine dof with no dofs to interpolate from"); 2529 pInd = &leafIndices[pIndOff]; 2530 offsets[0] = 0; 2531 offsetsCopy[0] = 0; 2532 newOffsets[0] = 0; 2533 newOffsetsCopy[0] = 0; 2534 if (numFields) { 2535 PetscInt f; 2536 for (f = 0; f < numFields; f++) { 2537 PetscInt rowDof; 2538 2539 PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof)); 2540 offsets[f + 1] = offsets[f] + rowDof; 2541 offsetsCopy[f + 1] = offsets[f + 1]; 2542 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2543 numD[f] = 0; 2544 numO[f] = 0; 2545 } 2546 PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices)); 2547 for (f = 0; f < numFields; f++) { 2548 PetscInt colOffset = newOffsets[f]; 2549 PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f]; 2550 2551 for (i = 0; i < numFieldCols; i++) { 2552 PetscInt gInd = pInd[i + colOffset]; 2553 2554 if (gInd >= colStart && gInd < colEnd) { 2555 numD[f]++; 2556 } else if (gInd >= 0) { /* negative means non-entry */ 2557 numO[f]++; 2558 } 2559 } 2560 } 2561 } else { 2562 PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices)); 2563 numD[0] = 0; 2564 numO[0] = 0; 2565 for (i = 0; i < numColIndices; i++) { 2566 PetscInt gInd = pInd[i]; 2567 2568 if (gInd >= colStart && gInd < colEnd) { 2569 numD[0]++; 2570 } else if (gInd >= 0) { /* negative means non-entry */ 2571 numO[0]++; 2572 } 2573 } 2574 } 2575 PetscCall(PetscSectionGetDof(leafMatricesSec, p, &matSize)); 2576 if (!matSize) { /* incoming matrix is identity */ 2577 PetscInt childId; 2578 2579 childId = childIds[p - pStartF]; 2580 if (childId < 0) { /* no child interpolation: one nnz per */ 2581 if (numFields) { 2582 PetscInt f; 2583 for (f = 0; f < numFields; f++) { 2584 PetscInt numRows = offsets[f + 1] - offsets[f], row; 2585 for (row = 0; row < numRows; row++) { 2586 PetscInt gIndCoarse = pInd[newOffsets[f] + row]; 2587 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2588 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2589 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs"); 2590 dnnz[gIndFine - rowStart] = 1; 2591 } else if (gIndCoarse >= 0) { /* remote */ 2592 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs"); 2593 onnz[gIndFine - rowStart] = 1; 2594 } else { /* constrained */ 2595 PetscCheck(gIndFine < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs"); 2596 } 2597 } 2598 } 2599 } else { 2600 PetscInt i; 2601 for (i = 0; i < gDof; i++) { 2602 PetscInt gIndCoarse = pInd[i]; 2603 PetscInt gIndFine = rowIndices[i]; 2604 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2605 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs"); 2606 dnnz[gIndFine - rowStart] = 1; 2607 } else if (gIndCoarse >= 0) { /* remote */ 2608 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs"); 2609 onnz[gIndFine - rowStart] = 1; 2610 } else { /* constrained */ 2611 PetscCheck(gIndFine < 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs"); 2612 } 2613 } 2614 } 2615 } else { /* interpolate from all */ 2616 if (numFields) { 2617 PetscInt f; 2618 for (f = 0; f < numFields; f++) { 2619 PetscInt numRows = offsets[f + 1] - offsets[f], row; 2620 for (row = 0; row < numRows; row++) { 2621 PetscInt gIndFine = rowIndices[offsets[f] + row]; 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[f]; 2625 onnz[gIndFine - rowStart] = numO[f]; 2626 } 2627 } 2628 } 2629 } else { 2630 PetscInt i; 2631 for (i = 0; i < gDof; i++) { 2632 PetscInt gIndFine = rowIndices[i]; 2633 if (gIndFine >= 0) { 2634 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs"); 2635 dnnz[gIndFine - rowStart] = numD[0]; 2636 onnz[gIndFine - rowStart] = numO[0]; 2637 } 2638 } 2639 } 2640 } 2641 } else { /* interpolate from all */ 2642 if (numFields) { 2643 PetscInt f; 2644 for (f = 0; f < numFields; f++) { 2645 PetscInt numRows = offsets[f + 1] - offsets[f], row; 2646 for (row = 0; row < numRows; row++) { 2647 PetscInt gIndFine = rowIndices[offsets[f] + row]; 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[f]; 2651 onnz[gIndFine - rowStart] = numO[f]; 2652 } 2653 } 2654 } 2655 } else { /* every dof get a full row */ 2656 PetscInt i; 2657 for (i = 0; i < gDof; i++) { 2658 PetscInt gIndFine = rowIndices[i]; 2659 if (gIndFine >= 0) { 2660 PetscCheck(gIndFine >= rowStart && gIndFine < rowEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mismatched number of constrained dofs"); 2661 dnnz[gIndFine - rowStart] = numD[0]; 2662 onnz[gIndFine - rowStart] = numO[0]; 2663 } 2664 } 2665 } 2666 } 2667 } 2668 PetscCall(MatXAIJSetPreallocation(mat, 1, dnnz, onnz, NULL, NULL)); 2669 PetscCall(PetscFree2(dnnz, onnz)); 2670 2671 PetscCall(DMPlexGetReferenceTree(fine, &refTree)); 2672 PetscCall(DMCopyDisc(fine, refTree)); 2673 PetscCall(DMSetLocalSection(refTree, NULL)); 2674 PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL)); 2675 PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN)); 2676 PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL)); 2677 PetscCall(DMPlexGetAnchors(refTree, &refAnSec, NULL)); 2678 PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd)); 2679 PetscCall(PetscSectionGetMaxDof(refConSec, &maxConDof)); 2680 PetscCall(PetscSectionGetMaxDof(leafIndicesSec, &maxColumns)); 2681 PetscCall(PetscMalloc1(maxConDof * maxColumns, &pointWork)); 2682 for (p = leafStart; p < leafEnd; p++) { 2683 PetscInt gDof, gcDof, gOff; 2684 PetscInt numColIndices, pIndOff, *pInd; 2685 PetscInt matSize; 2686 PetscInt childId; 2687 2688 PetscCall(PetscSectionGetDof(globalFine, p, &gDof)); 2689 PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof)); 2690 if ((gDof - gcDof) <= 0) continue; 2691 childId = childIds[p - pStartF]; 2692 PetscCall(PetscSectionGetOffset(globalFine, p, &gOff)); 2693 PetscCall(PetscSectionGetDof(leafIndicesSec, p, &numColIndices)); 2694 PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &pIndOff)); 2695 numColIndices -= 2 * numFields; 2696 pInd = &leafIndices[pIndOff]; 2697 offsets[0] = 0; 2698 offsetsCopy[0] = 0; 2699 newOffsets[0] = 0; 2700 newOffsetsCopy[0] = 0; 2701 rowOffsets[0] = 0; 2702 if (numFields) { 2703 PetscInt f; 2704 for (f = 0; f < numFields; f++) { 2705 PetscInt rowDof; 2706 2707 PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof)); 2708 offsets[f + 1] = offsets[f] + rowDof; 2709 offsetsCopy[f + 1] = offsets[f + 1]; 2710 rowOffsets[f + 1] = pInd[numColIndices + f]; 2711 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2712 } 2713 PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices)); 2714 } else { 2715 PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices)); 2716 } 2717 PetscCall(PetscSectionGetDof(leafMatricesSec, p, &matSize)); 2718 if (!matSize) { /* incoming matrix is identity */ 2719 if (childId < 0) { /* no child interpolation: scatter */ 2720 if (numFields) { 2721 PetscInt f; 2722 for (f = 0; f < numFields; f++) { 2723 PetscInt numRows = offsets[f + 1] - offsets[f], row; 2724 for (row = 0; row < numRows; row++) PetscCall(MatSetValue(mat, rowIndices[offsets[f] + row], pInd[newOffsets[f] + row], 1., INSERT_VALUES)); 2725 } 2726 } else { 2727 PetscInt numRows = gDof, row; 2728 for (row = 0; row < numRows; row++) PetscCall(MatSetValue(mat, rowIndices[row], pInd[row], 1., INSERT_VALUES)); 2729 } 2730 } else { /* interpolate from all */ 2731 if (numFields) { 2732 PetscInt f; 2733 for (f = 0; f < numFields; f++) { 2734 PetscInt numRows = offsets[f + 1] - offsets[f]; 2735 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 2736 PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], refPointFieldMats[childId - pRefStart][f], INSERT_VALUES)); 2737 } 2738 } else { 2739 PetscCall(MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, refPointFieldMats[childId - pRefStart][0], INSERT_VALUES)); 2740 } 2741 } 2742 } else { /* interpolate from all */ 2743 PetscInt pMatOff; 2744 PetscScalar *pMat; 2745 2746 PetscCall(PetscSectionGetOffset(leafMatricesSec, p, &pMatOff)); 2747 pMat = &leafMatrices[pMatOff]; 2748 if (childId < 0) { /* copy the incoming matrix */ 2749 if (numFields) { 2750 PetscInt f, count; 2751 for (f = 0, count = 0; f < numFields; f++) { 2752 PetscInt numRows = offsets[f + 1] - offsets[f]; 2753 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 2754 PetscInt numInRows = rowOffsets[f + 1] - rowOffsets[f]; 2755 PetscScalar *inMat = &pMat[count]; 2756 2757 PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], inMat, INSERT_VALUES)); 2758 count += numCols * numInRows; 2759 } 2760 } else { 2761 PetscCall(MatSetValues(mat, gDof, rowIndices, numColIndices, pInd, pMat, INSERT_VALUES)); 2762 } 2763 } else { /* multiply the incoming matrix by the child interpolation */ 2764 if (numFields) { 2765 PetscInt f, count; 2766 for (f = 0, count = 0; f < numFields; f++) { 2767 PetscInt numRows = offsets[f + 1] - offsets[f]; 2768 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 2769 PetscInt numInRows = rowOffsets[f + 1] - rowOffsets[f]; 2770 PetscScalar *inMat = &pMat[count]; 2771 PetscInt i, j, k; 2772 PetscCheck(refPointFieldN[childId - pRefStart][f] == numInRows, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point constraint matrix multiply dimension mismatch"); 2773 for (i = 0; i < numRows; i++) { 2774 for (j = 0; j < numCols; j++) { 2775 PetscScalar val = 0.; 2776 for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j]; 2777 pointWork[i * numCols + j] = val; 2778 } 2779 } 2780 PetscCall(MatSetValues(mat, numRows, &rowIndices[offsets[f]], numCols, &pInd[newOffsets[f]], pointWork, INSERT_VALUES)); 2781 count += numCols * numInRows; 2782 } 2783 } else { /* every dof gets a full row */ 2784 PetscInt numRows = gDof; 2785 PetscInt numCols = numColIndices; 2786 PetscInt numInRows = matSize / numColIndices; 2787 PetscInt i, j, k; 2788 PetscCheck(refPointFieldN[childId - pRefStart][0] == numInRows, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point constraint matrix multiply dimension mismatch"); 2789 for (i = 0; i < numRows; i++) { 2790 for (j = 0; j < numCols; j++) { 2791 PetscScalar val = 0.; 2792 for (k = 0; k < numInRows; k++) val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j]; 2793 pointWork[i * numCols + j] = val; 2794 } 2795 } 2796 PetscCall(MatSetValues(mat, numRows, rowIndices, numCols, pInd, pointWork, INSERT_VALUES)); 2797 } 2798 } 2799 } 2800 } 2801 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN)); 2802 PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices)); 2803 PetscCall(PetscFree(pointWork)); 2804 } 2805 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 2806 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 2807 PetscCall(PetscSectionDestroy(&leafIndicesSec)); 2808 PetscCall(PetscSectionDestroy(&leafMatricesSec)); 2809 PetscCall(PetscFree2(leafIndices, leafMatrices)); 2810 PetscCall(PetscFree2(*(PetscInt ****)&perms, *(PetscScalar ****)&flips)); 2811 PetscCall(PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO)); 2812 PetscCall(ISRestoreIndices(aIS, &anchors)); 2813 PetscFunctionReturn(PETSC_SUCCESS); 2814 } 2815 2816 /* 2817 * Assuming a nodal basis (w.r.t. the dual basis) basis: 2818 * 2819 * for each coarse dof \phi^c_i: 2820 * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i: 2821 * for each fine dof \phi^f_j; 2822 * a_{i,j} = 0; 2823 * for each fine dof \phi^f_k: 2824 * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l 2825 * [^^^ this is = \phi^c_i ^^^] 2826 */ 2827 PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj) 2828 { 2829 PetscDS ds; 2830 PetscSection section, cSection; 2831 DMLabel canonical, depth; 2832 Mat cMat, mat; 2833 PetscInt *nnz; 2834 PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd; 2835 PetscInt m, n; 2836 PetscScalar *pointScalar; 2837 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent; 2838 2839 PetscFunctionBegin; 2840 PetscCall(DMGetLocalSection(refTree, §ion)); 2841 PetscCall(DMGetDimension(refTree, &dim)); 2842 PetscCall(PetscMalloc6(dim, &v0, dim, &v0parent, dim, &vtmp, dim * dim, &J, dim * dim, &Jparent, dim * dim, &invJ)); 2843 PetscCall(PetscMalloc2(dim, &pointScalar, dim, &pointRef)); 2844 PetscCall(DMGetDS(refTree, &ds)); 2845 PetscCall(PetscDSGetNumFields(ds, &numFields)); 2846 PetscCall(PetscSectionGetNumFields(section, &numSecFields)); 2847 PetscCall(DMGetLabel(refTree, "canonical", &canonical)); 2848 PetscCall(DMGetLabel(refTree, "depth", &depth)); 2849 PetscCall(DMGetDefaultConstraints(refTree, &cSection, &cMat, NULL)); 2850 PetscCall(DMPlexGetChart(refTree, &pStart, &pEnd)); 2851 PetscCall(DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd)); 2852 PetscCall(MatGetSize(cMat, &n, &m)); /* the injector has transpose sizes from the constraint matrix */ 2853 /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */ 2854 PetscCall(PetscCalloc1(m, &nnz)); 2855 for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */ 2856 const PetscInt *children; 2857 PetscInt numChildren; 2858 PetscInt i, numChildDof, numSelfDof; 2859 2860 if (canonical) { 2861 PetscInt pCanonical; 2862 PetscCall(DMLabelGetValue(canonical, p, &pCanonical)); 2863 if (p != pCanonical) continue; 2864 } 2865 PetscCall(DMPlexGetTreeChildren(refTree, p, &numChildren, &children)); 2866 if (!numChildren) continue; 2867 for (i = 0, numChildDof = 0; i < numChildren; i++) { 2868 PetscInt child = children[i]; 2869 PetscInt dof; 2870 2871 PetscCall(PetscSectionGetDof(section, child, &dof)); 2872 numChildDof += dof; 2873 } 2874 PetscCall(PetscSectionGetDof(section, p, &numSelfDof)); 2875 if (!numChildDof || !numSelfDof) continue; 2876 for (f = 0; f < numFields; f++) { 2877 PetscInt selfOff; 2878 2879 if (numSecFields) { /* count the dofs for just this field */ 2880 for (i = 0, numChildDof = 0; i < numChildren; i++) { 2881 PetscInt child = children[i]; 2882 PetscInt dof; 2883 2884 PetscCall(PetscSectionGetFieldDof(section, child, f, &dof)); 2885 numChildDof += dof; 2886 } 2887 PetscCall(PetscSectionGetFieldDof(section, p, f, &numSelfDof)); 2888 PetscCall(PetscSectionGetFieldOffset(section, p, f, &selfOff)); 2889 } else { 2890 PetscCall(PetscSectionGetOffset(section, p, &selfOff)); 2891 } 2892 for (i = 0; i < numSelfDof; i++) nnz[selfOff + i] = numChildDof; 2893 } 2894 } 2895 PetscCall(MatCreateAIJ(PETSC_COMM_SELF, m, n, m, n, -1, nnz, -1, NULL, &mat)); 2896 PetscCall(PetscFree(nnz)); 2897 /* Setp 2: compute entries */ 2898 for (p = pStart; p < pEnd; p++) { 2899 const PetscInt *children; 2900 PetscInt numChildren; 2901 PetscInt i, numChildDof, numSelfDof; 2902 2903 /* same conditions about when entries occur */ 2904 if (canonical) { 2905 PetscInt pCanonical; 2906 PetscCall(DMLabelGetValue(canonical, p, &pCanonical)); 2907 if (p != pCanonical) continue; 2908 } 2909 PetscCall(DMPlexGetTreeChildren(refTree, p, &numChildren, &children)); 2910 if (!numChildren) continue; 2911 for (i = 0, numChildDof = 0; i < numChildren; i++) { 2912 PetscInt child = children[i]; 2913 PetscInt dof; 2914 2915 PetscCall(PetscSectionGetDof(section, child, &dof)); 2916 numChildDof += dof; 2917 } 2918 PetscCall(PetscSectionGetDof(section, p, &numSelfDof)); 2919 if (!numChildDof || !numSelfDof) continue; 2920 2921 for (f = 0; f < numFields; f++) { 2922 PetscInt pI = -1, cI = -1; 2923 PetscInt selfOff, Nc, parentCell; 2924 PetscInt cellShapeOff; 2925 PetscObject disc; 2926 PetscDualSpace dsp; 2927 PetscClassId classId; 2928 PetscScalar *pointMat; 2929 PetscInt *matRows, *matCols; 2930 PetscInt pO = PETSC_MIN_INT; 2931 const PetscInt *depthNumDof; 2932 2933 if (numSecFields) { 2934 for (i = 0, numChildDof = 0; i < numChildren; i++) { 2935 PetscInt child = children[i]; 2936 PetscInt dof; 2937 2938 PetscCall(PetscSectionGetFieldDof(section, child, f, &dof)); 2939 numChildDof += dof; 2940 } 2941 PetscCall(PetscSectionGetFieldDof(section, p, f, &numSelfDof)); 2942 PetscCall(PetscSectionGetFieldOffset(section, p, f, &selfOff)); 2943 } else { 2944 PetscCall(PetscSectionGetOffset(section, p, &selfOff)); 2945 } 2946 2947 /* find a cell whose closure contains p */ 2948 if (p >= cStart && p < cEnd) { 2949 parentCell = p; 2950 } else { 2951 PetscInt *star = NULL; 2952 PetscInt numStar; 2953 2954 parentCell = -1; 2955 PetscCall(DMPlexGetTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star)); 2956 for (i = numStar - 1; i >= 0; i--) { 2957 PetscInt c = star[2 * i]; 2958 2959 if (c >= cStart && c < cEnd) { 2960 parentCell = c; 2961 break; 2962 } 2963 } 2964 PetscCall(DMPlexRestoreTransitiveClosure(refTree, p, PETSC_FALSE, &numStar, &star)); 2965 } 2966 /* determine the offset of p's shape functions within parentCell's shape functions */ 2967 PetscCall(PetscDSGetDiscretization(ds, f, &disc)); 2968 PetscCall(PetscObjectGetClassId(disc, &classId)); 2969 if (classId == PETSCFE_CLASSID) { 2970 PetscCall(PetscFEGetDualSpace((PetscFE)disc, &dsp)); 2971 } else if (classId == PETSCFV_CLASSID) { 2972 PetscCall(PetscFVGetDualSpace((PetscFV)disc, &dsp)); 2973 } else { 2974 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported discretization object"); 2975 } 2976 PetscCall(PetscDualSpaceGetNumDof(dsp, &depthNumDof)); 2977 PetscCall(PetscDualSpaceGetNumComponents(dsp, &Nc)); 2978 { 2979 PetscInt *closure = NULL; 2980 PetscInt numClosure; 2981 2982 PetscCall(DMPlexGetTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure)); 2983 for (i = 0, pI = -1, cellShapeOff = 0; i < numClosure; i++) { 2984 PetscInt point = closure[2 * i], pointDepth; 2985 2986 pO = closure[2 * i + 1]; 2987 if (point == p) { 2988 pI = i; 2989 break; 2990 } 2991 PetscCall(DMLabelGetValue(depth, point, &pointDepth)); 2992 cellShapeOff += depthNumDof[pointDepth]; 2993 } 2994 PetscCall(DMPlexRestoreTransitiveClosure(refTree, parentCell, PETSC_TRUE, &numClosure, &closure)); 2995 } 2996 2997 PetscCall(DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat)); 2998 PetscCall(DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows)); 2999 matCols = matRows + numSelfDof; 3000 for (i = 0; i < numSelfDof; i++) matRows[i] = selfOff + i; 3001 for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.; 3002 { 3003 PetscInt colOff = 0; 3004 3005 for (i = 0; i < numChildren; i++) { 3006 PetscInt child = children[i]; 3007 PetscInt dof, off, j; 3008 3009 if (numSecFields) { 3010 PetscCall(PetscSectionGetFieldDof(cSection, child, f, &dof)); 3011 PetscCall(PetscSectionGetFieldOffset(cSection, child, f, &off)); 3012 } else { 3013 PetscCall(PetscSectionGetDof(cSection, child, &dof)); 3014 PetscCall(PetscSectionGetOffset(cSection, child, &off)); 3015 } 3016 3017 for (j = 0; j < dof; j++) matCols[colOff++] = off + j; 3018 } 3019 } 3020 if (classId == PETSCFE_CLASSID) { 3021 PetscFE fe = (PetscFE)disc; 3022 PetscInt fSize; 3023 const PetscInt ***perms; 3024 const PetscScalar ***flips; 3025 const PetscInt *pperms; 3026 3027 PetscCall(PetscFEGetDualSpace(fe, &dsp)); 3028 PetscCall(PetscDualSpaceGetDimension(dsp, &fSize)); 3029 PetscCall(PetscDualSpaceGetSymmetries(dsp, &perms, &flips)); 3030 pperms = perms ? perms[pI] ? perms[pI][pO] : NULL : NULL; 3031 for (i = 0; i < numSelfDof; i++) { /* for every shape function */ 3032 PetscQuadrature q; 3033 PetscInt dim, thisNc, numPoints, j, k; 3034 const PetscReal *points; 3035 const PetscReal *weights; 3036 PetscInt *closure = NULL; 3037 PetscInt numClosure; 3038 PetscInt iCell = pperms ? pperms[i] : i; 3039 PetscInt parentCellShapeDof = cellShapeOff + iCell; 3040 PetscTabulation Tparent; 3041 3042 PetscCall(PetscDualSpaceGetFunctional(dsp, parentCellShapeDof, &q)); 3043 PetscCall(PetscQuadratureGetData(q, &dim, &thisNc, &numPoints, &points, &weights)); 3044 PetscCheck(thisNc == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Functional dim %" PetscInt_FMT " does not much basis dim %" PetscInt_FMT, thisNc, Nc); 3045 PetscCall(PetscFECreateTabulation(fe, 1, numPoints, points, 0, &Tparent)); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */ 3046 for (j = 0; j < numPoints; j++) { 3047 PetscInt childCell = -1; 3048 PetscReal *parentValAtPoint; 3049 const PetscReal xi0[3] = {-1., -1., -1.}; 3050 const PetscReal *pointReal = &points[dim * j]; 3051 const PetscScalar *point; 3052 PetscTabulation Tchild; 3053 PetscInt childCellShapeOff, pointMatOff; 3054 #if defined(PETSC_USE_COMPLEX) 3055 PetscInt d; 3056 3057 for (d = 0; d < dim; d++) pointScalar[d] = points[dim * j + d]; 3058 point = pointScalar; 3059 #else 3060 point = pointReal; 3061 #endif 3062 3063 parentValAtPoint = &Tparent->T[0][(fSize * j + parentCellShapeDof) * Nc]; 3064 3065 for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/ 3066 PetscInt child = children[k]; 3067 PetscInt *star = NULL; 3068 PetscInt numStar, s; 3069 3070 PetscCall(DMPlexGetTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star)); 3071 for (s = numStar - 1; s >= 0; s--) { 3072 PetscInt c = star[2 * s]; 3073 3074 if (c < cStart || c >= cEnd) continue; 3075 PetscCall(DMPlexLocatePoint_Internal(refTree, dim, point, c, &childCell)); 3076 if (childCell >= 0) break; 3077 } 3078 PetscCall(DMPlexRestoreTransitiveClosure(refTree, child, PETSC_FALSE, &numStar, &star)); 3079 if (childCell >= 0) break; 3080 } 3081 PetscCheck(childCell >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not locate quadrature point"); 3082 PetscCall(DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ)); 3083 PetscCall(DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent)); 3084 CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp); 3085 CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef); 3086 3087 PetscCall(PetscFECreateTabulation(fe, 1, 1, pointRef, 0, &Tchild)); 3088 PetscCall(DMPlexGetTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure)); 3089 for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */ 3090 PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT; 3091 PetscInt l; 3092 const PetscInt *cperms; 3093 3094 PetscCall(DMLabelGetValue(depth, child, &childDepth)); 3095 childDof = depthNumDof[childDepth]; 3096 for (l = 0, cI = -1, childCellShapeOff = 0; l < numClosure; l++) { 3097 PetscInt point = closure[2 * l]; 3098 PetscInt pointDepth; 3099 3100 childO = closure[2 * l + 1]; 3101 if (point == child) { 3102 cI = l; 3103 break; 3104 } 3105 PetscCall(DMLabelGetValue(depth, point, &pointDepth)); 3106 childCellShapeOff += depthNumDof[pointDepth]; 3107 } 3108 if (l == numClosure) { 3109 pointMatOff += childDof; 3110 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */ 3111 } 3112 cperms = perms ? perms[cI] ? perms[cI][childO] : NULL : NULL; 3113 for (l = 0; l < childDof; l++) { 3114 PetscInt lCell = cperms ? cperms[l] : l; 3115 PetscInt childCellDof = childCellShapeOff + lCell; 3116 PetscReal *childValAtPoint; 3117 PetscReal val = 0.; 3118 3119 childValAtPoint = &Tchild->T[0][childCellDof * Nc]; 3120 for (m = 0; m < Nc; m++) val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m]; 3121 3122 pointMat[i * numChildDof + pointMatOff + l] += val; 3123 } 3124 pointMatOff += childDof; 3125 } 3126 PetscCall(DMPlexRestoreTransitiveClosure(refTree, childCell, PETSC_TRUE, &numClosure, &closure)); 3127 PetscCall(PetscTabulationDestroy(&Tchild)); 3128 } 3129 PetscCall(PetscTabulationDestroy(&Tparent)); 3130 } 3131 } else { /* just the volume-weighted averages of the children */ 3132 PetscReal parentVol; 3133 PetscInt childCell; 3134 3135 PetscCall(DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL)); 3136 for (i = 0, childCell = 0; i < numChildren; i++) { 3137 PetscInt child = children[i], j; 3138 PetscReal childVol; 3139 3140 if (child < cStart || child >= cEnd) continue; 3141 PetscCall(DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL)); 3142 for (j = 0; j < Nc; j++) pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol; 3143 childCell++; 3144 } 3145 } 3146 /* Insert pointMat into mat */ 3147 PetscCall(MatSetValues(mat, numSelfDof, matRows, numChildDof, matCols, pointMat, INSERT_VALUES)); 3148 PetscCall(DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT, &matRows)); 3149 PetscCall(DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR, &pointMat)); 3150 } 3151 } 3152 PetscCall(PetscFree6(v0, v0parent, vtmp, J, Jparent, invJ)); 3153 PetscCall(PetscFree2(pointScalar, pointRef)); 3154 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 3155 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 3156 *inj = mat; 3157 PetscFunctionReturn(PETSC_SUCCESS); 3158 } 3159 3160 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3161 { 3162 PetscDS ds; 3163 PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof; 3164 PetscScalar ***refPointFieldMats; 3165 PetscSection refConSec, refSection; 3166 3167 PetscFunctionBegin; 3168 PetscCall(DMGetDS(refTree, &ds)); 3169 PetscCall(PetscDSGetNumFields(ds, &numFields)); 3170 PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL)); 3171 PetscCall(DMGetLocalSection(refTree, &refSection)); 3172 PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd)); 3173 PetscCall(PetscMalloc1(pRefEnd - pRefStart, &refPointFieldMats)); 3174 PetscCall(PetscSectionGetMaxDof(refConSec, &maxDof)); 3175 PetscCall(PetscMalloc1(maxDof, &rows)); 3176 PetscCall(PetscMalloc1(maxDof * maxDof, &cols)); 3177 for (p = pRefStart; p < pRefEnd; p++) { 3178 PetscInt parent, pDof, parentDof; 3179 3180 PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL)); 3181 PetscCall(PetscSectionGetDof(refConSec, p, &pDof)); 3182 PetscCall(PetscSectionGetDof(refSection, parent, &parentDof)); 3183 if (!pDof || !parentDof || parent == p) continue; 3184 3185 PetscCall(PetscMalloc1(numFields, &refPointFieldMats[p - pRefStart])); 3186 for (f = 0; f < numFields; f++) { 3187 PetscInt cDof, cOff, numCols, r; 3188 3189 if (numFields > 1) { 3190 PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof)); 3191 PetscCall(PetscSectionGetFieldOffset(refConSec, p, f, &cOff)); 3192 } else { 3193 PetscCall(PetscSectionGetDof(refConSec, p, &cDof)); 3194 PetscCall(PetscSectionGetOffset(refConSec, p, &cOff)); 3195 } 3196 3197 for (r = 0; r < cDof; r++) rows[r] = cOff + r; 3198 numCols = 0; 3199 { 3200 PetscInt aDof, aOff, j; 3201 3202 if (numFields > 1) { 3203 PetscCall(PetscSectionGetFieldDof(refSection, parent, f, &aDof)); 3204 PetscCall(PetscSectionGetFieldOffset(refSection, parent, f, &aOff)); 3205 } else { 3206 PetscCall(PetscSectionGetDof(refSection, parent, &aDof)); 3207 PetscCall(PetscSectionGetOffset(refSection, parent, &aOff)); 3208 } 3209 3210 for (j = 0; j < aDof; j++) cols[numCols++] = aOff + j; 3211 } 3212 PetscCall(PetscMalloc1(cDof * numCols, &refPointFieldMats[p - pRefStart][f])); 3213 /* transpose of constraint matrix */ 3214 PetscCall(MatGetValues(inj, numCols, cols, cDof, rows, refPointFieldMats[p - pRefStart][f])); 3215 } 3216 } 3217 *childrenMats = refPointFieldMats; 3218 PetscCall(PetscFree(rows)); 3219 PetscCall(PetscFree(cols)); 3220 PetscFunctionReturn(PETSC_SUCCESS); 3221 } 3222 3223 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3224 { 3225 PetscDS ds; 3226 PetscScalar ***refPointFieldMats; 3227 PetscInt numFields, pRefStart, pRefEnd, p, f; 3228 PetscSection refConSec, refSection; 3229 3230 PetscFunctionBegin; 3231 refPointFieldMats = *childrenMats; 3232 *childrenMats = NULL; 3233 PetscCall(DMGetDS(refTree, &ds)); 3234 PetscCall(DMGetLocalSection(refTree, &refSection)); 3235 PetscCall(PetscDSGetNumFields(ds, &numFields)); 3236 PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL)); 3237 PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd)); 3238 for (p = pRefStart; p < pRefEnd; p++) { 3239 PetscInt parent, pDof, parentDof; 3240 3241 PetscCall(DMPlexGetTreeParent(refTree, p, &parent, NULL)); 3242 PetscCall(PetscSectionGetDof(refConSec, p, &pDof)); 3243 PetscCall(PetscSectionGetDof(refSection, parent, &parentDof)); 3244 if (!pDof || !parentDof || parent == p) continue; 3245 3246 for (f = 0; f < numFields; f++) { 3247 PetscInt cDof; 3248 3249 if (numFields > 1) { 3250 PetscCall(PetscSectionGetFieldDof(refConSec, p, f, &cDof)); 3251 } else { 3252 PetscCall(PetscSectionGetDof(refConSec, p, &cDof)); 3253 } 3254 3255 PetscCall(PetscFree(refPointFieldMats[p - pRefStart][f])); 3256 } 3257 PetscCall(PetscFree(refPointFieldMats[p - pRefStart])); 3258 } 3259 PetscCall(PetscFree(refPointFieldMats)); 3260 PetscFunctionReturn(PETSC_SUCCESS); 3261 } 3262 3263 static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree, Mat *injRef) 3264 { 3265 Mat cMatRef; 3266 PetscObject injRefObj; 3267 3268 PetscFunctionBegin; 3269 PetscCall(DMGetDefaultConstraints(refTree, NULL, &cMatRef, NULL)); 3270 PetscCall(PetscObjectQuery((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", &injRefObj)); 3271 *injRef = (Mat)injRefObj; 3272 if (!*injRef) { 3273 PetscCall(DMPlexComputeInjectorReferenceTree(refTree, injRef)); 3274 PetscCall(PetscObjectCompose((PetscObject)cMatRef, "DMPlexComputeInjectorTree_refTree", (PetscObject)*injRef)); 3275 /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */ 3276 PetscCall(PetscObjectDereference((PetscObject)*injRef)); 3277 } 3278 PetscFunctionReturn(PETSC_SUCCESS); 3279 } 3280 3281 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) 3282 { 3283 PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti; 3284 PetscSection globalCoarse, globalFine; 3285 PetscSection localCoarse, localFine, leafIndicesSec; 3286 PetscSection multiRootSec, rootIndicesSec; 3287 PetscInt *leafInds, *rootInds = NULL; 3288 const PetscInt *rootDegrees; 3289 PetscScalar *leafVals = NULL, *rootVals = NULL; 3290 PetscSF coarseToFineEmbedded; 3291 3292 PetscFunctionBegin; 3293 PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC)); 3294 PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF)); 3295 PetscCall(DMGetLocalSection(fine, &localFine)); 3296 PetscCall(DMGetGlobalSection(fine, &globalFine)); 3297 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafIndicesSec)); 3298 PetscCall(PetscSectionSetChart(leafIndicesSec, pStartF, pEndF)); 3299 PetscCall(PetscSectionGetMaxDof(localFine, &maxDof)); 3300 { /* winnow fine points that don't have global dofs out of the sf */ 3301 PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices; 3302 const PetscInt *leaves; 3303 3304 PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL)); 3305 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3306 p = leaves ? leaves[l] : l; 3307 PetscCall(PetscSectionGetDof(globalFine, p, &dof)); 3308 PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof)); 3309 if ((dof - cdof) > 0) { 3310 numPointsWithDofs++; 3311 3312 PetscCall(PetscSectionGetDof(localFine, p, &dof)); 3313 PetscCall(PetscSectionSetDof(leafIndicesSec, p, dof + 1)); 3314 } 3315 } 3316 PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs)); 3317 PetscCall(PetscSectionSetUp(leafIndicesSec)); 3318 PetscCall(PetscSectionGetStorageSize(leafIndicesSec, &numIndices)); 3319 PetscCall(PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1), &leafInds)); 3320 if (gatheredValues) PetscCall(PetscMalloc1(numIndices, &leafVals)); 3321 for (l = 0, offset = 0; l < nleaves; l++) { 3322 p = leaves ? leaves[l] : l; 3323 PetscCall(PetscSectionGetDof(globalFine, p, &dof)); 3324 PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof)); 3325 if ((dof - cdof) > 0) { 3326 PetscInt off, gOff; 3327 PetscInt *pInd; 3328 PetscScalar *pVal = NULL; 3329 3330 pointsWithDofs[offset++] = l; 3331 3332 PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &off)); 3333 3334 pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds; 3335 if (gatheredValues) { 3336 PetscInt i; 3337 3338 pVal = &leafVals[off + 1]; 3339 for (i = 0; i < dof; i++) pVal[i] = 0.; 3340 } 3341 PetscCall(PetscSectionGetOffset(globalFine, p, &gOff)); 3342 3343 offsets[0] = 0; 3344 if (numFields) { 3345 PetscInt f; 3346 3347 for (f = 0; f < numFields; f++) { 3348 PetscInt fDof; 3349 PetscCall(PetscSectionGetFieldDof(localFine, p, f, &fDof)); 3350 offsets[f + 1] = fDof + offsets[f]; 3351 } 3352 PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, -1, NULL, pInd)); 3353 } else { 3354 PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsets, PETSC_FALSE, NULL, NULL, pInd)); 3355 } 3356 if (gatheredValues) PetscCall(VecGetValues(fineVec, dof, pInd, pVal)); 3357 } 3358 } 3359 PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded)); 3360 PetscCall(PetscFree(pointsWithDofs)); 3361 } 3362 3363 PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC)); 3364 PetscCall(DMGetLocalSection(coarse, &localCoarse)); 3365 PetscCall(DMGetGlobalSection(coarse, &globalCoarse)); 3366 3367 { /* there may be the case where an sf root has a parent: broadcast parents back to children */ 3368 MPI_Datatype threeInt; 3369 PetscMPIInt rank; 3370 PetscInt(*parentNodeAndIdCoarse)[3]; 3371 PetscInt(*parentNodeAndIdFine)[3]; 3372 PetscInt p, nleaves, nleavesToParents; 3373 PetscSF pointSF, sfToParents; 3374 const PetscInt *ilocal; 3375 const PetscSFNode *iremote; 3376 PetscSFNode *iremoteToParents; 3377 PetscInt *ilocalToParents; 3378 3379 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank)); 3380 PetscCallMPI(MPI_Type_contiguous(3, MPIU_INT, &threeInt)); 3381 PetscCallMPI(MPI_Type_commit(&threeInt)); 3382 PetscCall(PetscMalloc2(pEndC - pStartC, &parentNodeAndIdCoarse, pEndF - pStartF, &parentNodeAndIdFine)); 3383 PetscCall(DMGetPointSF(coarse, &pointSF)); 3384 PetscCall(PetscSFGetGraph(pointSF, NULL, &nleaves, &ilocal, &iremote)); 3385 for (p = pStartC; p < pEndC; p++) { 3386 PetscInt parent, childId; 3387 PetscCall(DMPlexGetTreeParent(coarse, p, &parent, &childId)); 3388 parentNodeAndIdCoarse[p - pStartC][0] = rank; 3389 parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC; 3390 parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId; 3391 if (nleaves > 0) { 3392 PetscInt leaf = -1; 3393 3394 if (ilocal) { 3395 PetscCall(PetscFindInt(parent, nleaves, ilocal, &leaf)); 3396 } else { 3397 leaf = p - pStartC; 3398 } 3399 if (leaf >= 0) { 3400 parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank; 3401 parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index; 3402 } 3403 } 3404 } 3405 for (p = pStartF; p < pEndF; p++) { 3406 parentNodeAndIdFine[p - pStartF][0] = -1; 3407 parentNodeAndIdFine[p - pStartF][1] = -1; 3408 parentNodeAndIdFine[p - pStartF][2] = -1; 3409 } 3410 PetscCall(PetscSFBcastBegin(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE)); 3411 PetscCall(PetscSFBcastEnd(coarseToFineEmbedded, threeInt, parentNodeAndIdCoarse, parentNodeAndIdFine, MPI_REPLACE)); 3412 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3413 PetscInt dof; 3414 3415 PetscCall(PetscSectionGetDof(leafIndicesSec, p, &dof)); 3416 if (dof) { 3417 PetscInt off; 3418 3419 PetscCall(PetscSectionGetOffset(leafIndicesSec, p, &off)); 3420 if (gatheredIndices) { 3421 leafInds[off] = PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]); 3422 } else if (gatheredValues) { 3423 leafVals[off] = (PetscScalar)PetscMax(childIds[p - pStartF], parentNodeAndIdFine[p - pStartF][2]); 3424 } 3425 } 3426 if (parentNodeAndIdFine[p - pStartF][0] >= 0) nleavesToParents++; 3427 } 3428 PetscCall(PetscMalloc1(nleavesToParents, &ilocalToParents)); 3429 PetscCall(PetscMalloc1(nleavesToParents, &iremoteToParents)); 3430 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3431 if (parentNodeAndIdFine[p - pStartF][0] >= 0) { 3432 ilocalToParents[nleavesToParents] = p - pStartF; 3433 iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p - pStartF][0]; 3434 iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p - pStartF][1]; 3435 nleavesToParents++; 3436 } 3437 } 3438 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)coarse), &sfToParents)); 3439 PetscCall(PetscSFSetGraph(sfToParents, pEndC - pStartC, nleavesToParents, ilocalToParents, PETSC_OWN_POINTER, iremoteToParents, PETSC_OWN_POINTER)); 3440 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 3441 3442 coarseToFineEmbedded = sfToParents; 3443 3444 PetscCall(PetscFree2(parentNodeAndIdCoarse, parentNodeAndIdFine)); 3445 PetscCallMPI(MPI_Type_free(&threeInt)); 3446 } 3447 3448 { /* winnow out coarse points that don't have dofs */ 3449 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3450 PetscSF sfDofsOnly; 3451 3452 for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) { 3453 PetscCall(PetscSectionGetDof(globalCoarse, p, &dof)); 3454 PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof)); 3455 if ((dof - cdof) > 0) numPointsWithDofs++; 3456 } 3457 PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs)); 3458 for (p = pStartC, offset = 0; p < pEndC; p++) { 3459 PetscCall(PetscSectionGetDof(globalCoarse, p, &dof)); 3460 PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof)); 3461 if ((dof - cdof) > 0) pointsWithDofs[offset++] = p - pStartC; 3462 } 3463 PetscCall(PetscSFCreateEmbeddedRootSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly)); 3464 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 3465 PetscCall(PetscFree(pointsWithDofs)); 3466 coarseToFineEmbedded = sfDofsOnly; 3467 } 3468 3469 /* communicate back to the coarse mesh which coarse points have children (that may require injection) */ 3470 PetscCall(PetscSFComputeDegreeBegin(coarseToFineEmbedded, &rootDegrees)); 3471 PetscCall(PetscSFComputeDegreeEnd(coarseToFineEmbedded, &rootDegrees)); 3472 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &multiRootSec)); 3473 PetscCall(PetscSectionSetChart(multiRootSec, pStartC, pEndC)); 3474 for (p = pStartC; p < pEndC; p++) PetscCall(PetscSectionSetDof(multiRootSec, p, rootDegrees[p - pStartC])); 3475 PetscCall(PetscSectionSetUp(multiRootSec)); 3476 PetscCall(PetscSectionGetStorageSize(multiRootSec, &numMulti)); 3477 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootIndicesSec)); 3478 { /* distribute the leaf section */ 3479 PetscSF multi, multiInv, indicesSF; 3480 PetscInt *remoteOffsets, numRootIndices; 3481 3482 PetscCall(PetscSFGetMultiSF(coarseToFineEmbedded, &multi)); 3483 PetscCall(PetscSFCreateInverseSF(multi, &multiInv)); 3484 PetscCall(PetscSFDistributeSection(multiInv, leafIndicesSec, &remoteOffsets, rootIndicesSec)); 3485 PetscCall(PetscSFCreateSectionSF(multiInv, leafIndicesSec, remoteOffsets, rootIndicesSec, &indicesSF)); 3486 PetscCall(PetscFree(remoteOffsets)); 3487 PetscCall(PetscSFDestroy(&multiInv)); 3488 PetscCall(PetscSectionGetStorageSize(rootIndicesSec, &numRootIndices)); 3489 if (gatheredIndices) { 3490 PetscCall(PetscMalloc1(numRootIndices, &rootInds)); 3491 PetscCall(PetscSFBcastBegin(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE)); 3492 PetscCall(PetscSFBcastEnd(indicesSF, MPIU_INT, leafInds, rootInds, MPI_REPLACE)); 3493 } 3494 if (gatheredValues) { 3495 PetscCall(PetscMalloc1(numRootIndices, &rootVals)); 3496 PetscCall(PetscSFBcastBegin(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE)); 3497 PetscCall(PetscSFBcastEnd(indicesSF, MPIU_SCALAR, leafVals, rootVals, MPI_REPLACE)); 3498 } 3499 PetscCall(PetscSFDestroy(&indicesSF)); 3500 } 3501 PetscCall(PetscSectionDestroy(&leafIndicesSec)); 3502 PetscCall(PetscFree(leafInds)); 3503 PetscCall(PetscFree(leafVals)); 3504 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 3505 *rootMultiSec = multiRootSec; 3506 *multiLeafSec = rootIndicesSec; 3507 if (gatheredIndices) *gatheredIndices = rootInds; 3508 if (gatheredValues) *gatheredValues = rootVals; 3509 PetscFunctionReturn(PETSC_SUCCESS); 3510 } 3511 3512 PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 3513 { 3514 DM refTree; 3515 PetscSection multiRootSec, rootIndicesSec; 3516 PetscSection globalCoarse, globalFine; 3517 PetscSection localCoarse, localFine; 3518 PetscSection cSecRef; 3519 PetscInt *rootIndices = NULL, *parentIndices, pRefStart, pRefEnd; 3520 Mat injRef; 3521 PetscInt numFields, maxDof; 3522 PetscInt pStartC, pEndC, pStartF, pEndF, p; 3523 PetscInt *offsets, *offsetsCopy, *rowOffsets; 3524 PetscLayout rowMap, colMap; 3525 PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO; 3526 PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 3527 3528 PetscFunctionBegin; 3529 /* get the templates for the fine-to-coarse injection from the reference tree */ 3530 PetscCall(DMPlexGetReferenceTree(coarse, &refTree)); 3531 PetscCall(DMCopyDisc(coarse, refTree)); 3532 PetscCall(DMSetLocalSection(refTree, NULL)); 3533 PetscCall(DMSetDefaultConstraints(refTree, NULL, NULL, NULL)); 3534 PetscCall(DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL)); 3535 PetscCall(PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd)); 3536 PetscCall(DMPlexReferenceTreeGetInjector(refTree, &injRef)); 3537 3538 PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF)); 3539 PetscCall(DMGetLocalSection(fine, &localFine)); 3540 PetscCall(DMGetGlobalSection(fine, &globalFine)); 3541 PetscCall(PetscSectionGetNumFields(localFine, &numFields)); 3542 PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC)); 3543 PetscCall(DMGetLocalSection(coarse, &localCoarse)); 3544 PetscCall(DMGetGlobalSection(coarse, &globalCoarse)); 3545 PetscCall(PetscSectionGetMaxDof(localCoarse, &maxDof)); 3546 { 3547 PetscInt maxFields = PetscMax(1, numFields) + 1; 3548 PetscCall(PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets)); 3549 } 3550 3551 PetscCall(DMPlexTransferInjectorTree(coarse, fine, coarseToFine, childIds, NULL, numFields, offsets, &multiRootSec, &rootIndicesSec, &rootIndices, NULL)); 3552 3553 PetscCall(PetscMalloc1(maxDof, &parentIndices)); 3554 3555 /* count indices */ 3556 PetscCall(MatGetLayouts(mat, &rowMap, &colMap)); 3557 PetscCall(PetscLayoutSetUp(rowMap)); 3558 PetscCall(PetscLayoutSetUp(colMap)); 3559 PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd)); 3560 PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd)); 3561 PetscCall(PetscCalloc2(rowEnd - rowStart, &nnzD, rowEnd - rowStart, &nnzO)); 3562 for (p = pStartC; p < pEndC; p++) { 3563 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3564 3565 PetscCall(PetscSectionGetDof(globalCoarse, p, &dof)); 3566 PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof)); 3567 if ((dof - cdof) <= 0) continue; 3568 PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff)); 3569 3570 rowOffsets[0] = 0; 3571 offsetsCopy[0] = 0; 3572 if (numFields) { 3573 PetscInt f; 3574 3575 for (f = 0; f < numFields; f++) { 3576 PetscInt fDof; 3577 PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof)); 3578 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3579 } 3580 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices)); 3581 } else { 3582 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices)); 3583 rowOffsets[1] = offsetsCopy[0]; 3584 } 3585 3586 PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves)); 3587 PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart)); 3588 leafEnd = leafStart + numLeaves; 3589 for (l = leafStart; l < leafEnd; l++) { 3590 PetscInt numIndices, childId, offset; 3591 const PetscInt *childIndices; 3592 3593 PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices)); 3594 PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset)); 3595 childId = rootIndices[offset++]; 3596 childIndices = &rootIndices[offset]; 3597 numIndices--; 3598 3599 if (childId == -1) { /* equivalent points: scatter */ 3600 PetscInt i; 3601 3602 for (i = 0; i < numIndices; i++) { 3603 PetscInt colIndex = childIndices[i]; 3604 PetscInt rowIndex = parentIndices[i]; 3605 if (rowIndex < 0) continue; 3606 PetscCheck(colIndex >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unconstrained fine and constrained coarse"); 3607 if (colIndex >= colStart && colIndex < colEnd) { 3608 nnzD[rowIndex - rowStart] = 1; 3609 } else { 3610 nnzO[rowIndex - rowStart] = 1; 3611 } 3612 } 3613 } else { 3614 PetscInt parentId, f, lim; 3615 3616 PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL)); 3617 3618 lim = PetscMax(1, numFields); 3619 offsets[0] = 0; 3620 if (numFields) { 3621 PetscInt f; 3622 3623 for (f = 0; f < numFields; f++) { 3624 PetscInt fDof; 3625 PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof)); 3626 3627 offsets[f + 1] = fDof + offsets[f]; 3628 } 3629 } else { 3630 PetscInt cDof; 3631 3632 PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof)); 3633 offsets[1] = cDof; 3634 } 3635 for (f = 0; f < lim; f++) { 3636 PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1]; 3637 PetscInt childStart = offsets[f], childEnd = offsets[f + 1]; 3638 PetscInt i, numD = 0, numO = 0; 3639 3640 for (i = childStart; i < childEnd; i++) { 3641 PetscInt colIndex = childIndices[i]; 3642 3643 if (colIndex < 0) continue; 3644 if (colIndex >= colStart && colIndex < colEnd) { 3645 numD++; 3646 } else { 3647 numO++; 3648 } 3649 } 3650 for (i = parentStart; i < parentEnd; i++) { 3651 PetscInt rowIndex = parentIndices[i]; 3652 3653 if (rowIndex < 0) continue; 3654 nnzD[rowIndex - rowStart] += numD; 3655 nnzO[rowIndex - rowStart] += numO; 3656 } 3657 } 3658 } 3659 } 3660 } 3661 /* preallocate */ 3662 PetscCall(MatXAIJSetPreallocation(mat, 1, nnzD, nnzO, NULL, NULL)); 3663 PetscCall(PetscFree2(nnzD, nnzO)); 3664 /* insert values */ 3665 PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats)); 3666 for (p = pStartC; p < pEndC; p++) { 3667 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3668 3669 PetscCall(PetscSectionGetDof(globalCoarse, p, &dof)); 3670 PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof)); 3671 if ((dof - cdof) <= 0) continue; 3672 PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff)); 3673 3674 rowOffsets[0] = 0; 3675 offsetsCopy[0] = 0; 3676 if (numFields) { 3677 PetscInt f; 3678 3679 for (f = 0; f < numFields; f++) { 3680 PetscInt fDof; 3681 PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof)); 3682 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3683 } 3684 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices)); 3685 } else { 3686 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices)); 3687 rowOffsets[1] = offsetsCopy[0]; 3688 } 3689 3690 PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves)); 3691 PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart)); 3692 leafEnd = leafStart + numLeaves; 3693 for (l = leafStart; l < leafEnd; l++) { 3694 PetscInt numIndices, childId, offset; 3695 const PetscInt *childIndices; 3696 3697 PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices)); 3698 PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset)); 3699 childId = rootIndices[offset++]; 3700 childIndices = &rootIndices[offset]; 3701 numIndices--; 3702 3703 if (childId == -1) { /* equivalent points: scatter */ 3704 PetscInt i; 3705 3706 for (i = 0; i < numIndices; i++) PetscCall(MatSetValue(mat, parentIndices[i], childIndices[i], 1., INSERT_VALUES)); 3707 } else { 3708 PetscInt parentId, f, lim; 3709 3710 PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL)); 3711 3712 lim = PetscMax(1, numFields); 3713 offsets[0] = 0; 3714 if (numFields) { 3715 PetscInt f; 3716 3717 for (f = 0; f < numFields; f++) { 3718 PetscInt fDof; 3719 PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof)); 3720 3721 offsets[f + 1] = fDof + offsets[f]; 3722 } 3723 } else { 3724 PetscInt cDof; 3725 3726 PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof)); 3727 offsets[1] = cDof; 3728 } 3729 for (f = 0; f < lim; f++) { 3730 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 3731 PetscInt *rowIndices = &parentIndices[rowOffsets[f]]; 3732 const PetscInt *colIndices = &childIndices[offsets[f]]; 3733 3734 PetscCall(MatSetValues(mat, rowOffsets[f + 1] - rowOffsets[f], rowIndices, offsets[f + 1] - offsets[f], colIndices, childMat, INSERT_VALUES)); 3735 } 3736 } 3737 } 3738 } 3739 PetscCall(PetscSectionDestroy(&multiRootSec)); 3740 PetscCall(PetscSectionDestroy(&rootIndicesSec)); 3741 PetscCall(PetscFree(parentIndices)); 3742 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats)); 3743 PetscCall(PetscFree(rootIndices)); 3744 PetscCall(PetscFree3(offsets, offsetsCopy, rowOffsets)); 3745 3746 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 3747 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 3748 PetscFunctionReturn(PETSC_SUCCESS); 3749 } 3750 3751 static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom) 3752 { 3753 PetscSF coarseToFineEmbedded; 3754 PetscSection globalCoarse, globalFine; 3755 PetscSection localCoarse, localFine; 3756 PetscSection aSec, cSec; 3757 PetscSection rootValuesSec; 3758 PetscSection leafValuesSec; 3759 PetscScalar *rootValues, *leafValues; 3760 IS aIS; 3761 const PetscInt *anchors; 3762 Mat cMat; 3763 PetscInt numFields; 3764 PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd; 3765 PetscInt aStart, aEnd, cStart, cEnd; 3766 PetscInt *maxChildIds; 3767 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 3768 PetscFV fv = NULL; 3769 PetscInt dim, numFVcomps = -1, fvField = -1; 3770 DM cellDM = NULL, gradDM = NULL; 3771 const PetscScalar *cellGeomArray = NULL; 3772 const PetscScalar *gradArray = NULL; 3773 3774 PetscFunctionBegin; 3775 PetscCall(VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE)); 3776 PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC)); 3777 PetscCall(DMPlexGetSimplexOrBoxCells(coarse, 0, &cellStart, &cellEnd)); 3778 PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF)); 3779 PetscCall(DMGetGlobalSection(fine, &globalFine)); 3780 PetscCall(DMGetCoordinateDim(coarse, &dim)); 3781 { /* winnow fine points that don't have global dofs out of the sf */ 3782 PetscInt nleaves, l; 3783 const PetscInt *leaves; 3784 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3785 3786 PetscCall(PetscSFGetGraph(coarseToFine, NULL, &nleaves, &leaves, NULL)); 3787 3788 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3789 PetscInt p = leaves ? leaves[l] : l; 3790 3791 PetscCall(PetscSectionGetDof(globalFine, p, &dof)); 3792 PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof)); 3793 if ((dof - cdof) > 0) numPointsWithDofs++; 3794 } 3795 PetscCall(PetscMalloc1(numPointsWithDofs, &pointsWithDofs)); 3796 for (l = 0, offset = 0; l < nleaves; l++) { 3797 PetscInt p = leaves ? leaves[l] : l; 3798 3799 PetscCall(PetscSectionGetDof(globalFine, p, &dof)); 3800 PetscCall(PetscSectionGetConstraintDof(globalFine, p, &cdof)); 3801 if ((dof - cdof) > 0) pointsWithDofs[offset++] = l; 3802 } 3803 PetscCall(PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded)); 3804 PetscCall(PetscFree(pointsWithDofs)); 3805 } 3806 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 3807 PetscCall(PetscMalloc1(pEndC - pStartC, &maxChildIds)); 3808 for (p = pStartC; p < pEndC; p++) maxChildIds[p - pStartC] = -2; 3809 PetscCall(PetscSFReduceBegin(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX)); 3810 PetscCall(PetscSFReduceEnd(coarseToFineEmbedded, MPIU_INT, cids, maxChildIds, MPIU_MAX)); 3811 3812 PetscCall(DMGetLocalSection(coarse, &localCoarse)); 3813 PetscCall(DMGetGlobalSection(coarse, &globalCoarse)); 3814 3815 PetscCall(DMPlexGetAnchors(coarse, &aSec, &aIS)); 3816 PetscCall(ISGetIndices(aIS, &anchors)); 3817 PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd)); 3818 3819 PetscCall(DMGetDefaultConstraints(coarse, &cSec, &cMat, NULL)); 3820 PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd)); 3821 3822 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 3823 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)coarse), &rootValuesSec)); 3824 PetscCall(PetscSectionSetChart(rootValuesSec, pStartC, pEndC)); 3825 PetscCall(PetscSectionGetNumFields(localCoarse, &numFields)); 3826 { 3827 PetscInt maxFields = PetscMax(1, numFields) + 1; 3828 PetscCall(PetscMalloc7(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &newOffsets, maxFields, &newOffsetsCopy, maxFields, &rowOffsets, maxFields, &numD, maxFields, &numO)); 3829 } 3830 if (grad) { 3831 PetscInt i; 3832 3833 PetscCall(VecGetDM(cellGeom, &cellDM)); 3834 PetscCall(VecGetArrayRead(cellGeom, &cellGeomArray)); 3835 PetscCall(VecGetDM(grad, &gradDM)); 3836 PetscCall(VecGetArrayRead(grad, &gradArray)); 3837 for (i = 0; i < PetscMax(1, numFields); i++) { 3838 PetscObject obj; 3839 PetscClassId id; 3840 3841 PetscCall(DMGetField(coarse, i, NULL, &obj)); 3842 PetscCall(PetscObjectGetClassId(obj, &id)); 3843 if (id == PETSCFV_CLASSID) { 3844 fv = (PetscFV)obj; 3845 PetscCall(PetscFVGetNumComponents(fv, &numFVcomps)); 3846 fvField = i; 3847 break; 3848 } 3849 } 3850 } 3851 3852 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 3853 PetscInt dof; 3854 PetscInt maxChildId = maxChildIds[p - pStartC]; 3855 PetscInt numValues = 0; 3856 3857 PetscCall(PetscSectionGetDof(globalCoarse, p, &dof)); 3858 if (dof < 0) dof = -(dof + 1); 3859 offsets[0] = 0; 3860 newOffsets[0] = 0; 3861 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 3862 PetscInt *closure = NULL, closureSize, cl; 3863 3864 PetscCall(DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure)); 3865 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 3866 PetscInt c = closure[2 * cl], clDof; 3867 3868 PetscCall(PetscSectionGetDof(localCoarse, c, &clDof)); 3869 numValues += clDof; 3870 } 3871 PetscCall(DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure)); 3872 } else if (maxChildId == -1) { 3873 PetscCall(PetscSectionGetDof(localCoarse, p, &numValues)); 3874 } 3875 /* we will pack the column indices with the field offsets */ 3876 if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) { 3877 /* also send the centroid, and the gradient */ 3878 numValues += dim * (1 + numFVcomps); 3879 } 3880 PetscCall(PetscSectionSetDof(rootValuesSec, p, numValues)); 3881 } 3882 PetscCall(PetscSectionSetUp(rootValuesSec)); 3883 { 3884 PetscInt numRootValues; 3885 const PetscScalar *coarseArray; 3886 3887 PetscCall(PetscSectionGetStorageSize(rootValuesSec, &numRootValues)); 3888 PetscCall(PetscMalloc1(numRootValues, &rootValues)); 3889 PetscCall(VecGetArrayRead(vecCoarseLocal, &coarseArray)); 3890 for (p = pStartC; p < pEndC; p++) { 3891 PetscInt numValues; 3892 PetscInt pValOff; 3893 PetscScalar *pVal; 3894 PetscInt maxChildId = maxChildIds[p - pStartC]; 3895 3896 PetscCall(PetscSectionGetDof(rootValuesSec, p, &numValues)); 3897 if (!numValues) continue; 3898 PetscCall(PetscSectionGetOffset(rootValuesSec, p, &pValOff)); 3899 pVal = &rootValues[pValOff]; 3900 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 3901 PetscInt closureSize = numValues; 3902 PetscCall(DMPlexVecGetClosure(coarse, NULL, vecCoarseLocal, p, &closureSize, &pVal)); 3903 if (grad && p >= cellStart && p < cellEnd) { 3904 PetscFVCellGeom *cg; 3905 PetscScalar *gradVals = NULL; 3906 PetscInt i; 3907 3908 pVal += (numValues - dim * (1 + numFVcomps)); 3909 3910 PetscCall(DMPlexPointLocalRead(cellDM, p, cellGeomArray, (void *)&cg)); 3911 for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i]; 3912 pVal += dim; 3913 PetscCall(DMPlexPointGlobalRead(gradDM, p, gradArray, (void *)&gradVals)); 3914 for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i]; 3915 } 3916 } else if (maxChildId == -1) { 3917 PetscInt lDof, lOff, i; 3918 3919 PetscCall(PetscSectionGetDof(localCoarse, p, &lDof)); 3920 PetscCall(PetscSectionGetOffset(localCoarse, p, &lOff)); 3921 for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i]; 3922 } 3923 } 3924 PetscCall(VecRestoreArrayRead(vecCoarseLocal, &coarseArray)); 3925 PetscCall(PetscFree(maxChildIds)); 3926 } 3927 { 3928 PetscSF valuesSF; 3929 PetscInt *remoteOffsetsValues, numLeafValues; 3930 3931 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)fine), &leafValuesSec)); 3932 PetscCall(PetscSFDistributeSection(coarseToFineEmbedded, rootValuesSec, &remoteOffsetsValues, leafValuesSec)); 3933 PetscCall(PetscSFCreateSectionSF(coarseToFineEmbedded, rootValuesSec, remoteOffsetsValues, leafValuesSec, &valuesSF)); 3934 PetscCall(PetscSFDestroy(&coarseToFineEmbedded)); 3935 PetscCall(PetscFree(remoteOffsetsValues)); 3936 PetscCall(PetscSectionGetStorageSize(leafValuesSec, &numLeafValues)); 3937 PetscCall(PetscMalloc1(numLeafValues, &leafValues)); 3938 PetscCall(PetscSFBcastBegin(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE)); 3939 PetscCall(PetscSFBcastEnd(valuesSF, MPIU_SCALAR, rootValues, leafValues, MPI_REPLACE)); 3940 PetscCall(PetscSFDestroy(&valuesSF)); 3941 PetscCall(PetscFree(rootValues)); 3942 PetscCall(PetscSectionDestroy(&rootValuesSec)); 3943 } 3944 PetscCall(DMGetLocalSection(fine, &localFine)); 3945 { 3946 PetscInt maxDof; 3947 PetscInt *rowIndices; 3948 DM refTree; 3949 PetscInt **refPointFieldN; 3950 PetscScalar ***refPointFieldMats; 3951 PetscSection refConSec, refAnSec; 3952 PetscInt pRefStart, pRefEnd, leafStart, leafEnd; 3953 PetscScalar *pointWork; 3954 3955 PetscCall(PetscSectionGetMaxDof(localFine, &maxDof)); 3956 PetscCall(DMGetWorkArray(fine, maxDof, MPIU_INT, &rowIndices)); 3957 PetscCall(DMGetWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork)); 3958 PetscCall(DMPlexGetReferenceTree(fine, &refTree)); 3959 PetscCall(DMCopyDisc(fine, refTree)); 3960 PetscCall(DMPlexReferenceTreeGetChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN)); 3961 PetscCall(DMGetDefaultConstraints(refTree, &refConSec, NULL, NULL)); 3962 PetscCall(DMPlexGetAnchors(refTree, &refAnSec, NULL)); 3963 PetscCall(PetscSectionGetChart(refConSec, &pRefStart, &pRefEnd)); 3964 PetscCall(PetscSectionGetChart(leafValuesSec, &leafStart, &leafEnd)); 3965 PetscCall(DMPlexGetSimplexOrBoxCells(fine, 0, &cellStart, &cellEnd)); 3966 for (p = leafStart; p < leafEnd; p++) { 3967 PetscInt gDof, gcDof, gOff, lDof; 3968 PetscInt numValues, pValOff; 3969 PetscInt childId; 3970 const PetscScalar *pVal; 3971 const PetscScalar *fvGradData = NULL; 3972 3973 PetscCall(PetscSectionGetDof(globalFine, p, &gDof)); 3974 PetscCall(PetscSectionGetDof(localFine, p, &lDof)); 3975 PetscCall(PetscSectionGetConstraintDof(globalFine, p, &gcDof)); 3976 if ((gDof - gcDof) <= 0) continue; 3977 PetscCall(PetscSectionGetOffset(globalFine, p, &gOff)); 3978 PetscCall(PetscSectionGetDof(leafValuesSec, p, &numValues)); 3979 if (!numValues) continue; 3980 PetscCall(PetscSectionGetOffset(leafValuesSec, p, &pValOff)); 3981 pVal = &leafValues[pValOff]; 3982 offsets[0] = 0; 3983 offsetsCopy[0] = 0; 3984 newOffsets[0] = 0; 3985 newOffsetsCopy[0] = 0; 3986 childId = cids[p - pStartF]; 3987 if (numFields) { 3988 PetscInt f; 3989 for (f = 0; f < numFields; f++) { 3990 PetscInt rowDof; 3991 3992 PetscCall(PetscSectionGetFieldDof(localFine, p, f, &rowDof)); 3993 offsets[f + 1] = offsets[f] + rowDof; 3994 offsetsCopy[f + 1] = offsets[f + 1]; 3995 /* TODO: closure indices */ 3996 newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]); 3997 } 3998 PetscCall(DMPlexGetIndicesPointFields_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, rowIndices)); 3999 } else { 4000 offsets[0] = 0; 4001 offsets[1] = lDof; 4002 newOffsets[0] = 0; 4003 newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0]; 4004 PetscCall(DMPlexGetIndicesPoint_Internal(localFine, PETSC_FALSE, p, gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, rowIndices)); 4005 } 4006 if (childId == -1) { /* no child interpolation: one nnz per */ 4007 PetscCall(VecSetValues(vecFine, numValues, rowIndices, pVal, INSERT_VALUES)); 4008 } else { 4009 PetscInt f; 4010 4011 if (grad && p >= cellStart && p < cellEnd) { 4012 numValues -= (dim * (1 + numFVcomps)); 4013 fvGradData = &pVal[numValues]; 4014 } 4015 for (f = 0; f < PetscMax(1, numFields); f++) { 4016 const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f]; 4017 PetscInt numRows = offsets[f + 1] - offsets[f]; 4018 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 4019 const PetscScalar *cVal = &pVal[newOffsets[f]]; 4020 PetscScalar *rVal = &pointWork[offsets[f]]; 4021 PetscInt i, j; 4022 4023 #if 0 4024 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)); 4025 #endif 4026 for (i = 0; i < numRows; i++) { 4027 PetscScalar val = 0.; 4028 for (j = 0; j < numCols; j++) val += childMat[i * numCols + j] * cVal[j]; 4029 rVal[i] = val; 4030 } 4031 if (f == fvField && p >= cellStart && p < cellEnd) { 4032 PetscReal centroid[3]; 4033 PetscScalar diff[3]; 4034 const PetscScalar *parentCentroid = &fvGradData[0]; 4035 const PetscScalar *gradient = &fvGradData[dim]; 4036 4037 PetscCall(DMPlexComputeCellGeometryFVM(fine, p, NULL, centroid, NULL)); 4038 for (i = 0; i < dim; i++) diff[i] = centroid[i] - parentCentroid[i]; 4039 for (i = 0; i < numFVcomps; i++) { 4040 PetscScalar val = 0.; 4041 4042 for (j = 0; j < dim; j++) val += gradient[dim * i + j] * diff[j]; 4043 rVal[i] += val; 4044 } 4045 } 4046 PetscCall(VecSetValues(vecFine, numRows, &rowIndices[offsets[f]], rVal, INSERT_VALUES)); 4047 } 4048 } 4049 } 4050 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices(refTree, &refPointFieldMats, &refPointFieldN)); 4051 PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_SCALAR, &pointWork)); 4052 PetscCall(DMRestoreWorkArray(fine, maxDof, MPIU_INT, &rowIndices)); 4053 } 4054 PetscCall(PetscFree(leafValues)); 4055 PetscCall(PetscSectionDestroy(&leafValuesSec)); 4056 PetscCall(PetscFree7(offsets, offsetsCopy, newOffsets, newOffsetsCopy, rowOffsets, numD, numO)); 4057 PetscCall(ISRestoreIndices(aIS, &anchors)); 4058 PetscFunctionReturn(PETSC_SUCCESS); 4059 } 4060 4061 static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids) 4062 { 4063 DM refTree; 4064 PetscSection multiRootSec, rootIndicesSec; 4065 PetscSection globalCoarse, globalFine; 4066 PetscSection localCoarse, localFine; 4067 PetscSection cSecRef; 4068 PetscInt *parentIndices, pRefStart, pRefEnd; 4069 PetscScalar *rootValues, *parentValues; 4070 Mat injRef; 4071 PetscInt numFields, maxDof; 4072 PetscInt pStartC, pEndC, pStartF, pEndF, p; 4073 PetscInt *offsets, *offsetsCopy, *rowOffsets; 4074 PetscLayout rowMap, colMap; 4075 PetscInt rowStart, rowEnd, colStart, colEnd; 4076 PetscScalar ***childrenMats = NULL; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 4077 4078 PetscFunctionBegin; 4079 /* get the templates for the fine-to-coarse injection from the reference tree */ 4080 PetscCall(VecSetOption(vecFine, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE)); 4081 PetscCall(VecSetOption(vecCoarse, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE)); 4082 PetscCall(DMPlexGetReferenceTree(coarse, &refTree)); 4083 PetscCall(DMCopyDisc(coarse, refTree)); 4084 PetscCall(DMGetDefaultConstraints(refTree, &cSecRef, NULL, NULL)); 4085 PetscCall(PetscSectionGetChart(cSecRef, &pRefStart, &pRefEnd)); 4086 PetscCall(DMPlexReferenceTreeGetInjector(refTree, &injRef)); 4087 4088 PetscCall(DMPlexGetChart(fine, &pStartF, &pEndF)); 4089 PetscCall(DMGetLocalSection(fine, &localFine)); 4090 PetscCall(DMGetGlobalSection(fine, &globalFine)); 4091 PetscCall(PetscSectionGetNumFields(localFine, &numFields)); 4092 PetscCall(DMPlexGetChart(coarse, &pStartC, &pEndC)); 4093 PetscCall(DMGetLocalSection(coarse, &localCoarse)); 4094 PetscCall(DMGetGlobalSection(coarse, &globalCoarse)); 4095 PetscCall(PetscSectionGetMaxDof(localCoarse, &maxDof)); 4096 { 4097 PetscInt maxFields = PetscMax(1, numFields) + 1; 4098 PetscCall(PetscMalloc3(maxFields, &offsets, maxFields, &offsetsCopy, maxFields, &rowOffsets)); 4099 } 4100 4101 PetscCall(DMPlexTransferInjectorTree(coarse, fine, coarseToFine, cids, vecFine, numFields, offsets, &multiRootSec, &rootIndicesSec, NULL, &rootValues)); 4102 4103 PetscCall(PetscMalloc2(maxDof, &parentIndices, maxDof, &parentValues)); 4104 4105 /* count indices */ 4106 PetscCall(VecGetLayout(vecFine, &colMap)); 4107 PetscCall(VecGetLayout(vecCoarse, &rowMap)); 4108 PetscCall(PetscLayoutSetUp(rowMap)); 4109 PetscCall(PetscLayoutSetUp(colMap)); 4110 PetscCall(PetscLayoutGetRange(rowMap, &rowStart, &rowEnd)); 4111 PetscCall(PetscLayoutGetRange(colMap, &colStart, &colEnd)); 4112 /* insert values */ 4113 PetscCall(DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree, injRef, &childrenMats)); 4114 for (p = pStartC; p < pEndC; p++) { 4115 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 4116 PetscBool contribute = PETSC_FALSE; 4117 4118 PetscCall(PetscSectionGetDof(globalCoarse, p, &dof)); 4119 PetscCall(PetscSectionGetConstraintDof(globalCoarse, p, &cdof)); 4120 if ((dof - cdof) <= 0) continue; 4121 PetscCall(PetscSectionGetDof(localCoarse, p, &dof)); 4122 PetscCall(PetscSectionGetOffset(globalCoarse, p, &gOff)); 4123 4124 rowOffsets[0] = 0; 4125 offsetsCopy[0] = 0; 4126 if (numFields) { 4127 PetscInt f; 4128 4129 for (f = 0; f < numFields; f++) { 4130 PetscInt fDof; 4131 PetscCall(PetscSectionGetFieldDof(localCoarse, p, f, &fDof)); 4132 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 4133 } 4134 PetscCall(DMPlexGetIndicesPointFields_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, -1, NULL, parentIndices)); 4135 } else { 4136 PetscCall(DMPlexGetIndicesPoint_Internal(localCoarse, PETSC_FALSE, p, gOff < 0 ? -(gOff + 1) : gOff, offsetsCopy, PETSC_FALSE, NULL, NULL, parentIndices)); 4137 rowOffsets[1] = offsetsCopy[0]; 4138 } 4139 4140 PetscCall(PetscSectionGetDof(multiRootSec, p, &numLeaves)); 4141 PetscCall(PetscSectionGetOffset(multiRootSec, p, &leafStart)); 4142 leafEnd = leafStart + numLeaves; 4143 for (l = 0; l < dof; l++) parentValues[l] = 0.; 4144 for (l = leafStart; l < leafEnd; l++) { 4145 PetscInt numIndices, childId, offset; 4146 const PetscScalar *childValues; 4147 4148 PetscCall(PetscSectionGetDof(rootIndicesSec, l, &numIndices)); 4149 PetscCall(PetscSectionGetOffset(rootIndicesSec, l, &offset)); 4150 childId = (PetscInt)PetscRealPart(rootValues[offset++]); 4151 childValues = &rootValues[offset]; 4152 numIndices--; 4153 4154 if (childId == -2) { /* skip */ 4155 continue; 4156 } else if (childId == -1) { /* equivalent points: scatter */ 4157 PetscInt m; 4158 4159 contribute = PETSC_TRUE; 4160 for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m]; 4161 } else { /* contributions from children: sum with injectors from reference tree */ 4162 PetscInt parentId, f, lim; 4163 4164 contribute = PETSC_TRUE; 4165 PetscCall(DMPlexGetTreeParent(refTree, childId, &parentId, NULL)); 4166 4167 lim = PetscMax(1, numFields); 4168 offsets[0] = 0; 4169 if (numFields) { 4170 PetscInt f; 4171 4172 for (f = 0; f < numFields; f++) { 4173 PetscInt fDof; 4174 PetscCall(PetscSectionGetFieldDof(cSecRef, childId, f, &fDof)); 4175 4176 offsets[f + 1] = fDof + offsets[f]; 4177 } 4178 } else { 4179 PetscInt cDof; 4180 4181 PetscCall(PetscSectionGetDof(cSecRef, childId, &cDof)); 4182 offsets[1] = cDof; 4183 } 4184 for (f = 0; f < lim; f++) { 4185 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 4186 PetscInt n = offsets[f + 1] - offsets[f]; 4187 PetscInt m = rowOffsets[f + 1] - rowOffsets[f]; 4188 PetscInt i, j; 4189 const PetscScalar *colValues = &childValues[offsets[f]]; 4190 4191 for (i = 0; i < m; i++) { 4192 PetscScalar val = 0.; 4193 for (j = 0; j < n; j++) val += childMat[n * i + j] * colValues[j]; 4194 parentValues[rowOffsets[f] + i] += val; 4195 } 4196 } 4197 } 4198 } 4199 if (contribute) PetscCall(VecSetValues(vecCoarse, dof, parentIndices, parentValues, INSERT_VALUES)); 4200 } 4201 PetscCall(PetscSectionDestroy(&multiRootSec)); 4202 PetscCall(PetscSectionDestroy(&rootIndicesSec)); 4203 PetscCall(PetscFree2(parentIndices, parentValues)); 4204 PetscCall(DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree, injRef, &childrenMats)); 4205 PetscCall(PetscFree(rootValues)); 4206 PetscCall(PetscFree3(offsets, offsetsCopy, rowOffsets)); 4207 PetscFunctionReturn(PETSC_SUCCESS); 4208 } 4209 4210 /*@ 4211 DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening 4212 that can be represented by a common reference tree used by both. This routine can be used for a combination of 4213 coarsening and refinement at the same time. 4214 4215 Collective 4216 4217 Input Parameters: 4218 + dmIn - The `DMPLEX` mesh for the input vector 4219 . dmOut - The second `DMPLEX` mesh 4220 . vecIn - The input vector 4221 . sfRefine - A star forest indicating points in the mesh `dmIn` (roots in the star forest) that are parents to points in 4222 the mesh `dmOut` (leaves in the star forest), i.e. where `dmOut` is more refined than `dmIn` 4223 . sfCoarsen - A star forest indicating points in the mesh `dmOut` (roots in the star forest) that are parents to points in 4224 the mesh `dmIn` (leaves in the star forest), i.e. where `dmOut` is more coarsened than `dmIn` 4225 . cidsRefine - The childIds of the points in `dmOut`. These childIds relate back to the reference tree: childid[j] = k implies 4226 that mesh point j of `dmOut` was refined from a point in `dmIn` just as the mesh point k in the reference 4227 tree was refined from its parent. childid[j] = -1 indicates that the point j in `dmOut` is exactly 4228 equivalent to its root in `dmIn`, so no interpolation is necessary. childid[j] = -2 indicates that this 4229 point j in `dmOut` is not a leaf of `sfRefine`. 4230 . cidsCoarsen - The childIds of the points in `dmIn`. These childIds relate back to the reference tree: childid[j] = k implies 4231 that mesh point j of dmIn coarsens to a point in `dmOut` just as the mesh point k in the reference 4232 tree coarsens to its parent. childid[j] = -2 indicates that point j in `dmOut` is not a leaf in `sfCoarsen`. 4233 . useBCs - `PETSC_TRUE` indicates that boundary values should be inserted into `vecIn` before transfer. 4234 - time - Used if boundary values are time dependent. 4235 4236 Output Parameter: 4237 . vecOut - Using interpolation and injection operators calculated on the reference tree, the transferred 4238 projection of `vecIn` from `dmIn` to `dmOut`. Note that any field discretized with a `PetscFV` finite volume 4239 method that uses gradient reconstruction will use reconstructed gradients when interpolating from 4240 coarse points to fine points. 4241 4242 Level: developer 4243 4244 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `Vec`, `PetscFV`, `DMPlexSetReferenceTree()`, `DMPlexGetReferenceTree()`, `PetscFVGetComputeGradients()` 4245 @*/ 4246 PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time) 4247 { 4248 PetscFunctionBegin; 4249 PetscCall(VecSet(vecOut, 0.0)); 4250 if (sfRefine) { 4251 Vec vecInLocal; 4252 DM dmGrad = NULL; 4253 Vec faceGeom = NULL, cellGeom = NULL, grad = NULL; 4254 4255 PetscCall(DMGetLocalVector(dmIn, &vecInLocal)); 4256 PetscCall(VecSet(vecInLocal, 0.0)); 4257 { 4258 PetscInt numFields, i; 4259 4260 PetscCall(DMGetNumFields(dmIn, &numFields)); 4261 for (i = 0; i < numFields; i++) { 4262 PetscObject obj; 4263 PetscClassId classid; 4264 4265 PetscCall(DMGetField(dmIn, i, NULL, &obj)); 4266 PetscCall(PetscObjectGetClassId(obj, &classid)); 4267 if (classid == PETSCFV_CLASSID) { 4268 PetscCall(DMPlexGetDataFVM(dmIn, (PetscFV)obj, &cellGeom, &faceGeom, &dmGrad)); 4269 break; 4270 } 4271 } 4272 } 4273 if (useBCs) PetscCall(DMPlexInsertBoundaryValues(dmIn, PETSC_TRUE, vecInLocal, time, faceGeom, cellGeom, NULL)); 4274 PetscCall(DMGlobalToLocalBegin(dmIn, vecIn, INSERT_VALUES, vecInLocal)); 4275 PetscCall(DMGlobalToLocalEnd(dmIn, vecIn, INSERT_VALUES, vecInLocal)); 4276 if (dmGrad) { 4277 PetscCall(DMGetGlobalVector(dmGrad, &grad)); 4278 PetscCall(DMPlexReconstructGradientsFVM(dmIn, vecInLocal, grad)); 4279 } 4280 PetscCall(DMPlexTransferVecTree_Interpolate(dmIn, vecInLocal, dmOut, vecOut, sfRefine, cidsRefine, grad, cellGeom)); 4281 PetscCall(DMRestoreLocalVector(dmIn, &vecInLocal)); 4282 if (dmGrad) PetscCall(DMRestoreGlobalVector(dmGrad, &grad)); 4283 } 4284 if (sfCoarsen) PetscCall(DMPlexTransferVecTree_Inject(dmIn, vecIn, dmOut, vecOut, sfCoarsen, cidsCoarsen)); 4285 PetscCall(VecAssemblyBegin(vecOut)); 4286 PetscCall(VecAssemblyEnd(vecOut)); 4287 PetscFunctionReturn(PETSC_SUCCESS); 4288 } 4289