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