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