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