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