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