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