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