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