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