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