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