1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/isimpl.h> 3 #include <petsc/private/petscfeimpl.h> 4 #include <petscsf.h> 5 #include <petscds.h> 6 7 /** hierarchy routines */ 8 9 /*@ 10 DMPlexSetReferenceTree - set the reference tree for hierarchically non-conforming meshes. 11 12 Not collective 13 14 Input Parameters: 15 + dm - The DMPlex object 16 - ref - The reference tree DMPlex object 17 18 Level: intermediate 19 20 .seealso: DMPlexGetReferenceTree(), DMPlexCreateDefaultReferenceTree() 21 @*/ 22 PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref) 23 { 24 DM_Plex *mesh = (DM_Plex *)dm->data; 25 PetscErrorCode ierr; 26 27 PetscFunctionBegin; 28 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 29 if (ref) {PetscValidHeaderSpecific(ref, DM_CLASSID, 2);} 30 ierr = PetscObjectReference((PetscObject)ref);CHKERRQ(ierr); 31 ierr = DMDestroy(&mesh->referenceTree);CHKERRQ(ierr); 32 mesh->referenceTree = ref; 33 PetscFunctionReturn(0); 34 } 35 36 /*@ 37 DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes. 38 39 Not collective 40 41 Input Parameters: 42 . dm - The DMPlex object 43 44 Output Parameters 45 . ref - The reference tree DMPlex object 46 47 Level: intermediate 48 49 .seealso: DMPlexSetReferenceTree(), DMPlexCreateDefaultReferenceTree() 50 @*/ 51 PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref) 52 { 53 DM_Plex *mesh = (DM_Plex *)dm->data; 54 55 PetscFunctionBegin; 56 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 57 PetscValidPointer(ref,2); 58 *ref = mesh->referenceTree; 59 PetscFunctionReturn(0); 60 } 61 62 static PetscErrorCode DMPlexReferenceTreeGetChildSymmetry_Default(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB) 63 { 64 PetscInt coneSize, dStart, dEnd, dim, ABswap, oAvert, oBvert, ABswapVert; 65 PetscErrorCode ierr; 66 67 PetscFunctionBegin; 68 if (parentOrientA == parentOrientB) { 69 if (childOrientB) *childOrientB = childOrientA; 70 if (childB) *childB = childA; 71 PetscFunctionReturn(0); 72 } 73 for (dim = 0; dim < 3; dim++) { 74 ierr = DMPlexGetDepthStratum(dm,dim,&dStart,&dEnd);CHKERRQ(ierr); 75 if (parent >= dStart && parent <= dEnd) { 76 break; 77 } 78 } 79 if (dim > 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot perform child symmetry for %d-cells",dim); 80 if (!dim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A vertex has no children"); 81 if (childA < dStart || childA >= dEnd) { 82 /* this is a lower-dimensional child: bootstrap */ 83 PetscInt size, i, sA = -1, sB, sOrientB, sConeSize; 84 const PetscInt *supp, *coneA, *coneB, *oA, *oB; 85 86 ierr = DMPlexGetSupportSize(dm,childA,&size);CHKERRQ(ierr); 87 ierr = DMPlexGetSupport(dm,childA,&supp);CHKERRQ(ierr); 88 89 /* find a point sA in supp(childA) that has the same parent */ 90 for (i = 0; i < size; i++) { 91 PetscInt sParent; 92 93 sA = supp[i]; 94 if (sA == parent) continue; 95 ierr = DMPlexGetTreeParent(dm,sA,&sParent,NULL);CHKERRQ(ierr); 96 if (sParent == parent) { 97 break; 98 } 99 } 100 if (i == size) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"could not find support in children"); 101 /* find out which point sB is in an equivalent position to sA under 102 * parentOrientB */ 103 ierr = DMPlexReferenceTreeGetChildSymmetry_Default(dm,parent,parentOrientA,0,sA,parentOrientB,&sOrientB,&sB);CHKERRQ(ierr); 104 ierr = DMPlexGetConeSize(dm,sA,&sConeSize);CHKERRQ(ierr); 105 ierr = DMPlexGetCone(dm,sA,&coneA);CHKERRQ(ierr); 106 ierr = DMPlexGetCone(dm,sB,&coneB);CHKERRQ(ierr); 107 ierr = DMPlexGetConeOrientation(dm,sA,&oA);CHKERRQ(ierr); 108 ierr = DMPlexGetConeOrientation(dm,sB,&oB);CHKERRQ(ierr); 109 /* step through the cone of sA in natural order */ 110 for (i = 0; i < sConeSize; i++) { 111 if (coneA[i] == childA) { 112 /* if childA is at position i in coneA, 113 * then we want the point that is at sOrientB*i in coneB */ 114 PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize -(sOrientB+1) - i) % sConeSize); 115 if (childB) *childB = coneB[j]; 116 if (childOrientB) { 117 PetscInt oBtrue; 118 119 ierr = DMPlexGetConeSize(dm,childA,&coneSize);CHKERRQ(ierr); 120 /* compose sOrientB and oB[j] */ 121 if (coneSize != 0 && coneSize != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected a vertex or an edge"); 122 /* we may have to flip an edge */ 123 oBtrue = coneSize ? ((sOrientB >= 0) ? oB[j] : -(oB[j] + 2)) : 0; 124 ABswap = DihedralSwap(coneSize,oA[i],oBtrue);CHKERRQ(ierr); 125 *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap); 126 } 127 break; 128 } 129 } 130 if (i == sConeSize) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"support cone mismatch"); 131 PetscFunctionReturn(0); 132 } 133 /* get the cone size and symmetry swap */ 134 ierr = DMPlexGetConeSize(dm,parent,&coneSize);CHKERRQ(ierr); 135 ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB); 136 if (dim == 2) { 137 /* orientations refer to cones: we want them to refer to vertices: 138 * if it's a rotation, they are the same, but if the order is reversed, a 139 * permutation that puts side i first does *not* put vertex i first */ 140 oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1); 141 oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1); 142 ABswapVert = DihedralSwap(coneSize, oAvert, oBvert); 143 } else { 144 ABswapVert = ABswap; 145 } 146 if (childB) { 147 /* assume that each child corresponds to a vertex, in the same order */ 148 PetscInt p, posA = -1, numChildren, i; 149 const PetscInt *children; 150 151 /* count which position the child is in */ 152 ierr = DMPlexGetTreeChildren(dm,parent,&numChildren,&children);CHKERRQ(ierr); 153 for (i = 0; i < numChildren; i++) { 154 p = children[i]; 155 if (p == childA) { 156 posA = i; 157 break; 158 } 159 } 160 if (posA >= coneSize) { 161 /* this is the triangle in the middle of a uniformly refined triangle: it 162 * is invariant */ 163 if (dim != 2 || posA != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Expected a middle triangle, got something else"); 164 *childB = childA; 165 } 166 else { 167 /* figure out position B by applying ABswapVert */ 168 PetscInt posB; 169 170 posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize -(ABswapVert + 1) - posA) % coneSize); 171 if (childB) *childB = children[posB]; 172 } 173 } 174 if (childOrientB) *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap); 175 PetscFunctionReturn(0); 176 } 177 178 /*@ 179 DMPlexReferenceTreeGetChildSymmetry - Given a reference tree, transform a childid and orientation from one parent frame to another 180 181 Input Parameters: 182 + dm - the reference tree DMPlex object 183 . parent - the parent point 184 . parentOrientA - the reference orientation for describing the parent 185 . childOrientA - the reference orientation for describing the child 186 . childA - the reference childID for describing the child 187 - parentOrientB - the new orientation for describing the parent 188 189 Output Parameters: 190 + childOrientB - if not NULL, set to the new oreintation for describing the child 191 - childB - if not NULL, the new childID for describing the child 192 193 Level: developer 194 195 .seealso: DMPlexGetReferenceTree(), DMPlexSetReferenceTree(), DMPlexSetTree() 196 @*/ 197 PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB) 198 { 199 DM_Plex *mesh = (DM_Plex *)dm->data; 200 PetscErrorCode ierr; 201 202 PetscFunctionBegin; 203 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 204 if (!mesh->getchildsymmetry) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlexReferenceTreeGetChildSymmetry not implemented"); 205 ierr = mesh->getchildsymmetry(dm,parent,parentOrientA,childOrientA,childA,parentOrientB,childOrientB,childB);CHKERRQ(ierr); 206 PetscFunctionReturn(0); 207 } 208 209 static PetscErrorCode DMPlexSetTree_Internal(DM,PetscSection,PetscInt*,PetscInt*,PetscBool,PetscBool); 210 211 PetscErrorCode DMPlexCreateReferenceTree_Union(DM K, DM Kref, const char *labelName, DM *ref) 212 { 213 MPI_Comm comm; 214 PetscInt dim, p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs; 215 PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts; 216 DMLabel identity, identityRef; 217 PetscSection unionSection, unionConeSection, parentSection; 218 PetscScalar *unionCoords; 219 IS perm; 220 PetscErrorCode ierr; 221 222 PetscFunctionBegin; 223 comm = PetscObjectComm((PetscObject)K); 224 ierr = DMGetDimension(K, &dim);CHKERRQ(ierr); 225 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr); 226 ierr = DMGetLabel(K, labelName, &identity);CHKERRQ(ierr); 227 ierr = DMGetLabel(Kref, labelName, &identityRef);CHKERRQ(ierr); 228 ierr = DMPlexGetChart(Kref, &pRefStart, &pRefEnd);CHKERRQ(ierr); 229 ierr = PetscSectionCreate(comm, &unionSection);CHKERRQ(ierr); 230 ierr = PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart));CHKERRQ(ierr); 231 /* count points that will go in the union */ 232 for (p = pStart; p < pEnd; p++) { 233 ierr = PetscSectionSetDof(unionSection, p - pStart, 1);CHKERRQ(ierr); 234 } 235 for (p = pRefStart; p < pRefEnd; p++) { 236 PetscInt q, qSize; 237 ierr = DMLabelGetValue(identityRef, p, &q);CHKERRQ(ierr); 238 ierr = DMLabelGetStratumSize(identityRef, q, &qSize);CHKERRQ(ierr); 239 if (qSize > 1) { 240 ierr = PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1);CHKERRQ(ierr); 241 } 242 } 243 ierr = PetscMalloc1(pEnd - pStart + pRefEnd - pRefStart,&permvals);CHKERRQ(ierr); 244 offset = 0; 245 /* stratify points in the union by topological dimension */ 246 for (d = 0; d <= dim; d++) { 247 PetscInt cStart, cEnd, c; 248 249 ierr = DMPlexGetHeightStratum(K, d, &cStart, &cEnd);CHKERRQ(ierr); 250 for (c = cStart; c < cEnd; c++) { 251 permvals[offset++] = c; 252 } 253 254 ierr = DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);CHKERRQ(ierr); 255 for (c = cStart; c < cEnd; c++) { 256 permvals[offset++] = c + (pEnd - pStart); 257 } 258 } 259 ierr = ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);CHKERRQ(ierr); 260 ierr = PetscSectionSetPermutation(unionSection,perm);CHKERRQ(ierr); 261 ierr = PetscSectionSetUp(unionSection);CHKERRQ(ierr); 262 ierr = PetscSectionGetStorageSize(unionSection,&numUnionPoints);CHKERRQ(ierr); 263 ierr = PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints);CHKERRQ(ierr); 264 /* count dimension points */ 265 for (d = 0; d <= dim; d++) { 266 PetscInt cStart, cOff, cOff2; 267 ierr = DMPlexGetHeightStratum(K,d,&cStart,NULL);CHKERRQ(ierr); 268 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff);CHKERRQ(ierr); 269 if (d < dim) { 270 ierr = DMPlexGetHeightStratum(K,d+1,&cStart,NULL);CHKERRQ(ierr); 271 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2);CHKERRQ(ierr); 272 } 273 else { 274 cOff2 = numUnionPoints; 275 } 276 numDimPoints[dim - d] = cOff2 - cOff; 277 } 278 ierr = PetscSectionCreate(comm, &unionConeSection);CHKERRQ(ierr); 279 ierr = PetscSectionSetChart(unionConeSection, 0, numUnionPoints);CHKERRQ(ierr); 280 /* count the cones in the union */ 281 for (p = pStart; p < pEnd; p++) { 282 PetscInt dof, uOff; 283 284 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 285 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 286 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 287 coneSizes[uOff] = dof; 288 } 289 for (p = pRefStart; p < pRefEnd; p++) { 290 PetscInt dof, uDof, uOff; 291 292 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 293 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 294 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 295 if (uDof) { 296 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 297 coneSizes[uOff] = dof; 298 } 299 } 300 ierr = PetscSectionSetUp(unionConeSection);CHKERRQ(ierr); 301 ierr = PetscSectionGetStorageSize(unionConeSection,&numCones);CHKERRQ(ierr); 302 ierr = PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations);CHKERRQ(ierr); 303 /* write the cones in the union */ 304 for (p = pStart; p < pEnd; p++) { 305 PetscInt dof, uOff, c, cOff; 306 const PetscInt *cone, *orientation; 307 308 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 309 ierr = DMPlexGetCone(K, p, &cone);CHKERRQ(ierr); 310 ierr = DMPlexGetConeOrientation(K, p, &orientation);CHKERRQ(ierr); 311 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 312 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 313 for (c = 0; c < dof; c++) { 314 PetscInt e, eOff; 315 e = cone[c]; 316 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 317 unionCones[cOff + c] = eOff; 318 unionOrientations[cOff + c] = orientation[c]; 319 } 320 } 321 for (p = pRefStart; p < pRefEnd; p++) { 322 PetscInt dof, uDof, uOff, c, cOff; 323 const PetscInt *cone, *orientation; 324 325 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 326 ierr = DMPlexGetCone(Kref, p, &cone);CHKERRQ(ierr); 327 ierr = DMPlexGetConeOrientation(Kref, p, &orientation);CHKERRQ(ierr); 328 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 329 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 330 if (uDof) { 331 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 332 for (c = 0; c < dof; c++) { 333 PetscInt e, eOff, eDof; 334 335 e = cone[c]; 336 ierr = PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof);CHKERRQ(ierr); 337 if (eDof) { 338 ierr = PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);CHKERRQ(ierr); 339 } 340 else { 341 ierr = DMLabelGetValue(identityRef, e, &e);CHKERRQ(ierr); 342 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 343 } 344 unionCones[cOff + c] = eOff; 345 unionOrientations[cOff + c] = orientation[c]; 346 } 347 } 348 } 349 /* get the coordinates */ 350 { 351 PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff; 352 PetscSection KcoordsSec, KrefCoordsSec; 353 Vec KcoordsVec, KrefCoordsVec; 354 PetscScalar *Kcoords; 355 356 DMGetCoordinateSection(K, &KcoordsSec);CHKERRQ(ierr); 357 DMGetCoordinatesLocal(K, &KcoordsVec);CHKERRQ(ierr); 358 DMGetCoordinateSection(Kref, &KrefCoordsSec);CHKERRQ(ierr); 359 DMGetCoordinatesLocal(Kref, &KrefCoordsVec);CHKERRQ(ierr); 360 361 numVerts = numDimPoints[0]; 362 ierr = PetscMalloc1(numVerts * dim,&unionCoords);CHKERRQ(ierr); 363 ierr = DMPlexGetDepthStratum(K,0,&vStart,&vEnd);CHKERRQ(ierr); 364 365 offset = 0; 366 for (v = vStart; v < vEnd; v++) { 367 ierr = PetscSectionGetOffset(unionSection,v - pStart,&vOff);CHKERRQ(ierr); 368 ierr = VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);CHKERRQ(ierr); 369 for (d = 0; d < dim; d++) { 370 unionCoords[offset * dim + d] = Kcoords[d]; 371 } 372 offset++; 373 } 374 ierr = DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd);CHKERRQ(ierr); 375 for (v = vRefStart; v < vRefEnd; v++) { 376 ierr = PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof);CHKERRQ(ierr); 377 ierr = PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff);CHKERRQ(ierr); 378 ierr = VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);CHKERRQ(ierr); 379 if (vDof) { 380 for (d = 0; d < dim; d++) { 381 unionCoords[offset * dim + d] = Kcoords[d]; 382 } 383 offset++; 384 } 385 } 386 } 387 ierr = DMCreate(comm,ref);CHKERRQ(ierr); 388 ierr = DMSetType(*ref,DMPLEX);CHKERRQ(ierr); 389 ierr = DMSetDimension(*ref,dim);CHKERRQ(ierr); 390 ierr = DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords);CHKERRQ(ierr); 391 /* set the tree */ 392 ierr = PetscSectionCreate(comm,&parentSection);CHKERRQ(ierr); 393 ierr = PetscSectionSetChart(parentSection,0,numUnionPoints);CHKERRQ(ierr); 394 for (p = pRefStart; p < pRefEnd; p++) { 395 PetscInt uDof, uOff; 396 397 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 398 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 399 if (uDof) { 400 PetscSectionSetDof(parentSection,uOff,1);CHKERRQ(ierr); 401 } 402 } 403 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 404 ierr = PetscSectionGetStorageSize(parentSection,&parentSize);CHKERRQ(ierr); 405 ierr = PetscMalloc2(parentSize,&parents,parentSize,&childIDs);CHKERRQ(ierr); 406 for (p = pRefStart; p < pRefEnd; p++) { 407 PetscInt uDof, uOff; 408 409 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 410 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 411 if (uDof) { 412 PetscInt pOff, parent, parentU; 413 PetscSectionGetOffset(parentSection,uOff,&pOff);CHKERRQ(ierr); 414 DMLabelGetValue(identityRef,p,&parent);CHKERRQ(ierr); 415 ierr = PetscSectionGetOffset(unionSection, parent - pStart,&parentU);CHKERRQ(ierr); 416 parents[pOff] = parentU; 417 childIDs[pOff] = uOff; 418 } 419 } 420 ierr = DMPlexSetTree_Internal(*ref,parentSection,parents,childIDs,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 421 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 422 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 423 424 /* clean up */ 425 ierr = PetscSectionDestroy(&unionSection);CHKERRQ(ierr); 426 ierr = PetscSectionDestroy(&unionConeSection);CHKERRQ(ierr); 427 ierr = ISDestroy(&perm);CHKERRQ(ierr); 428 ierr = PetscFree(unionCoords);CHKERRQ(ierr); 429 ierr = PetscFree2(unionCones,unionOrientations);CHKERRQ(ierr); 430 ierr = PetscFree2(coneSizes,numDimPoints);CHKERRQ(ierr); 431 PetscFunctionReturn(0); 432 } 433 434 /*@ 435 DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement. 436 437 Collective on comm 438 439 Input Parameters: 440 + comm - the MPI communicator 441 . dim - the spatial dimension 442 - simplex - Flag for simplex, otherwise use a tensor-product cell 443 444 Output Parameters: 445 . ref - the reference tree DMPlex object 446 447 Level: intermediate 448 449 .keywords: reference cell 450 .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree() 451 @*/ 452 PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref) 453 { 454 DM_Plex *mesh; 455 DM K, Kref; 456 PetscInt p, pStart, pEnd; 457 DMLabel identity; 458 PetscErrorCode ierr; 459 460 PetscFunctionBegin; 461 #if 1 462 comm = PETSC_COMM_SELF; 463 #endif 464 /* create a reference element */ 465 ierr = DMPlexCreateReferenceCell(comm, dim, simplex, &K);CHKERRQ(ierr); 466 ierr = DMCreateLabel(K, "identity");CHKERRQ(ierr); 467 ierr = DMGetLabel(K, "identity", &identity);CHKERRQ(ierr); 468 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr); 469 for (p = pStart; p < pEnd; p++) { 470 ierr = DMLabelSetValue(identity, p, p);CHKERRQ(ierr); 471 } 472 /* refine it */ 473 ierr = DMRefine(K,comm,&Kref);CHKERRQ(ierr); 474 475 /* the reference tree is the union of these two, without duplicating 476 * points that appear in both */ 477 ierr = DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref);CHKERRQ(ierr); 478 mesh = (DM_Plex *) (*ref)->data; 479 mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default; 480 ierr = DMDestroy(&K);CHKERRQ(ierr); 481 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 482 PetscFunctionReturn(0); 483 } 484 485 static PetscErrorCode DMPlexTreeSymmetrize(DM dm) 486 { 487 DM_Plex *mesh = (DM_Plex *)dm->data; 488 PetscSection childSec, pSec; 489 PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT; 490 PetscInt *offsets, *children, pStart, pEnd; 491 PetscErrorCode ierr; 492 493 PetscFunctionBegin; 494 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 495 ierr = PetscSectionDestroy(&mesh->childSection);CHKERRQ(ierr); 496 ierr = PetscFree(mesh->children);CHKERRQ(ierr); 497 pSec = mesh->parentSection; 498 if (!pSec) PetscFunctionReturn(0); 499 ierr = PetscSectionGetStorageSize(pSec,&pSize);CHKERRQ(ierr); 500 for (p = 0; p < pSize; p++) { 501 PetscInt par = mesh->parents[p]; 502 503 parMax = PetscMax(parMax,par+1); 504 parMin = PetscMin(parMin,par); 505 } 506 if (parMin > parMax) { 507 parMin = -1; 508 parMax = -1; 509 } 510 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);CHKERRQ(ierr); 511 ierr = PetscSectionSetChart(childSec,parMin,parMax);CHKERRQ(ierr); 512 for (p = 0; p < pSize; p++) { 513 PetscInt par = mesh->parents[p]; 514 515 ierr = PetscSectionAddDof(childSec,par,1);CHKERRQ(ierr); 516 } 517 ierr = PetscSectionSetUp(childSec);CHKERRQ(ierr); 518 ierr = PetscSectionGetStorageSize(childSec,&cSize);CHKERRQ(ierr); 519 ierr = PetscMalloc1(cSize,&children);CHKERRQ(ierr); 520 ierr = PetscCalloc1(parMax-parMin,&offsets);CHKERRQ(ierr); 521 ierr = PetscSectionGetChart(pSec,&pStart,&pEnd);CHKERRQ(ierr); 522 for (p = pStart; p < pEnd; p++) { 523 PetscInt dof, off, i; 524 525 ierr = PetscSectionGetDof(pSec,p,&dof);CHKERRQ(ierr); 526 ierr = PetscSectionGetOffset(pSec,p,&off);CHKERRQ(ierr); 527 for (i = 0; i < dof; i++) { 528 PetscInt par = mesh->parents[off + i], cOff; 529 530 ierr = PetscSectionGetOffset(childSec,par,&cOff);CHKERRQ(ierr); 531 children[cOff + offsets[par-parMin]++] = p; 532 } 533 } 534 mesh->childSection = childSec; 535 mesh->children = children; 536 ierr = PetscFree(offsets);CHKERRQ(ierr); 537 PetscFunctionReturn(0); 538 } 539 540 static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew) 541 { 542 PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL; 543 const PetscInt *vals; 544 PetscSection secNew; 545 PetscBool anyNew, globalAnyNew; 546 PetscBool compress; 547 PetscErrorCode ierr; 548 549 PetscFunctionBegin; 550 ierr = PetscSectionGetChart(section,&pStart,&pEnd);CHKERRQ(ierr); 551 ierr = ISGetLocalSize(is,&size);CHKERRQ(ierr); 552 ierr = ISGetIndices(is,&vals);CHKERRQ(ierr); 553 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);CHKERRQ(ierr); 554 ierr = PetscSectionSetChart(secNew,pStart,pEnd);CHKERRQ(ierr); 555 for (i = 0; i < size; i++) { 556 PetscInt dof; 557 558 p = vals[i]; 559 if (p < pStart || p >= pEnd) continue; 560 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 561 if (dof) break; 562 } 563 if (i == size) { 564 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 565 anyNew = PETSC_FALSE; 566 compress = PETSC_FALSE; 567 sizeNew = 0; 568 } 569 else { 570 anyNew = PETSC_TRUE; 571 for (p = pStart; p < pEnd; p++) { 572 PetscInt dof, off; 573 574 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 575 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 576 for (i = 0; i < dof; i++) { 577 PetscInt q = vals[off + i], qDof = 0; 578 579 if (q >= pStart && q < pEnd) { 580 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 581 } 582 if (qDof) { 583 ierr = PetscSectionAddDof(secNew, p, qDof);CHKERRQ(ierr); 584 } 585 else { 586 ierr = PetscSectionAddDof(secNew, p, 1);CHKERRQ(ierr); 587 } 588 } 589 } 590 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 591 ierr = PetscSectionGetStorageSize(secNew,&sizeNew);CHKERRQ(ierr); 592 ierr = PetscMalloc1(sizeNew,&valsNew);CHKERRQ(ierr); 593 compress = PETSC_FALSE; 594 for (p = pStart; p < pEnd; p++) { 595 PetscInt dof, off, count, offNew, dofNew; 596 597 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 598 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 599 ierr = PetscSectionGetDof(secNew, p, &dofNew);CHKERRQ(ierr); 600 ierr = PetscSectionGetOffset(secNew, p, &offNew);CHKERRQ(ierr); 601 count = 0; 602 for (i = 0; i < dof; i++) { 603 PetscInt q = vals[off + i], qDof = 0, qOff = 0, j; 604 605 if (q >= pStart && q < pEnd) { 606 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 607 ierr = PetscSectionGetOffset(section, q, &qOff);CHKERRQ(ierr); 608 } 609 if (qDof) { 610 PetscInt oldCount = count; 611 612 for (j = 0; j < qDof; j++) { 613 PetscInt k, r = vals[qOff + j]; 614 615 for (k = 0; k < oldCount; k++) { 616 if (valsNew[offNew + k] == r) { 617 break; 618 } 619 } 620 if (k == oldCount) { 621 valsNew[offNew + count++] = r; 622 } 623 } 624 } 625 else { 626 PetscInt k, oldCount = count; 627 628 for (k = 0; k < oldCount; k++) { 629 if (valsNew[offNew + k] == q) { 630 break; 631 } 632 } 633 if (k == oldCount) { 634 valsNew[offNew + count++] = q; 635 } 636 } 637 } 638 if (count < dofNew) { 639 ierr = PetscSectionSetDof(secNew, p, count);CHKERRQ(ierr); 640 compress = PETSC_TRUE; 641 } 642 } 643 } 644 ierr = ISRestoreIndices(is,&vals);CHKERRQ(ierr); 645 ierr = MPIU_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 646 if (!globalAnyNew) { 647 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 648 *sectionNew = NULL; 649 *isNew = NULL; 650 } 651 else { 652 PetscBool globalCompress; 653 654 ierr = MPIU_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 655 if (compress) { 656 PetscSection secComp; 657 PetscInt *valsComp = NULL; 658 659 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);CHKERRQ(ierr); 660 ierr = PetscSectionSetChart(secComp,pStart,pEnd);CHKERRQ(ierr); 661 for (p = pStart; p < pEnd; p++) { 662 PetscInt dof; 663 664 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 665 ierr = PetscSectionSetDof(secComp, p, dof);CHKERRQ(ierr); 666 } 667 ierr = PetscSectionSetUp(secComp);CHKERRQ(ierr); 668 ierr = PetscSectionGetStorageSize(secComp,&sizeNew);CHKERRQ(ierr); 669 ierr = PetscMalloc1(sizeNew,&valsComp);CHKERRQ(ierr); 670 for (p = pStart; p < pEnd; p++) { 671 PetscInt dof, off, offNew, j; 672 673 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 674 ierr = PetscSectionGetOffset(secNew, p, &off);CHKERRQ(ierr); 675 ierr = PetscSectionGetOffset(secComp, p, &offNew);CHKERRQ(ierr); 676 for (j = 0; j < dof; j++) { 677 valsComp[offNew + j] = valsNew[off + j]; 678 } 679 } 680 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 681 secNew = secComp; 682 ierr = PetscFree(valsNew);CHKERRQ(ierr); 683 valsNew = valsComp; 684 } 685 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);CHKERRQ(ierr); 686 } 687 PetscFunctionReturn(0); 688 } 689 690 static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm) 691 { 692 PetscInt p, pStart, pEnd, *anchors, size; 693 PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT; 694 PetscSection aSec; 695 DMLabel canonLabel; 696 IS aIS; 697 PetscErrorCode ierr; 698 699 PetscFunctionBegin; 700 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 701 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 702 ierr = DMGetLabel(dm,"canonical",&canonLabel);CHKERRQ(ierr); 703 for (p = pStart; p < pEnd; p++) { 704 PetscInt parent; 705 706 if (canonLabel) { 707 PetscInt canon; 708 709 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 710 if (p != canon) continue; 711 } 712 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 713 if (parent != p) { 714 aMin = PetscMin(aMin,p); 715 aMax = PetscMax(aMax,p+1); 716 } 717 } 718 if (aMin > aMax) { 719 aMin = -1; 720 aMax = -1; 721 } 722 ierr = PetscSectionCreate(PETSC_COMM_SELF,&aSec);CHKERRQ(ierr); 723 ierr = PetscSectionSetChart(aSec,aMin,aMax);CHKERRQ(ierr); 724 for (p = aMin; p < aMax; p++) { 725 PetscInt parent, ancestor = p; 726 727 if (canonLabel) { 728 PetscInt canon; 729 730 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 731 if (p != canon) continue; 732 } 733 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 734 while (parent != ancestor) { 735 ancestor = parent; 736 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 737 } 738 if (ancestor != p) { 739 PetscInt closureSize, *closure = NULL; 740 741 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 742 ierr = PetscSectionSetDof(aSec,p,closureSize);CHKERRQ(ierr); 743 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 744 } 745 } 746 ierr = PetscSectionSetUp(aSec);CHKERRQ(ierr); 747 ierr = PetscSectionGetStorageSize(aSec,&size);CHKERRQ(ierr); 748 ierr = PetscMalloc1(size,&anchors);CHKERRQ(ierr); 749 for (p = aMin; p < aMax; p++) { 750 PetscInt parent, ancestor = p; 751 752 if (canonLabel) { 753 PetscInt canon; 754 755 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 756 if (p != canon) continue; 757 } 758 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 759 while (parent != ancestor) { 760 ancestor = parent; 761 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 762 } 763 if (ancestor != p) { 764 PetscInt j, closureSize, *closure = NULL, aOff; 765 766 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 767 768 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 769 for (j = 0; j < closureSize; j++) { 770 anchors[aOff + j] = closure[2*j]; 771 } 772 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 773 } 774 } 775 ierr = ISCreateGeneral(PETSC_COMM_SELF,size,anchors,PETSC_OWN_POINTER,&aIS);CHKERRQ(ierr); 776 { 777 PetscSection aSecNew = aSec; 778 IS aISNew = aIS; 779 780 ierr = PetscObjectReference((PetscObject)aSec);CHKERRQ(ierr); 781 ierr = PetscObjectReference((PetscObject)aIS);CHKERRQ(ierr); 782 while (aSecNew) { 783 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 784 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 785 aSec = aSecNew; 786 aIS = aISNew; 787 aSecNew = NULL; 788 aISNew = NULL; 789 ierr = AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);CHKERRQ(ierr); 790 } 791 } 792 ierr = DMPlexSetAnchors(dm,aSec,aIS);CHKERRQ(ierr); 793 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 794 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 795 PetscFunctionReturn(0); 796 } 797 798 static PetscErrorCode DMPlexGetTrueSupportSize(DM dm,PetscInt p,PetscInt *dof,PetscInt *numTrueSupp) 799 { 800 PetscErrorCode ierr; 801 802 PetscFunctionBegin; 803 if (numTrueSupp[p] == -1) { 804 PetscInt i, alldof; 805 const PetscInt *supp; 806 PetscInt count = 0; 807 808 ierr = DMPlexGetSupportSize(dm,p,&alldof);CHKERRQ(ierr); 809 ierr = DMPlexGetSupport(dm,p,&supp);CHKERRQ(ierr); 810 for (i = 0; i < alldof; i++) { 811 PetscInt q = supp[i], numCones, j; 812 const PetscInt *cone; 813 814 ierr = DMPlexGetConeSize(dm,q,&numCones);CHKERRQ(ierr); 815 ierr = DMPlexGetCone(dm,q,&cone);CHKERRQ(ierr); 816 for (j = 0; j < numCones; j++) { 817 if (cone[j] == p) break; 818 } 819 if (j < numCones) count++; 820 } 821 numTrueSupp[p] = count; 822 } 823 *dof = numTrueSupp[p]; 824 PetscFunctionReturn(0); 825 } 826 827 static PetscErrorCode DMPlexTreeExchangeSupports(DM dm) 828 { 829 DM_Plex *mesh = (DM_Plex *)dm->data; 830 PetscSection newSupportSection; 831 PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth; 832 PetscInt *numTrueSupp; 833 PetscInt *offsets; 834 PetscErrorCode ierr; 835 836 PetscFunctionBegin; 837 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 838 /* symmetrize the hierarchy */ 839 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 840 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)),&newSupportSection);CHKERRQ(ierr); 841 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 842 ierr = PetscSectionSetChart(newSupportSection,pStart,pEnd);CHKERRQ(ierr); 843 ierr = PetscCalloc1(pEnd,&offsets);CHKERRQ(ierr); 844 ierr = PetscMalloc1(pEnd,&numTrueSupp);CHKERRQ(ierr); 845 for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1; 846 /* if a point is in the (true) support of q, it should be in the support of 847 * parent(q) */ 848 for (d = 0; d <= depth; d++) { 849 ierr = DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);CHKERRQ(ierr); 850 for (p = pStart; p < pEnd; ++p) { 851 PetscInt dof, q, qdof, parent; 852 853 ierr = DMPlexGetTrueSupportSize(dm,p,&dof,numTrueSupp);CHKERRQ(ierr); 854 ierr = PetscSectionAddDof(newSupportSection, p, dof);CHKERRQ(ierr); 855 q = p; 856 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 857 while (parent != q && parent >= pStart && parent < pEnd) { 858 q = parent; 859 860 ierr = DMPlexGetTrueSupportSize(dm,q,&qdof,numTrueSupp);CHKERRQ(ierr); 861 ierr = PetscSectionAddDof(newSupportSection,p,qdof);CHKERRQ(ierr); 862 ierr = PetscSectionAddDof(newSupportSection,q,dof);CHKERRQ(ierr); 863 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 864 } 865 } 866 } 867 ierr = PetscSectionSetUp(newSupportSection);CHKERRQ(ierr); 868 ierr = PetscSectionGetStorageSize(newSupportSection,&newSize);CHKERRQ(ierr); 869 ierr = PetscMalloc1(newSize,&newSupports);CHKERRQ(ierr); 870 for (d = 0; d <= depth; d++) { 871 ierr = DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);CHKERRQ(ierr); 872 for (p = pStart; p < pEnd; p++) { 873 PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent; 874 875 ierr = PetscSectionGetDof(mesh->supportSection, p, &dof);CHKERRQ(ierr); 876 ierr = PetscSectionGetOffset(mesh->supportSection, p, &off);CHKERRQ(ierr); 877 ierr = PetscSectionGetDof(newSupportSection, p, &newDof);CHKERRQ(ierr); 878 ierr = PetscSectionGetOffset(newSupportSection, p, &newOff);CHKERRQ(ierr); 879 for (i = 0; i < dof; i++) { 880 PetscInt numCones, j; 881 const PetscInt *cone; 882 PetscInt q = mesh->supports[off + i]; 883 884 ierr = DMPlexGetConeSize(dm,q,&numCones);CHKERRQ(ierr); 885 ierr = DMPlexGetCone(dm,q,&cone);CHKERRQ(ierr); 886 for (j = 0; j < numCones; j++) { 887 if (cone[j] == p) break; 888 } 889 if (j < numCones) newSupports[newOff+offsets[p]++] = q; 890 } 891 mesh->maxSupportSize = PetscMax(mesh->maxSupportSize,newDof); 892 893 q = p; 894 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 895 while (parent != q && parent >= pStart && parent < pEnd) { 896 q = parent; 897 ierr = PetscSectionGetDof(mesh->supportSection, q, &qdof);CHKERRQ(ierr); 898 ierr = PetscSectionGetOffset(mesh->supportSection, q, &qoff);CHKERRQ(ierr); 899 ierr = PetscSectionGetOffset(newSupportSection, q, &newqOff);CHKERRQ(ierr); 900 for (i = 0; i < qdof; i++) { 901 PetscInt numCones, j; 902 const PetscInt *cone; 903 PetscInt r = mesh->supports[qoff + i]; 904 905 ierr = DMPlexGetConeSize(dm,r,&numCones);CHKERRQ(ierr); 906 ierr = DMPlexGetCone(dm,r,&cone);CHKERRQ(ierr); 907 for (j = 0; j < numCones; j++) { 908 if (cone[j] == q) break; 909 } 910 if (j < numCones) newSupports[newOff+offsets[p]++] = r; 911 } 912 for (i = 0; i < dof; i++) { 913 PetscInt numCones, j; 914 const PetscInt *cone; 915 PetscInt r = mesh->supports[off + i]; 916 917 ierr = DMPlexGetConeSize(dm,r,&numCones);CHKERRQ(ierr); 918 ierr = DMPlexGetCone(dm,r,&cone);CHKERRQ(ierr); 919 for (j = 0; j < numCones; j++) { 920 if (cone[j] == p) break; 921 } 922 if (j < numCones) newSupports[newqOff+offsets[q]++] = r; 923 } 924 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 925 } 926 } 927 } 928 ierr = PetscSectionDestroy(&mesh->supportSection);CHKERRQ(ierr); 929 mesh->supportSection = newSupportSection; 930 ierr = PetscFree(mesh->supports);CHKERRQ(ierr); 931 mesh->supports = newSupports; 932 ierr = PetscFree(offsets);CHKERRQ(ierr); 933 ierr = PetscFree(numTrueSupp);CHKERRQ(ierr); 934 935 PetscFunctionReturn(0); 936 } 937 938 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM,PetscSection,PetscSection,Mat); 939 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM,PetscSection,PetscSection,Mat); 940 941 static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports) 942 { 943 DM_Plex *mesh = (DM_Plex *)dm->data; 944 DM refTree; 945 PetscInt size; 946 PetscErrorCode ierr; 947 948 PetscFunctionBegin; 949 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 950 PetscValidHeaderSpecific(parentSection, PETSC_SECTION_CLASSID, 2); 951 ierr = PetscObjectReference((PetscObject)parentSection);CHKERRQ(ierr); 952 ierr = PetscSectionDestroy(&mesh->parentSection);CHKERRQ(ierr); 953 mesh->parentSection = parentSection; 954 ierr = PetscSectionGetStorageSize(parentSection,&size);CHKERRQ(ierr); 955 if (parents != mesh->parents) { 956 ierr = PetscFree(mesh->parents);CHKERRQ(ierr); 957 ierr = PetscMalloc1(size,&mesh->parents);CHKERRQ(ierr); 958 ierr = PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));CHKERRQ(ierr); 959 } 960 if (childIDs != mesh->childIDs) { 961 ierr = PetscFree(mesh->childIDs);CHKERRQ(ierr); 962 ierr = PetscMalloc1(size,&mesh->childIDs);CHKERRQ(ierr); 963 ierr = PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));CHKERRQ(ierr); 964 } 965 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr); 966 if (refTree) { 967 DMLabel canonLabel; 968 969 ierr = DMGetLabel(refTree,"canonical",&canonLabel);CHKERRQ(ierr); 970 if (canonLabel) { 971 PetscInt i; 972 973 for (i = 0; i < size; i++) { 974 PetscInt canon; 975 ierr = DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);CHKERRQ(ierr); 976 if (canon >= 0) { 977 mesh->childIDs[i] = canon; 978 } 979 } 980 } 981 mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference; 982 } 983 else { 984 mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct; 985 } 986 ierr = DMPlexTreeSymmetrize(dm);CHKERRQ(ierr); 987 if (computeCanonical) { 988 PetscInt d, dim; 989 990 /* add the canonical label */ 991 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 992 ierr = DMCreateLabel(dm,"canonical");CHKERRQ(ierr); 993 for (d = 0; d <= dim; d++) { 994 PetscInt p, dStart, dEnd, canon = -1, cNumChildren; 995 const PetscInt *cChildren; 996 997 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr); 998 for (p = dStart; p < dEnd; p++) { 999 ierr = DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);CHKERRQ(ierr); 1000 if (cNumChildren) { 1001 canon = p; 1002 break; 1003 } 1004 } 1005 if (canon == -1) continue; 1006 for (p = dStart; p < dEnd; p++) { 1007 PetscInt numChildren, i; 1008 const PetscInt *children; 1009 1010 ierr = DMPlexGetTreeChildren(dm,p,&numChildren,&children);CHKERRQ(ierr); 1011 if (numChildren) { 1012 if (numChildren != cNumChildren) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"All parent points in a stratum should have the same number of children: %d != %d", numChildren, cNumChildren); 1013 ierr = DMSetLabelValue(dm,"canonical",p,canon);CHKERRQ(ierr); 1014 for (i = 0; i < numChildren; i++) { 1015 ierr = DMSetLabelValue(dm,"canonical",children[i],cChildren[i]);CHKERRQ(ierr); 1016 } 1017 } 1018 } 1019 } 1020 } 1021 if (exchangeSupports) { 1022 ierr = DMPlexTreeExchangeSupports(dm);CHKERRQ(ierr); 1023 } 1024 mesh->createanchors = DMPlexCreateAnchors_Tree; 1025 /* reset anchors */ 1026 ierr = DMPlexSetAnchors(dm,NULL,NULL);CHKERRQ(ierr); 1027 PetscFunctionReturn(0); 1028 } 1029 1030 /*@ 1031 DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates 1032 the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its 1033 tree root. 1034 1035 Collective on dm 1036 1037 Input Parameters: 1038 + dm - the DMPlex object 1039 . parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 1040 offset indexes the parent and childID list; the reference count of parentSection is incremented 1041 . parents - a list of the point parents; copied, can be destroyed 1042 - childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 1043 the child corresponds to the point in the reference tree with index childIDs; copied, can be destroyed 1044 1045 Level: intermediate 1046 1047 .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 1048 @*/ 1049 PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[]) 1050 { 1051 PetscErrorCode ierr; 1052 1053 PetscFunctionBegin; 1054 ierr = DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE);CHKERRQ(ierr); 1055 PetscFunctionReturn(0); 1056 } 1057 1058 /*@ 1059 DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points. 1060 Collective on dm 1061 1062 Input Parameters: 1063 . dm - the DMPlex object 1064 1065 Output Parameters: 1066 + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 1067 offset indexes the parent and childID list 1068 . parents - a list of the point parents 1069 . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 1070 the child corresponds to the point in the reference tree with index childID 1071 . childSection - the inverse of the parent section 1072 - children - a list of the point children 1073 1074 Level: intermediate 1075 1076 .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 1077 @*/ 1078 PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[]) 1079 { 1080 DM_Plex *mesh = (DM_Plex *)dm->data; 1081 1082 PetscFunctionBegin; 1083 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1084 if (parentSection) *parentSection = mesh->parentSection; 1085 if (parents) *parents = mesh->parents; 1086 if (childIDs) *childIDs = mesh->childIDs; 1087 if (childSection) *childSection = mesh->childSection; 1088 if (children) *children = mesh->children; 1089 PetscFunctionReturn(0); 1090 } 1091 1092 /*@ 1093 DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the 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 = PetscSpaceSetOrder(bspace,0);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 = DMGetDefaultSection(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);CHKERRQ(ierr); 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);CHKERRQ(ierr); 2119 for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l]; 2120 offset += dim; 2121 } 2122 } 2123 ierr = VecRestoreArray(coords,&coordvals);CHKERRQ(ierr); 2124 } 2125 2126 /* need to reverse the order of pNewCount: vertices first, cells last */ 2127 for (d = 0; d < (dim + 1) / 2; d++) { 2128 PetscInt tmp; 2129 2130 tmp = pNewCount[d]; 2131 pNewCount[d] = pNewCount[dim - d]; 2132 pNewCount[dim - d] = tmp; 2133 } 2134 2135 ierr = DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);CHKERRQ(ierr); 2136 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr); 2137 ierr = DMPlexSetTree(*ncdm,parentSection,parents,childIDs);CHKERRQ(ierr); 2138 2139 /* clean up */ 2140 ierr = DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);CHKERRQ(ierr); 2141 ierr = PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);CHKERRQ(ierr); 2142 ierr = PetscFree(newConeSizes);CHKERRQ(ierr); 2143 ierr = PetscFree2(newCones,newOrientations);CHKERRQ(ierr); 2144 ierr = PetscFree(newVertexCoords);CHKERRQ(ierr); 2145 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 2146 ierr = PetscFree4(Kembedding,perm,iperm,preOrient);CHKERRQ(ierr); 2147 } 2148 else { 2149 PetscInt p, counts[4]; 2150 PetscInt *coneSizes, *cones, *orientations; 2151 Vec coordVec; 2152 PetscScalar *coords; 2153 2154 for (d = 0; d <= dim; d++) { 2155 PetscInt dStart, dEnd; 2156 2157 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr); 2158 counts[d] = dEnd - dStart; 2159 } 2160 ierr = PetscMalloc1(pEnd-pStart,&coneSizes);CHKERRQ(ierr); 2161 for (p = pStart; p < pEnd; p++) { 2162 ierr = DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);CHKERRQ(ierr); 2163 } 2164 ierr = DMPlexGetCones(dm, &cones);CHKERRQ(ierr); 2165 ierr = DMPlexGetConeOrientations(dm, &orientations);CHKERRQ(ierr); 2166 ierr = DMGetCoordinatesLocal(dm,&coordVec);CHKERRQ(ierr); 2167 ierr = VecGetArray(coordVec,&coords);CHKERRQ(ierr); 2168 2169 ierr = PetscSectionSetChart(parentSection,pStart,pEnd);CHKERRQ(ierr); 2170 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 2171 ierr = DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);CHKERRQ(ierr); 2172 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr); 2173 ierr = DMPlexSetTree(*ncdm,parentSection,NULL,NULL);CHKERRQ(ierr); 2174 ierr = VecRestoreArray(coordVec,&coords);CHKERRQ(ierr); 2175 } 2176 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 2177 2178 PetscFunctionReturn(0); 2179 } 2180 2181 PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 2182 { 2183 PetscSF coarseToFineEmbedded; 2184 PetscSection globalCoarse, globalFine; 2185 PetscSection localCoarse, localFine; 2186 PetscSection aSec, cSec; 2187 PetscSection rootIndicesSec, rootMatricesSec; 2188 PetscSection leafIndicesSec, leafMatricesSec; 2189 PetscInt *rootIndices, *leafIndices; 2190 PetscScalar *rootMatrices, *leafMatrices; 2191 IS aIS; 2192 const PetscInt *anchors; 2193 Mat cMat; 2194 PetscInt numFields, maxFields; 2195 PetscInt pStartC, pEndC, pStartF, pEndF, p; 2196 PetscInt aStart, aEnd, cStart, cEnd; 2197 PetscInt *maxChildIds; 2198 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 2199 const PetscInt ***perms; 2200 const PetscScalar ***flips; 2201 PetscErrorCode ierr; 2202 2203 PetscFunctionBegin; 2204 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 2205 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 2206 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 2207 { /* winnow fine points that don't have global dofs out of the sf */ 2208 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l; 2209 const PetscInt *leaves; 2210 2211 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 2212 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 2213 p = leaves ? leaves[l] : l; 2214 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 2215 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 2216 if ((dof - cdof) > 0) { 2217 numPointsWithDofs++; 2218 } 2219 } 2220 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 2221 for (l = 0, offset = 0; l < nleaves; l++) { 2222 p = leaves ? leaves[l] : l; 2223 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 2224 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 2225 if ((dof - cdof) > 0) { 2226 pointsWithDofs[offset++] = l; 2227 } 2228 } 2229 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 2230 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 2231 } 2232 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 2233 ierr = PetscMalloc1(pEndC-pStartC,&maxChildIds);CHKERRQ(ierr); 2234 for (p = pStartC; p < pEndC; p++) { 2235 maxChildIds[p - pStartC] = -2; 2236 } 2237 ierr = PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 2238 ierr = PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 2239 2240 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 2241 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 2242 2243 ierr = DMPlexGetAnchors(coarse,&aSec,&aIS);CHKERRQ(ierr); 2244 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 2245 ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr); 2246 2247 ierr = DMGetDefaultConstraints(coarse,&cSec,&cMat);CHKERRQ(ierr); 2248 ierr = PetscSectionGetChart(cSec,&cStart,&cEnd);CHKERRQ(ierr); 2249 2250 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 2251 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);CHKERRQ(ierr); 2252 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec);CHKERRQ(ierr); 2253 ierr = PetscSectionSetChart(rootIndicesSec,pStartC,pEndC);CHKERRQ(ierr); 2254 ierr = PetscSectionSetChart(rootMatricesSec,pStartC,pEndC);CHKERRQ(ierr); 2255 ierr = PetscSectionGetNumFields(localCoarse,&numFields);CHKERRQ(ierr); 2256 maxFields = PetscMax(1,numFields); 2257 ierr = PetscMalloc7(maxFields+1,&offsets,maxFields+1,&offsetsCopy,maxFields+1,&newOffsets,maxFields+1,&newOffsetsCopy,maxFields+1,&rowOffsets,maxFields+1,&numD,maxFields+1,&numO);CHKERRQ(ierr); 2258 ierr = PetscMalloc2(maxFields+1,(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 DMPlexGetIndicesPointFields_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,-1,rowIndices);CHKERRQ(ierr); 2534 for (a = 0; a < aDof; a++) { 2535 PetscInt anchor = anchors[a + aOff], lOff; 2536 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr); 2537 DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,-1,colIndices);CHKERRQ(ierr); 2538 } 2539 } 2540 else { 2541 DMPlexGetIndicesPoint_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,rowIndices);CHKERRQ(ierr); 2542 for (a = 0; a < aDof; a++) { 2543 PetscInt anchor = anchors[a + aOff], lOff; 2544 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr); 2545 DMPlexGetIndicesPoint_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,colIndices);CHKERRQ(ierr); 2546 } 2547 } 2548 if (numFields) { 2549 PetscInt count, a; 2550 2551 for (f = 0, count = 0; f < numFields; f++) { 2552 PetscInt iSize = offsets[f + 1] - offsets[f]; 2553 PetscInt jSize = newOffsets[f + 1] - newOffsets[f]; 2554 ierr = MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);CHKERRQ(ierr); 2555 count += iSize * jSize; 2556 pInd[numColIndices + f] = offsets[f+1]; 2557 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2558 } 2559 for (a = 0; a < aDof; a++) { 2560 PetscInt anchor = anchors[a + aOff]; 2561 PetscInt gOff; 2562 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr); 2563 DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 2564 } 2565 } 2566 else { 2567 PetscInt a; 2568 ierr = MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);CHKERRQ(ierr); 2569 for (a = 0; a < aDof; a++) { 2570 PetscInt anchor = anchors[a + aOff]; 2571 PetscInt gOff; 2572 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr); 2573 DMPlexGetIndicesPoint_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 2574 } 2575 } 2576 ierr = DMRestoreWorkArray(coarse,numColIndices,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 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 2594 } 2595 else { 2596 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 2597 } 2598 } 2599 } 2600 ierr = PetscFree(maxChildIds);CHKERRQ(ierr); 2601 } 2602 { 2603 PetscSF indicesSF, matricesSF; 2604 PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices; 2605 2606 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);CHKERRQ(ierr); 2607 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafMatricesSec);CHKERRQ(ierr); 2608 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootIndicesSec,&remoteOffsetsIndices,leafIndicesSec);CHKERRQ(ierr); 2609 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootMatricesSec,&remoteOffsetsMatrices,leafMatricesSec);CHKERRQ(ierr); 2610 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootIndicesSec,remoteOffsetsIndices,leafIndicesSec,&indicesSF);CHKERRQ(ierr); 2611 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootMatricesSec,remoteOffsetsMatrices,leafMatricesSec,&matricesSF);CHKERRQ(ierr); 2612 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 2613 ierr = PetscFree(remoteOffsetsIndices);CHKERRQ(ierr); 2614 ierr = PetscFree(remoteOffsetsMatrices);CHKERRQ(ierr); 2615 ierr = PetscSectionGetStorageSize(leafIndicesSec,&numLeafIndices);CHKERRQ(ierr); 2616 ierr = PetscSectionGetStorageSize(leafMatricesSec,&numLeafMatrices);CHKERRQ(ierr); 2617 ierr = PetscMalloc2(numLeafIndices,&leafIndices,numLeafMatrices,&leafMatrices);CHKERRQ(ierr); 2618 ierr = PetscSFBcastBegin(indicesSF,MPIU_INT,rootIndices,leafIndices);CHKERRQ(ierr); 2619 ierr = PetscSFBcastBegin(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);CHKERRQ(ierr); 2620 ierr = PetscSFBcastEnd(indicesSF,MPIU_INT,rootIndices,leafIndices);CHKERRQ(ierr); 2621 ierr = PetscSFBcastEnd(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);CHKERRQ(ierr); 2622 ierr = PetscSFDestroy(&matricesSF);CHKERRQ(ierr); 2623 ierr = PetscSFDestroy(&indicesSF);CHKERRQ(ierr); 2624 ierr = PetscFree2(rootIndices,rootMatrices);CHKERRQ(ierr); 2625 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 2626 ierr = PetscSectionDestroy(&rootMatricesSec);CHKERRQ(ierr); 2627 } 2628 /* count to preallocate */ 2629 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 2630 { 2631 PetscInt nGlobal; 2632 PetscInt *dnnz, *onnz; 2633 PetscLayout rowMap, colMap; 2634 PetscInt rowStart, rowEnd, colStart, colEnd; 2635 PetscInt maxDof; 2636 PetscInt *rowIndices; 2637 DM refTree; 2638 PetscInt **refPointFieldN; 2639 PetscScalar ***refPointFieldMats; 2640 PetscSection refConSec, refAnSec; 2641 PetscInt pRefStart,pRefEnd,maxConDof,maxColumns,leafStart,leafEnd; 2642 PetscScalar *pointWork; 2643 2644 ierr = PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal);CHKERRQ(ierr); 2645 ierr = PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz);CHKERRQ(ierr); 2646 ierr = MatGetLayouts(mat,&rowMap,&colMap);CHKERRQ(ierr); 2647 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 2648 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 2649 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 2650 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 2651 ierr = PetscSectionGetMaxDof(globalFine,&maxDof);CHKERRQ(ierr); 2652 ierr = PetscSectionGetChart(leafIndicesSec,&leafStart,&leafEnd);CHKERRQ(ierr); 2653 ierr = DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);CHKERRQ(ierr); 2654 for (p = leafStart; p < leafEnd; p++) { 2655 PetscInt gDof, gcDof, gOff; 2656 PetscInt numColIndices, pIndOff, *pInd; 2657 PetscInt matSize; 2658 PetscInt i; 2659 2660 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 2661 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 2662 if ((gDof - gcDof) <= 0) { 2663 continue; 2664 } 2665 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 2666 if (gOff < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I though having global dofs meant a non-negative offset"); 2667 if ((gOff < rowStart) || ((gOff + gDof - gcDof) > rowEnd)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I thought the row map would constrain the global dofs"); 2668 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2669 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2670 numColIndices -= 2 * numFields; 2671 if (numColIndices <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"global fine dof with no dofs to interpolate from"); 2672 pInd = &leafIndices[pIndOff]; 2673 offsets[0] = 0; 2674 offsetsCopy[0] = 0; 2675 newOffsets[0] = 0; 2676 newOffsetsCopy[0] = 0; 2677 if (numFields) { 2678 PetscInt f; 2679 for (f = 0; f < numFields; f++) { 2680 PetscInt rowDof; 2681 2682 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 2683 offsets[f + 1] = offsets[f] + rowDof; 2684 offsetsCopy[f + 1] = offsets[f + 1]; 2685 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2686 numD[f] = 0; 2687 numO[f] = 0; 2688 } 2689 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 2690 for (f = 0; f < numFields; f++) { 2691 PetscInt colOffset = newOffsets[f]; 2692 PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f]; 2693 2694 for (i = 0; i < numFieldCols; i++) { 2695 PetscInt gInd = pInd[i + colOffset]; 2696 2697 if (gInd >= colStart && gInd < colEnd) { 2698 numD[f]++; 2699 } 2700 else if (gInd >= 0) { /* negative means non-entry */ 2701 numO[f]++; 2702 } 2703 } 2704 } 2705 } 2706 else { 2707 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 2708 numD[0] = 0; 2709 numO[0] = 0; 2710 for (i = 0; i < numColIndices; i++) { 2711 PetscInt gInd = pInd[i]; 2712 2713 if (gInd >= colStart && gInd < colEnd) { 2714 numD[0]++; 2715 } 2716 else if (gInd >= 0) { /* negative means non-entry */ 2717 numO[0]++; 2718 } 2719 } 2720 } 2721 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr); 2722 if (!matSize) { /* incoming matrix is identity */ 2723 PetscInt childId; 2724 2725 childId = childIds[p-pStartF]; 2726 if (childId < 0) { /* no child interpolation: one nnz per */ 2727 if (numFields) { 2728 PetscInt f; 2729 for (f = 0; f < numFields; f++) { 2730 PetscInt numRows = offsets[f+1] - offsets[f], row; 2731 for (row = 0; row < numRows; row++) { 2732 PetscInt gIndCoarse = pInd[newOffsets[f] + row]; 2733 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2734 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2735 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2736 dnnz[gIndFine - rowStart] = 1; 2737 } 2738 else if (gIndCoarse >= 0) { /* remote */ 2739 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2740 onnz[gIndFine - rowStart] = 1; 2741 } 2742 else { /* constrained */ 2743 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2744 } 2745 } 2746 } 2747 } 2748 else { 2749 PetscInt i; 2750 for (i = 0; i < gDof; i++) { 2751 PetscInt gIndCoarse = pInd[i]; 2752 PetscInt gIndFine = rowIndices[i]; 2753 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2754 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2755 dnnz[gIndFine - rowStart] = 1; 2756 } 2757 else if (gIndCoarse >= 0) { /* remote */ 2758 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2759 onnz[gIndFine - rowStart] = 1; 2760 } 2761 else { /* constrained */ 2762 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2763 } 2764 } 2765 } 2766 } 2767 else { /* interpolate from all */ 2768 if (numFields) { 2769 PetscInt f; 2770 for (f = 0; f < numFields; f++) { 2771 PetscInt numRows = offsets[f+1] - offsets[f], row; 2772 for (row = 0; row < numRows; row++) { 2773 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2774 if (gIndFine >= 0) { 2775 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2776 dnnz[gIndFine - rowStart] = numD[f]; 2777 onnz[gIndFine - rowStart] = numO[f]; 2778 } 2779 } 2780 } 2781 } 2782 else { 2783 PetscInt i; 2784 for (i = 0; i < gDof; i++) { 2785 PetscInt gIndFine = rowIndices[i]; 2786 if (gIndFine >= 0) { 2787 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2788 dnnz[gIndFine - rowStart] = numD[0]; 2789 onnz[gIndFine - rowStart] = numO[0]; 2790 } 2791 } 2792 } 2793 } 2794 } 2795 else { /* interpolate from all */ 2796 if (numFields) { 2797 PetscInt f; 2798 for (f = 0; f < numFields; f++) { 2799 PetscInt numRows = offsets[f+1] - offsets[f], row; 2800 for (row = 0; row < numRows; row++) { 2801 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2802 if (gIndFine >= 0) { 2803 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2804 dnnz[gIndFine - rowStart] = numD[f]; 2805 onnz[gIndFine - rowStart] = numO[f]; 2806 } 2807 } 2808 } 2809 } 2810 else { /* every dof get a full row */ 2811 PetscInt i; 2812 for (i = 0; i < gDof; i++) { 2813 PetscInt gIndFine = rowIndices[i]; 2814 if (gIndFine >= 0) { 2815 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2816 dnnz[gIndFine - rowStart] = numD[0]; 2817 onnz[gIndFine - rowStart] = numO[0]; 2818 } 2819 } 2820 } 2821 } 2822 } 2823 ierr = MatXAIJSetPreallocation(mat,1,dnnz,onnz,NULL,NULL);CHKERRQ(ierr); 2824 ierr = PetscFree2(dnnz,onnz);CHKERRQ(ierr); 2825 2826 ierr = DMPlexGetReferenceTree(fine,&refTree);CHKERRQ(ierr); 2827 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 2828 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 2829 ierr = DMPlexGetAnchors(refTree,&refAnSec,NULL);CHKERRQ(ierr); 2830 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 2831 ierr = PetscSectionGetMaxDof(refConSec,&maxConDof);CHKERRQ(ierr); 2832 ierr = PetscSectionGetMaxDof(leafIndicesSec,&maxColumns);CHKERRQ(ierr); 2833 ierr = PetscMalloc1(maxConDof*maxColumns,&pointWork);CHKERRQ(ierr); 2834 for (p = leafStart; p < leafEnd; p++) { 2835 PetscInt gDof, gcDof, gOff; 2836 PetscInt numColIndices, pIndOff, *pInd; 2837 PetscInt matSize; 2838 PetscInt childId; 2839 2840 2841 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 2842 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 2843 if ((gDof - gcDof) <= 0) { 2844 continue; 2845 } 2846 childId = childIds[p-pStartF]; 2847 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 2848 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2849 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2850 numColIndices -= 2 * numFields; 2851 pInd = &leafIndices[pIndOff]; 2852 offsets[0] = 0; 2853 offsetsCopy[0] = 0; 2854 newOffsets[0] = 0; 2855 newOffsetsCopy[0] = 0; 2856 rowOffsets[0] = 0; 2857 if (numFields) { 2858 PetscInt f; 2859 for (f = 0; f < numFields; f++) { 2860 PetscInt rowDof; 2861 2862 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 2863 offsets[f + 1] = offsets[f] + rowDof; 2864 offsetsCopy[f + 1] = offsets[f + 1]; 2865 rowOffsets[f + 1] = pInd[numColIndices + f]; 2866 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2867 } 2868 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 2869 } 2870 else { 2871 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 2872 } 2873 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr); 2874 if (!matSize) { /* incoming matrix is identity */ 2875 if (childId < 0) { /* no child interpolation: scatter */ 2876 if (numFields) { 2877 PetscInt f; 2878 for (f = 0; f < numFields; f++) { 2879 PetscInt numRows = offsets[f+1] - offsets[f], row; 2880 for (row = 0; row < numRows; row++) { 2881 ierr = MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES);CHKERRQ(ierr); 2882 } 2883 } 2884 } 2885 else { 2886 PetscInt numRows = gDof, row; 2887 for (row = 0; row < numRows; row++) { 2888 ierr = MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES);CHKERRQ(ierr); 2889 } 2890 } 2891 } 2892 else { /* interpolate from all */ 2893 if (numFields) { 2894 PetscInt f; 2895 for (f = 0; f < numFields; f++) { 2896 PetscInt numRows = offsets[f+1] - offsets[f]; 2897 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 2898 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES);CHKERRQ(ierr); 2899 } 2900 } 2901 else { 2902 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES);CHKERRQ(ierr); 2903 } 2904 } 2905 } 2906 else { /* interpolate from all */ 2907 PetscInt pMatOff; 2908 PetscScalar *pMat; 2909 2910 ierr = PetscSectionGetOffset(leafMatricesSec,p,&pMatOff);CHKERRQ(ierr); 2911 pMat = &leafMatrices[pMatOff]; 2912 if (childId < 0) { /* copy the incoming matrix */ 2913 if (numFields) { 2914 PetscInt f, count; 2915 for (f = 0, count = 0; f < numFields; f++) { 2916 PetscInt numRows = offsets[f+1]-offsets[f]; 2917 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2918 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2919 PetscScalar *inMat = &pMat[count]; 2920 2921 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES);CHKERRQ(ierr); 2922 count += numCols * numInRows; 2923 } 2924 } 2925 else { 2926 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES);CHKERRQ(ierr); 2927 } 2928 } 2929 else { /* multiply the incoming matrix by the child interpolation */ 2930 if (numFields) { 2931 PetscInt f, count; 2932 for (f = 0, count = 0; f < numFields; f++) { 2933 PetscInt numRows = offsets[f+1]-offsets[f]; 2934 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2935 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2936 PetscScalar *inMat = &pMat[count]; 2937 PetscInt i, j, k; 2938 if (refPointFieldN[childId - pRefStart][f] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2939 for (i = 0; i < numRows; i++) { 2940 for (j = 0; j < numCols; j++) { 2941 PetscScalar val = 0.; 2942 for (k = 0; k < numInRows; k++) { 2943 val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j]; 2944 } 2945 pointWork[i * numCols + j] = val; 2946 } 2947 } 2948 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES);CHKERRQ(ierr); 2949 count += numCols * numInRows; 2950 } 2951 } 2952 else { /* every dof gets a full row */ 2953 PetscInt numRows = gDof; 2954 PetscInt numCols = numColIndices; 2955 PetscInt numInRows = matSize / numColIndices; 2956 PetscInt i, j, k; 2957 if (refPointFieldN[childId - pRefStart][0] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2958 for (i = 0; i < numRows; i++) { 2959 for (j = 0; j < numCols; j++) { 2960 PetscScalar val = 0.; 2961 for (k = 0; k < numInRows; k++) { 2962 val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j]; 2963 } 2964 pointWork[i * numCols + j] = val; 2965 } 2966 } 2967 ierr = MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES);CHKERRQ(ierr); 2968 } 2969 } 2970 } 2971 } 2972 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 2973 ierr = DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);CHKERRQ(ierr); 2974 ierr = PetscFree(pointWork);CHKERRQ(ierr); 2975 } 2976 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2977 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2978 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr); 2979 ierr = PetscSectionDestroy(&leafMatricesSec);CHKERRQ(ierr); 2980 ierr = PetscFree2(leafIndices,leafMatrices);CHKERRQ(ierr); 2981 ierr = PetscFree2(*(PetscInt****)&perms,*(PetscScalar****)&flips);CHKERRQ(ierr); 2982 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr); 2983 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 2984 PetscFunctionReturn(0); 2985 } 2986 2987 /* 2988 * Assuming a nodal basis (w.r.t. the dual basis) basis: 2989 * 2990 * for each coarse dof \phi^c_i: 2991 * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i: 2992 * for each fine dof \phi^f_j; 2993 * a_{i,j} = 0; 2994 * for each fine dof \phi^f_k: 2995 * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l 2996 * [^^^ this is = \phi^c_i ^^^] 2997 */ 2998 PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj) 2999 { 3000 PetscDS ds; 3001 PetscSection section, cSection; 3002 DMLabel canonical, depth; 3003 Mat cMat, mat; 3004 PetscInt *nnz; 3005 PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd; 3006 PetscInt m, n; 3007 PetscScalar *pointScalar; 3008 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent; 3009 PetscErrorCode ierr; 3010 3011 PetscFunctionBegin; 3012 ierr = DMGetDefaultSection(refTree,§ion);CHKERRQ(ierr); 3013 ierr = DMGetDimension(refTree, &dim);CHKERRQ(ierr); 3014 ierr = PetscMalloc6(dim,&v0,dim,&v0parent,dim,&vtmp,dim*dim,&J,dim*dim,&Jparent,dim*dim,&invJ);CHKERRQ(ierr); 3015 ierr = PetscMalloc2(dim,&pointScalar,dim,&pointRef);CHKERRQ(ierr); 3016 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3017 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3018 ierr = PetscSectionGetNumFields(section,&numSecFields);CHKERRQ(ierr); 3019 ierr = DMGetLabel(refTree,"canonical",&canonical);CHKERRQ(ierr); 3020 ierr = DMGetLabel(refTree,"depth",&depth);CHKERRQ(ierr); 3021 ierr = DMGetDefaultConstraints(refTree,&cSection,&cMat);CHKERRQ(ierr); 3022 ierr = DMPlexGetChart(refTree, &pStart, &pEnd);CHKERRQ(ierr); 3023 ierr = DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd);CHKERRQ(ierr); 3024 ierr = MatGetSize(cMat,&n,&m);CHKERRQ(ierr); /* the injector has transpose sizes from the constraint matrix */ 3025 /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */ 3026 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 3027 for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */ 3028 const PetscInt *children; 3029 PetscInt numChildren; 3030 PetscInt i, numChildDof, numSelfDof; 3031 3032 if (canonical) { 3033 PetscInt pCanonical; 3034 ierr = DMLabelGetValue(canonical,p,&pCanonical);CHKERRQ(ierr); 3035 if (p != pCanonical) continue; 3036 } 3037 ierr = DMPlexGetTreeChildren(refTree,p,&numChildren,&children);CHKERRQ(ierr); 3038 if (!numChildren) continue; 3039 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3040 PetscInt child = children[i]; 3041 PetscInt dof; 3042 3043 ierr = PetscSectionGetDof(section,child,&dof);CHKERRQ(ierr); 3044 numChildDof += dof; 3045 } 3046 ierr = PetscSectionGetDof(section,p,&numSelfDof);CHKERRQ(ierr); 3047 if (!numChildDof || !numSelfDof) continue; 3048 for (f = 0; f < numFields; f++) { 3049 PetscInt selfOff; 3050 3051 if (numSecFields) { /* count the dofs for just this field */ 3052 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3053 PetscInt child = children[i]; 3054 PetscInt dof; 3055 3056 ierr = PetscSectionGetFieldDof(section,child,f,&dof);CHKERRQ(ierr); 3057 numChildDof += dof; 3058 } 3059 ierr = PetscSectionGetFieldDof(section,p,f,&numSelfDof);CHKERRQ(ierr); 3060 ierr = PetscSectionGetFieldOffset(section,p,f,&selfOff);CHKERRQ(ierr); 3061 } 3062 else { 3063 ierr = PetscSectionGetOffset(section,p,&selfOff);CHKERRQ(ierr); 3064 } 3065 for (i = 0; i < numSelfDof; i++) { 3066 nnz[selfOff + i] = numChildDof; 3067 } 3068 } 3069 } 3070 ierr = MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat);CHKERRQ(ierr); 3071 ierr = PetscFree(nnz);CHKERRQ(ierr); 3072 /* Setp 2: compute entries */ 3073 for (p = pStart; p < pEnd; p++) { 3074 const PetscInt *children; 3075 PetscInt numChildren; 3076 PetscInt i, numChildDof, numSelfDof; 3077 3078 /* same conditions about when entries occur */ 3079 if (canonical) { 3080 PetscInt pCanonical; 3081 ierr = DMLabelGetValue(canonical,p,&pCanonical);CHKERRQ(ierr); 3082 if (p != pCanonical) continue; 3083 } 3084 ierr = DMPlexGetTreeChildren(refTree,p,&numChildren,&children);CHKERRQ(ierr); 3085 if (!numChildren) continue; 3086 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3087 PetscInt child = children[i]; 3088 PetscInt dof; 3089 3090 ierr = PetscSectionGetDof(section,child,&dof);CHKERRQ(ierr); 3091 numChildDof += dof; 3092 } 3093 ierr = PetscSectionGetDof(section,p,&numSelfDof);CHKERRQ(ierr); 3094 if (!numChildDof || !numSelfDof) continue; 3095 3096 for (f = 0; f < numFields; f++) { 3097 PetscInt selfOff, Nc, parentCell; 3098 PetscInt cellShapeOff; 3099 PetscObject disc; 3100 PetscDualSpace dsp; 3101 PetscClassId classId; 3102 PetscScalar *pointMat; 3103 PetscInt *matRows, *matCols; 3104 PetscInt pO = PETSC_MIN_INT; 3105 const PetscInt *depthNumDof; 3106 3107 if (numSecFields) { 3108 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3109 PetscInt child = children[i]; 3110 PetscInt dof; 3111 3112 ierr = PetscSectionGetFieldDof(section,child,f,&dof);CHKERRQ(ierr); 3113 numChildDof += dof; 3114 } 3115 ierr = PetscSectionGetFieldDof(section,p,f,&numSelfDof);CHKERRQ(ierr); 3116 ierr = PetscSectionGetFieldOffset(section,p,f,&selfOff);CHKERRQ(ierr); 3117 } 3118 else { 3119 ierr = PetscSectionGetOffset(section,p,&selfOff);CHKERRQ(ierr); 3120 } 3121 3122 /* find a cell whose closure contains p */ 3123 if (p >= cStart && p < cEnd) { 3124 parentCell = p; 3125 } 3126 else { 3127 PetscInt *star = NULL; 3128 PetscInt numStar; 3129 3130 parentCell = -1; 3131 ierr = DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3132 for (i = numStar - 1; i >= 0; i--) { 3133 PetscInt c = star[2 * i]; 3134 3135 if (c >= cStart && c < cEnd) { 3136 parentCell = c; 3137 break; 3138 } 3139 } 3140 ierr = DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3141 } 3142 /* determine the offset of p's shape functions withing parentCell's shape functions */ 3143 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 3144 ierr = PetscObjectGetClassId(disc,&classId);CHKERRQ(ierr); 3145 if (classId == PETSCFE_CLASSID) { 3146 ierr = PetscFEGetDualSpace((PetscFE)disc,&dsp);CHKERRQ(ierr); 3147 } 3148 else if (classId == PETSCFV_CLASSID) { 3149 ierr = PetscFVGetDualSpace((PetscFV)disc,&dsp);CHKERRQ(ierr); 3150 } 3151 else { 3152 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object"); 3153 } 3154 ierr = PetscDualSpaceGetNumDof(dsp,&depthNumDof);CHKERRQ(ierr); 3155 ierr = PetscDualSpaceGetNumComponents(dsp,&Nc);CHKERRQ(ierr); 3156 { 3157 PetscInt *closure = NULL; 3158 PetscInt numClosure; 3159 3160 ierr = DMPlexGetTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3161 for (i = 0, cellShapeOff = 0; i < numClosure; i++) { 3162 PetscInt point = closure[2 * i], pointDepth; 3163 3164 pO = closure[2 * i + 1]; 3165 if (point == p) break; 3166 ierr = DMLabelGetValue(depth,point,&pointDepth);CHKERRQ(ierr); 3167 cellShapeOff += depthNumDof[pointDepth]; 3168 } 3169 ierr = DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3170 } 3171 3172 ierr = DMGetWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);CHKERRQ(ierr); 3173 ierr = DMGetWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);CHKERRQ(ierr); 3174 matCols = matRows + numSelfDof; 3175 for (i = 0; i < numSelfDof; i++) { 3176 matRows[i] = selfOff + i; 3177 } 3178 for (i = 0; i < numSelfDof * numChildDof; i++) pointMat[i] = 0.; 3179 { 3180 PetscInt colOff = 0; 3181 3182 for (i = 0; i < numChildren; i++) { 3183 PetscInt child = children[i]; 3184 PetscInt dof, off, j; 3185 3186 if (numSecFields) { 3187 ierr = PetscSectionGetFieldDof(cSection,child,f,&dof);CHKERRQ(ierr); 3188 ierr = PetscSectionGetFieldOffset(cSection,child,f,&off);CHKERRQ(ierr); 3189 } 3190 else { 3191 ierr = PetscSectionGetDof(cSection,child,&dof);CHKERRQ(ierr); 3192 ierr = PetscSectionGetOffset(cSection,child,&off);CHKERRQ(ierr); 3193 } 3194 3195 for (j = 0; j < dof; j++) { 3196 matCols[colOff++] = off + j; 3197 } 3198 } 3199 } 3200 if (classId == PETSCFE_CLASSID) { 3201 PetscFE fe = (PetscFE) disc; 3202 PetscInt fSize; 3203 3204 ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr); 3205 ierr = PetscDualSpaceGetDimension(dsp,&fSize);CHKERRQ(ierr); 3206 for (i = 0; i < numSelfDof; i++) { /* for every shape function */ 3207 PetscQuadrature q; 3208 PetscInt dim, thisNc, numPoints, j, k; 3209 const PetscReal *points; 3210 const PetscReal *weights; 3211 PetscInt *closure = NULL; 3212 PetscInt numClosure; 3213 PetscInt parentCellShapeDof = cellShapeOff + (pO < 0 ? (numSelfDof - 1 - i) : i); 3214 PetscReal *Bparent; 3215 3216 ierr = PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q);CHKERRQ(ierr); 3217 ierr = PetscQuadratureGetData(q,&dim,&thisNc,&numPoints,&points,&weights);CHKERRQ(ierr); 3218 if (thisNc != Nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Functional dim %D does not much basis dim %D\n",thisNc,Nc); 3219 ierr = PetscFEGetTabulation(fe,numPoints,points,&Bparent,NULL,NULL);CHKERRQ(ierr); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */ 3220 for (j = 0; j < numPoints; j++) { 3221 PetscInt childCell = -1; 3222 PetscReal *parentValAtPoint; 3223 const PetscReal xi0[3] = {-1.,-1.,-1.}; 3224 const PetscReal *pointReal = &points[dim * j]; 3225 const PetscScalar *point; 3226 PetscReal *Bchild; 3227 PetscInt childCellShapeOff, pointMatOff; 3228 #if defined(PETSC_USE_COMPLEX) 3229 PetscInt d; 3230 3231 for (d = 0; d < dim; d++) { 3232 pointScalar[d] = points[dim * j + d]; 3233 } 3234 point = pointScalar; 3235 #else 3236 point = pointReal; 3237 #endif 3238 3239 parentValAtPoint = &Bparent[(fSize * j + parentCellShapeDof) * Nc]; 3240 3241 for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/ 3242 PetscInt child = children[k]; 3243 PetscInt *star = NULL; 3244 PetscInt numStar, s; 3245 3246 ierr = DMPlexGetTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3247 for (s = numStar - 1; s >= 0; s--) { 3248 PetscInt c = star[2 * s]; 3249 3250 if (c < cStart || c >= cEnd) continue; 3251 ierr = DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell);CHKERRQ(ierr); 3252 if (childCell >= 0) break; 3253 } 3254 ierr = DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3255 if (childCell >= 0) break; 3256 } 3257 if (childCell < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not locate quadrature point"); 3258 ierr = DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 3259 ierr = DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent);CHKERRQ(ierr); 3260 CoordinatesRefToReal(dim, dim, xi0, v0parent, Jparent, pointReal, vtmp); 3261 CoordinatesRealToRef(dim, dim, xi0, v0, invJ, vtmp, pointRef); 3262 3263 ierr = PetscFEGetTabulation(fe,1,pointRef,&Bchild,NULL,NULL);CHKERRQ(ierr); 3264 ierr = DMPlexGetTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3265 for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */ 3266 PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT; 3267 PetscInt l; 3268 3269 ierr = DMLabelGetValue(depth,child,&childDepth);CHKERRQ(ierr); 3270 childDof = depthNumDof[childDepth]; 3271 for (l = 0, childCellShapeOff = 0; l < numClosure; l++) { 3272 PetscInt point = closure[2 * l]; 3273 PetscInt pointDepth; 3274 3275 childO = closure[2 * l + 1]; 3276 if (point == child) break; 3277 ierr = DMLabelGetValue(depth,point,&pointDepth);CHKERRQ(ierr); 3278 childCellShapeOff += depthNumDof[pointDepth]; 3279 } 3280 if (l == numClosure) { 3281 pointMatOff += childDof; 3282 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */ 3283 } 3284 for (l = 0; l < childDof; l++) { 3285 PetscInt childCellDof = childCellShapeOff + (childO ? (childDof - 1 - l) : l), m; 3286 PetscReal *childValAtPoint; 3287 PetscReal val = 0.; 3288 3289 childValAtPoint = &Bchild[childCellDof * Nc]; 3290 for (m = 0; m < Nc; m++) { 3291 val += weights[j * Nc + m] * parentValAtPoint[m] * childValAtPoint[m]; 3292 } 3293 3294 pointMat[i * numChildDof + pointMatOff + l] += val; 3295 } 3296 pointMatOff += childDof; 3297 } 3298 ierr = DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3299 ierr = PetscFERestoreTabulation(fe,1,pointRef,&Bchild,NULL,NULL);CHKERRQ(ierr); 3300 } 3301 ierr = PetscFERestoreTabulation(fe,numPoints,points,&Bparent,NULL,NULL);CHKERRQ(ierr); 3302 } 3303 } 3304 else { /* just the volume-weighted averages of the children */ 3305 PetscReal parentVol; 3306 PetscInt childCell; 3307 3308 ierr = DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);CHKERRQ(ierr); 3309 for (i = 0, childCell = 0; i < numChildren; i++) { 3310 PetscInt child = children[i], j; 3311 PetscReal childVol; 3312 3313 if (child < cStart || child >= cEnd) continue; 3314 ierr = DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);CHKERRQ(ierr); 3315 for (j = 0; j < Nc; j++) { 3316 pointMat[j * numChildDof + Nc * childCell + j] = childVol / parentVol; 3317 } 3318 childCell++; 3319 } 3320 } 3321 /* Insert pointMat into mat */ 3322 ierr = MatSetValues(mat,numSelfDof,matRows,numChildDof,matCols,pointMat,INSERT_VALUES);CHKERRQ(ierr); 3323 ierr = DMRestoreWorkArray(refTree, numSelfDof + numChildDof, MPIU_INT,&matRows);CHKERRQ(ierr); 3324 ierr = DMRestoreWorkArray(refTree, numSelfDof * numChildDof, MPIU_SCALAR,&pointMat);CHKERRQ(ierr); 3325 } 3326 } 3327 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ);CHKERRQ(ierr); 3328 ierr = PetscFree2(pointScalar,pointRef);CHKERRQ(ierr); 3329 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3330 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3331 *inj = mat; 3332 PetscFunctionReturn(0); 3333 } 3334 3335 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3336 { 3337 PetscDS ds; 3338 PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof; 3339 PetscScalar ***refPointFieldMats; 3340 PetscSection refConSec, refSection; 3341 PetscErrorCode ierr; 3342 3343 PetscFunctionBegin; 3344 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3345 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3346 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 3347 ierr = DMGetDefaultSection(refTree,&refSection);CHKERRQ(ierr); 3348 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3349 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);CHKERRQ(ierr); 3350 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 3351 ierr = PetscMalloc1(maxDof,&rows);CHKERRQ(ierr); 3352 ierr = PetscMalloc1(maxDof*maxDof,&cols);CHKERRQ(ierr); 3353 for (p = pRefStart; p < pRefEnd; p++) { 3354 PetscInt parent, pDof, parentDof; 3355 3356 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 3357 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 3358 ierr = PetscSectionGetDof(refSection,parent,&parentDof);CHKERRQ(ierr); 3359 if (!pDof || !parentDof || parent == p) continue; 3360 3361 ierr = PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);CHKERRQ(ierr); 3362 for (f = 0; f < numFields; f++) { 3363 PetscInt cDof, cOff, numCols, r; 3364 3365 if (numFields > 1) { 3366 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 3367 ierr = PetscSectionGetFieldOffset(refConSec,p,f,&cOff);CHKERRQ(ierr); 3368 } 3369 else { 3370 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 3371 ierr = PetscSectionGetOffset(refConSec,p,&cOff);CHKERRQ(ierr); 3372 } 3373 3374 for (r = 0; r < cDof; r++) { 3375 rows[r] = cOff + r; 3376 } 3377 numCols = 0; 3378 { 3379 PetscInt aDof, aOff, j; 3380 3381 if (numFields > 1) { 3382 ierr = PetscSectionGetFieldDof(refSection,parent,f,&aDof);CHKERRQ(ierr); 3383 ierr = PetscSectionGetFieldOffset(refSection,parent,f,&aOff);CHKERRQ(ierr); 3384 } 3385 else { 3386 ierr = PetscSectionGetDof(refSection,parent,&aDof);CHKERRQ(ierr); 3387 ierr = PetscSectionGetOffset(refSection,parent,&aOff);CHKERRQ(ierr); 3388 } 3389 3390 for (j = 0; j < aDof; j++) { 3391 cols[numCols++] = aOff + j; 3392 } 3393 } 3394 ierr = PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 3395 /* transpose of constraint matrix */ 3396 ierr = MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 3397 } 3398 } 3399 *childrenMats = refPointFieldMats; 3400 ierr = PetscFree(rows);CHKERRQ(ierr); 3401 ierr = PetscFree(cols);CHKERRQ(ierr); 3402 PetscFunctionReturn(0); 3403 } 3404 3405 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3406 { 3407 PetscDS ds; 3408 PetscScalar ***refPointFieldMats; 3409 PetscInt numFields, pRefStart, pRefEnd, p, f; 3410 PetscSection refConSec, refSection; 3411 PetscErrorCode ierr; 3412 3413 PetscFunctionBegin; 3414 refPointFieldMats = *childrenMats; 3415 *childrenMats = NULL; 3416 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3417 ierr = DMGetDefaultSection(refTree,&refSection);CHKERRQ(ierr); 3418 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3419 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 3420 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3421 for (p = pRefStart; p < pRefEnd; p++) { 3422 PetscInt parent, pDof, parentDof; 3423 3424 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 3425 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 3426 ierr = PetscSectionGetDof(refSection,parent,&parentDof);CHKERRQ(ierr); 3427 if (!pDof || !parentDof || parent == p) continue; 3428 3429 for (f = 0; f < numFields; f++) { 3430 PetscInt cDof; 3431 3432 if (numFields > 1) { 3433 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 3434 } 3435 else { 3436 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 3437 } 3438 3439 ierr = PetscFree(refPointFieldMats[p - pRefStart][f]);CHKERRQ(ierr); 3440 } 3441 ierr = PetscFree(refPointFieldMats[p - pRefStart]);CHKERRQ(ierr); 3442 } 3443 ierr = PetscFree(refPointFieldMats);CHKERRQ(ierr); 3444 PetscFunctionReturn(0); 3445 } 3446 3447 static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree,Mat *injRef) 3448 { 3449 Mat cMatRef; 3450 PetscObject injRefObj; 3451 PetscErrorCode ierr; 3452 3453 PetscFunctionBegin; 3454 ierr = DMGetDefaultConstraints(refTree,NULL,&cMatRef);CHKERRQ(ierr); 3455 ierr = PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj);CHKERRQ(ierr); 3456 *injRef = (Mat) injRefObj; 3457 if (!*injRef) { 3458 ierr = DMPlexComputeInjectorReferenceTree(refTree,injRef);CHKERRQ(ierr); 3459 ierr = PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)*injRef);CHKERRQ(ierr); 3460 /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */ 3461 ierr = PetscObjectDereference((PetscObject)*injRef);CHKERRQ(ierr); 3462 } 3463 PetscFunctionReturn(0); 3464 } 3465 3466 static PetscErrorCode DMPlexTransferInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, const PetscInt *childIds, Vec fineVec, PetscInt numFields, PetscInt *offsets, PetscSection *rootMultiSec, PetscSection *multiLeafSec, PetscInt **gatheredIndices, PetscScalar **gatheredValues) 3467 { 3468 PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti; 3469 PetscSection globalCoarse, globalFine; 3470 PetscSection localCoarse, localFine, leafIndicesSec; 3471 PetscSection multiRootSec, rootIndicesSec; 3472 PetscInt *leafInds, *rootInds = NULL; 3473 const PetscInt *rootDegrees; 3474 PetscScalar *leafVals = NULL, *rootVals = NULL; 3475 PetscSF coarseToFineEmbedded; 3476 PetscErrorCode ierr; 3477 3478 PetscFunctionBegin; 3479 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3480 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 3481 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 3482 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 3483 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);CHKERRQ(ierr); 3484 ierr = PetscSectionSetChart(leafIndicesSec,pStartF, pEndF);CHKERRQ(ierr); 3485 ierr = PetscSectionGetMaxDof(localFine,&maxDof);CHKERRQ(ierr); 3486 { /* winnow fine points that don't have global dofs out of the sf */ 3487 PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices; 3488 const PetscInt *leaves; 3489 3490 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 3491 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3492 p = leaves ? leaves[l] : l; 3493 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 3494 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 3495 if ((dof - cdof) > 0) { 3496 numPointsWithDofs++; 3497 3498 ierr = PetscSectionGetDof(localFine,p,&dof);CHKERRQ(ierr); 3499 ierr = PetscSectionSetDof(leafIndicesSec,p,dof + 1);CHKERRQ(ierr); 3500 } 3501 } 3502 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 3503 ierr = PetscSectionSetUp(leafIndicesSec);CHKERRQ(ierr); 3504 ierr = PetscSectionGetStorageSize(leafIndicesSec,&numIndices);CHKERRQ(ierr); 3505 ierr = PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1),&leafInds);CHKERRQ(ierr); 3506 if (gatheredValues) {ierr = PetscMalloc1(numIndices,&leafVals);CHKERRQ(ierr);} 3507 for (l = 0, offset = 0; l < nleaves; l++) { 3508 p = leaves ? leaves[l] : l; 3509 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 3510 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 3511 if ((dof - cdof) > 0) { 3512 PetscInt off, gOff; 3513 PetscInt *pInd; 3514 PetscScalar *pVal = NULL; 3515 3516 pointsWithDofs[offset++] = l; 3517 3518 ierr = PetscSectionGetOffset(leafIndicesSec,p,&off);CHKERRQ(ierr); 3519 3520 pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds; 3521 if (gatheredValues) { 3522 PetscInt i; 3523 3524 pVal = &leafVals[off + 1]; 3525 for (i = 0; i < dof; i++) pVal[i] = 0.; 3526 } 3527 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 3528 3529 offsets[0] = 0; 3530 if (numFields) { 3531 PetscInt f; 3532 3533 for (f = 0; f < numFields; f++) { 3534 PetscInt fDof; 3535 ierr = PetscSectionGetFieldDof(localFine,p,f,&fDof);CHKERRQ(ierr); 3536 offsets[f + 1] = fDof + offsets[f]; 3537 } 3538 DMPlexGetIndicesPointFields_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 3539 } 3540 else { 3541 DMPlexGetIndicesPoint_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 3542 } 3543 if (gatheredValues) {ierr = VecGetValues(fineVec,dof,pInd,pVal);CHKERRQ(ierr);} 3544 } 3545 } 3546 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 3547 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 3548 } 3549 3550 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3551 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 3552 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 3553 3554 { /* there may be the case where an sf root has a parent: broadcast parents back to children */ 3555 MPI_Datatype threeInt; 3556 PetscMPIInt rank; 3557 PetscInt (*parentNodeAndIdCoarse)[3]; 3558 PetscInt (*parentNodeAndIdFine)[3]; 3559 PetscInt p, nleaves, nleavesToParents; 3560 PetscSF pointSF, sfToParents; 3561 const PetscInt *ilocal; 3562 const PetscSFNode *iremote; 3563 PetscSFNode *iremoteToParents; 3564 PetscInt *ilocalToParents; 3565 3566 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);CHKERRQ(ierr); 3567 ierr = MPI_Type_contiguous(3,MPIU_INT,&threeInt);CHKERRQ(ierr); 3568 ierr = MPI_Type_commit(&threeInt);CHKERRQ(ierr); 3569 ierr = PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine);CHKERRQ(ierr); 3570 ierr = DMGetPointSF(coarse,&pointSF);CHKERRQ(ierr); 3571 ierr = PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote);CHKERRQ(ierr); 3572 for (p = pStartC; p < pEndC; p++) { 3573 PetscInt parent, childId; 3574 ierr = DMPlexGetTreeParent(coarse,p,&parent,&childId);CHKERRQ(ierr); 3575 parentNodeAndIdCoarse[p - pStartC][0] = rank; 3576 parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC; 3577 parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId; 3578 if (nleaves > 0) { 3579 PetscInt leaf = -1; 3580 3581 if (ilocal) { 3582 ierr = PetscFindInt(parent,nleaves,ilocal,&leaf);CHKERRQ(ierr); 3583 } 3584 else { 3585 leaf = p - pStartC; 3586 } 3587 if (leaf >= 0) { 3588 parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank; 3589 parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index; 3590 } 3591 } 3592 } 3593 for (p = pStartF; p < pEndF; p++) { 3594 parentNodeAndIdFine[p - pStartF][0] = -1; 3595 parentNodeAndIdFine[p - pStartF][1] = -1; 3596 parentNodeAndIdFine[p - pStartF][2] = -1; 3597 } 3598 ierr = PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3599 ierr = PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3600 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3601 PetscInt dof; 3602 3603 ierr = PetscSectionGetDof(leafIndicesSec,p,&dof);CHKERRQ(ierr); 3604 if (dof) { 3605 PetscInt off; 3606 3607 ierr = PetscSectionGetOffset(leafIndicesSec,p,&off);CHKERRQ(ierr); 3608 if (gatheredIndices) { 3609 leafInds[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3610 } else if (gatheredValues) { 3611 leafVals[off] = (PetscScalar) PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3612 } 3613 } 3614 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3615 nleavesToParents++; 3616 } 3617 } 3618 ierr = PetscMalloc1(nleavesToParents,&ilocalToParents);CHKERRQ(ierr); 3619 ierr = PetscMalloc1(nleavesToParents,&iremoteToParents);CHKERRQ(ierr); 3620 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3621 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3622 ilocalToParents[nleavesToParents] = p - pStartF; 3623 iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p-pStartF][0]; 3624 iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1]; 3625 nleavesToParents++; 3626 } 3627 } 3628 ierr = PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents);CHKERRQ(ierr); 3629 ierr = PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER);CHKERRQ(ierr); 3630 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3631 3632 coarseToFineEmbedded = sfToParents; 3633 3634 ierr = PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3635 ierr = MPI_Type_free(&threeInt);CHKERRQ(ierr); 3636 } 3637 3638 { /* winnow out coarse points that don't have dofs */ 3639 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3640 PetscSF sfDofsOnly; 3641 3642 for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) { 3643 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3644 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3645 if ((dof - cdof) > 0) { 3646 numPointsWithDofs++; 3647 } 3648 } 3649 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 3650 for (p = pStartC, offset = 0; p < pEndC; p++) { 3651 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3652 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3653 if ((dof - cdof) > 0) { 3654 pointsWithDofs[offset++] = p - pStartC; 3655 } 3656 } 3657 ierr = PetscSFCreateEmbeddedSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);CHKERRQ(ierr); 3658 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3659 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 3660 coarseToFineEmbedded = sfDofsOnly; 3661 } 3662 3663 /* communicate back to the coarse mesh which coarse points have children (that may require injection) */ 3664 ierr = PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees);CHKERRQ(ierr); 3665 ierr = PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees);CHKERRQ(ierr); 3666 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec);CHKERRQ(ierr); 3667 ierr = PetscSectionSetChart(multiRootSec,pStartC,pEndC);CHKERRQ(ierr); 3668 for (p = pStartC; p < pEndC; p++) { 3669 ierr = PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC]);CHKERRQ(ierr); 3670 } 3671 ierr = PetscSectionSetUp(multiRootSec);CHKERRQ(ierr); 3672 ierr = PetscSectionGetStorageSize(multiRootSec,&numMulti);CHKERRQ(ierr); 3673 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);CHKERRQ(ierr); 3674 { /* distribute the leaf section */ 3675 PetscSF multi, multiInv, indicesSF; 3676 PetscInt *remoteOffsets, numRootIndices; 3677 3678 ierr = PetscSFGetMultiSF(coarseToFineEmbedded,&multi);CHKERRQ(ierr); 3679 ierr = PetscSFCreateInverseSF(multi,&multiInv);CHKERRQ(ierr); 3680 ierr = PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec);CHKERRQ(ierr); 3681 ierr = PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF);CHKERRQ(ierr); 3682 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 3683 ierr = PetscSFDestroy(&multiInv);CHKERRQ(ierr); 3684 ierr = PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);CHKERRQ(ierr); 3685 if (gatheredIndices) { 3686 ierr = PetscMalloc1(numRootIndices,&rootInds);CHKERRQ(ierr); 3687 ierr = PetscSFBcastBegin(indicesSF,MPIU_INT,leafInds,rootInds);CHKERRQ(ierr); 3688 ierr = PetscSFBcastEnd(indicesSF,MPIU_INT,leafInds,rootInds);CHKERRQ(ierr); 3689 } 3690 if (gatheredValues) { 3691 ierr = PetscMalloc1(numRootIndices,&rootVals);CHKERRQ(ierr); 3692 ierr = PetscSFBcastBegin(indicesSF,MPIU_SCALAR,leafVals,rootVals);CHKERRQ(ierr); 3693 ierr = PetscSFBcastEnd(indicesSF,MPIU_SCALAR,leafVals,rootVals);CHKERRQ(ierr); 3694 } 3695 ierr = PetscSFDestroy(&indicesSF);CHKERRQ(ierr); 3696 } 3697 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr); 3698 ierr = PetscFree(leafInds);CHKERRQ(ierr); 3699 ierr = PetscFree(leafVals);CHKERRQ(ierr); 3700 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3701 *rootMultiSec = multiRootSec; 3702 *multiLeafSec = rootIndicesSec; 3703 if (gatheredIndices) *gatheredIndices = rootInds; 3704 if (gatheredValues) *gatheredValues = rootVals; 3705 PetscFunctionReturn(0); 3706 } 3707 3708 PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 3709 { 3710 DM refTree; 3711 PetscSection multiRootSec, rootIndicesSec; 3712 PetscSection globalCoarse, globalFine; 3713 PetscSection localCoarse, localFine; 3714 PetscSection cSecRef; 3715 PetscInt *rootIndices, *parentIndices, pRefStart, pRefEnd; 3716 Mat injRef; 3717 PetscInt numFields, maxDof; 3718 PetscInt pStartC, pEndC, pStartF, pEndF, p; 3719 PetscInt *offsets, *offsetsCopy, *rowOffsets; 3720 PetscLayout rowMap, colMap; 3721 PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO; 3722 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 3723 PetscErrorCode ierr; 3724 3725 PetscFunctionBegin; 3726 3727 /* get the templates for the fine-to-coarse injection from the reference tree */ 3728 ierr = DMPlexGetReferenceTree(coarse,&refTree);CHKERRQ(ierr); 3729 ierr = DMGetDefaultConstraints(refTree,&cSecRef,NULL);CHKERRQ(ierr); 3730 ierr = PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3731 ierr = DMPlexReferenceTreeGetInjector(refTree,&injRef);CHKERRQ(ierr); 3732 3733 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 3734 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 3735 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 3736 ierr = PetscSectionGetNumFields(localFine,&numFields);CHKERRQ(ierr); 3737 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3738 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 3739 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 3740 ierr = PetscSectionGetMaxDof(localCoarse,&maxDof);CHKERRQ(ierr); 3741 { 3742 PetscInt maxFields = PetscMax(1,numFields) + 1; 3743 ierr = PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);CHKERRQ(ierr); 3744 } 3745 3746 ierr = DMPlexTransferInjectorTree(coarse,fine,coarseToFine,childIds,NULL,numFields,offsets,&multiRootSec,&rootIndicesSec,&rootIndices,NULL);CHKERRQ(ierr); 3747 3748 ierr = PetscMalloc1(maxDof,&parentIndices);CHKERRQ(ierr); 3749 3750 /* count indices */ 3751 ierr = MatGetLayouts(mat,&rowMap,&colMap);CHKERRQ(ierr); 3752 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 3753 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 3754 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 3755 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 3756 ierr = PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO);CHKERRQ(ierr); 3757 for (p = pStartC; p < pEndC; p++) { 3758 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3759 3760 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3761 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3762 if ((dof - cdof) <= 0) continue; 3763 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 3764 3765 rowOffsets[0] = 0; 3766 offsetsCopy[0] = 0; 3767 if (numFields) { 3768 PetscInt f; 3769 3770 for (f = 0; f < numFields; f++) { 3771 PetscInt fDof; 3772 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 3773 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3774 } 3775 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 3776 } 3777 else { 3778 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 3779 rowOffsets[1] = offsetsCopy[0]; 3780 } 3781 3782 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 3783 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 3784 leafEnd = leafStart + numLeaves; 3785 for (l = leafStart; l < leafEnd; l++) { 3786 PetscInt numIndices, childId, offset; 3787 const PetscInt *childIndices; 3788 3789 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 3790 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 3791 childId = rootIndices[offset++]; 3792 childIndices = &rootIndices[offset]; 3793 numIndices--; 3794 3795 if (childId == -1) { /* equivalent points: scatter */ 3796 PetscInt i; 3797 3798 for (i = 0; i < numIndices; i++) { 3799 PetscInt colIndex = childIndices[i]; 3800 PetscInt rowIndex = parentIndices[i]; 3801 if (rowIndex < 0) continue; 3802 if (colIndex < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unconstrained fine and constrained coarse"); 3803 if (colIndex >= colStart && colIndex < colEnd) { 3804 nnzD[rowIndex - rowStart] = 1; 3805 } 3806 else { 3807 nnzO[rowIndex - rowStart] = 1; 3808 } 3809 } 3810 } 3811 else { 3812 PetscInt parentId, f, lim; 3813 3814 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 3815 3816 lim = PetscMax(1,numFields); 3817 offsets[0] = 0; 3818 if (numFields) { 3819 PetscInt f; 3820 3821 for (f = 0; f < numFields; f++) { 3822 PetscInt fDof; 3823 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 3824 3825 offsets[f + 1] = fDof + offsets[f]; 3826 } 3827 } 3828 else { 3829 PetscInt cDof; 3830 3831 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 3832 offsets[1] = cDof; 3833 } 3834 for (f = 0; f < lim; f++) { 3835 PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1]; 3836 PetscInt childStart = offsets[f], childEnd = offsets[f + 1]; 3837 PetscInt i, numD = 0, numO = 0; 3838 3839 for (i = childStart; i < childEnd; i++) { 3840 PetscInt colIndex = childIndices[i]; 3841 3842 if (colIndex < 0) continue; 3843 if (colIndex >= colStart && colIndex < colEnd) { 3844 numD++; 3845 } 3846 else { 3847 numO++; 3848 } 3849 } 3850 for (i = parentStart; i < parentEnd; i++) { 3851 PetscInt rowIndex = parentIndices[i]; 3852 3853 if (rowIndex < 0) continue; 3854 nnzD[rowIndex - rowStart] += numD; 3855 nnzO[rowIndex - rowStart] += numO; 3856 } 3857 } 3858 } 3859 } 3860 } 3861 /* preallocate */ 3862 ierr = MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL);CHKERRQ(ierr); 3863 ierr = PetscFree2(nnzD,nnzO);CHKERRQ(ierr); 3864 /* insert values */ 3865 ierr = DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 3866 for (p = pStartC; p < pEndC; p++) { 3867 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3868 3869 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3870 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3871 if ((dof - cdof) <= 0) continue; 3872 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 3873 3874 rowOffsets[0] = 0; 3875 offsetsCopy[0] = 0; 3876 if (numFields) { 3877 PetscInt f; 3878 3879 for (f = 0; f < numFields; f++) { 3880 PetscInt fDof; 3881 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 3882 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3883 } 3884 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 3885 } 3886 else { 3887 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 3888 rowOffsets[1] = offsetsCopy[0]; 3889 } 3890 3891 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 3892 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 3893 leafEnd = leafStart + numLeaves; 3894 for (l = leafStart; l < leafEnd; l++) { 3895 PetscInt numIndices, childId, offset; 3896 const PetscInt *childIndices; 3897 3898 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 3899 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 3900 childId = rootIndices[offset++]; 3901 childIndices = &rootIndices[offset]; 3902 numIndices--; 3903 3904 if (childId == -1) { /* equivalent points: scatter */ 3905 PetscInt i; 3906 3907 for (i = 0; i < numIndices; i++) { 3908 ierr = MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES);CHKERRQ(ierr); 3909 } 3910 } 3911 else { 3912 PetscInt parentId, f, lim; 3913 3914 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 3915 3916 lim = PetscMax(1,numFields); 3917 offsets[0] = 0; 3918 if (numFields) { 3919 PetscInt f; 3920 3921 for (f = 0; f < numFields; f++) { 3922 PetscInt fDof; 3923 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 3924 3925 offsets[f + 1] = fDof + offsets[f]; 3926 } 3927 } 3928 else { 3929 PetscInt cDof; 3930 3931 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 3932 offsets[1] = cDof; 3933 } 3934 for (f = 0; f < lim; f++) { 3935 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 3936 PetscInt *rowIndices = &parentIndices[rowOffsets[f]]; 3937 const PetscInt *colIndices = &childIndices[offsets[f]]; 3938 3939 ierr = MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES);CHKERRQ(ierr); 3940 } 3941 } 3942 } 3943 } 3944 ierr = PetscSectionDestroy(&multiRootSec);CHKERRQ(ierr); 3945 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 3946 ierr = PetscFree(parentIndices);CHKERRQ(ierr); 3947 ierr = DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 3948 ierr = PetscFree(rootIndices);CHKERRQ(ierr); 3949 ierr = PetscFree3(offsets,offsetsCopy,rowOffsets);CHKERRQ(ierr); 3950 3951 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3952 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3953 PetscFunctionReturn(0); 3954 } 3955 3956 static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom) 3957 { 3958 PetscDS ds; 3959 PetscSF coarseToFineEmbedded; 3960 PetscSection globalCoarse, globalFine; 3961 PetscSection localCoarse, localFine; 3962 PetscSection aSec, cSec; 3963 PetscSection rootValuesSec; 3964 PetscSection leafValuesSec; 3965 PetscScalar *rootValues, *leafValues; 3966 IS aIS; 3967 const PetscInt *anchors; 3968 Mat cMat; 3969 PetscInt numFields; 3970 PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd, cellEndInterior; 3971 PetscInt aStart, aEnd, cStart, cEnd; 3972 PetscInt *maxChildIds; 3973 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 3974 PetscFV fv = NULL; 3975 PetscInt dim, numFVcomps = -1, fvField = -1; 3976 DM cellDM = NULL, gradDM = NULL; 3977 const PetscScalar *cellGeomArray = NULL; 3978 const PetscScalar *gradArray = NULL; 3979 PetscErrorCode ierr; 3980 3981 PetscFunctionBegin; 3982 ierr = VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 3983 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3984 ierr = DMPlexGetHeightStratum(coarse,0,&cellStart,&cellEnd);CHKERRQ(ierr); 3985 ierr = DMPlexGetHybridBounds(coarse,&cellEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 3986 cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior; 3987 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 3988 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 3989 ierr = DMGetCoordinateDim(coarse,&dim);CHKERRQ(ierr); 3990 { /* winnow fine points that don't have global dofs out of the sf */ 3991 PetscInt nleaves, l; 3992 const PetscInt *leaves; 3993 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3994 3995 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 3996 3997 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3998 PetscInt p = leaves ? leaves[l] : l; 3999 4000 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 4001 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 4002 if ((dof - cdof) > 0) { 4003 numPointsWithDofs++; 4004 } 4005 } 4006 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 4007 for (l = 0, offset = 0; l < nleaves; l++) { 4008 PetscInt p = leaves ? leaves[l] : l; 4009 4010 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 4011 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 4012 if ((dof - cdof) > 0) { 4013 pointsWithDofs[offset++] = l; 4014 } 4015 } 4016 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 4017 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 4018 } 4019 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 4020 ierr = PetscMalloc1(pEndC-pStartC,&maxChildIds);CHKERRQ(ierr); 4021 for (p = pStartC; p < pEndC; p++) { 4022 maxChildIds[p - pStartC] = -2; 4023 } 4024 ierr = PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 4025 ierr = PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 4026 4027 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 4028 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 4029 4030 ierr = DMPlexGetAnchors(coarse,&aSec,&aIS);CHKERRQ(ierr); 4031 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 4032 ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr); 4033 4034 ierr = DMGetDefaultConstraints(coarse,&cSec,&cMat);CHKERRQ(ierr); 4035 ierr = PetscSectionGetChart(cSec,&cStart,&cEnd);CHKERRQ(ierr); 4036 4037 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 4038 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootValuesSec);CHKERRQ(ierr); 4039 ierr = PetscSectionSetChart(rootValuesSec,pStartC,pEndC);CHKERRQ(ierr); 4040 ierr = PetscSectionGetNumFields(localCoarse,&numFields);CHKERRQ(ierr); 4041 { 4042 PetscInt maxFields = PetscMax(1,numFields) + 1; 4043 ierr = PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);CHKERRQ(ierr); 4044 } 4045 if (grad) { 4046 PetscInt i; 4047 4048 ierr = VecGetDM(cellGeom,&cellDM);CHKERRQ(ierr); 4049 ierr = VecGetArrayRead(cellGeom,&cellGeomArray);CHKERRQ(ierr); 4050 ierr = VecGetDM(grad,&gradDM);CHKERRQ(ierr); 4051 ierr = VecGetArrayRead(grad,&gradArray);CHKERRQ(ierr); 4052 for (i = 0; i < PetscMax(1,numFields); i++) { 4053 PetscObject obj; 4054 PetscClassId id; 4055 4056 ierr = DMGetField(coarse, i, &obj);CHKERRQ(ierr); 4057 ierr = PetscObjectGetClassId(obj,&id);CHKERRQ(ierr); 4058 if (id == PETSCFV_CLASSID) { 4059 fv = (PetscFV) obj; 4060 ierr = PetscFVGetNumComponents(fv,&numFVcomps);CHKERRQ(ierr); 4061 fvField = i; 4062 break; 4063 } 4064 } 4065 } 4066 4067 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 4068 PetscInt dof; 4069 PetscInt maxChildId = maxChildIds[p - pStartC]; 4070 PetscInt numValues = 0; 4071 4072 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 4073 if (dof < 0) { 4074 dof = -(dof + 1); 4075 } 4076 offsets[0] = 0; 4077 newOffsets[0] = 0; 4078 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 4079 PetscInt *closure = NULL, closureSize, cl; 4080 4081 ierr = DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 4082 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 4083 PetscInt c = closure[2 * cl], clDof; 4084 4085 ierr = PetscSectionGetDof(localCoarse,c,&clDof);CHKERRQ(ierr); 4086 numValues += clDof; 4087 } 4088 ierr = DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 4089 } 4090 else if (maxChildId == -1) { 4091 ierr = PetscSectionGetDof(localCoarse,p,&numValues);CHKERRQ(ierr); 4092 } 4093 /* we will pack the column indices with the field offsets */ 4094 if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) { 4095 /* also send the centroid, and the gradient */ 4096 numValues += dim * (1 + numFVcomps); 4097 } 4098 ierr = PetscSectionSetDof(rootValuesSec,p,numValues);CHKERRQ(ierr); 4099 } 4100 ierr = PetscSectionSetUp(rootValuesSec);CHKERRQ(ierr); 4101 { 4102 PetscInt numRootValues; 4103 const PetscScalar *coarseArray; 4104 4105 ierr = PetscSectionGetStorageSize(rootValuesSec,&numRootValues);CHKERRQ(ierr); 4106 ierr = PetscMalloc1(numRootValues,&rootValues);CHKERRQ(ierr); 4107 ierr = VecGetArrayRead(vecCoarseLocal,&coarseArray);CHKERRQ(ierr); 4108 for (p = pStartC; p < pEndC; p++) { 4109 PetscInt numValues; 4110 PetscInt pValOff; 4111 PetscScalar *pVal; 4112 PetscInt maxChildId = maxChildIds[p - pStartC]; 4113 4114 ierr = PetscSectionGetDof(rootValuesSec,p,&numValues);CHKERRQ(ierr); 4115 if (!numValues) { 4116 continue; 4117 } 4118 ierr = PetscSectionGetOffset(rootValuesSec,p,&pValOff);CHKERRQ(ierr); 4119 pVal = &(rootValues[pValOff]); 4120 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 4121 PetscInt closureSize = numValues; 4122 ierr = DMPlexVecGetClosure(coarse,NULL,vecCoarseLocal,p,&closureSize,&pVal);CHKERRQ(ierr); 4123 if (grad && p >= cellStart && p < cellEnd) { 4124 PetscFVCellGeom *cg; 4125 PetscScalar *gradVals = NULL; 4126 PetscInt i; 4127 4128 pVal += (numValues - dim * (1 + numFVcomps)); 4129 4130 ierr = DMPlexPointLocalRead(cellDM,p,cellGeomArray,(void *) &cg);CHKERRQ(ierr); 4131 for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i]; 4132 pVal += dim; 4133 ierr = DMPlexPointGlobalRead(gradDM,p,gradArray,(void *) &gradVals);CHKERRQ(ierr); 4134 for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i]; 4135 } 4136 } 4137 else if (maxChildId == -1) { 4138 PetscInt lDof, lOff, i; 4139 4140 ierr = PetscSectionGetDof(localCoarse,p,&lDof);CHKERRQ(ierr); 4141 ierr = PetscSectionGetOffset(localCoarse,p,&lOff);CHKERRQ(ierr); 4142 for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i]; 4143 } 4144 } 4145 ierr = VecRestoreArrayRead(vecCoarseLocal,&coarseArray);CHKERRQ(ierr); 4146 ierr = PetscFree(maxChildIds);CHKERRQ(ierr); 4147 } 4148 { 4149 PetscSF valuesSF; 4150 PetscInt *remoteOffsetsValues, numLeafValues; 4151 4152 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafValuesSec);CHKERRQ(ierr); 4153 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootValuesSec,&remoteOffsetsValues,leafValuesSec);CHKERRQ(ierr); 4154 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootValuesSec,remoteOffsetsValues,leafValuesSec,&valuesSF);CHKERRQ(ierr); 4155 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 4156 ierr = PetscFree(remoteOffsetsValues);CHKERRQ(ierr); 4157 ierr = PetscSectionGetStorageSize(leafValuesSec,&numLeafValues);CHKERRQ(ierr); 4158 ierr = PetscMalloc1(numLeafValues,&leafValues);CHKERRQ(ierr); 4159 ierr = PetscSFBcastBegin(valuesSF,MPIU_SCALAR,rootValues,leafValues);CHKERRQ(ierr); 4160 ierr = PetscSFBcastEnd(valuesSF,MPIU_SCALAR,rootValues,leafValues);CHKERRQ(ierr); 4161 ierr = PetscSFDestroy(&valuesSF);CHKERRQ(ierr); 4162 ierr = PetscFree(rootValues);CHKERRQ(ierr); 4163 ierr = PetscSectionDestroy(&rootValuesSec);CHKERRQ(ierr); 4164 } 4165 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 4166 { 4167 PetscInt maxDof; 4168 PetscInt *rowIndices; 4169 DM refTree; 4170 PetscInt **refPointFieldN; 4171 PetscScalar ***refPointFieldMats; 4172 PetscSection refConSec, refAnSec; 4173 PetscInt pRefStart,pRefEnd,leafStart,leafEnd; 4174 PetscScalar *pointWork; 4175 4176 ierr = PetscSectionGetMaxDof(globalFine,&maxDof);CHKERRQ(ierr); 4177 ierr = DMGetWorkArray(fine,maxDof,MPIU_INT,&rowIndices);CHKERRQ(ierr); 4178 ierr = DMGetWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);CHKERRQ(ierr); 4179 ierr = DMPlexGetReferenceTree(fine,&refTree);CHKERRQ(ierr); 4180 ierr = DMGetDS(fine,&ds);CHKERRQ(ierr); 4181 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 4182 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 4183 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 4184 ierr = DMPlexGetAnchors(refTree,&refAnSec,NULL);CHKERRQ(ierr); 4185 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 4186 ierr = PetscSectionGetChart(leafValuesSec,&leafStart,&leafEnd);CHKERRQ(ierr); 4187 ierr = DMPlexGetHeightStratum(fine,0,&cellStart,&cellEnd);CHKERRQ(ierr); 4188 ierr = DMPlexGetHybridBounds(fine,&cellEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 4189 cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior; 4190 for (p = leafStart; p < leafEnd; p++) { 4191 PetscInt gDof, gcDof, gOff, lDof; 4192 PetscInt numValues, pValOff; 4193 PetscInt childId; 4194 const PetscScalar *pVal; 4195 const PetscScalar *fvGradData = NULL; 4196 4197 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 4198 ierr = PetscSectionGetDof(localFine,p,&lDof);CHKERRQ(ierr); 4199 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 4200 if ((gDof - gcDof) <= 0) { 4201 continue; 4202 } 4203 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 4204 ierr = PetscSectionGetDof(leafValuesSec,p,&numValues);CHKERRQ(ierr); 4205 if (!numValues) continue; 4206 ierr = PetscSectionGetOffset(leafValuesSec,p,&pValOff);CHKERRQ(ierr); 4207 pVal = &leafValues[pValOff]; 4208 offsets[0] = 0; 4209 offsetsCopy[0] = 0; 4210 newOffsets[0] = 0; 4211 newOffsetsCopy[0] = 0; 4212 childId = cids[p - pStartF]; 4213 if (numFields) { 4214 PetscInt f; 4215 for (f = 0; f < numFields; f++) { 4216 PetscInt rowDof; 4217 4218 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 4219 offsets[f + 1] = offsets[f] + rowDof; 4220 offsetsCopy[f + 1] = offsets[f + 1]; 4221 /* TODO: closure indices */ 4222 newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]); 4223 } 4224 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 4225 } 4226 else { 4227 offsets[0] = 0; 4228 offsets[1] = lDof; 4229 newOffsets[0] = 0; 4230 newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0]; 4231 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 4232 } 4233 if (childId == -1) { /* no child interpolation: one nnz per */ 4234 ierr = VecSetValues(vecFine,numValues,rowIndices,pVal,INSERT_VALUES);CHKERRQ(ierr); 4235 } else { 4236 PetscInt f; 4237 4238 if (grad && p >= cellStart && p < cellEnd) { 4239 numValues -= (dim * (1 + numFVcomps)); 4240 fvGradData = &pVal[numValues]; 4241 } 4242 for (f = 0; f < PetscMax(1,numFields); f++) { 4243 const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f]; 4244 PetscInt numRows = offsets[f+1] - offsets[f]; 4245 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 4246 const PetscScalar *cVal = &pVal[newOffsets[f]]; 4247 PetscScalar *rVal = &pointWork[offsets[f]]; 4248 PetscInt i, j; 4249 4250 #if 0 4251 ierr = PetscInfo6(coarse,"field %D, numFields %D, childId %D, numRows %D, numCols %D, refPointFieldN %D\n",f,numFields,childId,numRows,numCols,refPointFieldN[childId - pRefStart][f]);CHKERRQ(ierr); 4252 #endif 4253 for (i = 0; i < numRows; i++) { 4254 PetscScalar val = 0.; 4255 for (j = 0; j < numCols; j++) { 4256 val += childMat[i * numCols + j] * cVal[j]; 4257 } 4258 rVal[i] = val; 4259 } 4260 if (f == fvField && p >= cellStart && p < cellEnd) { 4261 PetscReal centroid[3]; 4262 PetscScalar diff[3]; 4263 const PetscScalar *parentCentroid = &fvGradData[0]; 4264 const PetscScalar *gradient = &fvGradData[dim]; 4265 4266 ierr = DMPlexComputeCellGeometryFVM(fine,p,NULL,centroid,NULL);CHKERRQ(ierr); 4267 for (i = 0; i < dim; i++) { 4268 diff[i] = centroid[i] - parentCentroid[i]; 4269 } 4270 for (i = 0; i < numFVcomps; i++) { 4271 PetscScalar val = 0.; 4272 4273 for (j = 0; j < dim; j++) { 4274 val += gradient[dim * i + j] * diff[j]; 4275 } 4276 rVal[i] += val; 4277 } 4278 } 4279 ierr = VecSetValues(vecFine,numRows,&rowIndices[offsets[f]],rVal,INSERT_VALUES);CHKERRQ(ierr); 4280 } 4281 } 4282 } 4283 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 4284 ierr = DMRestoreWorkArray(fine,maxDof,MPIU_INT,&rowIndices);CHKERRQ(ierr); 4285 ierr = DMRestoreWorkArray(fine,maxDof,MPIU_SCALAR,&pointWork);CHKERRQ(ierr); 4286 } 4287 ierr = PetscFree(leafValues);CHKERRQ(ierr); 4288 ierr = PetscSectionDestroy(&leafValuesSec);CHKERRQ(ierr); 4289 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr); 4290 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 4291 PetscFunctionReturn(0); 4292 } 4293 4294 static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids) 4295 { 4296 PetscDS ds; 4297 DM refTree; 4298 PetscSection multiRootSec, rootIndicesSec; 4299 PetscSection globalCoarse, globalFine; 4300 PetscSection localCoarse, localFine; 4301 PetscSection cSecRef; 4302 PetscInt *parentIndices, pRefStart, pRefEnd; 4303 PetscScalar *rootValues, *parentValues; 4304 Mat injRef; 4305 PetscInt numFields, maxDof; 4306 PetscInt pStartC, pEndC, pStartF, pEndF, p; 4307 PetscInt *offsets, *offsetsCopy, *rowOffsets; 4308 PetscLayout rowMap, colMap; 4309 PetscInt rowStart, rowEnd, colStart, colEnd; 4310 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 4311 PetscErrorCode ierr; 4312 4313 PetscFunctionBegin; 4314 4315 /* get the templates for the fine-to-coarse injection from the reference tree */ 4316 ierr = VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4317 ierr = VecSetOption(vecCoarse,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4318 ierr = DMPlexGetReferenceTree(coarse,&refTree);CHKERRQ(ierr); 4319 ierr = DMGetDS(coarse,&ds);CHKERRQ(ierr); 4320 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 4321 ierr = DMGetDefaultConstraints(refTree,&cSecRef,NULL);CHKERRQ(ierr); 4322 ierr = PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);CHKERRQ(ierr); 4323 ierr = DMPlexReferenceTreeGetInjector(refTree,&injRef);CHKERRQ(ierr); 4324 4325 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 4326 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 4327 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 4328 ierr = PetscSectionGetNumFields(localFine,&numFields);CHKERRQ(ierr); 4329 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 4330 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 4331 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 4332 ierr = PetscSectionGetMaxDof(localCoarse,&maxDof);CHKERRQ(ierr); 4333 { 4334 PetscInt maxFields = PetscMax(1,numFields) + 1; 4335 ierr = PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);CHKERRQ(ierr); 4336 } 4337 4338 ierr = DMPlexTransferInjectorTree(coarse,fine,coarseToFine,cids,vecFine,numFields,offsets,&multiRootSec,&rootIndicesSec,NULL,&rootValues);CHKERRQ(ierr); 4339 4340 ierr = PetscMalloc2(maxDof,&parentIndices,maxDof,&parentValues);CHKERRQ(ierr); 4341 4342 /* count indices */ 4343 ierr = VecGetLayout(vecFine,&colMap);CHKERRQ(ierr); 4344 ierr = VecGetLayout(vecCoarse,&rowMap);CHKERRQ(ierr); 4345 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 4346 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 4347 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 4348 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 4349 /* insert values */ 4350 ierr = DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 4351 for (p = pStartC; p < pEndC; p++) { 4352 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 4353 PetscBool contribute = PETSC_FALSE; 4354 4355 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 4356 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 4357 if ((dof - cdof) <= 0) continue; 4358 ierr = PetscSectionGetDof(localCoarse,p,&dof);CHKERRQ(ierr); 4359 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 4360 4361 rowOffsets[0] = 0; 4362 offsetsCopy[0] = 0; 4363 if (numFields) { 4364 PetscInt f; 4365 4366 for (f = 0; f < numFields; f++) { 4367 PetscInt fDof; 4368 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 4369 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 4370 } 4371 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 4372 } 4373 else { 4374 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 4375 rowOffsets[1] = offsetsCopy[0]; 4376 } 4377 4378 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 4379 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 4380 leafEnd = leafStart + numLeaves; 4381 for (l = 0; l < dof; l++) parentValues[l] = 0.; 4382 for (l = leafStart; l < leafEnd; l++) { 4383 PetscInt numIndices, childId, offset; 4384 const PetscScalar *childValues; 4385 4386 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 4387 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 4388 childId = (PetscInt) PetscRealPart(rootValues[offset++]); 4389 childValues = &rootValues[offset]; 4390 numIndices--; 4391 4392 if (childId == -2) { /* skip */ 4393 continue; 4394 } else if (childId == -1) { /* equivalent points: scatter */ 4395 PetscInt m; 4396 4397 contribute = PETSC_TRUE; 4398 for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m]; 4399 } else { /* contributions from children: sum with injectors from reference tree */ 4400 PetscInt parentId, f, lim; 4401 4402 contribute = PETSC_TRUE; 4403 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 4404 4405 lim = PetscMax(1,numFields); 4406 offsets[0] = 0; 4407 if (numFields) { 4408 PetscInt f; 4409 4410 for (f = 0; f < numFields; f++) { 4411 PetscInt fDof; 4412 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 4413 4414 offsets[f + 1] = fDof + offsets[f]; 4415 } 4416 } 4417 else { 4418 PetscInt cDof; 4419 4420 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 4421 offsets[1] = cDof; 4422 } 4423 for (f = 0; f < lim; f++) { 4424 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 4425 PetscInt n = offsets[f+1]-offsets[f]; 4426 PetscInt m = rowOffsets[f+1]-rowOffsets[f]; 4427 PetscInt i, j; 4428 const PetscScalar *colValues = &childValues[offsets[f]]; 4429 4430 for (i = 0; i < m; i++) { 4431 PetscScalar val = 0.; 4432 for (j = 0; j < n; j++) { 4433 val += childMat[n * i + j] * colValues[j]; 4434 } 4435 parentValues[rowOffsets[f] + i] += val; 4436 } 4437 } 4438 } 4439 } 4440 if (contribute) {ierr = VecSetValues(vecCoarse,dof,parentIndices,parentValues,INSERT_VALUES);CHKERRQ(ierr);} 4441 } 4442 ierr = PetscSectionDestroy(&multiRootSec);CHKERRQ(ierr); 4443 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 4444 ierr = PetscFree2(parentIndices,parentValues);CHKERRQ(ierr); 4445 ierr = DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 4446 ierr = PetscFree(rootValues);CHKERRQ(ierr); 4447 ierr = PetscFree3(offsets,offsetsCopy,rowOffsets);CHKERRQ(ierr); 4448 PetscFunctionReturn(0); 4449 } 4450 4451 /*@ 4452 DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening 4453 that can be represented by a common reference tree used by both. This routine can be used for a combination of 4454 coarsening and refinement at the same time. 4455 4456 collective 4457 4458 Input Parameters: 4459 + dmIn - The DMPlex mesh for the input vector 4460 . vecIn - The input vector 4461 . sfRefine - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in 4462 the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn 4463 . sfCoarsen - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in 4464 the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn 4465 . cidsRefine - The childIds of the points in dmOut. These childIds relate back to the reference tree: childid[j] = k implies 4466 that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference 4467 tree was refined from its parent. childid[j] = -1 indicates that the point j in dmOut is exactly 4468 equivalent to its root in dmIn, so no interpolation is necessary. childid[j] = -2 indicates that this 4469 point j in dmOut is not a leaf of sfRefine. 4470 . cidsCoarsen - The childIds of the points in dmIn. These childIds relate back to the reference tree: childid[j] = k implies 4471 that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference 4472 tree coarsens to its parent. childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen. 4473 . useBCs - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer. 4474 - time - Used if boundary values are time dependent. 4475 4476 Output Parameters: 4477 . vecOut - Using interpolation and injection operators calculated on the reference tree, the transfered 4478 projection of vecIn from dmIn to dmOut. Note that any field discretized with a PetscFV finite volume 4479 method that uses gradient reconstruction will use reconstructed gradients when interpolating from 4480 coarse points to fine points. 4481 4482 Level: developer 4483 4484 .seealso(): DMPlexSetReferenceTree(), DMPlexGetReferenceTree(), PetscFVGetComputeGradients() 4485 @*/ 4486 PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time) 4487 { 4488 PetscErrorCode ierr; 4489 4490 PetscFunctionBegin; 4491 ierr = VecSet(vecOut,0.0);CHKERRQ(ierr); 4492 if (sfRefine) { 4493 Vec vecInLocal; 4494 DM dmGrad = NULL; 4495 Vec faceGeom = NULL, cellGeom = NULL, grad = NULL; 4496 4497 ierr = DMGetLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr); 4498 ierr = VecSet(vecInLocal,0.0);CHKERRQ(ierr); 4499 { 4500 PetscInt numFields, i; 4501 4502 ierr = DMGetNumFields(dmIn, &numFields);CHKERRQ(ierr); 4503 for (i = 0; i < numFields; i++) { 4504 PetscObject obj; 4505 PetscClassId classid; 4506 4507 ierr = DMGetField(dmIn, i, &obj);CHKERRQ(ierr); 4508 ierr = PetscObjectGetClassId(obj, &classid);CHKERRQ(ierr); 4509 if (classid == PETSCFV_CLASSID) { 4510 ierr = DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad);CHKERRQ(ierr); 4511 break; 4512 } 4513 } 4514 } 4515 if (useBCs) { 4516 ierr = DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL);CHKERRQ(ierr); 4517 } 4518 ierr = DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr); 4519 ierr = DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr); 4520 if (dmGrad) { 4521 ierr = DMGetGlobalVector(dmGrad,&grad);CHKERRQ(ierr); 4522 ierr = DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad);CHKERRQ(ierr); 4523 } 4524 ierr = DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfRefine,cidsRefine,grad,cellGeom);CHKERRQ(ierr); 4525 ierr = DMRestoreLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr); 4526 if (dmGrad) { 4527 ierr = DMRestoreGlobalVector(dmGrad,&grad);CHKERRQ(ierr); 4528 } 4529 } 4530 if (sfCoarsen) { 4531 ierr = DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfCoarsen,cidsCoarsen);CHKERRQ(ierr); 4532 } 4533 ierr = VecAssemblyBegin(vecOut);CHKERRQ(ierr); 4534 ierr = VecAssemblyEnd(vecOut);CHKERRQ(ierr); 4535 PetscFunctionReturn(0); 4536 } 4537