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