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