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