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 #undef __FUNCT__ 11 #define __FUNCT__ "DMPlexSetReferenceTree" 12 /*@ 13 DMPlexSetReferenceTree - set the reference tree for hierarchically non-conforming meshes. 14 15 Not collective 16 17 Input Parameters: 18 + dm - The DMPlex object 19 - ref - The reference tree DMPlex object 20 21 Level: intermediate 22 23 .seealso: DMPlexGetReferenceTree(), DMPlexCreateDefaultReferenceTree() 24 @*/ 25 PetscErrorCode DMPlexSetReferenceTree(DM dm, DM ref) 26 { 27 DM_Plex *mesh = (DM_Plex *)dm->data; 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 32 if (ref) {PetscValidHeaderSpecific(ref, DM_CLASSID, 2);} 33 ierr = PetscObjectReference((PetscObject)ref);CHKERRQ(ierr); 34 ierr = DMDestroy(&mesh->referenceTree);CHKERRQ(ierr); 35 mesh->referenceTree = ref; 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "DMPlexGetReferenceTree" 41 /*@ 42 DMPlexGetReferenceTree - get the reference tree for hierarchically non-conforming meshes. 43 44 Not collective 45 46 Input Parameters: 47 . dm - The DMPlex object 48 49 Output Parameters 50 . ref - The reference tree DMPlex object 51 52 Level: intermediate 53 54 .seealso: DMPlexSetReferenceTree(), DMPlexCreateDefaultReferenceTree() 55 @*/ 56 PetscErrorCode DMPlexGetReferenceTree(DM dm, DM *ref) 57 { 58 DM_Plex *mesh = (DM_Plex *)dm->data; 59 60 PetscFunctionBegin; 61 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 62 PetscValidPointer(ref,2); 63 *ref = mesh->referenceTree; 64 PetscFunctionReturn(0); 65 } 66 67 #undef __FUNCT__ 68 #define __FUNCT__ "DMPlexReferenceTreeGetChildSymmetry_Default" 69 static PetscErrorCode DMPlexReferenceTreeGetChildSymmetry_Default(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB) 70 { 71 PetscInt coneSize, dStart, dEnd, dim, ABswap, oAvert, oBvert, ABswapVert; 72 PetscErrorCode ierr; 73 74 PetscFunctionBegin; 75 if (parentOrientA == parentOrientB) { 76 if (childOrientB) *childOrientB = childOrientA; 77 if (childB) *childB = childA; 78 PetscFunctionReturn(0); 79 } 80 for (dim = 0; dim < 3; dim++) { 81 ierr = DMPlexGetDepthStratum(dm,dim,&dStart,&dEnd);CHKERRQ(ierr); 82 if (parent >= dStart && parent <= dEnd) { 83 break; 84 } 85 } 86 if (dim > 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot perform child symmetry for %d-cells",dim); 87 if (!dim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A vertex has no children"); 88 if (childA < dStart || childA >= dEnd) { 89 /* this is a lower-dimensional child: bootstrap */ 90 PetscInt size, i, sA = -1, sB, sOrientB, sConeSize; 91 const PetscInt *supp, *coneA, *coneB, *oA, *oB; 92 93 ierr = DMPlexGetSupportSize(dm,childA,&size);CHKERRQ(ierr); 94 ierr = DMPlexGetSupport(dm,childA,&supp);CHKERRQ(ierr); 95 96 /* find a point sA in supp(childA) that has the same parent */ 97 for (i = 0; i < size; i++) { 98 PetscInt sParent; 99 100 sA = supp[i]; 101 if (sA == parent) continue; 102 ierr = DMPlexGetTreeParent(dm,sA,&sParent,NULL);CHKERRQ(ierr); 103 if (sParent == parent) { 104 break; 105 } 106 } 107 if (i == size) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"could not find support in children"); 108 /* find out which point sB is in an equivalent position to sA under 109 * parentOrientB */ 110 ierr = DMPlexReferenceTreeGetChildSymmetry_Default(dm,parent,parentOrientA,0,sA,parentOrientB,&sOrientB,&sB);CHKERRQ(ierr); 111 ierr = DMPlexGetConeSize(dm,sA,&sConeSize);CHKERRQ(ierr); 112 ierr = DMPlexGetCone(dm,sA,&coneA);CHKERRQ(ierr); 113 ierr = DMPlexGetCone(dm,sB,&coneB);CHKERRQ(ierr); 114 ierr = DMPlexGetConeOrientation(dm,sA,&oA);CHKERRQ(ierr); 115 ierr = DMPlexGetConeOrientation(dm,sB,&oB);CHKERRQ(ierr); 116 /* step through the cone of sA in natural order */ 117 for (i = 0; i < sConeSize; i++) { 118 if (coneA[i] == childA) { 119 /* if childA is at position i in coneA, 120 * then we want the point that is at sOrientB*i in coneB */ 121 PetscInt j = (sOrientB >= 0) ? ((sOrientB + i) % sConeSize) : ((sConeSize -(sOrientB+1) - i) % sConeSize); 122 if (childB) *childB = coneB[j]; 123 if (childOrientB) { 124 PetscInt oBtrue; 125 126 ierr = DMPlexGetConeSize(dm,childA,&coneSize);CHKERRQ(ierr); 127 /* compose sOrientB and oB[j] */ 128 if (coneSize != 0 && coneSize != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected a vertex or an edge"); 129 /* we may have to flip an edge */ 130 oBtrue = coneSize ? ((sOrientB >= 0) ? oB[j] : -(oB[j] + 2)) : 0; 131 ABswap = DihedralSwap(coneSize,oA[i],oBtrue);CHKERRQ(ierr); 132 *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap); 133 } 134 break; 135 } 136 } 137 if (i == sConeSize) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"support cone mismatch"); 138 PetscFunctionReturn(0); 139 } 140 /* get the cone size and symmetry swap */ 141 ierr = DMPlexGetConeSize(dm,parent,&coneSize);CHKERRQ(ierr); 142 ABswap = DihedralSwap(coneSize, parentOrientA, parentOrientB); 143 if (dim == 2) { 144 /* orientations refer to cones: we want them to refer to vertices: 145 * if it's a rotation, they are the same, but if the order is reversed, a 146 * permutation that puts side i first does *not* put vertex i first */ 147 oAvert = (parentOrientA >= 0) ? parentOrientA : -((-parentOrientA % coneSize) + 1); 148 oBvert = (parentOrientB >= 0) ? parentOrientB : -((-parentOrientB % coneSize) + 1); 149 ABswapVert = DihedralSwap(coneSize, oAvert, oBvert); 150 } 151 else { 152 ABswapVert = ABswap; 153 } 154 if (childB) { 155 /* assume that each child corresponds to a vertex, in the same order */ 156 PetscInt p, posA = -1, numChildren, i; 157 const PetscInt *children; 158 159 /* count which position the child is in */ 160 ierr = DMPlexGetTreeChildren(dm,parent,&numChildren,&children);CHKERRQ(ierr); 161 for (i = 0; i < numChildren; i++) { 162 p = children[i]; 163 if (p == childA) { 164 posA = i; 165 break; 166 } 167 } 168 if (posA >= coneSize) { 169 /* this is the triangle in the middle of a uniformly refined triangle: it 170 * is invariant */ 171 if (dim != 2 || posA != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Expected a middle triangle, got something else"); 172 *childB = childA; 173 } 174 else { 175 /* figure out position B by applying ABswapVert */ 176 PetscInt posB; 177 178 posB = (ABswapVert >= 0) ? ((ABswapVert + posA) % coneSize) : ((coneSize -(ABswapVert + 1) - posA) % coneSize); 179 if (childB) *childB = children[posB]; 180 } 181 } 182 if (childOrientB) *childOrientB = DihedralCompose(coneSize,childOrientA,ABswap); 183 PetscFunctionReturn(0); 184 } 185 186 #undef __FUNCT__ 187 #define __FUNCT__ "DMPlexReferenceTreeGetChildSymmetry" 188 /*@ 189 DMPlexReferenceTreeGetChildSymmetry - Given a reference tree, transform a childid and orientation from one parent frame to another 190 191 Input Parameters: 192 + dm - the reference tree DMPlex object 193 . parent - the parent point 194 . parentOrientA - the reference orientation for describing the parent 195 . childOrientA - the reference orientation for describing the child 196 . childA - the reference childID for describing the child 197 - parentOrientB - the new orientation for describing the parent 198 199 Output Parameters: 200 + childOrientB - if not NULL, set to the new oreintation for describing the child 201 . childB - if not NULL, the new childID for describing the child 202 203 Level: developer 204 205 .seealso: DMPlexGetReferenceTree(), DMPlexSetReferenceTree(), DMPlexSetTree() 206 @*/ 207 PetscErrorCode DMPlexReferenceTreeGetChildSymmetry(DM dm, PetscInt parent, PetscInt parentOrientA, PetscInt childOrientA, PetscInt childA, PetscInt parentOrientB, PetscInt *childOrientB, PetscInt *childB) 208 { 209 DM_Plex *mesh = (DM_Plex *)dm->data; 210 PetscErrorCode ierr; 211 212 PetscFunctionBegin; 213 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 214 if (!mesh->getchildsymmetry) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DMPlexReferenceTreeGetChildSymmetry not implemented"); 215 ierr = mesh->getchildsymmetry(dm,parent,parentOrientA,childOrientA,childA,parentOrientB,childOrientB,childB);CHKERRQ(ierr); 216 PetscFunctionReturn(0); 217 } 218 219 static PetscErrorCode DMPlexSetTree_Internal(DM,PetscSection,PetscInt*,PetscInt*,PetscBool,PetscBool); 220 221 #undef __FUNCT__ 222 #define __FUNCT__ "DMPlexCreateReferenceTree_Union" 223 PetscErrorCode DMPlexCreateReferenceTree_Union(DM K, DM Kref, const char *labelName, DM *ref) 224 { 225 MPI_Comm comm; 226 PetscInt dim, p, pStart, pEnd, pRefStart, pRefEnd, d, offset, parentSize, *parents, *childIDs; 227 PetscInt *permvals, *unionCones, *coneSizes, *unionOrientations, numUnionPoints, *numDimPoints, numCones, numVerts; 228 DMLabel identity, identityRef; 229 PetscSection unionSection, unionConeSection, parentSection; 230 PetscScalar *unionCoords; 231 IS perm; 232 PetscErrorCode ierr; 233 234 PetscFunctionBegin; 235 comm = PetscObjectComm((PetscObject)K); 236 ierr = DMGetDimension(K, &dim);CHKERRQ(ierr); 237 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr); 238 ierr = DMGetLabel(K, labelName, &identity);CHKERRQ(ierr); 239 ierr = DMGetLabel(Kref, labelName, &identityRef);CHKERRQ(ierr); 240 ierr = DMPlexGetChart(Kref, &pRefStart, &pRefEnd);CHKERRQ(ierr); 241 ierr = PetscSectionCreate(comm, &unionSection);CHKERRQ(ierr); 242 ierr = PetscSectionSetChart(unionSection, 0, (pEnd - pStart) + (pRefEnd - pRefStart));CHKERRQ(ierr); 243 /* count points that will go in the union */ 244 for (p = pStart; p < pEnd; p++) { 245 ierr = PetscSectionSetDof(unionSection, p - pStart, 1);CHKERRQ(ierr); 246 } 247 for (p = pRefStart; p < pRefEnd; p++) { 248 PetscInt q, qSize; 249 ierr = DMLabelGetValue(identityRef, p, &q);CHKERRQ(ierr); 250 ierr = DMLabelGetStratumSize(identityRef, q, &qSize);CHKERRQ(ierr); 251 if (qSize > 1) { 252 ierr = PetscSectionSetDof(unionSection, p - pRefStart + (pEnd - pStart), 1);CHKERRQ(ierr); 253 } 254 } 255 ierr = PetscMalloc1(pEnd - pStart + pRefEnd - pRefStart,&permvals);CHKERRQ(ierr); 256 offset = 0; 257 /* stratify points in the union by topological dimension */ 258 for (d = 0; d <= dim; d++) { 259 PetscInt cStart, cEnd, c; 260 261 ierr = DMPlexGetHeightStratum(K, d, &cStart, &cEnd);CHKERRQ(ierr); 262 for (c = cStart; c < cEnd; c++) { 263 permvals[offset++] = c; 264 } 265 266 ierr = DMPlexGetHeightStratum(Kref, d, &cStart, &cEnd);CHKERRQ(ierr); 267 for (c = cStart; c < cEnd; c++) { 268 permvals[offset++] = c + (pEnd - pStart); 269 } 270 } 271 ierr = ISCreateGeneral(comm, (pEnd - pStart) + (pRefEnd - pRefStart), permvals, PETSC_OWN_POINTER, &perm);CHKERRQ(ierr); 272 ierr = PetscSectionSetPermutation(unionSection,perm);CHKERRQ(ierr); 273 ierr = PetscSectionSetUp(unionSection);CHKERRQ(ierr); 274 ierr = PetscSectionGetStorageSize(unionSection,&numUnionPoints);CHKERRQ(ierr); 275 ierr = PetscMalloc2(numUnionPoints,&coneSizes,dim+1,&numDimPoints);CHKERRQ(ierr); 276 /* count dimension points */ 277 for (d = 0; d <= dim; d++) { 278 PetscInt cStart, cOff, cOff2; 279 ierr = DMPlexGetHeightStratum(K,d,&cStart,NULL);CHKERRQ(ierr); 280 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff);CHKERRQ(ierr); 281 if (d < dim) { 282 ierr = DMPlexGetHeightStratum(K,d+1,&cStart,NULL);CHKERRQ(ierr); 283 ierr = PetscSectionGetOffset(unionSection,cStart-pStart,&cOff2);CHKERRQ(ierr); 284 } 285 else { 286 cOff2 = numUnionPoints; 287 } 288 numDimPoints[dim - d] = cOff2 - cOff; 289 } 290 ierr = PetscSectionCreate(comm, &unionConeSection);CHKERRQ(ierr); 291 ierr = PetscSectionSetChart(unionConeSection, 0, numUnionPoints);CHKERRQ(ierr); 292 /* count the cones in the union */ 293 for (p = pStart; p < pEnd; p++) { 294 PetscInt dof, uOff; 295 296 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 297 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 298 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 299 coneSizes[uOff] = dof; 300 } 301 for (p = pRefStart; p < pRefEnd; p++) { 302 PetscInt dof, uDof, uOff; 303 304 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 305 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 306 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 307 if (uDof) { 308 ierr = PetscSectionSetDof(unionConeSection, uOff, dof);CHKERRQ(ierr); 309 coneSizes[uOff] = dof; 310 } 311 } 312 ierr = PetscSectionSetUp(unionConeSection);CHKERRQ(ierr); 313 ierr = PetscSectionGetStorageSize(unionConeSection,&numCones);CHKERRQ(ierr); 314 ierr = PetscMalloc2(numCones,&unionCones,numCones,&unionOrientations);CHKERRQ(ierr); 315 /* write the cones in the union */ 316 for (p = pStart; p < pEnd; p++) { 317 PetscInt dof, uOff, c, cOff; 318 const PetscInt *cone, *orientation; 319 320 ierr = DMPlexGetConeSize(K, p, &dof);CHKERRQ(ierr); 321 ierr = DMPlexGetCone(K, p, &cone);CHKERRQ(ierr); 322 ierr = DMPlexGetConeOrientation(K, p, &orientation);CHKERRQ(ierr); 323 ierr = PetscSectionGetOffset(unionSection, p - pStart,&uOff);CHKERRQ(ierr); 324 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 325 for (c = 0; c < dof; c++) { 326 PetscInt e, eOff; 327 e = cone[c]; 328 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 329 unionCones[cOff + c] = eOff; 330 unionOrientations[cOff + c] = orientation[c]; 331 } 332 } 333 for (p = pRefStart; p < pRefEnd; p++) { 334 PetscInt dof, uDof, uOff, c, cOff; 335 const PetscInt *cone, *orientation; 336 337 ierr = DMPlexGetConeSize(Kref, p, &dof);CHKERRQ(ierr); 338 ierr = DMPlexGetCone(Kref, p, &cone);CHKERRQ(ierr); 339 ierr = DMPlexGetConeOrientation(Kref, p, &orientation);CHKERRQ(ierr); 340 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 341 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 342 if (uDof) { 343 ierr = PetscSectionGetOffset(unionConeSection,uOff,&cOff);CHKERRQ(ierr); 344 for (c = 0; c < dof; c++) { 345 PetscInt e, eOff, eDof; 346 347 e = cone[c]; 348 ierr = PetscSectionGetDof(unionSection, e - pRefStart + (pEnd - pStart),&eDof);CHKERRQ(ierr); 349 if (eDof) { 350 ierr = PetscSectionGetOffset(unionSection, e - pRefStart + (pEnd - pStart), &eOff);CHKERRQ(ierr); 351 } 352 else { 353 ierr = DMLabelGetValue(identityRef, e, &e);CHKERRQ(ierr); 354 ierr = PetscSectionGetOffset(unionSection, e - pStart, &eOff);CHKERRQ(ierr); 355 } 356 unionCones[cOff + c] = eOff; 357 unionOrientations[cOff + c] = orientation[c]; 358 } 359 } 360 } 361 /* get the coordinates */ 362 { 363 PetscInt vStart, vEnd, vRefStart, vRefEnd, v, vDof, vOff; 364 PetscSection KcoordsSec, KrefCoordsSec; 365 Vec KcoordsVec, KrefCoordsVec; 366 PetscScalar *Kcoords; 367 368 DMGetCoordinateSection(K, &KcoordsSec);CHKERRQ(ierr); 369 DMGetCoordinatesLocal(K, &KcoordsVec);CHKERRQ(ierr); 370 DMGetCoordinateSection(Kref, &KrefCoordsSec);CHKERRQ(ierr); 371 DMGetCoordinatesLocal(Kref, &KrefCoordsVec);CHKERRQ(ierr); 372 373 numVerts = numDimPoints[0]; 374 ierr = PetscMalloc1(numVerts * dim,&unionCoords);CHKERRQ(ierr); 375 ierr = DMPlexGetDepthStratum(K,0,&vStart,&vEnd);CHKERRQ(ierr); 376 377 offset = 0; 378 for (v = vStart; v < vEnd; v++) { 379 ierr = PetscSectionGetOffset(unionSection,v - pStart,&vOff);CHKERRQ(ierr); 380 ierr = VecGetValuesSection(KcoordsVec, KcoordsSec, v, &Kcoords);CHKERRQ(ierr); 381 for (d = 0; d < dim; d++) { 382 unionCoords[offset * dim + d] = Kcoords[d]; 383 } 384 offset++; 385 } 386 ierr = DMPlexGetDepthStratum(Kref,0,&vRefStart,&vRefEnd);CHKERRQ(ierr); 387 for (v = vRefStart; v < vRefEnd; v++) { 388 ierr = PetscSectionGetDof(unionSection,v - pRefStart + (pEnd - pStart),&vDof);CHKERRQ(ierr); 389 ierr = PetscSectionGetOffset(unionSection,v - pRefStart + (pEnd - pStart),&vOff);CHKERRQ(ierr); 390 ierr = VecGetValuesSection(KrefCoordsVec, KrefCoordsSec, v, &Kcoords);CHKERRQ(ierr); 391 if (vDof) { 392 for (d = 0; d < dim; d++) { 393 unionCoords[offset * dim + d] = Kcoords[d]; 394 } 395 offset++; 396 } 397 } 398 } 399 ierr = DMCreate(comm,ref);CHKERRQ(ierr); 400 ierr = DMSetType(*ref,DMPLEX);CHKERRQ(ierr); 401 ierr = DMSetDimension(*ref,dim);CHKERRQ(ierr); 402 ierr = DMPlexCreateFromDAG(*ref,dim,numDimPoints,coneSizes,unionCones,unionOrientations,unionCoords);CHKERRQ(ierr); 403 /* set the tree */ 404 ierr = PetscSectionCreate(comm,&parentSection);CHKERRQ(ierr); 405 ierr = PetscSectionSetChart(parentSection,0,numUnionPoints);CHKERRQ(ierr); 406 for (p = pRefStart; p < pRefEnd; p++) { 407 PetscInt uDof, uOff; 408 409 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 410 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 411 if (uDof) { 412 PetscSectionSetDof(parentSection,uOff,1);CHKERRQ(ierr); 413 } 414 } 415 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 416 ierr = PetscSectionGetStorageSize(parentSection,&parentSize);CHKERRQ(ierr); 417 ierr = PetscMalloc2(parentSize,&parents,parentSize,&childIDs);CHKERRQ(ierr); 418 for (p = pRefStart; p < pRefEnd; p++) { 419 PetscInt uDof, uOff; 420 421 ierr = PetscSectionGetDof(unionSection, p - pRefStart + (pEnd - pStart),&uDof);CHKERRQ(ierr); 422 ierr = PetscSectionGetOffset(unionSection, p - pRefStart + (pEnd - pStart),&uOff);CHKERRQ(ierr); 423 if (uDof) { 424 PetscInt pOff, parent, parentU; 425 PetscSectionGetOffset(parentSection,uOff,&pOff);CHKERRQ(ierr); 426 DMLabelGetValue(identityRef,p,&parent);CHKERRQ(ierr); 427 ierr = PetscSectionGetOffset(unionSection, parent - pStart,&parentU);CHKERRQ(ierr); 428 parents[pOff] = parentU; 429 childIDs[pOff] = uOff; 430 } 431 } 432 ierr = DMPlexSetTree_Internal(*ref,parentSection,parents,childIDs,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 433 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 434 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 435 436 /* clean up */ 437 ierr = PetscSectionDestroy(&unionSection);CHKERRQ(ierr); 438 ierr = PetscSectionDestroy(&unionConeSection);CHKERRQ(ierr); 439 ierr = ISDestroy(&perm);CHKERRQ(ierr); 440 ierr = PetscFree(unionCoords);CHKERRQ(ierr); 441 ierr = PetscFree2(unionCones,unionOrientations);CHKERRQ(ierr); 442 ierr = PetscFree2(coneSizes,numDimPoints);CHKERRQ(ierr); 443 PetscFunctionReturn(0); 444 } 445 446 #undef __FUNCT__ 447 #define __FUNCT__ "DMPlexCreateDefaultReferenceTree" 448 /*@ 449 DMPlexCreateDefaultReferenceTree - create a reference tree for isotropic hierarchical mesh refinement. 450 451 Collective on comm 452 453 Input Parameters: 454 + comm - the MPI communicator 455 . dim - the spatial dimension 456 - simplex - Flag for simplex, otherwise use a tensor-product cell 457 458 Output Parameters: 459 . ref - the reference tree DMPlex object 460 461 Level: intermediate 462 463 .keywords: reference cell 464 .seealso: DMPlexSetReferenceTree(), DMPlexGetReferenceTree() 465 @*/ 466 PetscErrorCode DMPlexCreateDefaultReferenceTree(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *ref) 467 { 468 DM_Plex *mesh; 469 DM K, Kref; 470 PetscInt p, pStart, pEnd; 471 DMLabel identity; 472 PetscErrorCode ierr; 473 474 PetscFunctionBegin; 475 #if 1 476 comm = PETSC_COMM_SELF; 477 #endif 478 /* create a reference element */ 479 ierr = DMPlexCreateReferenceCell(comm, dim, simplex, &K);CHKERRQ(ierr); 480 ierr = DMCreateLabel(K, "identity");CHKERRQ(ierr); 481 ierr = DMGetLabel(K, "identity", &identity);CHKERRQ(ierr); 482 ierr = DMPlexGetChart(K, &pStart, &pEnd);CHKERRQ(ierr); 483 for (p = pStart; p < pEnd; p++) { 484 ierr = DMLabelSetValue(identity, p, p);CHKERRQ(ierr); 485 } 486 /* refine it */ 487 ierr = DMRefine(K,comm,&Kref);CHKERRQ(ierr); 488 489 /* the reference tree is the union of these two, without duplicating 490 * points that appear in both */ 491 ierr = DMPlexCreateReferenceTree_Union(K, Kref, "identity", ref);CHKERRQ(ierr); 492 mesh = (DM_Plex *) (*ref)->data; 493 mesh->getchildsymmetry = DMPlexReferenceTreeGetChildSymmetry_Default; 494 ierr = DMDestroy(&K);CHKERRQ(ierr); 495 ierr = DMDestroy(&Kref);CHKERRQ(ierr); 496 PetscFunctionReturn(0); 497 } 498 499 #undef __FUNCT__ 500 #define __FUNCT__ "DMPlexTreeSymmetrize" 501 static PetscErrorCode DMPlexTreeSymmetrize(DM dm) 502 { 503 DM_Plex *mesh = (DM_Plex *)dm->data; 504 PetscSection childSec, pSec; 505 PetscInt p, pSize, cSize, parMax = PETSC_MIN_INT, parMin = PETSC_MAX_INT; 506 PetscInt *offsets, *children, pStart, pEnd; 507 PetscErrorCode ierr; 508 509 PetscFunctionBegin; 510 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 511 ierr = PetscSectionDestroy(&mesh->childSection);CHKERRQ(ierr); 512 ierr = PetscFree(mesh->children);CHKERRQ(ierr); 513 pSec = mesh->parentSection; 514 if (!pSec) PetscFunctionReturn(0); 515 ierr = PetscSectionGetStorageSize(pSec,&pSize);CHKERRQ(ierr); 516 for (p = 0; p < pSize; p++) { 517 PetscInt par = mesh->parents[p]; 518 519 parMax = PetscMax(parMax,par+1); 520 parMin = PetscMin(parMin,par); 521 } 522 if (parMin > parMax) { 523 parMin = -1; 524 parMax = -1; 525 } 526 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)pSec),&childSec);CHKERRQ(ierr); 527 ierr = PetscSectionSetChart(childSec,parMin,parMax);CHKERRQ(ierr); 528 for (p = 0; p < pSize; p++) { 529 PetscInt par = mesh->parents[p]; 530 531 ierr = PetscSectionAddDof(childSec,par,1);CHKERRQ(ierr); 532 } 533 ierr = PetscSectionSetUp(childSec);CHKERRQ(ierr); 534 ierr = PetscSectionGetStorageSize(childSec,&cSize);CHKERRQ(ierr); 535 ierr = PetscMalloc1(cSize,&children);CHKERRQ(ierr); 536 ierr = PetscCalloc1(parMax-parMin,&offsets);CHKERRQ(ierr); 537 ierr = PetscSectionGetChart(pSec,&pStart,&pEnd);CHKERRQ(ierr); 538 for (p = pStart; p < pEnd; p++) { 539 PetscInt dof, off, i; 540 541 ierr = PetscSectionGetDof(pSec,p,&dof);CHKERRQ(ierr); 542 ierr = PetscSectionGetOffset(pSec,p,&off);CHKERRQ(ierr); 543 for (i = 0; i < dof; i++) { 544 PetscInt par = mesh->parents[off + i], cOff; 545 546 ierr = PetscSectionGetOffset(childSec,par,&cOff);CHKERRQ(ierr); 547 children[cOff + offsets[par-parMin]++] = p; 548 } 549 } 550 mesh->childSection = childSec; 551 mesh->children = children; 552 ierr = PetscFree(offsets);CHKERRQ(ierr); 553 PetscFunctionReturn(0); 554 } 555 556 #undef __FUNCT__ 557 #define __FUNCT__ "AnchorsFlatten" 558 static PetscErrorCode AnchorsFlatten (PetscSection section, IS is, PetscSection *sectionNew, IS *isNew) 559 { 560 PetscInt pStart, pEnd, size, sizeNew, i, p, *valsNew = NULL; 561 const PetscInt *vals; 562 PetscSection secNew; 563 PetscBool anyNew, globalAnyNew; 564 PetscBool compress; 565 PetscErrorCode ierr; 566 567 PetscFunctionBegin; 568 ierr = PetscSectionGetChart(section,&pStart,&pEnd);CHKERRQ(ierr); 569 ierr = ISGetLocalSize(is,&size);CHKERRQ(ierr); 570 ierr = ISGetIndices(is,&vals);CHKERRQ(ierr); 571 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secNew);CHKERRQ(ierr); 572 ierr = PetscSectionSetChart(secNew,pStart,pEnd);CHKERRQ(ierr); 573 for (i = 0; i < size; i++) { 574 PetscInt dof; 575 576 p = vals[i]; 577 if (p < pStart || p >= pEnd) continue; 578 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 579 if (dof) break; 580 } 581 if (i == size) { 582 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 583 anyNew = PETSC_FALSE; 584 compress = PETSC_FALSE; 585 sizeNew = 0; 586 } 587 else { 588 anyNew = PETSC_TRUE; 589 for (p = pStart; p < pEnd; p++) { 590 PetscInt dof, off; 591 592 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 593 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 594 for (i = 0; i < dof; i++) { 595 PetscInt q = vals[off + i], qDof = 0; 596 597 if (q >= pStart && q < pEnd) { 598 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 599 } 600 if (qDof) { 601 ierr = PetscSectionAddDof(secNew, p, qDof);CHKERRQ(ierr); 602 } 603 else { 604 ierr = PetscSectionAddDof(secNew, p, 1);CHKERRQ(ierr); 605 } 606 } 607 } 608 ierr = PetscSectionSetUp(secNew);CHKERRQ(ierr); 609 ierr = PetscSectionGetStorageSize(secNew,&sizeNew);CHKERRQ(ierr); 610 ierr = PetscMalloc1(sizeNew,&valsNew);CHKERRQ(ierr); 611 compress = PETSC_FALSE; 612 for (p = pStart; p < pEnd; p++) { 613 PetscInt dof, off, count, offNew, dofNew; 614 615 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 616 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 617 ierr = PetscSectionGetDof(secNew, p, &dofNew);CHKERRQ(ierr); 618 ierr = PetscSectionGetOffset(secNew, p, &offNew);CHKERRQ(ierr); 619 count = 0; 620 for (i = 0; i < dof; i++) { 621 PetscInt q = vals[off + i], qDof = 0, qOff = 0, j; 622 623 if (q >= pStart && q < pEnd) { 624 ierr = PetscSectionGetDof(section, q, &qDof);CHKERRQ(ierr); 625 ierr = PetscSectionGetOffset(section, q, &qOff);CHKERRQ(ierr); 626 } 627 if (qDof) { 628 PetscInt oldCount = count; 629 630 for (j = 0; j < qDof; j++) { 631 PetscInt k, r = vals[qOff + j]; 632 633 for (k = 0; k < oldCount; k++) { 634 if (valsNew[offNew + k] == r) { 635 break; 636 } 637 } 638 if (k == oldCount) { 639 valsNew[offNew + count++] = r; 640 } 641 } 642 } 643 else { 644 PetscInt k, oldCount = count; 645 646 for (k = 0; k < oldCount; k++) { 647 if (valsNew[offNew + k] == q) { 648 break; 649 } 650 } 651 if (k == oldCount) { 652 valsNew[offNew + count++] = q; 653 } 654 } 655 } 656 if (count < dofNew) { 657 ierr = PetscSectionSetDof(secNew, p, count);CHKERRQ(ierr); 658 compress = PETSC_TRUE; 659 } 660 } 661 } 662 ierr = ISRestoreIndices(is,&vals);CHKERRQ(ierr); 663 ierr = MPIU_Allreduce(&anyNew,&globalAnyNew,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 664 if (!globalAnyNew) { 665 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 666 *sectionNew = NULL; 667 *isNew = NULL; 668 } 669 else { 670 PetscBool globalCompress; 671 672 ierr = MPIU_Allreduce(&compress,&globalCompress,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)secNew));CHKERRQ(ierr); 673 if (compress) { 674 PetscSection secComp; 675 PetscInt *valsComp = NULL; 676 677 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&secComp);CHKERRQ(ierr); 678 ierr = PetscSectionSetChart(secComp,pStart,pEnd);CHKERRQ(ierr); 679 for (p = pStart; p < pEnd; p++) { 680 PetscInt dof; 681 682 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 683 ierr = PetscSectionSetDof(secComp, p, dof);CHKERRQ(ierr); 684 } 685 ierr = PetscSectionSetUp(secComp);CHKERRQ(ierr); 686 ierr = PetscSectionGetStorageSize(secComp,&sizeNew);CHKERRQ(ierr); 687 ierr = PetscMalloc1(sizeNew,&valsComp);CHKERRQ(ierr); 688 for (p = pStart; p < pEnd; p++) { 689 PetscInt dof, off, offNew, j; 690 691 ierr = PetscSectionGetDof(secNew, p, &dof);CHKERRQ(ierr); 692 ierr = PetscSectionGetOffset(secNew, p, &off);CHKERRQ(ierr); 693 ierr = PetscSectionGetOffset(secComp, p, &offNew);CHKERRQ(ierr); 694 for (j = 0; j < dof; j++) { 695 valsComp[offNew + j] = valsNew[off + j]; 696 } 697 } 698 ierr = PetscSectionDestroy(&secNew);CHKERRQ(ierr); 699 secNew = secComp; 700 ierr = PetscFree(valsNew);CHKERRQ(ierr); 701 valsNew = valsComp; 702 } 703 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is),sizeNew,valsNew,PETSC_OWN_POINTER,isNew);CHKERRQ(ierr); 704 } 705 PetscFunctionReturn(0); 706 } 707 708 #undef __FUNCT__ 709 #define __FUNCT__ "DMPlexCreateAnchors_Tree" 710 static PetscErrorCode DMPlexCreateAnchors_Tree(DM dm) 711 { 712 PetscInt p, pStart, pEnd, *anchors, size; 713 PetscInt aMin = PETSC_MAX_INT, aMax = PETSC_MIN_INT; 714 PetscSection aSec; 715 DMLabel canonLabel; 716 IS aIS; 717 PetscErrorCode ierr; 718 719 PetscFunctionBegin; 720 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 721 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 722 ierr = DMGetLabel(dm,"canonical",&canonLabel);CHKERRQ(ierr); 723 for (p = pStart; p < pEnd; p++) { 724 PetscInt parent; 725 726 if (canonLabel) { 727 PetscInt canon; 728 729 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 730 if (p != canon) continue; 731 } 732 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 733 if (parent != p) { 734 aMin = PetscMin(aMin,p); 735 aMax = PetscMax(aMax,p+1); 736 } 737 } 738 if (aMin > aMax) { 739 aMin = -1; 740 aMax = -1; 741 } 742 ierr = PetscSectionCreate(PETSC_COMM_SELF,&aSec);CHKERRQ(ierr); 743 ierr = PetscSectionSetChart(aSec,aMin,aMax);CHKERRQ(ierr); 744 for (p = aMin; p < aMax; p++) { 745 PetscInt parent, ancestor = p; 746 747 if (canonLabel) { 748 PetscInt canon; 749 750 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 751 if (p != canon) continue; 752 } 753 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 754 while (parent != ancestor) { 755 ancestor = parent; 756 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 757 } 758 if (ancestor != p) { 759 PetscInt closureSize, *closure = NULL; 760 761 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 762 ierr = PetscSectionSetDof(aSec,p,closureSize);CHKERRQ(ierr); 763 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 764 } 765 } 766 ierr = PetscSectionSetUp(aSec);CHKERRQ(ierr); 767 ierr = PetscSectionGetStorageSize(aSec,&size);CHKERRQ(ierr); 768 ierr = PetscMalloc1(size,&anchors);CHKERRQ(ierr); 769 for (p = aMin; p < aMax; p++) { 770 PetscInt parent, ancestor = p; 771 772 if (canonLabel) { 773 PetscInt canon; 774 775 ierr = DMLabelGetValue(canonLabel,p,&canon);CHKERRQ(ierr); 776 if (p != canon) continue; 777 } 778 ierr = DMPlexGetTreeParent(dm,p,&parent,NULL);CHKERRQ(ierr); 779 while (parent != ancestor) { 780 ancestor = parent; 781 ierr = DMPlexGetTreeParent(dm,ancestor,&parent,NULL);CHKERRQ(ierr); 782 } 783 if (ancestor != p) { 784 PetscInt j, closureSize, *closure = NULL, aOff; 785 786 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 787 788 ierr = DMPlexGetTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 789 for (j = 0; j < closureSize; j++) { 790 anchors[aOff + j] = closure[2*j]; 791 } 792 ierr = DMPlexRestoreTransitiveClosure(dm,ancestor,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 793 } 794 } 795 ierr = ISCreateGeneral(PETSC_COMM_SELF,size,anchors,PETSC_OWN_POINTER,&aIS);CHKERRQ(ierr); 796 { 797 PetscSection aSecNew = aSec; 798 IS aISNew = aIS; 799 800 ierr = PetscObjectReference((PetscObject)aSec);CHKERRQ(ierr); 801 ierr = PetscObjectReference((PetscObject)aIS);CHKERRQ(ierr); 802 while (aSecNew) { 803 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 804 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 805 aSec = aSecNew; 806 aIS = aISNew; 807 aSecNew = NULL; 808 aISNew = NULL; 809 ierr = AnchorsFlatten(aSec,aIS,&aSecNew,&aISNew);CHKERRQ(ierr); 810 } 811 } 812 ierr = DMPlexSetAnchors(dm,aSec,aIS);CHKERRQ(ierr); 813 ierr = PetscSectionDestroy(&aSec);CHKERRQ(ierr); 814 ierr = ISDestroy(&aIS);CHKERRQ(ierr); 815 PetscFunctionReturn(0); 816 } 817 818 #undef __FUNCT__ 819 #define __FUNCT__ "DMPlexGetTrueSupportSize" 820 static PetscErrorCode DMPlexGetTrueSupportSize(DM dm,PetscInt p,PetscInt *dof,PetscInt *numTrueSupp) 821 { 822 PetscErrorCode ierr; 823 824 PetscFunctionBegin; 825 if (numTrueSupp[p] == -1) { 826 PetscInt i, alldof; 827 const PetscInt *supp; 828 PetscInt count = 0; 829 830 ierr = DMPlexGetSupportSize(dm,p,&alldof);CHKERRQ(ierr); 831 ierr = DMPlexGetSupport(dm,p,&supp);CHKERRQ(ierr); 832 for (i = 0; i < alldof; i++) { 833 PetscInt q = supp[i], numCones, j; 834 const PetscInt *cone; 835 836 ierr = DMPlexGetConeSize(dm,q,&numCones);CHKERRQ(ierr); 837 ierr = DMPlexGetCone(dm,q,&cone);CHKERRQ(ierr); 838 for (j = 0; j < numCones; j++) { 839 if (cone[j] == p) break; 840 } 841 if (j < numCones) count++; 842 } 843 numTrueSupp[p] = count; 844 } 845 *dof = numTrueSupp[p]; 846 PetscFunctionReturn(0); 847 } 848 849 #undef __FUNCT__ 850 #define __FUNCT__ "DMPlexTreeExchangeSupports" 851 static PetscErrorCode DMPlexTreeExchangeSupports(DM dm) 852 { 853 DM_Plex *mesh = (DM_Plex *)dm->data; 854 PetscSection newSupportSection; 855 PetscInt newSize, *newSupports, pStart, pEnd, p, d, depth; 856 PetscInt *numTrueSupp; 857 PetscInt *offsets; 858 PetscErrorCode ierr; 859 860 PetscFunctionBegin; 861 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 862 /* symmetrize the hierarchy */ 863 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 864 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)(mesh->supportSection)),&newSupportSection);CHKERRQ(ierr); 865 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 866 ierr = PetscSectionSetChart(newSupportSection,pStart,pEnd);CHKERRQ(ierr); 867 ierr = PetscCalloc1(pEnd,&offsets);CHKERRQ(ierr); 868 ierr = PetscMalloc1(pEnd,&numTrueSupp);CHKERRQ(ierr); 869 for (p = 0; p < pEnd; p++) numTrueSupp[p] = -1; 870 /* if a point is in the (true) support of q, it should be in the support of 871 * parent(q) */ 872 for (d = 0; d <= depth; d++) { 873 ierr = DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);CHKERRQ(ierr); 874 for (p = pStart; p < pEnd; ++p) { 875 PetscInt dof, q, qdof, parent; 876 877 ierr = DMPlexGetTrueSupportSize(dm,p,&dof,numTrueSupp);CHKERRQ(ierr); 878 ierr = PetscSectionAddDof(newSupportSection, p, dof);CHKERRQ(ierr); 879 q = p; 880 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 881 while (parent != q && parent >= pStart && parent < pEnd) { 882 q = parent; 883 884 ierr = DMPlexGetTrueSupportSize(dm,q,&qdof,numTrueSupp);CHKERRQ(ierr); 885 ierr = PetscSectionAddDof(newSupportSection,p,qdof);CHKERRQ(ierr); 886 ierr = PetscSectionAddDof(newSupportSection,q,dof);CHKERRQ(ierr); 887 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 888 } 889 } 890 } 891 ierr = PetscSectionSetUp(newSupportSection);CHKERRQ(ierr); 892 ierr = PetscSectionGetStorageSize(newSupportSection,&newSize);CHKERRQ(ierr); 893 ierr = PetscMalloc1(newSize,&newSupports);CHKERRQ(ierr); 894 for (d = 0; d <= depth; d++) { 895 ierr = DMPlexGetHeightStratum(dm,d,&pStart,&pEnd);CHKERRQ(ierr); 896 for (p = pStart; p < pEnd; p++) { 897 PetscInt dof, off, q, qdof, qoff, newDof, newOff, newqOff, i, parent; 898 899 ierr = PetscSectionGetDof(mesh->supportSection, p, &dof);CHKERRQ(ierr); 900 ierr = PetscSectionGetOffset(mesh->supportSection, p, &off);CHKERRQ(ierr); 901 ierr = PetscSectionGetDof(newSupportSection, p, &newDof);CHKERRQ(ierr); 902 ierr = PetscSectionGetOffset(newSupportSection, p, &newOff);CHKERRQ(ierr); 903 for (i = 0; i < dof; i++) { 904 PetscInt numCones, j; 905 const PetscInt *cone; 906 PetscInt q = mesh->supports[off + i]; 907 908 ierr = DMPlexGetConeSize(dm,q,&numCones);CHKERRQ(ierr); 909 ierr = DMPlexGetCone(dm,q,&cone);CHKERRQ(ierr); 910 for (j = 0; j < numCones; j++) { 911 if (cone[j] == p) break; 912 } 913 if (j < numCones) newSupports[newOff+offsets[p]++] = q; 914 } 915 mesh->maxSupportSize = PetscMax(mesh->maxSupportSize,newDof); 916 917 q = p; 918 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 919 while (parent != q && parent >= pStart && parent < pEnd) { 920 q = parent; 921 ierr = PetscSectionGetDof(mesh->supportSection, q, &qdof);CHKERRQ(ierr); 922 ierr = PetscSectionGetOffset(mesh->supportSection, q, &qoff);CHKERRQ(ierr); 923 ierr = PetscSectionGetOffset(newSupportSection, q, &newqOff);CHKERRQ(ierr); 924 for (i = 0; i < qdof; i++) { 925 PetscInt numCones, j; 926 const PetscInt *cone; 927 PetscInt r = mesh->supports[qoff + i]; 928 929 ierr = DMPlexGetConeSize(dm,r,&numCones);CHKERRQ(ierr); 930 ierr = DMPlexGetCone(dm,r,&cone);CHKERRQ(ierr); 931 for (j = 0; j < numCones; j++) { 932 if (cone[j] == q) break; 933 } 934 if (j < numCones) newSupports[newOff+offsets[p]++] = r; 935 } 936 for (i = 0; i < dof; i++) { 937 PetscInt numCones, j; 938 const PetscInt *cone; 939 PetscInt r = mesh->supports[off + i]; 940 941 ierr = DMPlexGetConeSize(dm,r,&numCones);CHKERRQ(ierr); 942 ierr = DMPlexGetCone(dm,r,&cone);CHKERRQ(ierr); 943 for (j = 0; j < numCones; j++) { 944 if (cone[j] == p) break; 945 } 946 if (j < numCones) newSupports[newqOff+offsets[q]++] = r; 947 } 948 ierr = DMPlexGetTreeParent(dm,q,&parent,NULL);CHKERRQ(ierr); 949 } 950 } 951 } 952 ierr = PetscSectionDestroy(&mesh->supportSection);CHKERRQ(ierr); 953 mesh->supportSection = newSupportSection; 954 ierr = PetscFree(mesh->supports);CHKERRQ(ierr); 955 mesh->supports = newSupports; 956 ierr = PetscFree(offsets);CHKERRQ(ierr); 957 ierr = PetscFree(numTrueSupp);CHKERRQ(ierr); 958 959 PetscFunctionReturn(0); 960 } 961 962 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM,PetscSection,PetscSection,Mat); 963 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM,PetscSection,PetscSection,Mat); 964 965 #undef __FUNCT__ 966 #define __FUNCT__ "DMPlexSetTree_Internal" 967 static PetscErrorCode DMPlexSetTree_Internal(DM dm, PetscSection parentSection, PetscInt *parents, PetscInt *childIDs, PetscBool computeCanonical, PetscBool exchangeSupports) 968 { 969 DM_Plex *mesh = (DM_Plex *)dm->data; 970 DM refTree; 971 PetscInt size; 972 PetscErrorCode ierr; 973 974 PetscFunctionBegin; 975 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 976 PetscValidHeaderSpecific(parentSection, PETSC_SECTION_CLASSID, 2); 977 ierr = PetscObjectReference((PetscObject)parentSection);CHKERRQ(ierr); 978 ierr = PetscSectionDestroy(&mesh->parentSection);CHKERRQ(ierr); 979 mesh->parentSection = parentSection; 980 ierr = PetscSectionGetStorageSize(parentSection,&size);CHKERRQ(ierr); 981 if (parents != mesh->parents) { 982 ierr = PetscFree(mesh->parents);CHKERRQ(ierr); 983 ierr = PetscMalloc1(size,&mesh->parents);CHKERRQ(ierr); 984 ierr = PetscMemcpy(mesh->parents, parents, size * sizeof(*parents));CHKERRQ(ierr); 985 } 986 if (childIDs != mesh->childIDs) { 987 ierr = PetscFree(mesh->childIDs);CHKERRQ(ierr); 988 ierr = PetscMalloc1(size,&mesh->childIDs);CHKERRQ(ierr); 989 ierr = PetscMemcpy(mesh->childIDs, childIDs, size * sizeof(*childIDs));CHKERRQ(ierr); 990 } 991 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr); 992 if (refTree) { 993 DMLabel canonLabel; 994 995 ierr = DMGetLabel(refTree,"canonical",&canonLabel);CHKERRQ(ierr); 996 if (canonLabel) { 997 PetscInt i; 998 999 for (i = 0; i < size; i++) { 1000 PetscInt canon; 1001 ierr = DMLabelGetValue(canonLabel, mesh->childIDs[i], &canon);CHKERRQ(ierr); 1002 if (canon >= 0) { 1003 mesh->childIDs[i] = canon; 1004 } 1005 } 1006 } 1007 mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_FromReference; 1008 } 1009 else { 1010 mesh->computeanchormatrix = DMPlexComputeAnchorMatrix_Tree_Direct; 1011 } 1012 ierr = DMPlexTreeSymmetrize(dm);CHKERRQ(ierr); 1013 if (computeCanonical) { 1014 PetscInt d, dim; 1015 1016 /* add the canonical label */ 1017 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 1018 ierr = DMCreateLabel(dm,"canonical");CHKERRQ(ierr); 1019 for (d = 0; d <= dim; d++) { 1020 PetscInt p, dStart, dEnd, canon = -1, cNumChildren; 1021 const PetscInt *cChildren; 1022 1023 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr); 1024 for (p = dStart; p < dEnd; p++) { 1025 ierr = DMPlexGetTreeChildren(dm,p,&cNumChildren,&cChildren);CHKERRQ(ierr); 1026 if (cNumChildren) { 1027 canon = p; 1028 break; 1029 } 1030 } 1031 if (canon == -1) continue; 1032 for (p = dStart; p < dEnd; p++) { 1033 PetscInt numChildren, i; 1034 const PetscInt *children; 1035 1036 ierr = DMPlexGetTreeChildren(dm,p,&numChildren,&children);CHKERRQ(ierr); 1037 if (numChildren) { 1038 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); 1039 ierr = DMSetLabelValue(dm,"canonical",p,canon);CHKERRQ(ierr); 1040 for (i = 0; i < numChildren; i++) { 1041 ierr = DMSetLabelValue(dm,"canonical",children[i],cChildren[i]);CHKERRQ(ierr); 1042 } 1043 } 1044 } 1045 } 1046 } 1047 if (exchangeSupports) { 1048 ierr = DMPlexTreeExchangeSupports(dm);CHKERRQ(ierr); 1049 } 1050 mesh->createanchors = DMPlexCreateAnchors_Tree; 1051 /* reset anchors */ 1052 ierr = DMPlexSetAnchors(dm,NULL,NULL);CHKERRQ(ierr); 1053 PetscFunctionReturn(0); 1054 } 1055 1056 #undef __FUNCT__ 1057 #define __FUNCT__ "DMPlexSetTree" 1058 /*@ 1059 DMPlexSetTree - set the tree that describes the hierarchy of non-conforming mesh points. This routine also creates 1060 the point-to-point constraints determined by the tree: a point is constained to the points in the closure of its 1061 tree root. 1062 1063 Collective on dm 1064 1065 Input Parameters: 1066 + dm - the DMPlex object 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; the reference count of parentSection is incremented 1069 . parents - a list of the point parents; copied, can be destroyed 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 childIDs; copied, can be destroyed 1072 1073 Level: intermediate 1074 1075 .seealso: DMPlexGetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 1076 @*/ 1077 PetscErrorCode DMPlexSetTree(DM dm, PetscSection parentSection, PetscInt parents[], PetscInt childIDs[]) 1078 { 1079 PetscErrorCode ierr; 1080 1081 PetscFunctionBegin; 1082 ierr = DMPlexSetTree_Internal(dm,parentSection,parents,childIDs,PETSC_FALSE,PETSC_TRUE);CHKERRQ(ierr); 1083 PetscFunctionReturn(0); 1084 } 1085 1086 #undef __FUNCT__ 1087 #define __FUNCT__ "DMPlexGetTree" 1088 /*@ 1089 DMPlexGetTree - get the tree that describes the hierarchy of non-conforming mesh points. 1090 Collective on dm 1091 1092 Input Parameters: 1093 . dm - the DMPlex object 1094 1095 Output Parameters: 1096 + parentSection - a section describing the tree: a point has a parent if it has 1 dof in the section; the section 1097 offset indexes the parent and childID list 1098 . parents - a list of the point parents 1099 . childIDs - identifies the relationship of the child point to the parent point; if there is a reference tree, then 1100 the child corresponds to the point in the reference tree with index childID 1101 . childSection - the inverse of the parent section 1102 - children - a list of the point children 1103 1104 Level: intermediate 1105 1106 .seealso: DMPlexSetTree(), DMPlexSetReferenceTree(), DMPlexSetAnchors(), DMPlexGetTreeParent(), DMPlexGetTreeChildren() 1107 @*/ 1108 PetscErrorCode DMPlexGetTree(DM dm, PetscSection *parentSection, PetscInt *parents[], PetscInt *childIDs[], PetscSection *childSection, PetscInt *children[]) 1109 { 1110 DM_Plex *mesh = (DM_Plex *)dm->data; 1111 1112 PetscFunctionBegin; 1113 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1114 if (parentSection) *parentSection = mesh->parentSection; 1115 if (parents) *parents = mesh->parents; 1116 if (childIDs) *childIDs = mesh->childIDs; 1117 if (childSection) *childSection = mesh->childSection; 1118 if (children) *children = mesh->children; 1119 PetscFunctionReturn(0); 1120 } 1121 1122 #undef __FUNCT__ 1123 #define __FUNCT__ "DMPlexGetTreeParent" 1124 /*@ 1125 DMPlexGetTreeParent - get the parent of a point in the tree describing the point hierarchy (not the Sieve DAG) 1126 1127 Input Parameters: 1128 + dm - the DMPlex object 1129 - point - the query point 1130 1131 Output Parameters: 1132 + parent - if not NULL, set to the parent of the point, or the point itself if the point does not have a parent 1133 - childID - if not NULL, set to the child ID of the point with respect to its parent, or 0 if the point 1134 does not have a parent 1135 1136 Level: intermediate 1137 1138 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeChildren() 1139 @*/ 1140 PetscErrorCode DMPlexGetTreeParent(DM dm, PetscInt point, PetscInt *parent, PetscInt *childID) 1141 { 1142 DM_Plex *mesh = (DM_Plex *)dm->data; 1143 PetscSection pSec; 1144 PetscErrorCode ierr; 1145 1146 PetscFunctionBegin; 1147 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1148 pSec = mesh->parentSection; 1149 if (pSec && point >= pSec->pStart && point < pSec->pEnd) { 1150 PetscInt dof; 1151 1152 ierr = PetscSectionGetDof (pSec, point, &dof);CHKERRQ(ierr); 1153 if (dof) { 1154 PetscInt off; 1155 1156 ierr = PetscSectionGetOffset (pSec, point, &off);CHKERRQ(ierr); 1157 if (parent) *parent = mesh->parents[off]; 1158 if (childID) *childID = mesh->childIDs[off]; 1159 PetscFunctionReturn(0); 1160 } 1161 } 1162 if (parent) { 1163 *parent = point; 1164 } 1165 if (childID) { 1166 *childID = 0; 1167 } 1168 PetscFunctionReturn(0); 1169 } 1170 1171 #undef __FUNCT__ 1172 #define __FUNCT__ "DMPlexGetTreeChildren" 1173 /*@C 1174 DMPlexGetTreeChildren - get the children of a point in the tree describing the point hierarchy (not the Sieve DAG) 1175 1176 Input Parameters: 1177 + dm - the DMPlex object 1178 - point - the query point 1179 1180 Output Parameters: 1181 + numChildren - if not NULL, set to the number of children 1182 - children - if not NULL, set to a list children, or set to NULL if the point has no children 1183 1184 Level: intermediate 1185 1186 Fortran Notes: 1187 Since it returns an array, this routine is only available in Fortran 90, and you must 1188 include petsc.h90 in your code. 1189 1190 .seealso: DMPlexSetTree(), DMPlexGetTree(), DMPlexGetTreeParent() 1191 @*/ 1192 PetscErrorCode DMPlexGetTreeChildren(DM dm, PetscInt point, PetscInt *numChildren, const PetscInt *children[]) 1193 { 1194 DM_Plex *mesh = (DM_Plex *)dm->data; 1195 PetscSection childSec; 1196 PetscInt dof = 0; 1197 PetscErrorCode ierr; 1198 1199 PetscFunctionBegin; 1200 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1201 childSec = mesh->childSection; 1202 if (childSec && point >= childSec->pStart && point < childSec->pEnd) { 1203 ierr = PetscSectionGetDof (childSec, point, &dof);CHKERRQ(ierr); 1204 } 1205 if (numChildren) *numChildren = dof; 1206 if (children) { 1207 if (dof) { 1208 PetscInt off; 1209 1210 ierr = PetscSectionGetOffset (childSec, point, &off);CHKERRQ(ierr); 1211 *children = &mesh->children[off]; 1212 } 1213 else { 1214 *children = NULL; 1215 } 1216 } 1217 PetscFunctionReturn(0); 1218 } 1219 1220 #undef __FUNCT__ 1221 #define __FUNCT__ "DMPlexComputeAnchorMatrix_Tree_Direct" 1222 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat) 1223 { 1224 PetscDS ds; 1225 PetscInt spdim; 1226 PetscInt numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd; 1227 const PetscInt *anchors; 1228 PetscSection aSec; 1229 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent; 1230 IS aIS; 1231 PetscErrorCode ierr; 1232 1233 PetscFunctionBegin; 1234 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 1235 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 1236 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1237 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 1238 ierr = DMPlexGetAnchors(dm,&aSec,&aIS);CHKERRQ(ierr); 1239 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 1240 ierr = PetscSectionGetChart(cSec,&conStart,&conEnd);CHKERRQ(ierr); 1241 ierr = DMGetDimension(dm,&spdim);CHKERRQ(ierr); 1242 ierr = PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent);CHKERRQ(ierr); 1243 1244 for (f = 0; f < numFields; f++) { 1245 PetscObject disc; 1246 PetscFE fe = NULL; 1247 PetscFV fv = NULL; 1248 PetscClassId id; 1249 PetscDualSpace space; 1250 PetscInt i, j, k, nPoints, offset; 1251 PetscInt fSize, fComp; 1252 PetscReal *B = NULL; 1253 PetscReal *weights, *pointsRef, *pointsReal; 1254 Mat Amat, Bmat, Xmat; 1255 1256 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 1257 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 1258 if (id == PETSCFE_CLASSID) { 1259 fe = (PetscFE) disc; 1260 ierr = PetscFEGetDualSpace(fe,&space);CHKERRQ(ierr); 1261 ierr = PetscDualSpaceGetDimension(space,&fSize);CHKERRQ(ierr); 1262 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 1263 } 1264 else if (id == PETSCFV_CLASSID) { 1265 fv = (PetscFV) disc; 1266 ierr = PetscFVGetDualSpace(fv,&space);CHKERRQ(ierr); 1267 ierr = PetscDualSpaceGetDimension(space,&fSize);CHKERRQ(ierr); 1268 ierr = PetscFVGetNumComponents(fv,&fComp);CHKERRQ(ierr); 1269 } 1270 else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id); 1271 1272 ierr = MatCreate(PETSC_COMM_SELF,&Amat);CHKERRQ(ierr); 1273 ierr = MatSetSizes(Amat,fSize,fSize,fSize,fSize);CHKERRQ(ierr); 1274 ierr = MatSetType(Amat,MATSEQDENSE);CHKERRQ(ierr); 1275 ierr = MatSetUp(Amat);CHKERRQ(ierr); 1276 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);CHKERRQ(ierr); 1277 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);CHKERRQ(ierr); 1278 nPoints = 0; 1279 for (i = 0; i < fSize; i++) { 1280 PetscInt qPoints; 1281 PetscQuadrature quad; 1282 1283 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 1284 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 1285 nPoints += qPoints; 1286 } 1287 ierr = PetscMalloc3(nPoints,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal);CHKERRQ(ierr); 1288 offset = 0; 1289 for (i = 0; i < fSize; i++) { 1290 PetscInt qPoints; 1291 const PetscReal *p, *w; 1292 PetscQuadrature quad; 1293 1294 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 1295 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,&p,&w);CHKERRQ(ierr); 1296 ierr = PetscMemcpy(weights+offset,w,qPoints*sizeof(*w));CHKERRQ(ierr); 1297 ierr = PetscMemcpy(pointsRef+spdim*offset,p,spdim*qPoints*sizeof(*p));CHKERRQ(ierr); 1298 offset += qPoints; 1299 } 1300 if (id == PETSCFE_CLASSID) { 1301 ierr = PetscFEGetTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);CHKERRQ(ierr); 1302 } 1303 else { 1304 ierr = PetscFVGetTabulation(fv,nPoints,pointsRef,&B,NULL,NULL);CHKERRQ(ierr); 1305 } 1306 offset = 0; 1307 for (i = 0; i < fSize; i++) { 1308 PetscInt qPoints; 1309 PetscQuadrature quad; 1310 1311 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 1312 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 1313 for (j = 0; j < fSize; j++) { 1314 PetscScalar val = 0.; 1315 1316 for (k = 0; k < qPoints; k++) { 1317 val += B[((offset + k) * fSize + j) * fComp] * weights[k]; 1318 } 1319 ierr = MatSetValue(Amat,i,j,val,INSERT_VALUES);CHKERRQ(ierr); 1320 } 1321 offset += qPoints; 1322 } 1323 if (id == PETSCFE_CLASSID) { 1324 ierr = PetscFERestoreTabulation(fe,nPoints,pointsRef,&B,NULL,NULL);CHKERRQ(ierr); 1325 } 1326 else { 1327 ierr = PetscFVRestoreTabulation(fv,nPoints,pointsRef,&B,NULL,NULL);CHKERRQ(ierr); 1328 } 1329 ierr = MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1330 ierr = MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1331 ierr = MatLUFactor(Amat,NULL,NULL,NULL);CHKERRQ(ierr); 1332 for (c = cStart; c < cEnd; c++) { 1333 PetscInt parent; 1334 PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL; 1335 PetscInt *childOffsets, *parentOffsets; 1336 1337 ierr = DMPlexGetTreeParent(dm,c,&parent,NULL);CHKERRQ(ierr); 1338 if (parent == c) continue; 1339 ierr = DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1340 for (i = 0; i < closureSize; i++) { 1341 PetscInt p = closure[2*i]; 1342 PetscInt conDof; 1343 1344 if (p < conStart || p >= conEnd) continue; 1345 if (numFields > 1) { 1346 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 1347 } 1348 else { 1349 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 1350 } 1351 if (conDof) break; 1352 } 1353 if (i == closureSize) { 1354 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1355 continue; 1356 } 1357 1358 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ);CHKERRQ(ierr); 1359 ierr = DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent);CHKERRQ(ierr); 1360 for (i = 0; i < nPoints; i++) { 1361 CoordinatesRefToReal(spdim, spdim, v0, J, &pointsRef[i*spdim],vtmp);CHKERRQ(ierr); 1362 CoordinatesRealToRef(spdim, spdim, v0parent, invJparent, vtmp, &pointsReal[i*spdim]);CHKERRQ(ierr); 1363 } 1364 if (id == PETSCFE_CLASSID) { 1365 ierr = PetscFEGetTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);CHKERRQ(ierr); 1366 } 1367 else { 1368 ierr = PetscFVGetTabulation(fv,nPoints,pointsReal,&B,NULL,NULL);CHKERRQ(ierr); 1369 } 1370 offset = 0; 1371 for (i = 0; i < fSize; i++) { 1372 PetscInt qPoints; 1373 PetscQuadrature quad; 1374 1375 ierr = PetscDualSpaceGetFunctional(space,i,&quad);CHKERRQ(ierr); 1376 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 1377 for (j = 0; j < fSize; j++) { 1378 PetscScalar val = 0.; 1379 1380 for (k = 0; k < qPoints; k++) { 1381 val += B[((offset + k) * fSize + j) * fComp] * weights[k]; 1382 } 1383 MatSetValue(Bmat,i,j,val,INSERT_VALUES);CHKERRQ(ierr); 1384 } 1385 offset += qPoints; 1386 } 1387 if (id == PETSCFE_CLASSID) { 1388 ierr = PetscFERestoreTabulation(fe,nPoints,pointsReal,&B,NULL,NULL);CHKERRQ(ierr); 1389 } 1390 else { 1391 ierr = PetscFVRestoreTabulation(fv,nPoints,pointsReal,&B,NULL,NULL);CHKERRQ(ierr); 1392 } 1393 ierr = MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1394 ierr = MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1395 ierr = MatMatSolve(Amat,Bmat,Xmat);CHKERRQ(ierr); 1396 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 1397 ierr = PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);CHKERRQ(ierr); 1398 childOffsets[0] = 0; 1399 for (i = 0; i < closureSize; i++) { 1400 PetscInt p = closure[2*i]; 1401 PetscInt dof; 1402 1403 if (numFields > 1) { 1404 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 1405 } 1406 else { 1407 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 1408 } 1409 childOffsets[i+1]=childOffsets[i]+dof / fComp; 1410 } 1411 parentOffsets[0] = 0; 1412 for (i = 0; i < closureSizeP; i++) { 1413 PetscInt p = closureP[2*i]; 1414 PetscInt dof; 1415 1416 if (numFields > 1) { 1417 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 1418 } 1419 else { 1420 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 1421 } 1422 parentOffsets[i+1]=parentOffsets[i]+dof / fComp; 1423 } 1424 for (i = 0; i < closureSize; i++) { 1425 PetscInt conDof, conOff, aDof, aOff; 1426 PetscInt p = closure[2*i]; 1427 PetscInt o = closure[2*i+1]; 1428 1429 if (p < conStart || p >= conEnd) continue; 1430 if (numFields > 1) { 1431 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 1432 ierr = PetscSectionGetFieldOffset(cSec,p,f,&conOff);CHKERRQ(ierr); 1433 } 1434 else { 1435 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 1436 ierr = PetscSectionGetOffset(cSec,p,&conOff);CHKERRQ(ierr); 1437 } 1438 if (!conDof) continue; 1439 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 1440 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 1441 for (k = 0; k < aDof; k++) { 1442 PetscInt a = anchors[aOff + k]; 1443 PetscInt aSecDof, aSecOff; 1444 1445 if (numFields > 1) { 1446 ierr = PetscSectionGetFieldDof(section,a,f,&aSecDof);CHKERRQ(ierr); 1447 ierr = PetscSectionGetFieldOffset(section,a,f,&aSecOff);CHKERRQ(ierr); 1448 } 1449 else { 1450 ierr = PetscSectionGetDof(section,a,&aSecDof);CHKERRQ(ierr); 1451 ierr = PetscSectionGetOffset(section,a,&aSecOff);CHKERRQ(ierr); 1452 } 1453 if (!aSecDof) continue; 1454 1455 for (j = 0; j < closureSizeP; j++) { 1456 PetscInt q = closureP[2*j]; 1457 PetscInt oq = closureP[2*j+1]; 1458 1459 if (q == a) { 1460 PetscInt r, s, t; 1461 1462 for (r = childOffsets[i]; r < childOffsets[i+1]; r++) { 1463 for (s = parentOffsets[j]; s < parentOffsets[j+1]; s++) { 1464 PetscScalar val; 1465 PetscInt insertCol, insertRow; 1466 1467 ierr = MatGetValue(Xmat,r,s,&val);CHKERRQ(ierr); 1468 if (o >= 0) { 1469 insertRow = conOff + fComp * (r - childOffsets[i]); 1470 } 1471 else { 1472 insertRow = conOff + fComp * (childOffsets[i + 1] - 1 - r); 1473 } 1474 if (oq >= 0) { 1475 insertCol = aSecOff + fComp * (s - parentOffsets[j]); 1476 } 1477 else { 1478 insertCol = aSecOff + fComp * (parentOffsets[j + 1] - 1 - s); 1479 } 1480 for (t = 0; t < fComp; t++) { 1481 ierr = MatSetValue(cMat,insertRow + t,insertCol + t,val,INSERT_VALUES);CHKERRQ(ierr); 1482 } 1483 } 1484 } 1485 } 1486 } 1487 } 1488 } 1489 ierr = PetscFree2(childOffsets,parentOffsets);CHKERRQ(ierr); 1490 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1491 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 1492 } 1493 ierr = MatDestroy(&Amat);CHKERRQ(ierr); 1494 ierr = MatDestroy(&Bmat);CHKERRQ(ierr); 1495 ierr = MatDestroy(&Xmat);CHKERRQ(ierr); 1496 ierr = PetscFree3(weights,pointsRef,pointsReal);CHKERRQ(ierr); 1497 } 1498 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1499 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1500 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);CHKERRQ(ierr); 1501 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 1502 1503 PetscFunctionReturn(0); 1504 } 1505 1506 #undef __FUNCT__ 1507 #define __FUNCT__ "DMPlexReferenceTreeGetChildrenMatrices" 1508 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN) 1509 { 1510 Mat refCmat; 1511 PetscDS ds; 1512 PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN; 1513 PetscScalar ***refPointFieldMats; 1514 PetscSection refConSec, refAnSec, refSection; 1515 IS refAnIS; 1516 const PetscInt *refAnchors; 1517 PetscErrorCode ierr; 1518 1519 PetscFunctionBegin; 1520 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 1521 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1522 ierr = DMGetDefaultConstraints(refTree,&refConSec,&refCmat);CHKERRQ(ierr); 1523 ierr = DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);CHKERRQ(ierr); 1524 ierr = ISGetIndices(refAnIS,&refAnchors);CHKERRQ(ierr); 1525 ierr = DMGetDefaultSection(refTree,&refSection);CHKERRQ(ierr); 1526 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1527 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);CHKERRQ(ierr); 1528 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);CHKERRQ(ierr); 1529 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 1530 ierr = PetscSectionGetMaxDof(refAnSec,&maxAnDof);CHKERRQ(ierr); 1531 ierr = PetscMalloc1(maxDof,&rows);CHKERRQ(ierr); 1532 ierr = PetscMalloc1(maxDof*maxAnDof,&cols);CHKERRQ(ierr); 1533 for (p = pRefStart; p < pRefEnd; p++) { 1534 PetscInt parent, closureSize, *closure = NULL, pDof; 1535 1536 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 1537 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 1538 if (!pDof || parent == p) continue; 1539 1540 ierr = PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);CHKERRQ(ierr); 1541 ierr = PetscCalloc1(numFields,&refPointFieldN[p-pRefStart]);CHKERRQ(ierr); 1542 ierr = DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1543 for (f = 0; f < numFields; f++) { 1544 PetscInt cDof, cOff, numCols, r, i, fComp; 1545 PetscObject disc; 1546 PetscClassId id; 1547 PetscFE fe = NULL; 1548 PetscFV fv = NULL; 1549 1550 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 1551 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 1552 if (id == PETSCFE_CLASSID) { 1553 fe = (PetscFE) disc; 1554 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 1555 } 1556 else if (id == PETSCFV_CLASSID) { 1557 fv = (PetscFV) disc; 1558 ierr = PetscFVGetNumComponents(fv,&fComp);CHKERRQ(ierr); 1559 } 1560 else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id); 1561 1562 if (numFields > 1) { 1563 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 1564 ierr = PetscSectionGetFieldOffset(refConSec,p,f,&cOff);CHKERRQ(ierr); 1565 } 1566 else { 1567 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 1568 ierr = PetscSectionGetOffset(refConSec,p,&cOff);CHKERRQ(ierr); 1569 } 1570 1571 if (!cDof) continue; 1572 for (r = 0; r < cDof; r++) { 1573 rows[r] = cOff + r; 1574 } 1575 numCols = 0; 1576 for (i = 0; i < closureSize; i++) { 1577 PetscInt q = closure[2*i]; 1578 PetscInt o = closure[2*i+1]; 1579 PetscInt aDof, aOff, j; 1580 1581 if (numFields > 1) { 1582 ierr = PetscSectionGetFieldDof(refSection,q,f,&aDof);CHKERRQ(ierr); 1583 ierr = PetscSectionGetFieldOffset(refSection,q,f,&aOff);CHKERRQ(ierr); 1584 } 1585 else { 1586 ierr = PetscSectionGetDof(refSection,q,&aDof);CHKERRQ(ierr); 1587 ierr = PetscSectionGetOffset(refSection,q,&aOff);CHKERRQ(ierr); 1588 } 1589 1590 for (j = 0; j < aDof; j++) { 1591 PetscInt node = (o >= 0) ? (j / fComp) : ((aDof - 1 - j) / fComp); 1592 PetscInt comp = (j % fComp); 1593 1594 cols[numCols++] = aOff + node * fComp + comp; 1595 } 1596 } 1597 refPointFieldN[p-pRefStart][f] = numCols; 1598 ierr = PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 1599 ierr = MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 1600 } 1601 ierr = DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1602 } 1603 *childrenMats = refPointFieldMats; 1604 *childrenN = refPointFieldN; 1605 ierr = ISRestoreIndices(refAnIS,&refAnchors);CHKERRQ(ierr); 1606 ierr = PetscFree(rows);CHKERRQ(ierr); 1607 ierr = PetscFree(cols);CHKERRQ(ierr); 1608 PetscFunctionReturn(0); 1609 } 1610 1611 #undef __FUNCT__ 1612 #define __FUNCT__ "DMPlexReferenceTreeRestoreChildrenMatrices" 1613 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN) 1614 { 1615 PetscDS ds; 1616 PetscInt **refPointFieldN; 1617 PetscScalar ***refPointFieldMats; 1618 PetscInt numFields, pRefStart, pRefEnd, p, f; 1619 PetscSection refConSec; 1620 PetscErrorCode ierr; 1621 1622 PetscFunctionBegin; 1623 refPointFieldN = *childrenN; 1624 *childrenN = NULL; 1625 refPointFieldMats = *childrenMats; 1626 *childrenMats = NULL; 1627 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 1628 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1629 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 1630 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1631 for (p = pRefStart; p < pRefEnd; p++) { 1632 PetscInt parent, pDof; 1633 1634 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 1635 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 1636 if (!pDof || parent == p) continue; 1637 1638 for (f = 0; f < numFields; f++) { 1639 PetscInt cDof; 1640 1641 if (numFields > 1) { 1642 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 1643 } 1644 else { 1645 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 1646 } 1647 1648 if (!cDof) continue; 1649 ierr = PetscFree(refPointFieldMats[p - pRefStart][f]);CHKERRQ(ierr); 1650 } 1651 ierr = PetscFree(refPointFieldMats[p - pRefStart]);CHKERRQ(ierr); 1652 ierr = PetscFree(refPointFieldN[p - pRefStart]);CHKERRQ(ierr); 1653 } 1654 ierr = PetscFree(refPointFieldMats);CHKERRQ(ierr); 1655 ierr = PetscFree(refPointFieldN);CHKERRQ(ierr); 1656 PetscFunctionReturn(0); 1657 } 1658 1659 #undef __FUNCT__ 1660 #define __FUNCT__ "DMPlexComputeAnchorMatrix_Tree_FromReference" 1661 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat) 1662 { 1663 DM refTree; 1664 PetscDS ds; 1665 Mat refCmat; 1666 PetscInt numFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN; 1667 PetscScalar ***refPointFieldMats, *pointWork; 1668 PetscSection refConSec, refAnSec, anSec; 1669 IS refAnIS, anIS; 1670 const PetscInt *anchors; 1671 PetscErrorCode ierr; 1672 1673 PetscFunctionBegin; 1674 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1675 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 1676 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1677 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr); 1678 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 1679 ierr = DMGetDefaultConstraints(refTree,&refConSec,&refCmat);CHKERRQ(ierr); 1680 ierr = DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);CHKERRQ(ierr); 1681 ierr = DMPlexGetAnchors(dm,&anSec,&anIS);CHKERRQ(ierr); 1682 ierr = ISGetIndices(anIS,&anchors);CHKERRQ(ierr); 1683 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1684 ierr = PetscSectionGetChart(conSec,&conStart,&conEnd);CHKERRQ(ierr); 1685 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 1686 ierr = PetscSectionGetMaxDof(refAnSec,&maxAnDof);CHKERRQ(ierr); 1687 ierr = PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);CHKERRQ(ierr); 1688 1689 /* step 1: get submats for every constrained point in the reference tree */ 1690 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 1691 1692 /* step 2: compute the preorder */ 1693 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 1694 ierr = PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);CHKERRQ(ierr); 1695 for (p = pStart; p < pEnd; p++) { 1696 perm[p - pStart] = p; 1697 iperm[p - pStart] = p-pStart; 1698 } 1699 for (p = 0; p < pEnd - pStart;) { 1700 PetscInt point = perm[p]; 1701 PetscInt parent; 1702 1703 ierr = DMPlexGetTreeParent(dm,point,&parent,NULL);CHKERRQ(ierr); 1704 if (parent == point) { 1705 p++; 1706 } 1707 else { 1708 PetscInt size, closureSize, *closure = NULL, i; 1709 1710 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1711 for (i = 0; i < closureSize; i++) { 1712 PetscInt q = closure[2*i]; 1713 if (iperm[q-pStart] > iperm[point-pStart]) { 1714 /* swap */ 1715 perm[p] = q; 1716 perm[iperm[q-pStart]] = point; 1717 iperm[point-pStart] = iperm[q-pStart]; 1718 iperm[q-pStart] = p; 1719 break; 1720 } 1721 } 1722 size = closureSize; 1723 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1724 if (i == size) { 1725 p++; 1726 } 1727 } 1728 } 1729 1730 /* step 3: fill the constraint matrix */ 1731 /* we are going to use a preorder progressive fill strategy. Mat doesn't 1732 * allow progressive fill without assembly, so we are going to set up the 1733 * values outside of the Mat first. 1734 */ 1735 { 1736 PetscInt nRows, row, nnz; 1737 PetscBool done; 1738 const PetscInt *ia, *ja; 1739 PetscScalar *vals; 1740 1741 ierr = MatGetRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);CHKERRQ(ierr); 1742 if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not get RowIJ of constraint matrix"); 1743 nnz = ia[nRows]; 1744 /* malloc and then zero rows right before we fill them: this way valgrind 1745 * can tell if we are doing progressive fill in the wrong order */ 1746 ierr = PetscMalloc1(nnz,&vals);CHKERRQ(ierr); 1747 for (p = 0; p < pEnd - pStart; p++) { 1748 PetscInt parent, childid, closureSize, *closure = NULL; 1749 PetscInt point = perm[p], pointDof; 1750 1751 ierr = DMPlexGetTreeParent(dm,point,&parent,&childid);CHKERRQ(ierr); 1752 if ((point < conStart) || (point >= conEnd) || (parent == point)) continue; 1753 ierr = PetscSectionGetDof(conSec,point,&pointDof);CHKERRQ(ierr); 1754 if (!pointDof) continue; 1755 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1756 for (f = 0; f < numFields; f++) { 1757 PetscInt cDof, cOff, numCols, numFillCols, i, r, fComp = -1, matOffset, offset; 1758 PetscScalar *pointMat; 1759 PetscObject disc; 1760 PetscClassId id; 1761 PetscFE fe = NULL; 1762 PetscFV fv = NULL; 1763 1764 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 1765 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 1766 if (id == PETSCFE_CLASSID) { 1767 fe = (PetscFE) disc; 1768 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 1769 } 1770 else if (id == PETSCFV_CLASSID) { 1771 fv = (PetscFV) disc; 1772 ierr = PetscFVGetNumComponents(fv,&fComp);CHKERRQ(ierr); 1773 } 1774 else { 1775 SETERRQ(PetscObjectComm((PetscObject)disc),PETSC_ERR_SUP,"Unsupported discretization"); 1776 } 1777 1778 if (numFields > 1) { 1779 ierr = PetscSectionGetFieldDof(conSec,point,f,&cDof);CHKERRQ(ierr); 1780 ierr = PetscSectionGetFieldOffset(conSec,point,f,&cOff);CHKERRQ(ierr); 1781 } 1782 else { 1783 ierr = PetscSectionGetDof(conSec,point,&cDof);CHKERRQ(ierr); 1784 ierr = PetscSectionGetOffset(conSec,point,&cOff);CHKERRQ(ierr); 1785 } 1786 if (!cDof) continue; 1787 1788 /* make sure that every row for this point is the same size */ 1789 #if defined(PETSC_USE_DEBUG) 1790 for (r = 0; r < cDof; r++) { 1791 if (cDof > 1 && r) { 1792 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])); 1793 } 1794 } 1795 #endif 1796 /* zero rows */ 1797 for (i = ia[cOff] ; i< ia[cOff+cDof];i++) { 1798 vals[i] = 0.; 1799 } 1800 matOffset = ia[cOff]; 1801 numFillCols = ia[cOff+1] - matOffset; 1802 pointMat = refPointFieldMats[childid-pRefStart][f]; 1803 numCols = refPointFieldN[childid-pRefStart][f]; 1804 offset = 0; 1805 for (i = 0; i < closureSize; i++) { 1806 PetscInt q = closure[2*i]; 1807 PetscInt o = closure[2*i+1]; 1808 PetscInt aDof, aOff, j, k, qConDof, qConOff; 1809 1810 qConDof = qConOff = 0; 1811 if (numFields > 1) { 1812 ierr = PetscSectionGetFieldDof(section,q,f,&aDof);CHKERRQ(ierr); 1813 ierr = PetscSectionGetFieldOffset(section,q,f,&aOff);CHKERRQ(ierr); 1814 if (q >= conStart && q < conEnd) { 1815 ierr = PetscSectionGetFieldDof(conSec,q,f,&qConDof);CHKERRQ(ierr); 1816 ierr = PetscSectionGetFieldOffset(conSec,q,f,&qConOff);CHKERRQ(ierr); 1817 } 1818 } 1819 else { 1820 ierr = PetscSectionGetDof(section,q,&aDof);CHKERRQ(ierr); 1821 ierr = PetscSectionGetOffset(section,q,&aOff);CHKERRQ(ierr); 1822 if (q >= conStart && q < conEnd) { 1823 ierr = PetscSectionGetDof(conSec,q,&qConDof);CHKERRQ(ierr); 1824 ierr = PetscSectionGetOffset(conSec,q,&qConOff);CHKERRQ(ierr); 1825 } 1826 } 1827 if (!aDof) continue; 1828 if (qConDof) { 1829 /* this point has anchors: its rows of the matrix should already 1830 * be filled, thanks to preordering */ 1831 /* first multiply into pointWork, then set in matrix */ 1832 PetscInt aMatOffset = ia[qConOff]; 1833 PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset; 1834 for (r = 0; r < cDof; r++) { 1835 for (j = 0; j < aNumFillCols; j++) { 1836 PetscScalar inVal = 0; 1837 for (k = 0; k < aDof; k++) { 1838 PetscInt node = (o >= 0) ? (k / fComp) : ((aDof - 1 - k) / fComp); 1839 PetscInt comp = (k % fComp); 1840 PetscInt col = node * fComp + comp; 1841 1842 inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j]; 1843 } 1844 pointWork[r * aNumFillCols + j] = inVal; 1845 } 1846 } 1847 /* assume that the columns are sorted, spend less time searching */ 1848 for (j = 0, k = 0; j < aNumFillCols; j++) { 1849 PetscInt col = ja[aMatOffset + j]; 1850 for (;k < numFillCols; k++) { 1851 if (ja[matOffset + k] == col) { 1852 break; 1853 } 1854 } 1855 if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, col); 1856 for (r = 0; r < cDof; r++) { 1857 vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j]; 1858 } 1859 } 1860 } 1861 else { 1862 /* find where to put this portion of pointMat into the matrix */ 1863 for (k = 0; k < numFillCols; k++) { 1864 if (ja[matOffset + k] == aOff) { 1865 break; 1866 } 1867 } 1868 if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, aOff); 1869 for (j = 0; j < aDof; j++) { 1870 PetscInt node = (o >= 0) ? (j / fComp) : ((aDof - 1 - j) / fComp); 1871 PetscInt comp = (j % fComp); 1872 PetscInt col = node * fComp + comp; 1873 for (r = 0; r < cDof; r++) { 1874 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + col]; 1875 } 1876 } 1877 } 1878 offset += aDof; 1879 } 1880 } 1881 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1882 } 1883 for (row = 0; row < nRows; row++) { 1884 ierr = MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);CHKERRQ(ierr); 1885 } 1886 ierr = MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);CHKERRQ(ierr); 1887 if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix"); 1888 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1889 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1890 ierr = PetscFree(vals);CHKERRQ(ierr); 1891 } 1892 1893 /* clean up */ 1894 ierr = ISRestoreIndices(anIS,&anchors);CHKERRQ(ierr); 1895 ierr = PetscFree2(perm,iperm);CHKERRQ(ierr); 1896 ierr = PetscFree(pointWork);CHKERRQ(ierr); 1897 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 1898 PetscFunctionReturn(0); 1899 } 1900 1901 #undef __FUNCT__ 1902 #define __FUNCT__ "DMPlexTreeRefineCell" 1903 /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of 1904 * a non-conforming mesh. Local refinement comes later */ 1905 PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm) 1906 { 1907 DM K; 1908 PetscMPIInt rank; 1909 PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd; 1910 PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations; 1911 PetscInt *Kembedding; 1912 PetscInt *cellClosure=NULL, nc; 1913 PetscScalar *newVertexCoords; 1914 PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset; 1915 PetscSection parentSection; 1916 PetscErrorCode ierr; 1917 1918 PetscFunctionBegin; 1919 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 1920 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 1921 ierr = DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);CHKERRQ(ierr); 1922 ierr = DMSetDimension(*ncdm,dim);CHKERRQ(ierr); 1923 1924 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1925 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection);CHKERRQ(ierr); 1926 ierr = DMPlexGetReferenceTree(dm,&K);CHKERRQ(ierr); 1927 if (!rank) { 1928 /* compute the new charts */ 1929 ierr = PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd);CHKERRQ(ierr); 1930 offset = 0; 1931 for (d = 0; d <= dim; d++) { 1932 PetscInt pOldCount, kStart, kEnd, k; 1933 1934 pNewStart[d] = offset; 1935 ierr = DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d]);CHKERRQ(ierr); 1936 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr); 1937 pOldCount = pOldEnd[d] - pOldStart[d]; 1938 /* adding the new points */ 1939 pNewCount[d] = pOldCount + kEnd - kStart; 1940 if (!d) { 1941 /* removing the cell */ 1942 pNewCount[d]--; 1943 } 1944 for (k = kStart; k < kEnd; k++) { 1945 PetscInt parent; 1946 ierr = DMPlexGetTreeParent(K,k,&parent,NULL);CHKERRQ(ierr); 1947 if (parent == k) { 1948 /* avoid double counting points that won't actually be new */ 1949 pNewCount[d]--; 1950 } 1951 } 1952 pNewEnd[d] = pNewStart[d] + pNewCount[d]; 1953 offset = pNewEnd[d]; 1954 1955 } 1956 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]); 1957 /* get the current closure of the cell that we are removing */ 1958 ierr = DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);CHKERRQ(ierr); 1959 1960 ierr = PetscMalloc1(pNewEnd[dim],&newConeSizes);CHKERRQ(ierr); 1961 { 1962 PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j; 1963 1964 ierr = DMPlexGetChart(K,&kStart,&kEnd);CHKERRQ(ierr); 1965 ierr = PetscMalloc4(kEnd-kStart,&Kembedding,kEnd-kStart,&perm,kEnd-kStart,&iperm,kEnd-kStart,&preOrient);CHKERRQ(ierr); 1966 1967 for (k = kStart; k < kEnd; k++) { 1968 perm[k - kStart] = k; 1969 iperm [k - kStart] = k - kStart; 1970 preOrient[k - kStart] = 0; 1971 } 1972 1973 ierr = DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);CHKERRQ(ierr); 1974 for (j = 1; j < closureSizeK; j++) { 1975 PetscInt parentOrientA = closureK[2*j+1]; 1976 PetscInt parentOrientB = cellClosure[2*j+1]; 1977 PetscInt p, q; 1978 1979 p = closureK[2*j]; 1980 q = cellClosure[2*j]; 1981 for (d = 0; d <= dim; d++) { 1982 if (q >= pOldStart[d] && q < pOldEnd[d]) { 1983 Kembedding[p] = (q - pOldStart[d]) + pNewStart[d]; 1984 } 1985 } 1986 if (parentOrientA != parentOrientB) { 1987 PetscInt numChildren, i; 1988 const PetscInt *children; 1989 1990 ierr = DMPlexGetTreeChildren(K,p,&numChildren,&children);CHKERRQ(ierr); 1991 for (i = 0; i < numChildren; i++) { 1992 PetscInt kPerm, oPerm; 1993 1994 k = children[i]; 1995 ierr = DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm);CHKERRQ(ierr); 1996 /* perm = what refTree position I'm in */ 1997 perm[kPerm-kStart] = k; 1998 /* iperm = who is at this position */ 1999 iperm[k-kStart] = kPerm-kStart; 2000 preOrient[kPerm-kStart] = oPerm; 2001 } 2002 } 2003 } 2004 ierr = DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);CHKERRQ(ierr); 2005 } 2006 ierr = PetscSectionSetChart(parentSection,0,pNewEnd[dim]);CHKERRQ(ierr); 2007 offset = 0; 2008 numNewCones = 0; 2009 for (d = 0; d <= dim; d++) { 2010 PetscInt kStart, kEnd, k; 2011 PetscInt p; 2012 PetscInt size; 2013 2014 for (p = pOldStart[d]; p < pOldEnd[d]; p++) { 2015 /* skip cell 0 */ 2016 if (p == cell) continue; 2017 /* old cones to new cones */ 2018 ierr = DMPlexGetConeSize(dm,p,&size);CHKERRQ(ierr); 2019 newConeSizes[offset++] = size; 2020 numNewCones += size; 2021 } 2022 2023 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr); 2024 for (k = kStart; k < kEnd; k++) { 2025 PetscInt kParent; 2026 2027 ierr = DMPlexGetTreeParent(K,k,&kParent,NULL);CHKERRQ(ierr); 2028 if (kParent != k) { 2029 Kembedding[k] = offset; 2030 ierr = DMPlexGetConeSize(K,k,&size);CHKERRQ(ierr); 2031 newConeSizes[offset++] = size; 2032 numNewCones += size; 2033 if (kParent != 0) { 2034 ierr = PetscSectionSetDof(parentSection,Kembedding[k],1);CHKERRQ(ierr); 2035 } 2036 } 2037 } 2038 } 2039 2040 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 2041 ierr = PetscSectionGetStorageSize(parentSection,&numPointsWithParents);CHKERRQ(ierr); 2042 ierr = PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations);CHKERRQ(ierr); 2043 ierr = PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs);CHKERRQ(ierr); 2044 2045 /* fill new cones */ 2046 offset = 0; 2047 for (d = 0; d <= dim; d++) { 2048 PetscInt kStart, kEnd, k, l; 2049 PetscInt p; 2050 PetscInt size; 2051 const PetscInt *cone, *orientation; 2052 2053 for (p = pOldStart[d]; p < pOldEnd[d]; p++) { 2054 /* skip cell 0 */ 2055 if (p == cell) continue; 2056 /* old cones to new cones */ 2057 ierr = DMPlexGetConeSize(dm,p,&size);CHKERRQ(ierr); 2058 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 2059 ierr = DMPlexGetConeOrientation(dm,p,&orientation);CHKERRQ(ierr); 2060 for (l = 0; l < size; l++) { 2061 newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1]; 2062 newOrientations[offset++] = orientation[l]; 2063 } 2064 } 2065 2066 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr); 2067 for (k = kStart; k < kEnd; k++) { 2068 PetscInt kPerm = perm[k], kParent; 2069 PetscInt preO = preOrient[k]; 2070 2071 ierr = DMPlexGetTreeParent(K,k,&kParent,NULL);CHKERRQ(ierr); 2072 if (kParent != k) { 2073 /* embed new cones */ 2074 ierr = DMPlexGetConeSize(K,k,&size);CHKERRQ(ierr); 2075 ierr = DMPlexGetCone(K,kPerm,&cone);CHKERRQ(ierr); 2076 ierr = DMPlexGetConeOrientation(K,kPerm,&orientation);CHKERRQ(ierr); 2077 for (l = 0; l < size; l++) { 2078 PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size); 2079 PetscInt newO, lSize, oTrue; 2080 2081 q = iperm[cone[m]]; 2082 newCones[offset] = Kembedding[q]; 2083 ierr = DMPlexGetConeSize(K,q,&lSize);CHKERRQ(ierr); 2084 oTrue = orientation[m]; 2085 oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2); 2086 newO = DihedralCompose(lSize,oTrue,preOrient[q]); 2087 newOrientations[offset++] = newO; 2088 } 2089 if (kParent != 0) { 2090 PetscInt newPoint = Kembedding[kParent]; 2091 ierr = PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset);CHKERRQ(ierr); 2092 parents[pOffset] = newPoint; 2093 childIDs[pOffset] = k; 2094 } 2095 } 2096 } 2097 } 2098 2099 ierr = PetscMalloc1(dim*(pNewEnd[dim]-pNewStart[dim]),&newVertexCoords);CHKERRQ(ierr); 2100 2101 /* fill coordinates */ 2102 offset = 0; 2103 { 2104 PetscInt kStart, kEnd, l; 2105 PetscSection vSection; 2106 PetscInt v; 2107 Vec coords; 2108 PetscScalar *coordvals; 2109 PetscInt dof, off; 2110 PetscReal v0[3], J[9], detJ; 2111 2112 #if defined(PETSC_USE_DEBUG) 2113 { 2114 PetscInt k; 2115 ierr = DMPlexGetHeightStratum(K,0,&kStart,&kEnd);CHKERRQ(ierr); 2116 for (k = kStart; k < kEnd; k++) { 2117 ierr = DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ);CHKERRQ(ierr); 2118 if (detJ <= 0.) SETERRQ1 (PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %d has bad determinant",k); 2119 } 2120 } 2121 #endif 2122 ierr = DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ);CHKERRQ(ierr); 2123 ierr = DMGetCoordinateSection(dm,&vSection);CHKERRQ(ierr); 2124 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2125 ierr = VecGetArray(coords,&coordvals);CHKERRQ(ierr); 2126 for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) { 2127 2128 ierr = PetscSectionGetDof(vSection,v,&dof);CHKERRQ(ierr); 2129 ierr = PetscSectionGetOffset(vSection,v,&off);CHKERRQ(ierr); 2130 for (l = 0; l < dof; l++) { 2131 newVertexCoords[offset++] = coordvals[off + l]; 2132 } 2133 } 2134 ierr = VecRestoreArray(coords,&coordvals);CHKERRQ(ierr); 2135 2136 ierr = DMGetCoordinateSection(K,&vSection);CHKERRQ(ierr); 2137 ierr = DMGetCoordinatesLocal(K,&coords);CHKERRQ(ierr); 2138 ierr = VecGetArray(coords,&coordvals);CHKERRQ(ierr); 2139 ierr = DMPlexGetDepthStratum(K,0,&kStart,&kEnd);CHKERRQ(ierr); 2140 for (v = kStart; v < kEnd; v++) { 2141 PetscReal coord[3], newCoord[3]; 2142 PetscInt vPerm = perm[v]; 2143 PetscInt kParent; 2144 2145 ierr = DMPlexGetTreeParent(K,v,&kParent,NULL);CHKERRQ(ierr); 2146 if (kParent != v) { 2147 /* this is a new vertex */ 2148 ierr = PetscSectionGetOffset(vSection,vPerm,&off);CHKERRQ(ierr); 2149 for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]); 2150 CoordinatesRefToReal(dim, dim, v0, J, coord, newCoord);CHKERRQ(ierr); 2151 for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l]; 2152 offset += dim; 2153 } 2154 } 2155 ierr = VecRestoreArray(coords,&coordvals);CHKERRQ(ierr); 2156 } 2157 2158 /* need to reverse the order of pNewCount: vertices first, cells last */ 2159 for (d = 0; d < (dim + 1) / 2; d++) { 2160 PetscInt tmp; 2161 2162 tmp = pNewCount[d]; 2163 pNewCount[d] = pNewCount[dim - d]; 2164 pNewCount[dim - d] = tmp; 2165 } 2166 2167 ierr = DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);CHKERRQ(ierr); 2168 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr); 2169 ierr = DMPlexSetTree(*ncdm,parentSection,parents,childIDs);CHKERRQ(ierr); 2170 2171 /* clean up */ 2172 ierr = DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);CHKERRQ(ierr); 2173 ierr = PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);CHKERRQ(ierr); 2174 ierr = PetscFree(newConeSizes);CHKERRQ(ierr); 2175 ierr = PetscFree2(newCones,newOrientations);CHKERRQ(ierr); 2176 ierr = PetscFree(newVertexCoords);CHKERRQ(ierr); 2177 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 2178 ierr = PetscFree4(Kembedding,perm,iperm,preOrient);CHKERRQ(ierr); 2179 } 2180 else { 2181 PetscInt p, counts[4]; 2182 PetscInt *coneSizes, *cones, *orientations; 2183 Vec coordVec; 2184 PetscScalar *coords; 2185 2186 for (d = 0; d <= dim; d++) { 2187 PetscInt dStart, dEnd; 2188 2189 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr); 2190 counts[d] = dEnd - dStart; 2191 } 2192 ierr = PetscMalloc1(pEnd-pStart,&coneSizes);CHKERRQ(ierr); 2193 for (p = pStart; p < pEnd; p++) { 2194 ierr = DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);CHKERRQ(ierr); 2195 } 2196 ierr = DMPlexGetCones(dm, &cones);CHKERRQ(ierr); 2197 ierr = DMPlexGetConeOrientations(dm, &orientations);CHKERRQ(ierr); 2198 ierr = DMGetCoordinatesLocal(dm,&coordVec);CHKERRQ(ierr); 2199 ierr = VecGetArray(coordVec,&coords);CHKERRQ(ierr); 2200 2201 ierr = PetscSectionSetChart(parentSection,pStart,pEnd);CHKERRQ(ierr); 2202 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 2203 ierr = DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);CHKERRQ(ierr); 2204 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr); 2205 ierr = DMPlexSetTree(*ncdm,parentSection,NULL,NULL);CHKERRQ(ierr); 2206 ierr = VecRestoreArray(coordVec,&coords);CHKERRQ(ierr); 2207 } 2208 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 2209 2210 PetscFunctionReturn(0); 2211 } 2212 2213 PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt[],PetscBool,PetscInt,PetscInt[]); 2214 PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt[]); 2215 2216 #undef __FUNCT__ 2217 #define __FUNCT__ "DMPlexComputeInterpolatorTree" 2218 PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 2219 { 2220 PetscSF coarseToFineEmbedded; 2221 PetscSection globalCoarse, globalFine; 2222 PetscSection localCoarse, localFine; 2223 PetscSection aSec, cSec; 2224 PetscSection rootIndicesSec, rootMatricesSec; 2225 PetscSection leafIndicesSec, leafMatricesSec; 2226 PetscInt *rootIndices, *leafIndices; 2227 PetscScalar *rootMatrices, *leafMatrices; 2228 IS aIS; 2229 const PetscInt *anchors; 2230 Mat cMat; 2231 PetscInt numFields; 2232 PetscInt pStartC, pEndC, pStartF, pEndF, p; 2233 PetscInt aStart, aEnd, cStart, cEnd; 2234 PetscInt *maxChildIds; 2235 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 2236 PetscErrorCode ierr; 2237 2238 PetscFunctionBegin; 2239 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 2240 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 2241 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 2242 { /* winnow fine points that don't have global dofs out of the sf */ 2243 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 2244 2245 for (p = pStartF, numPointsWithDofs = 0; p < pEndF; p++) { 2246 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 2247 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 2248 if ((dof - cdof) > 0) { 2249 numPointsWithDofs++; 2250 } 2251 } 2252 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 2253 for (p = pStartF, offset = 0; p < pEndF; p++) { 2254 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 2255 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 2256 if ((dof - cdof) > 0) { 2257 pointsWithDofs[offset++] = p - pStartF; 2258 } 2259 } 2260 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 2261 } 2262 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 2263 ierr = PetscMalloc1(pEndC-pStartC,&maxChildIds);CHKERRQ(ierr); 2264 for (p = pStartC; p < pEndC; p++) { 2265 maxChildIds[p - pStartC] = -1; 2266 } 2267 ierr = PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 2268 ierr = PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 2269 2270 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 2271 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 2272 2273 ierr = DMPlexGetAnchors(coarse,&aSec,&aIS);CHKERRQ(ierr); 2274 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 2275 ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr); 2276 2277 ierr = DMGetDefaultConstraints(coarse,&cSec,&cMat);CHKERRQ(ierr); 2278 ierr = PetscSectionGetChart(cSec,&cStart,&cEnd);CHKERRQ(ierr); 2279 2280 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 2281 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);CHKERRQ(ierr); 2282 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec);CHKERRQ(ierr); 2283 ierr = PetscSectionSetChart(rootIndicesSec,pStartC,pEndC);CHKERRQ(ierr); 2284 ierr = PetscSectionSetChart(rootMatricesSec,pStartC,pEndC);CHKERRQ(ierr); 2285 ierr = PetscSectionGetNumFields(globalCoarse,&numFields);CHKERRQ(ierr); 2286 { 2287 PetscInt maxFields = PetscMax(1,numFields) + 1; 2288 ierr = PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);CHKERRQ(ierr); 2289 } 2290 2291 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 2292 PetscInt dof, matSize; 2293 PetscInt aDof = 0; 2294 PetscInt cDof = 0; 2295 PetscInt maxChildId = maxChildIds[p - pStartC]; 2296 PetscInt numRowIndices = 0; 2297 PetscInt numColIndices = 0; 2298 2299 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 2300 if (dof < 0) { 2301 dof = -(dof + 1); 2302 } 2303 if (p >= aStart && p < aEnd) { 2304 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 2305 } 2306 if (p >= cStart && p < cEnd) { 2307 ierr = PetscSectionGetDof(cSec,p,&cDof);CHKERRQ(ierr); 2308 } 2309 offsets[0] = 0; 2310 newOffsets[0] = 0; 2311 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 2312 PetscInt *closure = NULL, closureSize, cl, f; 2313 2314 ierr = DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2315 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 2316 PetscInt c = closure[2 * cl], clDof; 2317 2318 ierr = PetscSectionGetDof(localCoarse,c,&clDof);CHKERRQ(ierr); 2319 numRowIndices += clDof; 2320 for (f = 0; f < numFields; f++) { 2321 ierr = PetscSectionGetFieldDof(localCoarse,c,f,&clDof);CHKERRQ(ierr); 2322 offsets[f + 1] += clDof; 2323 } 2324 } 2325 for (f = 0; f < numFields; f++) { 2326 offsets[f + 1] += offsets[f]; 2327 newOffsets[f + 1] = offsets[f + 1]; 2328 } 2329 /* get the number of indices needed and their field offsets */ 2330 ierr = DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE);CHKERRQ(ierr); 2331 ierr = DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2332 if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */ 2333 numColIndices = numRowIndices; 2334 matSize = 0; 2335 } 2336 else if (numFields) { /* we send one submat for each field: sum their sizes */ 2337 matSize = 0; 2338 for (f = 0; f < numFields; f++) { 2339 PetscInt numRow, numCol; 2340 2341 numRow = offsets[f + 1] - offsets[f]; 2342 numCol = newOffsets[f + 1] - newOffsets[f + 1]; 2343 matSize += numRow * numCol; 2344 } 2345 } 2346 else { 2347 matSize = numRowIndices * numColIndices; 2348 } 2349 } 2350 else if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */ 2351 PetscInt aOff, a, f; 2352 2353 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 2354 for (f = 0; f < numFields; f++) { 2355 PetscInt fDof; 2356 2357 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 2358 offsets[f+1] = fDof; 2359 } 2360 for (a = 0; a < aDof; a++) { 2361 PetscInt anchor = anchors[a + aOff], aLocalDof; 2362 2363 ierr = PetscSectionGetDof(localCoarse,anchor,&aLocalDof);CHKERRQ(ierr); 2364 numColIndices += aLocalDof; 2365 for (f = 0; f < numFields; f++) { 2366 PetscInt fDof; 2367 2368 ierr = PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);CHKERRQ(ierr); 2369 newOffsets[f+1] += fDof; 2370 } 2371 } 2372 if (numFields) { 2373 matSize = 0; 2374 for (f = 0; f < numFields; f++) { 2375 matSize += offsets[f+1] * newOffsets[f+1]; 2376 } 2377 } 2378 else { 2379 matSize = numColIndices * dof; 2380 } 2381 } 2382 else { /* no children, and no constraints on dofs: just get the global indices */ 2383 numColIndices = dof; 2384 matSize = 0; 2385 } 2386 /* we will pack the column indices with the field offsets */ 2387 ierr = PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0);CHKERRQ(ierr); 2388 ierr = PetscSectionSetDof(rootMatricesSec,p,matSize);CHKERRQ(ierr); 2389 } 2390 ierr = PetscSectionSetUp(rootIndicesSec);CHKERRQ(ierr); 2391 ierr = PetscSectionSetUp(rootMatricesSec);CHKERRQ(ierr); 2392 { 2393 PetscInt numRootIndices, numRootMatrices; 2394 2395 ierr = PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);CHKERRQ(ierr); 2396 ierr = PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices);CHKERRQ(ierr); 2397 ierr = PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices);CHKERRQ(ierr); 2398 for (p = pStartC; p < pEndC; p++) { 2399 PetscInt numRowIndices, numColIndices, matSize, dof; 2400 PetscInt pIndOff, pMatOff; 2401 PetscInt *pInd; 2402 PetscInt maxChildId = maxChildIds[p - pStartC]; 2403 PetscScalar *pMat = NULL; 2404 2405 ierr = PetscSectionGetDof(rootIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2406 if (!numColIndices) { 2407 continue; 2408 } 2409 offsets[0] = 0; 2410 newOffsets[0] = 0; 2411 offsetsCopy[0] = 0; 2412 newOffsetsCopy[0] = 0; 2413 numColIndices -= 2 * numFields; 2414 ierr = PetscSectionGetOffset(rootIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2415 pInd = &(rootIndices[pIndOff]); 2416 ierr = PetscSectionGetDof(rootMatricesSec,p,&matSize);CHKERRQ(ierr); 2417 if (matSize) { 2418 ierr = PetscSectionGetOffset(rootMatricesSec,p,&pMatOff);CHKERRQ(ierr); 2419 pMat = &rootMatrices[pMatOff]; 2420 } 2421 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 2422 if (dof < 0) { 2423 dof = -(dof + 1); 2424 } 2425 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 2426 PetscInt i, j; 2427 PetscInt numRowIndices = matSize / numColIndices; 2428 2429 if (!numRowIndices) { /* don't need to calculate the mat, just the indices */ 2430 PetscInt numIndices, *indices; 2431 ierr = DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);CHKERRQ(ierr); 2432 if (numIndices != numColIndices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations"); 2433 for (i = 0; i < numColIndices; i++) { 2434 pInd[i] = indices[i]; 2435 } 2436 for (i = 0; i < numFields; i++) { 2437 pInd[numColIndices + i] = offsets[i+1]; 2438 pInd[numColIndices + numFields + i] = offsets[i+1]; 2439 } 2440 ierr = DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);CHKERRQ(ierr); 2441 } 2442 else { 2443 PetscInt closureSize, *closure = NULL, cl; 2444 PetscScalar *pMatIn, *pMatModified; 2445 PetscInt numPoints,*points; 2446 2447 ierr = DMGetWorkArray(coarse,numRowIndices * numRowIndices,PETSC_SCALAR,&pMatIn);CHKERRQ(ierr); 2448 for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */ 2449 for (j = 0; j < numRowIndices; j++) { 2450 pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.; 2451 } 2452 } 2453 ierr = DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2454 if (numFields) { 2455 PetscInt f; 2456 2457 for (cl = 0; cl < closureSize; cl++) { 2458 PetscInt c = closure[2 * cl]; 2459 2460 for (f = 0; f < numFields; f++) { 2461 PetscInt fDof; 2462 2463 ierr = PetscSectionGetFieldDof(localCoarse,c,f,&fDof);CHKERRQ(ierr); 2464 offsets[f + 1] += fDof; 2465 } 2466 } 2467 for (f = 0; f < numFields; f++) { 2468 offsets[f + 1] += offsets[f]; 2469 newOffsets[f + 1] = offsets[f + 1]; 2470 } 2471 } 2472 /* apply hanging node constraints on the right, get the new points and the new offsets */ 2473 ierr = DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE);CHKERRQ(ierr); 2474 if (!numFields) { 2475 for (i = 0; i < numRowIndices * numColIndices; i++) { 2476 pMat[i] = pMatModified[i]; 2477 } 2478 } 2479 else { 2480 PetscInt f, i, j, count; 2481 for (f = 0, count = 0; f < numFields; f++) { 2482 for (i = offsets[f]; i < offsets[f+1]; i++) { 2483 for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) { 2484 pMat[count] = pMatModified[i * numColIndices + j]; 2485 } 2486 } 2487 } 2488 } 2489 ierr = DMRestoreWorkArray(coarse,numRowIndices * numColIndices,PETSC_SCALAR,&pMatModified);CHKERRQ(ierr); 2490 ierr = DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2491 ierr = DMRestoreWorkArray(coarse,numRowIndices * numColIndices,PETSC_SCALAR,&pMatIn);CHKERRQ(ierr); 2492 if (numFields) { 2493 PetscInt f; 2494 for (f = 0; f < numFields; f++) { 2495 pInd[numColIndices + f] = offsets[f+1]; 2496 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2497 } 2498 for (cl = 0; cl < numPoints*2; cl += 2) { 2499 PetscInt o = points[cl+1], globalOff, c = points[cl]; 2500 ierr = PetscSectionGetOffset(globalCoarse, c, &globalOff);CHKERRQ(ierr); 2501 indicesPointFields_private(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, o, pInd); 2502 } 2503 } else { 2504 for (cl = 0; cl < numPoints*2; cl += 2) { 2505 PetscInt o = points[cl+1], c = points[cl], globalOff; 2506 ierr = PetscSectionGetOffset(globalCoarse, c, &globalOff);CHKERRQ(ierr); 2507 indicesPoint_private(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, o, pInd); 2508 } 2509 } 2510 } 2511 } 2512 else if (matSize) { 2513 PetscInt cOff; 2514 PetscInt *rowIndices, *colIndices, a, aDof, aOff; 2515 2516 numRowIndices = matSize / numColIndices; 2517 if (numRowIndices != dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Miscounted dofs"); 2518 ierr = DMGetWorkArray(coarse,numRowIndices,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2519 ierr = DMGetWorkArray(coarse,numColIndices,PETSC_INT,&colIndices);CHKERRQ(ierr); 2520 ierr = PetscSectionGetOffset(cSec,p,&cOff);CHKERRQ(ierr); 2521 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 2522 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 2523 if (numFields) { 2524 PetscInt f; 2525 2526 for (f = 0; f < numFields; f++) { 2527 PetscInt fDof; 2528 ierr = PetscSectionGetFieldDof(cSec,p,f,&fDof);CHKERRQ(ierr); 2529 offsets[f + 1] = fDof; 2530 for (a = 0; a < aDof; a++) { 2531 PetscInt anchor = anchors[a + aOff]; 2532 ierr = PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);CHKERRQ(ierr); 2533 newOffsets[f + 1] += fDof; 2534 } 2535 } 2536 for (f = 0; f < numFields; f++) { 2537 offsets[f + 1] += offsets[f]; 2538 offsetsCopy[f + 1] = offsets[f + 1]; 2539 newOffsets[f + 1] += newOffsets[f]; 2540 newOffsetsCopy[f + 1] = newOffsets[f + 1]; 2541 } 2542 indicesPointFields_private(cSec,p,cOff,offsetsCopy,PETSC_TRUE,0,rowIndices);CHKERRQ(ierr); 2543 for (a = 0; a < aDof; a++) { 2544 PetscInt anchor = anchors[a + aOff], lOff; 2545 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr); 2546 indicesPointFields_private(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,0,colIndices);CHKERRQ(ierr); 2547 } 2548 } 2549 else { 2550 indicesPoint_private(cSec,p,cOff,offsetsCopy,PETSC_TRUE,0,rowIndices);CHKERRQ(ierr); 2551 for (a = 0; a < aDof; a++) { 2552 PetscInt anchor = anchors[a + aOff], lOff; 2553 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr); 2554 indicesPoint_private(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,0,colIndices);CHKERRQ(ierr); 2555 } 2556 } 2557 if (numFields) { 2558 PetscInt count, f, a; 2559 for (f = 0, count = 0; f < numFields; f++) { 2560 PetscInt iSize = offsets[f + 1] - offsets[f]; 2561 PetscInt jSize = newOffsets[f + 1] - newOffsets[f]; 2562 ierr = MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);CHKERRQ(ierr); 2563 count += iSize * jSize; 2564 pInd[numColIndices + f] = offsets[f+1]; 2565 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2566 } 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 indicesPointFields_private(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,0,pInd);CHKERRQ(ierr); 2572 } 2573 } 2574 else { 2575 PetscInt a; 2576 ierr = MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);CHKERRQ(ierr); 2577 for (a = 0; a < aDof; a++) { 2578 PetscInt anchor = anchors[a + aOff]; 2579 PetscInt gOff; 2580 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr); 2581 indicesPoint_private(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,0,pInd);CHKERRQ(ierr); 2582 } 2583 } 2584 ierr = DMRestoreWorkArray(coarse,numColIndices,PETSC_INT,&colIndices);CHKERRQ(ierr); 2585 ierr = DMRestoreWorkArray(coarse,numRowIndices,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2586 } 2587 else { 2588 PetscInt gOff; 2589 2590 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 2591 if (numFields) { 2592 PetscInt f; 2593 2594 for (f = 0; f < numFields; f++) { 2595 PetscInt fDof; 2596 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 2597 offsets[f + 1] = fDof + offsets[f]; 2598 } 2599 for (f = 0; f < numFields; f++) { 2600 pInd[numColIndices + f] = offsets[f+1]; 2601 pInd[numColIndices + numFields + f] = offsets[f+1]; 2602 } 2603 indicesPointFields_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,0,pInd);CHKERRQ(ierr); 2604 } 2605 else { 2606 indicesPoint_private(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,0,pInd);CHKERRQ(ierr); 2607 } 2608 } 2609 } 2610 ierr = PetscFree(maxChildIds);CHKERRQ(ierr); 2611 } 2612 { 2613 PetscSF indicesSF, matricesSF; 2614 PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices; 2615 2616 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);CHKERRQ(ierr); 2617 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafMatricesSec);CHKERRQ(ierr); 2618 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootIndicesSec,&remoteOffsetsIndices,leafIndicesSec);CHKERRQ(ierr); 2619 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootMatricesSec,&remoteOffsetsMatrices,leafMatricesSec);CHKERRQ(ierr); 2620 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootIndicesSec,remoteOffsetsIndices,leafIndicesSec,&indicesSF);CHKERRQ(ierr); 2621 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootMatricesSec,remoteOffsetsMatrices,leafMatricesSec,&matricesSF);CHKERRQ(ierr); 2622 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 2623 ierr = PetscFree(remoteOffsetsIndices);CHKERRQ(ierr); 2624 ierr = PetscFree(remoteOffsetsMatrices);CHKERRQ(ierr); 2625 ierr = PetscSectionGetStorageSize(leafIndicesSec,&numLeafIndices);CHKERRQ(ierr); 2626 ierr = PetscSectionGetStorageSize(leafMatricesSec,&numLeafMatrices);CHKERRQ(ierr); 2627 ierr = PetscMalloc2(numLeafIndices,&leafIndices,numLeafMatrices,&leafMatrices);CHKERRQ(ierr); 2628 ierr = PetscSFBcastBegin(indicesSF,MPIU_INT,rootIndices,leafIndices);CHKERRQ(ierr); 2629 ierr = PetscSFBcastBegin(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);CHKERRQ(ierr); 2630 ierr = PetscSFBcastEnd(indicesSF,MPIU_INT,rootIndices,leafIndices);CHKERRQ(ierr); 2631 ierr = PetscSFBcastEnd(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);CHKERRQ(ierr); 2632 ierr = PetscSFDestroy(&matricesSF);CHKERRQ(ierr); 2633 ierr = PetscSFDestroy(&indicesSF);CHKERRQ(ierr); 2634 ierr = PetscFree2(rootIndices,rootMatrices);CHKERRQ(ierr); 2635 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 2636 ierr = PetscSectionDestroy(&rootMatricesSec);CHKERRQ(ierr); 2637 } 2638 /* count to preallocate */ 2639 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 2640 { 2641 PetscInt nGlobal; 2642 PetscInt *dnnz, *onnz; 2643 PetscLayout rowMap, colMap; 2644 PetscInt rowStart, rowEnd, colStart, colEnd; 2645 PetscInt maxDof; 2646 PetscInt *rowIndices; 2647 DM refTree; 2648 PetscInt **refPointFieldN; 2649 PetscScalar ***refPointFieldMats; 2650 PetscSection refConSec, refAnSec; 2651 PetscInt pRefStart,pRefEnd,maxConDof,maxColumns; 2652 PetscScalar *pointWork; 2653 2654 ierr = PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal);CHKERRQ(ierr); 2655 ierr = PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz);CHKERRQ(ierr); 2656 ierr = MatGetLayouts(mat,&rowMap,&colMap);CHKERRQ(ierr); 2657 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 2658 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 2659 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 2660 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 2661 ierr = PetscSectionGetMaxDof(globalFine,&maxDof);CHKERRQ(ierr); 2662 ierr = DMGetWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2663 for (p = pStartF; p < pEndF; p++) { 2664 PetscInt gDof, gcDof, gOff; 2665 PetscInt numColIndices, pIndOff, *pInd; 2666 PetscInt matSize; 2667 PetscInt i; 2668 2669 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 2670 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 2671 if ((gDof - gcDof) <= 0) { 2672 continue; 2673 } 2674 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 2675 if (gOff < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I though having global dofs meant a non-negative offset"); 2676 if ((gOff < rowStart) || ((gOff + gDof - gcDof) > rowEnd)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I thought the row map would constrain the global dofs"); 2677 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2678 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2679 numColIndices -= 2 * numFields; 2680 if (numColIndices <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"global fine dof with no dofs to interpolate from"); 2681 pInd = &leafIndices[pIndOff]; 2682 offsets[0] = 0; 2683 offsetsCopy[0] = 0; 2684 newOffsets[0] = 0; 2685 newOffsetsCopy[0] = 0; 2686 if (numFields) { 2687 PetscInt f; 2688 for (f = 0; f < numFields; f++) { 2689 PetscInt rowDof; 2690 2691 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 2692 offsets[f + 1] = offsets[f] + rowDof; 2693 offsetsCopy[f + 1] = offsets[f + 1]; 2694 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2695 numD[f] = 0; 2696 numO[f] = 0; 2697 } 2698 ierr = indicesPointFields_private(localFine,p,gOff,offsetsCopy,PETSC_FALSE,0,rowIndices);CHKERRQ(ierr); 2699 for (f = 0; f < numFields; f++) { 2700 PetscInt colOffset = newOffsets[f]; 2701 PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f]; 2702 2703 for (i = 0; i < numFieldCols; i++) { 2704 PetscInt gInd = pInd[i + colOffset]; 2705 2706 if (gInd >= colStart && gInd < colEnd) { 2707 numD[f]++; 2708 } 2709 else if (gInd >= 0) { /* negative means non-entry */ 2710 numO[f]++; 2711 } 2712 } 2713 } 2714 } 2715 else { 2716 ierr = indicesPoint_private(localFine,p,gOff,offsetsCopy,PETSC_FALSE,0,rowIndices);CHKERRQ(ierr); 2717 numD[0] = 0; 2718 numO[0] = 0; 2719 for (i = 0; i < numColIndices; i++) { 2720 PetscInt gInd = pInd[i]; 2721 2722 if (gInd >= colStart && gInd < colEnd) { 2723 numD[0]++; 2724 } 2725 else if (gInd >= 0) { /* negative means non-entry */ 2726 numO[0]++; 2727 } 2728 } 2729 } 2730 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr); 2731 if (!matSize) { /* incoming matrix is identity */ 2732 PetscInt childId; 2733 2734 childId = childIds[p-pStartF]; 2735 if (childId < 0) { /* no child interpolation: one nnz per */ 2736 if (numFields) { 2737 PetscInt f; 2738 for (f = 0; f < numFields; f++) { 2739 PetscInt numRows = offsets[f+1] - offsets[f], row; 2740 for (row = 0; row < numRows; row++) { 2741 PetscInt gIndCoarse = pInd[newOffsets[f] + row]; 2742 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2743 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2744 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2745 dnnz[gIndFine - rowStart] = 1; 2746 } 2747 else if (gIndCoarse >= 0) { /* remote */ 2748 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2749 onnz[gIndFine - rowStart] = 1; 2750 } 2751 else { /* constrained */ 2752 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2753 } 2754 } 2755 } 2756 } 2757 else { 2758 PetscInt i; 2759 for (i = 0; i < gDof; i++) { 2760 PetscInt gIndCoarse = pInd[i]; 2761 PetscInt gIndFine = rowIndices[i]; 2762 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2763 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2764 dnnz[gIndFine - rowStart] = 1; 2765 } 2766 else if (gIndCoarse >= 0) { /* remote */ 2767 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2768 onnz[gIndFine - rowStart] = 1; 2769 } 2770 else { /* constrained */ 2771 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2772 } 2773 } 2774 } 2775 } 2776 else { /* interpolate from all */ 2777 if (numFields) { 2778 PetscInt f; 2779 for (f = 0; f < numFields; f++) { 2780 PetscInt numRows = offsets[f+1] - offsets[f], row; 2781 for (row = 0; row < numRows; row++) { 2782 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2783 if (gIndFine >= 0) { 2784 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2785 dnnz[gIndFine - rowStart] = numD[f]; 2786 onnz[gIndFine - rowStart] = numO[f]; 2787 } 2788 } 2789 } 2790 } 2791 else { 2792 PetscInt i; 2793 for (i = 0; i < gDof; i++) { 2794 PetscInt gIndFine = rowIndices[i]; 2795 if (gIndFine >= 0) { 2796 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2797 dnnz[gIndFine - rowStart] = numD[0]; 2798 onnz[gIndFine - rowStart] = numO[0]; 2799 } 2800 } 2801 } 2802 } 2803 } 2804 else { /* interpolate from all */ 2805 if (numFields) { 2806 PetscInt f; 2807 for (f = 0; f < numFields; f++) { 2808 PetscInt numRows = offsets[f+1] - offsets[f], row; 2809 for (row = 0; row < numRows; row++) { 2810 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2811 if (gIndFine >= 0) { 2812 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2813 dnnz[gIndFine - rowStart] = numD[f]; 2814 onnz[gIndFine - rowStart] = numO[f]; 2815 } 2816 } 2817 } 2818 } 2819 else { /* every dof get a full row */ 2820 PetscInt i; 2821 for (i = 0; i < gDof; i++) { 2822 PetscInt gIndFine = rowIndices[i]; 2823 if (gIndFine >= 0) { 2824 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2825 dnnz[gIndFine - rowStart] = numD[0]; 2826 onnz[gIndFine - rowStart] = numO[0]; 2827 } 2828 } 2829 } 2830 } 2831 } 2832 ierr = MatXAIJSetPreallocation(mat,1,dnnz,onnz,NULL,NULL);CHKERRQ(ierr); 2833 ierr = PetscFree2(dnnz,onnz);CHKERRQ(ierr); 2834 2835 ierr = DMPlexGetReferenceTree(fine,&refTree);CHKERRQ(ierr); 2836 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 2837 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 2838 ierr = DMPlexGetAnchors(refTree,&refAnSec,NULL);CHKERRQ(ierr); 2839 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 2840 ierr = PetscSectionGetMaxDof(refConSec,&maxConDof);CHKERRQ(ierr); 2841 ierr = PetscSectionGetMaxDof(leafIndicesSec,&maxColumns);CHKERRQ(ierr); 2842 ierr = PetscMalloc1(maxConDof*maxColumns,&pointWork);CHKERRQ(ierr); 2843 for (p = pStartF; p < pEndF; p++) { 2844 PetscInt gDof, gcDof, gOff; 2845 PetscInt numColIndices, pIndOff, *pInd; 2846 PetscInt matSize; 2847 PetscInt childId; 2848 2849 2850 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 2851 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 2852 if ((gDof - gcDof) <= 0) { 2853 continue; 2854 } 2855 childId = childIds[p-pStartF]; 2856 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 2857 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2858 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2859 numColIndices -= 2 * numFields; 2860 pInd = &leafIndices[pIndOff]; 2861 offsets[0] = 0; 2862 offsetsCopy[0] = 0; 2863 newOffsets[0] = 0; 2864 newOffsetsCopy[0] = 0; 2865 rowOffsets[0] = 0; 2866 if (numFields) { 2867 PetscInt f; 2868 for (f = 0; f < numFields; f++) { 2869 PetscInt rowDof; 2870 2871 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 2872 offsets[f + 1] = offsets[f] + rowDof; 2873 offsetsCopy[f + 1] = offsets[f + 1]; 2874 rowOffsets[f + 1] = pInd[numColIndices + f]; 2875 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2876 } 2877 ierr = indicesPointFields_private(localFine,p,gOff,offsetsCopy,PETSC_FALSE,0,rowIndices);CHKERRQ(ierr); 2878 } 2879 else { 2880 ierr = indicesPoint_private(localFine,p,gOff,offsetsCopy,PETSC_FALSE,0,rowIndices);CHKERRQ(ierr); 2881 } 2882 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr); 2883 if (!matSize) { /* incoming matrix is identity */ 2884 if (childId < 0) { /* no child interpolation: scatter */ 2885 if (numFields) { 2886 PetscInt f; 2887 for (f = 0; f < numFields; f++) { 2888 PetscInt numRows = offsets[f+1] - offsets[f], row; 2889 for (row = 0; row < numRows; row++) { 2890 ierr = MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES);CHKERRQ(ierr); 2891 } 2892 } 2893 } 2894 else { 2895 PetscInt numRows = gDof, row; 2896 for (row = 0; row < numRows; row++) { 2897 ierr = MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES);CHKERRQ(ierr); 2898 } 2899 } 2900 } 2901 else { /* interpolate from all */ 2902 if (numFields) { 2903 PetscInt f; 2904 for (f = 0; f < numFields; f++) { 2905 PetscInt numRows = offsets[f+1] - offsets[f]; 2906 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 2907 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES);CHKERRQ(ierr); 2908 } 2909 } 2910 else { 2911 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES);CHKERRQ(ierr); 2912 } 2913 } 2914 } 2915 else { /* interpolate from all */ 2916 PetscInt pMatOff; 2917 PetscScalar *pMat; 2918 2919 ierr = PetscSectionGetOffset(leafMatricesSec,p,&pMatOff);CHKERRQ(ierr); 2920 pMat = &leafMatrices[pMatOff]; 2921 if (childId < 0) { /* copy the incoming matrix */ 2922 if (numFields) { 2923 PetscInt f, count; 2924 for (f = 0, count = 0; f < numFields; f++) { 2925 PetscInt numRows = offsets[f+1]-offsets[f]; 2926 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2927 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2928 PetscScalar *inMat = &pMat[count]; 2929 2930 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES);CHKERRQ(ierr); 2931 count += numCols * numInRows; 2932 } 2933 } 2934 else { 2935 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES);CHKERRQ(ierr); 2936 } 2937 } 2938 else { /* multiply the incoming matrix by the child interpolation */ 2939 if (numFields) { 2940 PetscInt f, count; 2941 for (f = 0, count = 0; f < numFields; f++) { 2942 PetscInt numRows = offsets[f+1]-offsets[f]; 2943 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2944 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2945 PetscScalar *inMat = &pMat[count]; 2946 PetscInt i, j, k; 2947 if (refPointFieldN[childId - pRefStart][f] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2948 for (i = 0; i < numRows; i++) { 2949 for (j = 0; j < numCols; j++) { 2950 PetscScalar val = 0.; 2951 for (k = 0; k < numInRows; k++) { 2952 val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j]; 2953 } 2954 pointWork[i * numCols + j] = val; 2955 } 2956 } 2957 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES);CHKERRQ(ierr); 2958 count += numCols * numInRows; 2959 } 2960 } 2961 else { /* every dof gets a full row */ 2962 PetscInt numRows = gDof; 2963 PetscInt numCols = numColIndices; 2964 PetscInt numInRows = matSize / numColIndices; 2965 PetscInt i, j, k; 2966 if (refPointFieldN[childId - pRefStart][0] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2967 for (i = 0; i < numRows; i++) { 2968 for (j = 0; j < numCols; j++) { 2969 PetscScalar val = 0.; 2970 for (k = 0; k < numInRows; k++) { 2971 val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j]; 2972 } 2973 pointWork[i * numCols + j] = val; 2974 } 2975 } 2976 ierr = MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES);CHKERRQ(ierr); 2977 } 2978 } 2979 } 2980 } 2981 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 2982 ierr = DMRestoreWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2983 ierr = PetscFree(pointWork);CHKERRQ(ierr); 2984 } 2985 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2986 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2987 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr); 2988 ierr = PetscSectionDestroy(&leafMatricesSec);CHKERRQ(ierr); 2989 ierr = PetscFree2(leafIndices,leafMatrices);CHKERRQ(ierr); 2990 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr); 2991 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 2992 PetscFunctionReturn(0); 2993 } 2994