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__ "EvaluateBasis" 1221 static PetscErrorCode EvaluateBasis(PetscSpace space, PetscInt nFunctionals, PetscInt nPoints, const PetscInt *pointsPerFn, const PetscReal *points, const PetscReal *weights, PetscReal *work, Mat basisAtPoints) 1222 { 1223 PetscInt i, j, k, offset, qPoints; 1224 PetscErrorCode ierr; 1225 1226 PetscFunctionBegin; 1227 ierr = PetscSpaceEvaluate(space,nPoints,points,work,NULL,NULL);CHKERRQ(ierr); 1228 for (i = 0, offset = 0; i < nFunctionals; i++) { 1229 qPoints = pointsPerFn[i]; 1230 for (j = 0; j < nFunctionals; j++) { 1231 PetscScalar val = 0.; 1232 1233 for (k = 0; k < qPoints; k++) { 1234 val += work[(offset + k) * nFunctionals + j] * weights[k]; 1235 } 1236 ierr = MatSetValue(basisAtPoints,j,i,val,INSERT_VALUES);CHKERRQ(ierr); 1237 } 1238 offset += qPoints; 1239 } 1240 ierr = MatAssemblyBegin(basisAtPoints,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1241 ierr = MatAssemblyEnd(basisAtPoints,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1242 PetscFunctionReturn(0); 1243 } 1244 1245 #undef __FUNCT__ 1246 #define __FUNCT__ "DMPlexComputeAnchorMatrix_Tree_Direct" 1247 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_Direct(DM dm, PetscSection section, PetscSection cSec, Mat cMat) 1248 { 1249 PetscDS ds; 1250 PetscInt spdim; 1251 PetscInt numFields, f, c, cStart, cEnd, pStart, pEnd, conStart, conEnd; 1252 const PetscInt *anchors; 1253 PetscSection aSec; 1254 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJparent, detJ, detJparent; 1255 IS aIS; 1256 PetscErrorCode ierr; 1257 1258 PetscFunctionBegin; 1259 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 1260 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 1261 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1262 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 1263 ierr = DMPlexGetAnchors(dm,&aSec,&aIS);CHKERRQ(ierr); 1264 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 1265 ierr = PetscSectionGetChart(cSec,&conStart,&conEnd);CHKERRQ(ierr); 1266 ierr = DMGetDimension(dm,&spdim);CHKERRQ(ierr); 1267 ierr = PetscMalloc6(spdim,&v0,spdim,&v0parent,spdim,&vtmp,spdim*spdim,&J,spdim*spdim,&Jparent,spdim*spdim,&invJparent);CHKERRQ(ierr); 1268 1269 for (f = 0; f < numFields; f++) { 1270 PetscObject disc; 1271 PetscClassId id; 1272 PetscSpace bspace; 1273 PetscDualSpace dspace; 1274 PetscInt i, j, k, nPoints, offset; 1275 PetscInt fSize, fComp, maxDof; 1276 PetscReal *weights, *pointsRef, *pointsReal, *work; 1277 PetscScalar *scwork, *X; 1278 PetscInt *sizes, *workIndRow, *workIndCol; 1279 Mat Amat, Bmat, Xmat; 1280 const PetscInt *numDof = NULL; 1281 const PetscInt ***perms = NULL; 1282 const PetscScalar ***flips = NULL; 1283 1284 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 1285 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 1286 if (id == PETSCFE_CLASSID) { 1287 PetscFE fe = (PetscFE) disc; 1288 1289 ierr = PetscFEGetBasisSpace(fe,&bspace);CHKERRQ(ierr); 1290 ierr = PetscFEGetDualSpace(fe,&dspace);CHKERRQ(ierr); 1291 ierr = PetscDualSpaceGetDimension(dspace,&fSize);CHKERRQ(ierr); 1292 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 1293 } 1294 else if (id == PETSCFV_CLASSID) { 1295 PetscFV fv = (PetscFV) disc; 1296 1297 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject)fv),&bspace);CHKERRQ(ierr); 1298 ierr = PetscSpaceSetType(bspace,PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 1299 ierr = PetscSpaceSetOrder(bspace,0);CHKERRQ(ierr); 1300 ierr = PetscSpacePolynomialSetNumVariables(bspace,spdim);CHKERRQ(ierr); 1301 ierr = PetscSpaceSetUp(bspace);CHKERRQ(ierr); 1302 ierr = PetscFVGetDualSpace(fv,&dspace);CHKERRQ(ierr); 1303 ierr = PetscDualSpaceGetDimension(dspace,&fSize);CHKERRQ(ierr); 1304 ierr = PetscFVGetNumComponents(fv,&fComp);CHKERRQ(ierr); 1305 } 1306 else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id); 1307 ierr = PetscDualSpaceGetNumDof(dspace,&numDof);CHKERRQ(ierr); 1308 for (i = 0, maxDof = 0; i <= spdim; i++) {maxDof = PetscMax(maxDof,numDof[i]);} 1309 ierr = PetscDualSpaceGetSymmetries(dspace,&perms,&flips);CHKERRQ(ierr); 1310 1311 ierr = MatCreate(PETSC_COMM_SELF,&Amat);CHKERRQ(ierr); 1312 ierr = MatSetSizes(Amat,fSize,fSize,fSize,fSize);CHKERRQ(ierr); 1313 ierr = MatSetType(Amat,MATSEQDENSE);CHKERRQ(ierr); 1314 ierr = MatSetUp(Amat);CHKERRQ(ierr); 1315 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Bmat);CHKERRQ(ierr); 1316 ierr = MatDuplicate(Amat,MAT_DO_NOT_COPY_VALUES,&Xmat);CHKERRQ(ierr); 1317 nPoints = 0; 1318 for (i = 0; i < fSize; i++) { 1319 PetscInt qPoints; 1320 PetscQuadrature quad; 1321 1322 ierr = PetscDualSpaceGetFunctional(dspace,i,&quad);CHKERRQ(ierr); 1323 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,NULL,NULL);CHKERRQ(ierr); 1324 nPoints += qPoints; 1325 } 1326 ierr = PetscMalloc7(fSize,&sizes,nPoints,&weights,spdim*nPoints,&pointsRef,spdim*nPoints,&pointsReal,fSize*nPoints,&work,maxDof,&workIndRow,maxDof,&workIndCol);CHKERRQ(ierr); 1327 ierr = PetscMalloc1(maxDof * maxDof,&scwork);CHKERRQ(ierr); 1328 offset = 0; 1329 for (i = 0; i < fSize; i++) { 1330 PetscInt qPoints; 1331 const PetscReal *p, *w; 1332 PetscQuadrature quad; 1333 1334 ierr = PetscDualSpaceGetFunctional(dspace,i,&quad);CHKERRQ(ierr); 1335 ierr = PetscQuadratureGetData(quad,NULL,&qPoints,&p,&w);CHKERRQ(ierr); 1336 ierr = PetscMemcpy(weights+offset,w,qPoints*sizeof(*w));CHKERRQ(ierr); 1337 ierr = PetscMemcpy(pointsRef+spdim*offset,p,spdim*qPoints*sizeof(*p));CHKERRQ(ierr); 1338 sizes[i] = qPoints; 1339 offset += qPoints; 1340 } 1341 ierr = EvaluateBasis(bspace,fSize,nPoints,sizes,pointsRef,weights,work,Amat);CHKERRQ(ierr); 1342 ierr = MatLUFactor(Amat,NULL,NULL,NULL);CHKERRQ(ierr); 1343 for (c = cStart; c < cEnd; c++) { 1344 PetscInt parent; 1345 PetscInt closureSize, closureSizeP, *closure = NULL, *closureP = NULL; 1346 PetscInt *childOffsets, *parentOffsets; 1347 1348 ierr = DMPlexGetTreeParent(dm,c,&parent,NULL);CHKERRQ(ierr); 1349 if (parent == c) continue; 1350 ierr = DMPlexGetTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1351 for (i = 0; i < closureSize; i++) { 1352 PetscInt p = closure[2*i]; 1353 PetscInt conDof; 1354 1355 if (p < conStart || p >= conEnd) continue; 1356 if (numFields) { 1357 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 1358 } 1359 else { 1360 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 1361 } 1362 if (conDof) break; 1363 } 1364 if (i == closureSize) { 1365 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1366 continue; 1367 } 1368 1369 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, NULL, &detJ);CHKERRQ(ierr); 1370 ierr = DMPlexComputeCellGeometryFEM(dm, parent, NULL, v0parent, Jparent, invJparent, &detJparent);CHKERRQ(ierr); 1371 for (i = 0; i < nPoints; i++) { 1372 CoordinatesRefToReal(spdim, spdim, v0, J, &pointsRef[i*spdim],vtmp); 1373 CoordinatesRealToRef(spdim, spdim, v0parent, invJparent, vtmp, &pointsReal[i*spdim]); 1374 } 1375 ierr = EvaluateBasis(bspace,fSize,nPoints,sizes,pointsReal,weights,work,Bmat);CHKERRQ(ierr); 1376 ierr = MatMatSolve(Amat,Bmat,Xmat);CHKERRQ(ierr); 1377 ierr = MatDenseGetArray(Xmat,&X);CHKERRQ(ierr); 1378 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 1379 ierr = PetscMalloc2(closureSize+1,&childOffsets,closureSizeP+1,&parentOffsets);CHKERRQ(ierr); 1380 childOffsets[0] = 0; 1381 for (i = 0; i < closureSize; i++) { 1382 PetscInt p = closure[2*i]; 1383 PetscInt dof; 1384 1385 if (numFields) { 1386 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 1387 } 1388 else { 1389 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 1390 } 1391 childOffsets[i+1]=childOffsets[i]+dof / fComp; 1392 } 1393 parentOffsets[0] = 0; 1394 for (i = 0; i < closureSizeP; i++) { 1395 PetscInt p = closureP[2*i]; 1396 PetscInt dof; 1397 1398 if (numFields) { 1399 ierr = PetscSectionGetFieldDof(section,p,f,&dof);CHKERRQ(ierr); 1400 } 1401 else { 1402 ierr = PetscSectionGetDof(section,p,&dof);CHKERRQ(ierr); 1403 } 1404 parentOffsets[i+1]=parentOffsets[i]+dof / fComp; 1405 } 1406 for (i = 0; i < closureSize; i++) { 1407 PetscInt conDof, conOff, aDof, aOff, nWork; 1408 PetscInt p = closure[2*i]; 1409 PetscInt o = closure[2*i+1]; 1410 const PetscInt *perm; 1411 const PetscScalar *flip; 1412 1413 if (p < conStart || p >= conEnd) continue; 1414 if (numFields) { 1415 ierr = PetscSectionGetFieldDof(cSec,p,f,&conDof);CHKERRQ(ierr); 1416 ierr = PetscSectionGetFieldOffset(cSec,p,f,&conOff);CHKERRQ(ierr); 1417 } 1418 else { 1419 ierr = PetscSectionGetDof(cSec,p,&conDof);CHKERRQ(ierr); 1420 ierr = PetscSectionGetOffset(cSec,p,&conOff);CHKERRQ(ierr); 1421 } 1422 if (!conDof) continue; 1423 perm = (perms && perms[i]) ? perms[i][o] : NULL; 1424 flip = (flips && flips[i]) ? flips[i][o] : NULL; 1425 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 1426 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 1427 nWork = childOffsets[i+1]-childOffsets[i]; 1428 for (k = 0; k < aDof; k++) { 1429 PetscInt a = anchors[aOff + k]; 1430 PetscInt aSecDof, aSecOff; 1431 1432 if (numFields) { 1433 ierr = PetscSectionGetFieldDof(section,a,f,&aSecDof);CHKERRQ(ierr); 1434 ierr = PetscSectionGetFieldOffset(section,a,f,&aSecOff);CHKERRQ(ierr); 1435 } 1436 else { 1437 ierr = PetscSectionGetDof(section,a,&aSecDof);CHKERRQ(ierr); 1438 ierr = PetscSectionGetOffset(section,a,&aSecOff);CHKERRQ(ierr); 1439 } 1440 if (!aSecDof) continue; 1441 1442 for (j = 0; j < closureSizeP; j++) { 1443 PetscInt q = closureP[2*j]; 1444 PetscInt oq = closureP[2*j+1]; 1445 1446 if (q == a) { 1447 PetscInt r, s, t, nWorkP; 1448 const PetscInt *permP; 1449 const PetscScalar *flipP; 1450 1451 permP = (perms && perms[j]) ? perms[j][oq] : NULL; 1452 flipP = (flips && flips[j]) ? flips[j][oq] : NULL; 1453 nWorkP = parentOffsets[j+1]-parentOffsets[j]; 1454 /* get a copy of the child-to-anchor portion of the matrix, and transpose so that rows correspond to the 1455 * child and columns correspond to the anchor: BUT the maxrix returned by MatDenseGetArray is 1456 * column-major, so transpose-transpose = do nothing */ 1457 for (r = 0; r < nWork; r++) { 1458 for (s = 0; s < nWorkP; s++) { 1459 scwork[r * nWorkP + s] = X[fSize * (r + childOffsets[i]) + (s + parentOffsets[j])]; 1460 } 1461 } 1462 for (r = 0; r < nWork; r++) {workIndRow[perm ? perm[r] : r] = conOff + fComp * r;} 1463 for (s = 0; s < nWorkP; s++) {workIndCol[permP ? permP[s] : s] = aSecOff + fComp * s;} 1464 if (flip) { 1465 for (r = 0; r < nWork; r++) { 1466 for (s = 0; s < nWorkP; s++) { 1467 scwork[r * nWorkP + s] *= flip[r]; 1468 } 1469 } 1470 } 1471 if (flipP) { 1472 for (r = 0; r < nWork; r++) { 1473 for (s = 0; s < nWorkP; s++) { 1474 scwork[r * nWorkP + s] *= flipP[s]; 1475 } 1476 } 1477 } 1478 for (t = 0; t < fComp; t++) { 1479 ierr = MatSetValues(cMat,nWork,workIndRow,nWorkP,workIndCol,scwork,INSERT_VALUES);CHKERRQ(ierr); 1480 if (t < fComp - 1) { 1481 for (r = 0; r < nWork; r++) {workIndRow[r]++;} 1482 for (s = 0; s < nWorkP; s++) {workIndCol[s]++;} 1483 } 1484 } 1485 break; 1486 } 1487 } 1488 } 1489 } 1490 ierr = MatDenseRestoreArray(Xmat,&X);CHKERRQ(ierr); 1491 ierr = PetscFree2(childOffsets,parentOffsets);CHKERRQ(ierr); 1492 ierr = DMPlexRestoreTransitiveClosure(dm,c,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1493 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSizeP,&closureP);CHKERRQ(ierr); 1494 } 1495 ierr = MatDestroy(&Amat);CHKERRQ(ierr); 1496 ierr = MatDestroy(&Bmat);CHKERRQ(ierr); 1497 ierr = MatDestroy(&Xmat);CHKERRQ(ierr); 1498 ierr = PetscFree(scwork);CHKERRQ(ierr); 1499 ierr = PetscFree7(sizes,weights,pointsRef,pointsReal,work,workIndRow,workIndCol);CHKERRQ(ierr); 1500 if (id == PETSCFV_CLASSID) { 1501 ierr = PetscSpaceDestroy(&bspace);CHKERRQ(ierr); 1502 } 1503 } 1504 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1505 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1506 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJparent);CHKERRQ(ierr); 1507 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 1508 1509 PetscFunctionReturn(0); 1510 } 1511 1512 #undef __FUNCT__ 1513 #define __FUNCT__ "DMPlexReferenceTreeGetChildrenMatrices" 1514 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN) 1515 { 1516 Mat refCmat; 1517 PetscDS ds; 1518 PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof, maxAnDof, **refPointFieldN; 1519 PetscScalar ***refPointFieldMats; 1520 PetscSection refConSec, refAnSec, refSection; 1521 IS refAnIS; 1522 const PetscInt *refAnchors; 1523 const PetscInt **perms; 1524 const PetscScalar **flips; 1525 PetscErrorCode ierr; 1526 1527 PetscFunctionBegin; 1528 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 1529 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1530 maxFields = PetscMax(1,numFields); 1531 ierr = DMGetDefaultConstraints(refTree,&refConSec,&refCmat);CHKERRQ(ierr); 1532 ierr = DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);CHKERRQ(ierr); 1533 ierr = ISGetIndices(refAnIS,&refAnchors);CHKERRQ(ierr); 1534 ierr = DMGetDefaultSection(refTree,&refSection);CHKERRQ(ierr); 1535 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1536 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);CHKERRQ(ierr); 1537 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldN);CHKERRQ(ierr); 1538 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 1539 ierr = PetscSectionGetMaxDof(refAnSec,&maxAnDof);CHKERRQ(ierr); 1540 ierr = PetscMalloc1(maxDof,&rows);CHKERRQ(ierr); 1541 ierr = PetscMalloc1(maxDof*maxAnDof,&cols);CHKERRQ(ierr); 1542 for (p = pRefStart; p < pRefEnd; p++) { 1543 PetscInt parent, closureSize, *closure = NULL, pDof; 1544 1545 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 1546 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 1547 if (!pDof || parent == p) continue; 1548 1549 ierr = PetscMalloc1(maxFields,&refPointFieldMats[p-pRefStart]);CHKERRQ(ierr); 1550 ierr = PetscCalloc1(maxFields,&refPointFieldN[p-pRefStart]);CHKERRQ(ierr); 1551 ierr = DMPlexGetTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1552 for (f = 0; f < maxFields; f++) { 1553 PetscInt cDof, cOff, numCols, r, i; 1554 1555 if (f < numFields) { 1556 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 1557 ierr = PetscSectionGetFieldOffset(refConSec,p,f,&cOff);CHKERRQ(ierr); 1558 ierr = PetscSectionGetFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);CHKERRQ(ierr); 1559 } else { 1560 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 1561 ierr = PetscSectionGetOffset(refConSec,p,&cOff);CHKERRQ(ierr); 1562 ierr = PetscSectionGetPointSyms(refSection,closureSize,closure,&perms,&flips);CHKERRQ(ierr); 1563 } 1564 1565 for (r = 0; r < cDof; r++) { 1566 rows[r] = cOff + r; 1567 } 1568 numCols = 0; 1569 for (i = 0; i < closureSize; i++) { 1570 PetscInt q = closure[2*i]; 1571 PetscInt aDof, aOff, j; 1572 const PetscInt *perm = perms ? perms[i] : NULL; 1573 1574 if (numFields) { 1575 ierr = PetscSectionGetFieldDof(refSection,q,f,&aDof);CHKERRQ(ierr); 1576 ierr = PetscSectionGetFieldOffset(refSection,q,f,&aOff);CHKERRQ(ierr); 1577 } 1578 else { 1579 ierr = PetscSectionGetDof(refSection,q,&aDof);CHKERRQ(ierr); 1580 ierr = PetscSectionGetOffset(refSection,q,&aOff);CHKERRQ(ierr); 1581 } 1582 1583 for (j = 0; j < aDof; j++) { 1584 cols[numCols++] = aOff + (perm ? perm[j] : j); 1585 } 1586 } 1587 refPointFieldN[p-pRefStart][f] = numCols; 1588 ierr = PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 1589 ierr = MatGetValues(refCmat,cDof,rows,numCols,cols,refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 1590 if (flips) { 1591 PetscInt colOff = 0; 1592 1593 for (i = 0; i < closureSize; i++) { 1594 PetscInt q = closure[2*i]; 1595 PetscInt aDof, aOff, j; 1596 const PetscScalar *flip = flips ? flips[i] : NULL; 1597 1598 if (numFields) { 1599 ierr = PetscSectionGetFieldDof(refSection,q,f,&aDof);CHKERRQ(ierr); 1600 ierr = PetscSectionGetFieldOffset(refSection,q,f,&aOff);CHKERRQ(ierr); 1601 } 1602 else { 1603 ierr = PetscSectionGetDof(refSection,q,&aDof);CHKERRQ(ierr); 1604 ierr = PetscSectionGetOffset(refSection,q,&aOff);CHKERRQ(ierr); 1605 } 1606 if (flip) { 1607 PetscInt k; 1608 for (k = 0; k < cDof; k++) { 1609 for (j = 0; j < aDof; j++) { 1610 refPointFieldMats[p-pRefStart][f][k * numCols + colOff + j] *= flip[j]; 1611 } 1612 } 1613 } 1614 colOff += aDof; 1615 } 1616 } 1617 if (numFields) { 1618 ierr = PetscSectionRestoreFieldPointSyms(refSection,f,closureSize,closure,&perms,&flips);CHKERRQ(ierr); 1619 } else { 1620 ierr = PetscSectionRestorePointSyms(refSection,closureSize,closure,&perms,&flips);CHKERRQ(ierr); 1621 } 1622 } 1623 ierr = DMPlexRestoreTransitiveClosure(refTree,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1624 } 1625 *childrenMats = refPointFieldMats; 1626 *childrenN = refPointFieldN; 1627 ierr = ISRestoreIndices(refAnIS,&refAnchors);CHKERRQ(ierr); 1628 ierr = PetscFree(rows);CHKERRQ(ierr); 1629 ierr = PetscFree(cols);CHKERRQ(ierr); 1630 PetscFunctionReturn(0); 1631 } 1632 1633 #undef __FUNCT__ 1634 #define __FUNCT__ "DMPlexReferenceTreeRestoreChildrenMatrices" 1635 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices(DM refTree, PetscScalar ****childrenMats, PetscInt ***childrenN) 1636 { 1637 PetscDS ds; 1638 PetscInt **refPointFieldN; 1639 PetscScalar ***refPointFieldMats; 1640 PetscInt numFields, maxFields, pRefStart, pRefEnd, p, f; 1641 PetscSection refConSec; 1642 PetscErrorCode ierr; 1643 1644 PetscFunctionBegin; 1645 refPointFieldN = *childrenN; 1646 *childrenN = NULL; 1647 refPointFieldMats = *childrenMats; 1648 *childrenMats = NULL; 1649 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 1650 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1651 maxFields = PetscMax(1,numFields);CHKERRQ(ierr); 1652 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 1653 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1654 for (p = pRefStart; p < pRefEnd; p++) { 1655 PetscInt parent, pDof; 1656 1657 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 1658 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 1659 if (!pDof || parent == p) continue; 1660 1661 for (f = 0; f < maxFields; f++) { 1662 PetscInt cDof; 1663 1664 if (numFields) { 1665 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 1666 } 1667 else { 1668 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 1669 } 1670 1671 ierr = PetscFree(refPointFieldMats[p - pRefStart][f]);CHKERRQ(ierr); 1672 } 1673 ierr = PetscFree(refPointFieldMats[p - pRefStart]);CHKERRQ(ierr); 1674 ierr = PetscFree(refPointFieldN[p - pRefStart]);CHKERRQ(ierr); 1675 } 1676 ierr = PetscFree(refPointFieldMats);CHKERRQ(ierr); 1677 ierr = PetscFree(refPointFieldN);CHKERRQ(ierr); 1678 PetscFunctionReturn(0); 1679 } 1680 1681 #undef __FUNCT__ 1682 #define __FUNCT__ "DMPlexComputeAnchorMatrix_Tree_FromReference" 1683 static PetscErrorCode DMPlexComputeAnchorMatrix_Tree_FromReference(DM dm, PetscSection section, PetscSection conSec, Mat cMat) 1684 { 1685 DM refTree; 1686 PetscDS ds; 1687 Mat refCmat; 1688 PetscInt numFields, maxFields, f, pRefStart, pRefEnd, p, maxDof, maxAnDof, *perm, *iperm, pStart, pEnd, conStart, conEnd, **refPointFieldN; 1689 PetscScalar ***refPointFieldMats, *pointWork; 1690 PetscSection refConSec, refAnSec, anSec; 1691 IS refAnIS, anIS; 1692 const PetscInt *anchors; 1693 PetscErrorCode ierr; 1694 1695 PetscFunctionBegin; 1696 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1697 ierr = DMGetDS(dm,&ds);CHKERRQ(ierr); 1698 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 1699 maxFields = PetscMax(1,numFields); 1700 ierr = DMPlexGetReferenceTree(dm,&refTree);CHKERRQ(ierr); 1701 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 1702 ierr = DMGetDefaultConstraints(refTree,&refConSec,&refCmat);CHKERRQ(ierr); 1703 ierr = DMPlexGetAnchors(refTree,&refAnSec,&refAnIS);CHKERRQ(ierr); 1704 ierr = DMPlexGetAnchors(dm,&anSec,&anIS);CHKERRQ(ierr); 1705 ierr = ISGetIndices(anIS,&anchors);CHKERRQ(ierr); 1706 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 1707 ierr = PetscSectionGetChart(conSec,&conStart,&conEnd);CHKERRQ(ierr); 1708 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 1709 ierr = PetscSectionGetMaxDof(refAnSec,&maxAnDof);CHKERRQ(ierr); 1710 ierr = PetscMalloc1(maxDof*maxDof*maxAnDof,&pointWork);CHKERRQ(ierr); 1711 1712 /* step 1: get submats for every constrained point in the reference tree */ 1713 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 1714 1715 /* step 2: compute the preorder */ 1716 ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr); 1717 ierr = PetscMalloc2(pEnd-pStart,&perm,pEnd-pStart,&iperm);CHKERRQ(ierr); 1718 for (p = pStart; p < pEnd; p++) { 1719 perm[p - pStart] = p; 1720 iperm[p - pStart] = p-pStart; 1721 } 1722 for (p = 0; p < pEnd - pStart;) { 1723 PetscInt point = perm[p]; 1724 PetscInt parent; 1725 1726 ierr = DMPlexGetTreeParent(dm,point,&parent,NULL);CHKERRQ(ierr); 1727 if (parent == point) { 1728 p++; 1729 } 1730 else { 1731 PetscInt size, closureSize, *closure = NULL, i; 1732 1733 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1734 for (i = 0; i < closureSize; i++) { 1735 PetscInt q = closure[2*i]; 1736 if (iperm[q-pStart] > iperm[point-pStart]) { 1737 /* swap */ 1738 perm[p] = q; 1739 perm[iperm[q-pStart]] = point; 1740 iperm[point-pStart] = iperm[q-pStart]; 1741 iperm[q-pStart] = p; 1742 break; 1743 } 1744 } 1745 size = closureSize; 1746 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1747 if (i == size) { 1748 p++; 1749 } 1750 } 1751 } 1752 1753 /* step 3: fill the constraint matrix */ 1754 /* we are going to use a preorder progressive fill strategy. Mat doesn't 1755 * allow progressive fill without assembly, so we are going to set up the 1756 * values outside of the Mat first. 1757 */ 1758 { 1759 PetscInt nRows, row, nnz; 1760 PetscBool done; 1761 const PetscInt *ia, *ja; 1762 PetscScalar *vals; 1763 1764 ierr = MatGetRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);CHKERRQ(ierr); 1765 if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not get RowIJ of constraint matrix"); 1766 nnz = ia[nRows]; 1767 /* malloc and then zero rows right before we fill them: this way valgrind 1768 * can tell if we are doing progressive fill in the wrong order */ 1769 ierr = PetscMalloc1(nnz,&vals);CHKERRQ(ierr); 1770 for (p = 0; p < pEnd - pStart; p++) { 1771 PetscInt parent, childid, closureSize, *closure = NULL; 1772 PetscInt point = perm[p], pointDof; 1773 1774 ierr = DMPlexGetTreeParent(dm,point,&parent,&childid);CHKERRQ(ierr); 1775 if ((point < conStart) || (point >= conEnd) || (parent == point)) continue; 1776 ierr = PetscSectionGetDof(conSec,point,&pointDof);CHKERRQ(ierr); 1777 if (!pointDof) continue; 1778 ierr = DMPlexGetTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1779 for (f = 0; f < maxFields; f++) { 1780 PetscInt cDof, cOff, numCols, numFillCols, i, r, matOffset, offset; 1781 PetscScalar *pointMat; 1782 const PetscInt **perms; 1783 const PetscScalar **flips; 1784 1785 if (numFields) { 1786 ierr = PetscSectionGetFieldDof(conSec,point,f,&cDof);CHKERRQ(ierr); 1787 ierr = PetscSectionGetFieldOffset(conSec,point,f,&cOff);CHKERRQ(ierr); 1788 } 1789 else { 1790 ierr = PetscSectionGetDof(conSec,point,&cDof);CHKERRQ(ierr); 1791 ierr = PetscSectionGetOffset(conSec,point,&cOff);CHKERRQ(ierr); 1792 } 1793 if (!cDof) continue; 1794 if (numFields) {ierr = PetscSectionGetFieldPointSyms(section,f,closureSize,closure,&perms,&flips);CHKERRQ(ierr);} 1795 else {ierr = PetscSectionGetPointSyms(section,closureSize,closure,&perms,&flips);CHKERRQ(ierr);} 1796 1797 /* make sure that every row for this point is the same size */ 1798 #if defined(PETSC_USE_DEBUG) 1799 for (r = 0; r < cDof; r++) { 1800 if (cDof > 1 && r) { 1801 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])); 1802 } 1803 } 1804 #endif 1805 /* zero rows */ 1806 for (i = ia[cOff] ; i< ia[cOff+cDof];i++) { 1807 vals[i] = 0.; 1808 } 1809 matOffset = ia[cOff]; 1810 numFillCols = ia[cOff+1] - matOffset; 1811 pointMat = refPointFieldMats[childid-pRefStart][f]; 1812 numCols = refPointFieldN[childid-pRefStart][f]; 1813 offset = 0; 1814 for (i = 0; i < closureSize; i++) { 1815 PetscInt q = closure[2*i]; 1816 PetscInt aDof, aOff, j, k, qConDof, qConOff; 1817 const PetscInt *perm = perms ? perms[i] : NULL; 1818 const PetscScalar *flip = flips ? flips[i] : NULL; 1819 1820 qConDof = qConOff = 0; 1821 if (numFields) { 1822 ierr = PetscSectionGetFieldDof(section,q,f,&aDof);CHKERRQ(ierr); 1823 ierr = PetscSectionGetFieldOffset(section,q,f,&aOff);CHKERRQ(ierr); 1824 if (q >= conStart && q < conEnd) { 1825 ierr = PetscSectionGetFieldDof(conSec,q,f,&qConDof);CHKERRQ(ierr); 1826 ierr = PetscSectionGetFieldOffset(conSec,q,f,&qConOff);CHKERRQ(ierr); 1827 } 1828 } 1829 else { 1830 ierr = PetscSectionGetDof(section,q,&aDof);CHKERRQ(ierr); 1831 ierr = PetscSectionGetOffset(section,q,&aOff);CHKERRQ(ierr); 1832 if (q >= conStart && q < conEnd) { 1833 ierr = PetscSectionGetDof(conSec,q,&qConDof);CHKERRQ(ierr); 1834 ierr = PetscSectionGetOffset(conSec,q,&qConOff);CHKERRQ(ierr); 1835 } 1836 } 1837 if (!aDof) continue; 1838 if (qConDof) { 1839 /* this point has anchors: its rows of the matrix should already 1840 * be filled, thanks to preordering */ 1841 /* first multiply into pointWork, then set in matrix */ 1842 PetscInt aMatOffset = ia[qConOff]; 1843 PetscInt aNumFillCols = ia[qConOff + 1] - aMatOffset; 1844 for (r = 0; r < cDof; r++) { 1845 for (j = 0; j < aNumFillCols; j++) { 1846 PetscScalar inVal = 0; 1847 for (k = 0; k < aDof; k++) { 1848 PetscInt col = perm ? perm[k] : k; 1849 1850 inVal += pointMat[r * numCols + offset + col] * vals[aMatOffset + aNumFillCols * k + j] * (flip ? flip[col] : 1.); 1851 } 1852 pointWork[r * aNumFillCols + j] = inVal; 1853 } 1854 } 1855 /* assume that the columns are sorted, spend less time searching */ 1856 for (j = 0, k = 0; j < aNumFillCols; j++) { 1857 PetscInt col = ja[aMatOffset + j]; 1858 for (;k < numFillCols; k++) { 1859 if (ja[matOffset + k] == col) { 1860 break; 1861 } 1862 } 1863 if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, col); 1864 for (r = 0; r < cDof; r++) { 1865 vals[matOffset + numFillCols * r + k] = pointWork[r * aNumFillCols + j]; 1866 } 1867 } 1868 } 1869 else { 1870 /* find where to put this portion of pointMat into the matrix */ 1871 for (k = 0; k < numFillCols; k++) { 1872 if (ja[matOffset + k] == aOff) { 1873 break; 1874 } 1875 } 1876 if (k == numFillCols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No nonzero space for (%d, %d)", cOff, aOff); 1877 for (r = 0; r < cDof; r++) { 1878 for (j = 0; j < aDof; j++) { 1879 PetscInt col = perm ? perm[j] : j; 1880 1881 vals[matOffset + numFillCols * r + k + col] += pointMat[r * numCols + offset + j] * (flip ? flip[col] : 1.); 1882 } 1883 } 1884 } 1885 offset += aDof; 1886 } 1887 if (numFields) { 1888 ierr = PetscSectionRestoreFieldPointSyms(section,f,closureSize,closure,&perms,&flips);CHKERRQ(ierr); 1889 } else { 1890 ierr = PetscSectionRestorePointSyms(section,closureSize,closure,&perms,&flips);CHKERRQ(ierr); 1891 } 1892 } 1893 ierr = DMPlexRestoreTransitiveClosure(dm,parent,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 1894 } 1895 for (row = 0; row < nRows; row++) { 1896 ierr = MatSetValues(cMat,1,&row,ia[row+1]-ia[row],&ja[ia[row]],&vals[ia[row]],INSERT_VALUES);CHKERRQ(ierr); 1897 } 1898 ierr = MatRestoreRowIJ(cMat,0,PETSC_FALSE,PETSC_FALSE,&nRows,&ia,&ja,&done);CHKERRQ(ierr); 1899 if (!done) SETERRQ(PetscObjectComm((PetscObject)cMat),PETSC_ERR_PLIB,"Could not restore RowIJ of constraint matrix"); 1900 ierr = MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1901 ierr = MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1902 ierr = PetscFree(vals);CHKERRQ(ierr); 1903 } 1904 1905 /* clean up */ 1906 ierr = ISRestoreIndices(anIS,&anchors);CHKERRQ(ierr); 1907 ierr = PetscFree2(perm,iperm);CHKERRQ(ierr); 1908 ierr = PetscFree(pointWork);CHKERRQ(ierr); 1909 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 1910 PetscFunctionReturn(0); 1911 } 1912 1913 #undef __FUNCT__ 1914 #define __FUNCT__ "DMPlexTreeRefineCell" 1915 /* refine a single cell on rank 0: this is not intended to provide good local refinement, only to create an example of 1916 * a non-conforming mesh. Local refinement comes later */ 1917 PetscErrorCode DMPlexTreeRefineCell (DM dm, PetscInt cell, DM *ncdm) 1918 { 1919 DM K; 1920 PetscMPIInt rank; 1921 PetscInt dim, *pNewStart, *pNewEnd, *pNewCount, *pOldStart, *pOldEnd, offset, d, pStart, pEnd; 1922 PetscInt numNewCones, *newConeSizes, *newCones, *newOrientations; 1923 PetscInt *Kembedding; 1924 PetscInt *cellClosure=NULL, nc; 1925 PetscScalar *newVertexCoords; 1926 PetscInt numPointsWithParents, *parents, *childIDs, *perm, *iperm, *preOrient, pOffset; 1927 PetscSection parentSection; 1928 PetscErrorCode ierr; 1929 1930 PetscFunctionBegin; 1931 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr); 1932 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 1933 ierr = DMPlexCreate(PetscObjectComm((PetscObject)dm), ncdm);CHKERRQ(ierr); 1934 ierr = DMSetDimension(*ncdm,dim);CHKERRQ(ierr); 1935 1936 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1937 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm),&parentSection);CHKERRQ(ierr); 1938 ierr = DMPlexGetReferenceTree(dm,&K);CHKERRQ(ierr); 1939 if (!rank) { 1940 /* compute the new charts */ 1941 ierr = PetscMalloc5(dim+1,&pNewCount,dim+1,&pNewStart,dim+1,&pNewEnd,dim+1,&pOldStart,dim+1,&pOldEnd);CHKERRQ(ierr); 1942 offset = 0; 1943 for (d = 0; d <= dim; d++) { 1944 PetscInt pOldCount, kStart, kEnd, k; 1945 1946 pNewStart[d] = offset; 1947 ierr = DMPlexGetHeightStratum(dm,d,&pOldStart[d],&pOldEnd[d]);CHKERRQ(ierr); 1948 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr); 1949 pOldCount = pOldEnd[d] - pOldStart[d]; 1950 /* adding the new points */ 1951 pNewCount[d] = pOldCount + kEnd - kStart; 1952 if (!d) { 1953 /* removing the cell */ 1954 pNewCount[d]--; 1955 } 1956 for (k = kStart; k < kEnd; k++) { 1957 PetscInt parent; 1958 ierr = DMPlexGetTreeParent(K,k,&parent,NULL);CHKERRQ(ierr); 1959 if (parent == k) { 1960 /* avoid double counting points that won't actually be new */ 1961 pNewCount[d]--; 1962 } 1963 } 1964 pNewEnd[d] = pNewStart[d] + pNewCount[d]; 1965 offset = pNewEnd[d]; 1966 1967 } 1968 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]); 1969 /* get the current closure of the cell that we are removing */ 1970 ierr = DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);CHKERRQ(ierr); 1971 1972 ierr = PetscMalloc1(pNewEnd[dim],&newConeSizes);CHKERRQ(ierr); 1973 { 1974 PetscInt kStart, kEnd, k, closureSizeK, *closureK = NULL, j; 1975 1976 ierr = DMPlexGetChart(K,&kStart,&kEnd);CHKERRQ(ierr); 1977 ierr = PetscMalloc4(kEnd-kStart,&Kembedding,kEnd-kStart,&perm,kEnd-kStart,&iperm,kEnd-kStart,&preOrient);CHKERRQ(ierr); 1978 1979 for (k = kStart; k < kEnd; k++) { 1980 perm[k - kStart] = k; 1981 iperm [k - kStart] = k - kStart; 1982 preOrient[k - kStart] = 0; 1983 } 1984 1985 ierr = DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);CHKERRQ(ierr); 1986 for (j = 1; j < closureSizeK; j++) { 1987 PetscInt parentOrientA = closureK[2*j+1]; 1988 PetscInt parentOrientB = cellClosure[2*j+1]; 1989 PetscInt p, q; 1990 1991 p = closureK[2*j]; 1992 q = cellClosure[2*j]; 1993 for (d = 0; d <= dim; d++) { 1994 if (q >= pOldStart[d] && q < pOldEnd[d]) { 1995 Kembedding[p] = (q - pOldStart[d]) + pNewStart[d]; 1996 } 1997 } 1998 if (parentOrientA != parentOrientB) { 1999 PetscInt numChildren, i; 2000 const PetscInt *children; 2001 2002 ierr = DMPlexGetTreeChildren(K,p,&numChildren,&children);CHKERRQ(ierr); 2003 for (i = 0; i < numChildren; i++) { 2004 PetscInt kPerm, oPerm; 2005 2006 k = children[i]; 2007 ierr = DMPlexReferenceTreeGetChildSymmetry(K,p,parentOrientA,0,k,parentOrientB,&oPerm,&kPerm);CHKERRQ(ierr); 2008 /* perm = what refTree position I'm in */ 2009 perm[kPerm-kStart] = k; 2010 /* iperm = who is at this position */ 2011 iperm[k-kStart] = kPerm-kStart; 2012 preOrient[kPerm-kStart] = oPerm; 2013 } 2014 } 2015 } 2016 ierr = DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&closureSizeK,&closureK);CHKERRQ(ierr); 2017 } 2018 ierr = PetscSectionSetChart(parentSection,0,pNewEnd[dim]);CHKERRQ(ierr); 2019 offset = 0; 2020 numNewCones = 0; 2021 for (d = 0; d <= dim; d++) { 2022 PetscInt kStart, kEnd, k; 2023 PetscInt p; 2024 PetscInt size; 2025 2026 for (p = pOldStart[d]; p < pOldEnd[d]; p++) { 2027 /* skip cell 0 */ 2028 if (p == cell) continue; 2029 /* old cones to new cones */ 2030 ierr = DMPlexGetConeSize(dm,p,&size);CHKERRQ(ierr); 2031 newConeSizes[offset++] = size; 2032 numNewCones += size; 2033 } 2034 2035 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr); 2036 for (k = kStart; k < kEnd; k++) { 2037 PetscInt kParent; 2038 2039 ierr = DMPlexGetTreeParent(K,k,&kParent,NULL);CHKERRQ(ierr); 2040 if (kParent != k) { 2041 Kembedding[k] = offset; 2042 ierr = DMPlexGetConeSize(K,k,&size);CHKERRQ(ierr); 2043 newConeSizes[offset++] = size; 2044 numNewCones += size; 2045 if (kParent != 0) { 2046 ierr = PetscSectionSetDof(parentSection,Kembedding[k],1);CHKERRQ(ierr); 2047 } 2048 } 2049 } 2050 } 2051 2052 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 2053 ierr = PetscSectionGetStorageSize(parentSection,&numPointsWithParents);CHKERRQ(ierr); 2054 ierr = PetscMalloc2(numNewCones,&newCones,numNewCones,&newOrientations);CHKERRQ(ierr); 2055 ierr = PetscMalloc2(numPointsWithParents,&parents,numPointsWithParents,&childIDs);CHKERRQ(ierr); 2056 2057 /* fill new cones */ 2058 offset = 0; 2059 for (d = 0; d <= dim; d++) { 2060 PetscInt kStart, kEnd, k, l; 2061 PetscInt p; 2062 PetscInt size; 2063 const PetscInt *cone, *orientation; 2064 2065 for (p = pOldStart[d]; p < pOldEnd[d]; p++) { 2066 /* skip cell 0 */ 2067 if (p == cell) continue; 2068 /* old cones to new cones */ 2069 ierr = DMPlexGetConeSize(dm,p,&size);CHKERRQ(ierr); 2070 ierr = DMPlexGetCone(dm,p,&cone);CHKERRQ(ierr); 2071 ierr = DMPlexGetConeOrientation(dm,p,&orientation);CHKERRQ(ierr); 2072 for (l = 0; l < size; l++) { 2073 newCones[offset] = (cone[l] - pOldStart[d + 1]) + pNewStart[d + 1]; 2074 newOrientations[offset++] = orientation[l]; 2075 } 2076 } 2077 2078 ierr = DMPlexGetHeightStratum(K,d,&kStart,&kEnd);CHKERRQ(ierr); 2079 for (k = kStart; k < kEnd; k++) { 2080 PetscInt kPerm = perm[k], kParent; 2081 PetscInt preO = preOrient[k]; 2082 2083 ierr = DMPlexGetTreeParent(K,k,&kParent,NULL);CHKERRQ(ierr); 2084 if (kParent != k) { 2085 /* embed new cones */ 2086 ierr = DMPlexGetConeSize(K,k,&size);CHKERRQ(ierr); 2087 ierr = DMPlexGetCone(K,kPerm,&cone);CHKERRQ(ierr); 2088 ierr = DMPlexGetConeOrientation(K,kPerm,&orientation);CHKERRQ(ierr); 2089 for (l = 0; l < size; l++) { 2090 PetscInt q, m = (preO >= 0) ? ((preO + l) % size) : ((size -(preO + 1) - l) % size); 2091 PetscInt newO, lSize, oTrue; 2092 2093 q = iperm[cone[m]]; 2094 newCones[offset] = Kembedding[q]; 2095 ierr = DMPlexGetConeSize(K,q,&lSize);CHKERRQ(ierr); 2096 oTrue = orientation[m]; 2097 oTrue = ((!lSize) || (preOrient[k] >= 0)) ? oTrue : -(oTrue + 2); 2098 newO = DihedralCompose(lSize,oTrue,preOrient[q]); 2099 newOrientations[offset++] = newO; 2100 } 2101 if (kParent != 0) { 2102 PetscInt newPoint = Kembedding[kParent]; 2103 ierr = PetscSectionGetOffset(parentSection,Kembedding[k],&pOffset);CHKERRQ(ierr); 2104 parents[pOffset] = newPoint; 2105 childIDs[pOffset] = k; 2106 } 2107 } 2108 } 2109 } 2110 2111 ierr = PetscMalloc1(dim*(pNewEnd[dim]-pNewStart[dim]),&newVertexCoords);CHKERRQ(ierr); 2112 2113 /* fill coordinates */ 2114 offset = 0; 2115 { 2116 PetscInt kStart, kEnd, l; 2117 PetscSection vSection; 2118 PetscInt v; 2119 Vec coords; 2120 PetscScalar *coordvals; 2121 PetscInt dof, off; 2122 PetscReal v0[3], J[9], detJ; 2123 2124 #if defined(PETSC_USE_DEBUG) 2125 { 2126 PetscInt k; 2127 ierr = DMPlexGetHeightStratum(K,0,&kStart,&kEnd);CHKERRQ(ierr); 2128 for (k = kStart; k < kEnd; k++) { 2129 ierr = DMPlexComputeCellGeometryFEM(K, k, NULL, v0, J, NULL, &detJ);CHKERRQ(ierr); 2130 if (detJ <= 0.) SETERRQ1 (PETSC_COMM_SELF,PETSC_ERR_PLIB,"reference tree cell %d has bad determinant",k); 2131 } 2132 } 2133 #endif 2134 ierr = DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, NULL, &detJ);CHKERRQ(ierr); 2135 ierr = DMGetCoordinateSection(dm,&vSection);CHKERRQ(ierr); 2136 ierr = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr); 2137 ierr = VecGetArray(coords,&coordvals);CHKERRQ(ierr); 2138 for (v = pOldStart[dim]; v < pOldEnd[dim]; v++) { 2139 2140 ierr = PetscSectionGetDof(vSection,v,&dof);CHKERRQ(ierr); 2141 ierr = PetscSectionGetOffset(vSection,v,&off);CHKERRQ(ierr); 2142 for (l = 0; l < dof; l++) { 2143 newVertexCoords[offset++] = coordvals[off + l]; 2144 } 2145 } 2146 ierr = VecRestoreArray(coords,&coordvals);CHKERRQ(ierr); 2147 2148 ierr = DMGetCoordinateSection(K,&vSection);CHKERRQ(ierr); 2149 ierr = DMGetCoordinatesLocal(K,&coords);CHKERRQ(ierr); 2150 ierr = VecGetArray(coords,&coordvals);CHKERRQ(ierr); 2151 ierr = DMPlexGetDepthStratum(K,0,&kStart,&kEnd);CHKERRQ(ierr); 2152 for (v = kStart; v < kEnd; v++) { 2153 PetscReal coord[3], newCoord[3]; 2154 PetscInt vPerm = perm[v]; 2155 PetscInt kParent; 2156 2157 ierr = DMPlexGetTreeParent(K,v,&kParent,NULL);CHKERRQ(ierr); 2158 if (kParent != v) { 2159 /* this is a new vertex */ 2160 ierr = PetscSectionGetOffset(vSection,vPerm,&off);CHKERRQ(ierr); 2161 for (l = 0; l < dim; ++l) coord[l] = PetscRealPart(coordvals[off+l]); 2162 CoordinatesRefToReal(dim, dim, v0, J, coord, newCoord);CHKERRQ(ierr); 2163 for (l = 0; l < dim; ++l) newVertexCoords[offset+l] = newCoord[l]; 2164 offset += dim; 2165 } 2166 } 2167 ierr = VecRestoreArray(coords,&coordvals);CHKERRQ(ierr); 2168 } 2169 2170 /* need to reverse the order of pNewCount: vertices first, cells last */ 2171 for (d = 0; d < (dim + 1) / 2; d++) { 2172 PetscInt tmp; 2173 2174 tmp = pNewCount[d]; 2175 pNewCount[d] = pNewCount[dim - d]; 2176 pNewCount[dim - d] = tmp; 2177 } 2178 2179 ierr = DMPlexCreateFromDAG(*ncdm,dim,pNewCount,newConeSizes,newCones,newOrientations,newVertexCoords);CHKERRQ(ierr); 2180 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr); 2181 ierr = DMPlexSetTree(*ncdm,parentSection,parents,childIDs);CHKERRQ(ierr); 2182 2183 /* clean up */ 2184 ierr = DMPlexRestoreTransitiveClosure(dm,cell,PETSC_TRUE,&nc,&cellClosure);CHKERRQ(ierr); 2185 ierr = PetscFree5(pNewCount,pNewStart,pNewEnd,pOldStart,pOldEnd);CHKERRQ(ierr); 2186 ierr = PetscFree(newConeSizes);CHKERRQ(ierr); 2187 ierr = PetscFree2(newCones,newOrientations);CHKERRQ(ierr); 2188 ierr = PetscFree(newVertexCoords);CHKERRQ(ierr); 2189 ierr = PetscFree2(parents,childIDs);CHKERRQ(ierr); 2190 ierr = PetscFree4(Kembedding,perm,iperm,preOrient);CHKERRQ(ierr); 2191 } 2192 else { 2193 PetscInt p, counts[4]; 2194 PetscInt *coneSizes, *cones, *orientations; 2195 Vec coordVec; 2196 PetscScalar *coords; 2197 2198 for (d = 0; d <= dim; d++) { 2199 PetscInt dStart, dEnd; 2200 2201 ierr = DMPlexGetDepthStratum(dm,d,&dStart,&dEnd);CHKERRQ(ierr); 2202 counts[d] = dEnd - dStart; 2203 } 2204 ierr = PetscMalloc1(pEnd-pStart,&coneSizes);CHKERRQ(ierr); 2205 for (p = pStart; p < pEnd; p++) { 2206 ierr = DMPlexGetConeSize(dm,p,&coneSizes[p-pStart]);CHKERRQ(ierr); 2207 } 2208 ierr = DMPlexGetCones(dm, &cones);CHKERRQ(ierr); 2209 ierr = DMPlexGetConeOrientations(dm, &orientations);CHKERRQ(ierr); 2210 ierr = DMGetCoordinatesLocal(dm,&coordVec);CHKERRQ(ierr); 2211 ierr = VecGetArray(coordVec,&coords);CHKERRQ(ierr); 2212 2213 ierr = PetscSectionSetChart(parentSection,pStart,pEnd);CHKERRQ(ierr); 2214 ierr = PetscSectionSetUp(parentSection);CHKERRQ(ierr); 2215 ierr = DMPlexCreateFromDAG(*ncdm,dim,counts,coneSizes,cones,orientations,NULL);CHKERRQ(ierr); 2216 ierr = DMPlexSetReferenceTree(*ncdm,K);CHKERRQ(ierr); 2217 ierr = DMPlexSetTree(*ncdm,parentSection,NULL,NULL);CHKERRQ(ierr); 2218 ierr = VecRestoreArray(coordVec,&coords);CHKERRQ(ierr); 2219 } 2220 ierr = PetscSectionDestroy(&parentSection);CHKERRQ(ierr); 2221 2222 PetscFunctionReturn(0); 2223 } 2224 2225 #undef __FUNCT__ 2226 #define __FUNCT__ "DMPlexComputeInterpolatorTree" 2227 PetscErrorCode DMPlexComputeInterpolatorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 2228 { 2229 PetscSF coarseToFineEmbedded; 2230 PetscSection globalCoarse, globalFine; 2231 PetscSection localCoarse, localFine; 2232 PetscSection aSec, cSec; 2233 PetscSection rootIndicesSec, rootMatricesSec; 2234 PetscSection leafIndicesSec, leafMatricesSec; 2235 PetscInt *rootIndices, *leafIndices; 2236 PetscScalar *rootMatrices, *leafMatrices; 2237 IS aIS; 2238 const PetscInt *anchors; 2239 Mat cMat; 2240 PetscInt numFields, maxFields; 2241 PetscInt pStartC, pEndC, pStartF, pEndF, p; 2242 PetscInt aStart, aEnd, cStart, cEnd; 2243 PetscInt *maxChildIds; 2244 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 2245 const PetscInt ***perms; 2246 const PetscScalar ***flips; 2247 PetscErrorCode ierr; 2248 2249 PetscFunctionBegin; 2250 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 2251 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 2252 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 2253 { /* winnow fine points that don't have global dofs out of the sf */ 2254 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, nleaves, l; 2255 const PetscInt *leaves; 2256 2257 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 2258 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 2259 p = leaves ? leaves[l] : l; 2260 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 2261 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 2262 if ((dof - cdof) > 0) { 2263 numPointsWithDofs++; 2264 } 2265 } 2266 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 2267 for (l = 0, offset = 0; l < nleaves; l++) { 2268 p = leaves ? leaves[l] : l; 2269 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 2270 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 2271 if ((dof - cdof) > 0) { 2272 pointsWithDofs[offset++] = l; 2273 } 2274 } 2275 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 2276 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 2277 } 2278 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 2279 ierr = PetscMalloc1(pEndC-pStartC,&maxChildIds);CHKERRQ(ierr); 2280 for (p = pStartC; p < pEndC; p++) { 2281 maxChildIds[p - pStartC] = -2; 2282 } 2283 ierr = PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 2284 ierr = PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,childIds,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 2285 2286 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 2287 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 2288 2289 ierr = DMPlexGetAnchors(coarse,&aSec,&aIS);CHKERRQ(ierr); 2290 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 2291 ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr); 2292 2293 ierr = DMGetDefaultConstraints(coarse,&cSec,&cMat);CHKERRQ(ierr); 2294 ierr = PetscSectionGetChart(cSec,&cStart,&cEnd);CHKERRQ(ierr); 2295 2296 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 2297 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);CHKERRQ(ierr); 2298 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootMatricesSec);CHKERRQ(ierr); 2299 ierr = PetscSectionSetChart(rootIndicesSec,pStartC,pEndC);CHKERRQ(ierr); 2300 ierr = PetscSectionSetChart(rootMatricesSec,pStartC,pEndC);CHKERRQ(ierr); 2301 ierr = PetscSectionGetNumFields(localCoarse,&numFields);CHKERRQ(ierr); 2302 maxFields = PetscMax(1,numFields); 2303 ierr = PetscMalloc7(maxFields+1,&offsets,maxFields+1,&offsetsCopy,maxFields+1,&newOffsets,maxFields+1,&newOffsetsCopy,maxFields+1,&rowOffsets,maxFields+1,&numD,maxFields+1,&numO);CHKERRQ(ierr); 2304 ierr = PetscMalloc2(maxFields+1,&perms,maxFields+1,&flips);CHKERRQ(ierr); 2305 ierr = PetscMemzero((void *) perms, (maxFields+1) * sizeof(const PetscInt **));CHKERRQ(ierr); 2306 ierr = PetscMemzero((void *) flips, (maxFields+1) * sizeof(const PetscScalar **));CHKERRQ(ierr); 2307 2308 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 2309 PetscInt dof, matSize = 0; 2310 PetscInt aDof = 0; 2311 PetscInt cDof = 0; 2312 PetscInt maxChildId = maxChildIds[p - pStartC]; 2313 PetscInt numRowIndices = 0; 2314 PetscInt numColIndices = 0; 2315 PetscInt f; 2316 2317 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 2318 if (dof < 0) { 2319 dof = -(dof + 1); 2320 } 2321 if (p >= aStart && p < aEnd) { 2322 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 2323 } 2324 if (p >= cStart && p < cEnd) { 2325 ierr = PetscSectionGetDof(cSec,p,&cDof);CHKERRQ(ierr); 2326 } 2327 for (f = 0; f <= numFields; f++) offsets[f] = 0; 2328 for (f = 0; f <= numFields; f++) newOffsets[f] = 0; 2329 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 2330 PetscInt *closure = NULL, closureSize, cl; 2331 2332 ierr = DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2333 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 2334 PetscInt c = closure[2 * cl], clDof; 2335 2336 ierr = PetscSectionGetDof(localCoarse,c,&clDof);CHKERRQ(ierr); 2337 numRowIndices += clDof; 2338 for (f = 0; f < numFields; f++) { 2339 ierr = PetscSectionGetFieldDof(localCoarse,c,f,&clDof);CHKERRQ(ierr); 2340 offsets[f + 1] += clDof; 2341 } 2342 } 2343 for (f = 0; f < numFields; f++) { 2344 offsets[f + 1] += offsets[f]; 2345 newOffsets[f + 1] = offsets[f + 1]; 2346 } 2347 /* get the number of indices needed and their field offsets */ 2348 ierr = DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,NULL,NULL,NULL,&numColIndices,NULL,NULL,newOffsets,PETSC_FALSE);CHKERRQ(ierr); 2349 ierr = DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2350 if (!numColIndices) { /* there are no hanging constraint modifications, so the matrix is just the identity: do not send it */ 2351 numColIndices = numRowIndices; 2352 matSize = 0; 2353 } 2354 else if (numFields) { /* we send one submat for each field: sum their sizes */ 2355 matSize = 0; 2356 for (f = 0; f < numFields; f++) { 2357 PetscInt numRow, numCol; 2358 2359 numRow = offsets[f + 1] - offsets[f]; 2360 numCol = newOffsets[f + 1] - newOffsets[f]; 2361 matSize += numRow * numCol; 2362 } 2363 } 2364 else { 2365 matSize = numRowIndices * numColIndices; 2366 } 2367 } else if (maxChildId == -1) { 2368 if (cDof > 0) { /* this point's dofs are interpolated via cMat: get the submatrix of cMat */ 2369 PetscInt aOff, a; 2370 2371 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 2372 for (f = 0; f < numFields; f++) { 2373 PetscInt fDof; 2374 2375 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 2376 offsets[f+1] = fDof; 2377 } 2378 for (a = 0; a < aDof; a++) { 2379 PetscInt anchor = anchors[a + aOff], aLocalDof; 2380 2381 ierr = PetscSectionGetDof(localCoarse,anchor,&aLocalDof);CHKERRQ(ierr); 2382 numColIndices += aLocalDof; 2383 for (f = 0; f < numFields; f++) { 2384 PetscInt fDof; 2385 2386 ierr = PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);CHKERRQ(ierr); 2387 newOffsets[f+1] += fDof; 2388 } 2389 } 2390 if (numFields) { 2391 matSize = 0; 2392 for (f = 0; f < numFields; f++) { 2393 matSize += offsets[f+1] * newOffsets[f+1]; 2394 } 2395 } 2396 else { 2397 matSize = numColIndices * dof; 2398 } 2399 } 2400 else { /* no children, and no constraints on dofs: just get the global indices */ 2401 numColIndices = dof; 2402 matSize = 0; 2403 } 2404 } 2405 /* we will pack the column indices with the field offsets */ 2406 ierr = PetscSectionSetDof(rootIndicesSec,p,numColIndices ? numColIndices+2*numFields : 0);CHKERRQ(ierr); 2407 ierr = PetscSectionSetDof(rootMatricesSec,p,matSize);CHKERRQ(ierr); 2408 } 2409 ierr = PetscSectionSetUp(rootIndicesSec);CHKERRQ(ierr); 2410 ierr = PetscSectionSetUp(rootMatricesSec);CHKERRQ(ierr); 2411 { 2412 PetscInt numRootIndices, numRootMatrices; 2413 2414 ierr = PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);CHKERRQ(ierr); 2415 ierr = PetscSectionGetStorageSize(rootMatricesSec,&numRootMatrices);CHKERRQ(ierr); 2416 ierr = PetscMalloc2(numRootIndices,&rootIndices,numRootMatrices,&rootMatrices);CHKERRQ(ierr); 2417 for (p = pStartC; p < pEndC; p++) { 2418 PetscInt numRowIndices, numColIndices, matSize, dof; 2419 PetscInt pIndOff, pMatOff, f; 2420 PetscInt *pInd; 2421 PetscInt maxChildId = maxChildIds[p - pStartC]; 2422 PetscScalar *pMat = NULL; 2423 2424 ierr = PetscSectionGetDof(rootIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2425 if (!numColIndices) { 2426 continue; 2427 } 2428 for (f = 0; f <= numFields; f++) { 2429 offsets[f] = 0; 2430 newOffsets[f] = 0; 2431 offsetsCopy[f] = 0; 2432 newOffsetsCopy[f] = 0; 2433 } 2434 numColIndices -= 2 * numFields; 2435 ierr = PetscSectionGetOffset(rootIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2436 pInd = &(rootIndices[pIndOff]); 2437 ierr = PetscSectionGetDof(rootMatricesSec,p,&matSize);CHKERRQ(ierr); 2438 if (matSize) { 2439 ierr = PetscSectionGetOffset(rootMatricesSec,p,&pMatOff);CHKERRQ(ierr); 2440 pMat = &rootMatrices[pMatOff]; 2441 } 2442 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 2443 if (dof < 0) { 2444 dof = -(dof + 1); 2445 } 2446 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 2447 PetscInt i, j; 2448 PetscInt numRowIndices = matSize / numColIndices; 2449 2450 if (!numRowIndices) { /* don't need to calculate the mat, just the indices */ 2451 PetscInt numIndices, *indices; 2452 ierr = DMPlexGetClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);CHKERRQ(ierr); 2453 if (numIndices != numColIndices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"mismatching constraint indices calculations"); 2454 for (i = 0; i < numColIndices; i++) { 2455 pInd[i] = indices[i]; 2456 } 2457 for (i = 0; i < numFields; i++) { 2458 pInd[numColIndices + i] = offsets[i+1]; 2459 pInd[numColIndices + numFields + i] = offsets[i+1]; 2460 } 2461 ierr = DMPlexRestoreClosureIndices(coarse,localCoarse,globalCoarse,p,&numIndices,&indices,offsets);CHKERRQ(ierr); 2462 } 2463 else { 2464 PetscInt closureSize, *closure = NULL, cl; 2465 PetscScalar *pMatIn, *pMatModified; 2466 PetscInt numPoints,*points; 2467 2468 ierr = DMGetWorkArray(coarse,numRowIndices * numRowIndices,PETSC_SCALAR,&pMatIn);CHKERRQ(ierr); 2469 for (i = 0; i < numRowIndices; i++) { /* initialize to the identity */ 2470 for (j = 0; j < numRowIndices; j++) { 2471 pMatIn[i * numRowIndices + j] = (i == j) ? 1. : 0.; 2472 } 2473 } 2474 ierr = DMPlexGetTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2475 for (f = 0; f < maxFields; f++) { 2476 if (numFields) {ierr = PetscSectionGetFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2477 else {ierr = PetscSectionGetPointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2478 } 2479 if (numFields) { 2480 for (cl = 0; cl < closureSize; cl++) { 2481 PetscInt c = closure[2 * cl]; 2482 2483 for (f = 0; f < numFields; f++) { 2484 PetscInt fDof; 2485 2486 ierr = PetscSectionGetFieldDof(localCoarse,c,f,&fDof);CHKERRQ(ierr); 2487 offsets[f + 1] += fDof; 2488 } 2489 } 2490 for (f = 0; f < numFields; f++) { 2491 offsets[f + 1] += offsets[f]; 2492 newOffsets[f + 1] = offsets[f + 1]; 2493 } 2494 } 2495 /* TODO : flips here ? */ 2496 /* apply hanging node constraints on the right, get the new points and the new offsets */ 2497 ierr = DMPlexAnchorsModifyMat(coarse,localCoarse,closureSize,numRowIndices,closure,perms,pMatIn,&numPoints,NULL,&points,&pMatModified,newOffsets,PETSC_FALSE);CHKERRQ(ierr); 2498 for (f = 0; f < maxFields; f++) { 2499 if (numFields) {ierr = PetscSectionRestoreFieldPointSyms(localCoarse,f,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2500 else {ierr = PetscSectionRestorePointSyms(localCoarse,closureSize,closure,&perms[f],&flips[f]);CHKERRQ(ierr);} 2501 } 2502 for (f = 0; f < maxFields; f++) { 2503 if (numFields) {ierr = PetscSectionGetFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2504 else {ierr = PetscSectionGetPointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2505 } 2506 if (!numFields) { 2507 for (i = 0; i < numRowIndices * numColIndices; i++) { 2508 pMat[i] = pMatModified[i]; 2509 } 2510 } 2511 else { 2512 PetscInt i, j, count; 2513 for (f = 0, count = 0; f < numFields; f++) { 2514 for (i = offsets[f]; i < offsets[f+1]; i++) { 2515 for (j = newOffsets[f]; j < newOffsets[f+1]; j++, count++) { 2516 pMat[count] = pMatModified[i * numColIndices + j]; 2517 } 2518 } 2519 } 2520 } 2521 ierr = DMRestoreWorkArray(coarse,numRowIndices * numColIndices,PETSC_SCALAR,&pMatModified);CHKERRQ(ierr); 2522 ierr = DMPlexRestoreTransitiveClosure(coarse, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 2523 ierr = DMRestoreWorkArray(coarse,numRowIndices * numColIndices,PETSC_SCALAR,&pMatIn);CHKERRQ(ierr); 2524 if (numFields) { 2525 for (f = 0; f < numFields; f++) { 2526 pInd[numColIndices + f] = offsets[f+1]; 2527 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2528 } 2529 for (cl = 0; cl < numPoints; cl++) { 2530 PetscInt globalOff, c = points[2*cl]; 2531 ierr = PetscSectionGetOffset(globalCoarse, c, &globalOff);CHKERRQ(ierr); 2532 DMPlexGetIndicesPointFields_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perms, cl, pInd); 2533 } 2534 } else { 2535 for (cl = 0; cl < numPoints; cl++) { 2536 PetscInt c = points[2*cl], globalOff; 2537 const PetscInt *perm = perms[0] ? perms[0][cl] : NULL; 2538 2539 ierr = PetscSectionGetOffset(globalCoarse, c, &globalOff);CHKERRQ(ierr); 2540 DMPlexGetIndicesPoint_Internal(localCoarse, c, globalOff < 0 ? -(globalOff+1) : globalOff, newOffsets, PETSC_FALSE, perm, pInd); 2541 } 2542 } 2543 for (f = 0; f < maxFields; f++) { 2544 if (numFields) {ierr = PetscSectionRestoreFieldPointSyms(localCoarse,f,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2545 else {ierr = PetscSectionRestorePointSyms(localCoarse,numPoints,points,&perms[f],&flips[f]);CHKERRQ(ierr);} 2546 } 2547 ierr = DMRestoreWorkArray(coarse,numPoints,PETSC_SCALAR,&points);CHKERRQ(ierr); 2548 } 2549 } 2550 else if (matSize) { 2551 PetscInt cOff; 2552 PetscInt *rowIndices, *colIndices, a, aDof, aOff; 2553 2554 numRowIndices = matSize / numColIndices; 2555 if (numRowIndices != dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Miscounted dofs"); 2556 ierr = DMGetWorkArray(coarse,numRowIndices,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2557 ierr = DMGetWorkArray(coarse,numColIndices,PETSC_INT,&colIndices);CHKERRQ(ierr); 2558 ierr = PetscSectionGetOffset(cSec,p,&cOff);CHKERRQ(ierr); 2559 ierr = PetscSectionGetDof(aSec,p,&aDof);CHKERRQ(ierr); 2560 ierr = PetscSectionGetOffset(aSec,p,&aOff);CHKERRQ(ierr); 2561 if (numFields) { 2562 for (f = 0; f < numFields; f++) { 2563 PetscInt fDof; 2564 2565 ierr = PetscSectionGetFieldDof(cSec,p,f,&fDof);CHKERRQ(ierr); 2566 offsets[f + 1] = fDof; 2567 for (a = 0; a < aDof; a++) { 2568 PetscInt anchor = anchors[a + aOff]; 2569 ierr = PetscSectionGetFieldDof(localCoarse,anchor,f,&fDof);CHKERRQ(ierr); 2570 newOffsets[f + 1] += fDof; 2571 } 2572 } 2573 for (f = 0; f < numFields; f++) { 2574 offsets[f + 1] += offsets[f]; 2575 offsetsCopy[f + 1] = offsets[f + 1]; 2576 newOffsets[f + 1] += newOffsets[f]; 2577 newOffsetsCopy[f + 1] = newOffsets[f + 1]; 2578 } 2579 DMPlexGetIndicesPointFields_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,-1,rowIndices);CHKERRQ(ierr); 2580 for (a = 0; a < aDof; a++) { 2581 PetscInt anchor = anchors[a + aOff], lOff; 2582 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr); 2583 DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,-1,colIndices);CHKERRQ(ierr); 2584 } 2585 } 2586 else { 2587 DMPlexGetIndicesPoint_Internal(cSec,p,cOff,offsetsCopy,PETSC_TRUE,NULL,rowIndices);CHKERRQ(ierr); 2588 for (a = 0; a < aDof; a++) { 2589 PetscInt anchor = anchors[a + aOff], lOff; 2590 ierr = PetscSectionGetOffset(localCoarse,anchor,&lOff);CHKERRQ(ierr); 2591 DMPlexGetIndicesPoint_Internal(localCoarse,anchor,lOff,newOffsetsCopy,PETSC_TRUE,NULL,colIndices);CHKERRQ(ierr); 2592 } 2593 } 2594 if (numFields) { 2595 PetscInt count, a; 2596 2597 for (f = 0, count = 0; f < numFields; f++) { 2598 PetscInt iSize = offsets[f + 1] - offsets[f]; 2599 PetscInt jSize = newOffsets[f + 1] - newOffsets[f]; 2600 ierr = MatGetValues(cMat,iSize,&rowIndices[offsets[f]],jSize,&colIndices[newOffsets[f]],&pMat[count]);CHKERRQ(ierr); 2601 count += iSize * jSize; 2602 pInd[numColIndices + f] = offsets[f+1]; 2603 pInd[numColIndices + numFields + f] = newOffsets[f+1]; 2604 } 2605 for (a = 0; a < aDof; a++) { 2606 PetscInt anchor = anchors[a + aOff]; 2607 PetscInt gOff; 2608 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr); 2609 DMPlexGetIndicesPointFields_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 2610 } 2611 } 2612 else { 2613 PetscInt a; 2614 ierr = MatGetValues(cMat,numRowIndices,rowIndices,numColIndices,colIndices,pMat);CHKERRQ(ierr); 2615 for (a = 0; a < aDof; a++) { 2616 PetscInt anchor = anchors[a + aOff]; 2617 PetscInt gOff; 2618 ierr = PetscSectionGetOffset(globalCoarse,anchor,&gOff);CHKERRQ(ierr); 2619 DMPlexGetIndicesPoint_Internal(localCoarse,anchor,gOff < 0 ? -(gOff + 1) : gOff,newOffsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 2620 } 2621 } 2622 ierr = DMRestoreWorkArray(coarse,numColIndices,PETSC_INT,&colIndices);CHKERRQ(ierr); 2623 ierr = DMRestoreWorkArray(coarse,numRowIndices,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2624 } 2625 else { 2626 PetscInt gOff; 2627 2628 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 2629 if (numFields) { 2630 for (f = 0; f < numFields; f++) { 2631 PetscInt fDof; 2632 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 2633 offsets[f + 1] = fDof + offsets[f]; 2634 } 2635 for (f = 0; f < numFields; f++) { 2636 pInd[numColIndices + f] = offsets[f+1]; 2637 pInd[numColIndices + numFields + f] = offsets[f+1]; 2638 } 2639 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 2640 } 2641 else { 2642 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 2643 } 2644 } 2645 } 2646 ierr = PetscFree(maxChildIds);CHKERRQ(ierr); 2647 } 2648 { 2649 PetscSF indicesSF, matricesSF; 2650 PetscInt *remoteOffsetsIndices, *remoteOffsetsMatrices, numLeafIndices, numLeafMatrices; 2651 2652 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);CHKERRQ(ierr); 2653 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafMatricesSec);CHKERRQ(ierr); 2654 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootIndicesSec,&remoteOffsetsIndices,leafIndicesSec);CHKERRQ(ierr); 2655 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootMatricesSec,&remoteOffsetsMatrices,leafMatricesSec);CHKERRQ(ierr); 2656 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootIndicesSec,remoteOffsetsIndices,leafIndicesSec,&indicesSF);CHKERRQ(ierr); 2657 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootMatricesSec,remoteOffsetsMatrices,leafMatricesSec,&matricesSF);CHKERRQ(ierr); 2658 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 2659 ierr = PetscFree(remoteOffsetsIndices);CHKERRQ(ierr); 2660 ierr = PetscFree(remoteOffsetsMatrices);CHKERRQ(ierr); 2661 ierr = PetscSectionGetStorageSize(leafIndicesSec,&numLeafIndices);CHKERRQ(ierr); 2662 ierr = PetscSectionGetStorageSize(leafMatricesSec,&numLeafMatrices);CHKERRQ(ierr); 2663 ierr = PetscMalloc2(numLeafIndices,&leafIndices,numLeafMatrices,&leafMatrices);CHKERRQ(ierr); 2664 ierr = PetscSFBcastBegin(indicesSF,MPIU_INT,rootIndices,leafIndices);CHKERRQ(ierr); 2665 ierr = PetscSFBcastBegin(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);CHKERRQ(ierr); 2666 ierr = PetscSFBcastEnd(indicesSF,MPIU_INT,rootIndices,leafIndices);CHKERRQ(ierr); 2667 ierr = PetscSFBcastEnd(matricesSF,MPIU_SCALAR,rootMatrices,leafMatrices);CHKERRQ(ierr); 2668 ierr = PetscSFDestroy(&matricesSF);CHKERRQ(ierr); 2669 ierr = PetscSFDestroy(&indicesSF);CHKERRQ(ierr); 2670 ierr = PetscFree2(rootIndices,rootMatrices);CHKERRQ(ierr); 2671 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 2672 ierr = PetscSectionDestroy(&rootMatricesSec);CHKERRQ(ierr); 2673 } 2674 /* count to preallocate */ 2675 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 2676 { 2677 PetscInt nGlobal; 2678 PetscInt *dnnz, *onnz; 2679 PetscLayout rowMap, colMap; 2680 PetscInt rowStart, rowEnd, colStart, colEnd; 2681 PetscInt maxDof; 2682 PetscInt *rowIndices; 2683 DM refTree; 2684 PetscInt **refPointFieldN; 2685 PetscScalar ***refPointFieldMats; 2686 PetscSection refConSec, refAnSec; 2687 PetscInt pRefStart,pRefEnd,maxConDof,maxColumns,leafStart,leafEnd; 2688 PetscScalar *pointWork; 2689 2690 ierr = PetscSectionGetConstrainedStorageSize(globalFine,&nGlobal);CHKERRQ(ierr); 2691 ierr = PetscCalloc2(nGlobal,&dnnz,nGlobal,&onnz);CHKERRQ(ierr); 2692 ierr = MatGetLayouts(mat,&rowMap,&colMap);CHKERRQ(ierr); 2693 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 2694 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 2695 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 2696 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 2697 ierr = PetscSectionGetMaxDof(globalFine,&maxDof);CHKERRQ(ierr); 2698 ierr = PetscSectionGetChart(leafIndicesSec,&leafStart,&leafEnd);CHKERRQ(ierr); 2699 ierr = DMGetWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 2700 for (p = leafStart; p < leafEnd; p++) { 2701 PetscInt gDof, gcDof, gOff; 2702 PetscInt numColIndices, pIndOff, *pInd; 2703 PetscInt matSize; 2704 PetscInt i; 2705 2706 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 2707 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 2708 if ((gDof - gcDof) <= 0) { 2709 continue; 2710 } 2711 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 2712 if (gOff < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I though having global dofs meant a non-negative offset"); 2713 if ((gOff < rowStart) || ((gOff + gDof - gcDof) > rowEnd)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"I thought the row map would constrain the global dofs"); 2714 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2715 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2716 numColIndices -= 2 * numFields; 2717 if (numColIndices <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"global fine dof with no dofs to interpolate from"); 2718 pInd = &leafIndices[pIndOff]; 2719 offsets[0] = 0; 2720 offsetsCopy[0] = 0; 2721 newOffsets[0] = 0; 2722 newOffsetsCopy[0] = 0; 2723 if (numFields) { 2724 PetscInt f; 2725 for (f = 0; f < numFields; f++) { 2726 PetscInt rowDof; 2727 2728 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 2729 offsets[f + 1] = offsets[f] + rowDof; 2730 offsetsCopy[f + 1] = offsets[f + 1]; 2731 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2732 numD[f] = 0; 2733 numO[f] = 0; 2734 } 2735 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 2736 for (f = 0; f < numFields; f++) { 2737 PetscInt colOffset = newOffsets[f]; 2738 PetscInt numFieldCols = newOffsets[f + 1] - newOffsets[f]; 2739 2740 for (i = 0; i < numFieldCols; i++) { 2741 PetscInt gInd = pInd[i + colOffset]; 2742 2743 if (gInd >= colStart && gInd < colEnd) { 2744 numD[f]++; 2745 } 2746 else if (gInd >= 0) { /* negative means non-entry */ 2747 numO[f]++; 2748 } 2749 } 2750 } 2751 } 2752 else { 2753 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 2754 numD[0] = 0; 2755 numO[0] = 0; 2756 for (i = 0; i < numColIndices; i++) { 2757 PetscInt gInd = pInd[i]; 2758 2759 if (gInd >= colStart && gInd < colEnd) { 2760 numD[0]++; 2761 } 2762 else if (gInd >= 0) { /* negative means non-entry */ 2763 numO[0]++; 2764 } 2765 } 2766 } 2767 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr); 2768 if (!matSize) { /* incoming matrix is identity */ 2769 PetscInt childId; 2770 2771 childId = childIds[p-pStartF]; 2772 if (childId < 0) { /* no child interpolation: one nnz per */ 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 gIndCoarse = pInd[newOffsets[f] + row]; 2779 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2780 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2781 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2782 dnnz[gIndFine - rowStart] = 1; 2783 } 2784 else if (gIndCoarse >= 0) { /* remote */ 2785 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2786 onnz[gIndFine - rowStart] = 1; 2787 } 2788 else { /* constrained */ 2789 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2790 } 2791 } 2792 } 2793 } 2794 else { 2795 PetscInt i; 2796 for (i = 0; i < gDof; i++) { 2797 PetscInt gIndCoarse = pInd[i]; 2798 PetscInt gIndFine = rowIndices[i]; 2799 if (gIndCoarse >= colStart && gIndCoarse < colEnd) { /* local */ 2800 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2801 dnnz[gIndFine - rowStart] = 1; 2802 } 2803 else if (gIndCoarse >= 0) { /* remote */ 2804 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2805 onnz[gIndFine - rowStart] = 1; 2806 } 2807 else { /* constrained */ 2808 if (gIndFine >= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2809 } 2810 } 2811 } 2812 } 2813 else { /* interpolate from all */ 2814 if (numFields) { 2815 PetscInt f; 2816 for (f = 0; f < numFields; f++) { 2817 PetscInt numRows = offsets[f+1] - offsets[f], row; 2818 for (row = 0; row < numRows; row++) { 2819 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2820 if (gIndFine >= 0) { 2821 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2822 dnnz[gIndFine - rowStart] = numD[f]; 2823 onnz[gIndFine - rowStart] = numO[f]; 2824 } 2825 } 2826 } 2827 } 2828 else { 2829 PetscInt i; 2830 for (i = 0; i < gDof; i++) { 2831 PetscInt gIndFine = rowIndices[i]; 2832 if (gIndFine >= 0) { 2833 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2834 dnnz[gIndFine - rowStart] = numD[0]; 2835 onnz[gIndFine - rowStart] = numO[0]; 2836 } 2837 } 2838 } 2839 } 2840 } 2841 else { /* interpolate from all */ 2842 if (numFields) { 2843 PetscInt f; 2844 for (f = 0; f < numFields; f++) { 2845 PetscInt numRows = offsets[f+1] - offsets[f], row; 2846 for (row = 0; row < numRows; row++) { 2847 PetscInt gIndFine = rowIndices[offsets[f] + row]; 2848 if (gIndFine >= 0) { 2849 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2850 dnnz[gIndFine - rowStart] = numD[f]; 2851 onnz[gIndFine - rowStart] = numO[f]; 2852 } 2853 } 2854 } 2855 } 2856 else { /* every dof get a full row */ 2857 PetscInt i; 2858 for (i = 0; i < gDof; i++) { 2859 PetscInt gIndFine = rowIndices[i]; 2860 if (gIndFine >= 0) { 2861 if (gIndFine < rowStart || gIndFine >= rowEnd) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatched number of constrained dofs"); 2862 dnnz[gIndFine - rowStart] = numD[0]; 2863 onnz[gIndFine - rowStart] = numO[0]; 2864 } 2865 } 2866 } 2867 } 2868 } 2869 ierr = MatXAIJSetPreallocation(mat,1,dnnz,onnz,NULL,NULL);CHKERRQ(ierr); 2870 ierr = PetscFree2(dnnz,onnz);CHKERRQ(ierr); 2871 2872 ierr = DMPlexGetReferenceTree(fine,&refTree);CHKERRQ(ierr); 2873 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 2874 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 2875 ierr = DMPlexGetAnchors(refTree,&refAnSec,NULL);CHKERRQ(ierr); 2876 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 2877 ierr = PetscSectionGetMaxDof(refConSec,&maxConDof);CHKERRQ(ierr); 2878 ierr = PetscSectionGetMaxDof(leafIndicesSec,&maxColumns);CHKERRQ(ierr); 2879 ierr = PetscMalloc1(maxConDof*maxColumns,&pointWork);CHKERRQ(ierr); 2880 for (p = leafStart; p < leafEnd; p++) { 2881 PetscInt gDof, gcDof, gOff; 2882 PetscInt numColIndices, pIndOff, *pInd; 2883 PetscInt matSize; 2884 PetscInt childId; 2885 2886 2887 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 2888 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 2889 if ((gDof - gcDof) <= 0) { 2890 continue; 2891 } 2892 childId = childIds[p-pStartF]; 2893 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 2894 ierr = PetscSectionGetDof(leafIndicesSec,p,&numColIndices);CHKERRQ(ierr); 2895 ierr = PetscSectionGetOffset(leafIndicesSec,p,&pIndOff);CHKERRQ(ierr); 2896 numColIndices -= 2 * numFields; 2897 pInd = &leafIndices[pIndOff]; 2898 offsets[0] = 0; 2899 offsetsCopy[0] = 0; 2900 newOffsets[0] = 0; 2901 newOffsetsCopy[0] = 0; 2902 rowOffsets[0] = 0; 2903 if (numFields) { 2904 PetscInt f; 2905 for (f = 0; f < numFields; f++) { 2906 PetscInt rowDof; 2907 2908 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 2909 offsets[f + 1] = offsets[f] + rowDof; 2910 offsetsCopy[f + 1] = offsets[f + 1]; 2911 rowOffsets[f + 1] = pInd[numColIndices + f]; 2912 newOffsets[f + 1] = pInd[numColIndices + numFields + f]; 2913 } 2914 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 2915 } 2916 else { 2917 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 2918 } 2919 ierr = PetscSectionGetDof(leafMatricesSec,p,&matSize);CHKERRQ(ierr); 2920 if (!matSize) { /* incoming matrix is identity */ 2921 if (childId < 0) { /* no child interpolation: scatter */ 2922 if (numFields) { 2923 PetscInt f; 2924 for (f = 0; f < numFields; f++) { 2925 PetscInt numRows = offsets[f+1] - offsets[f], row; 2926 for (row = 0; row < numRows; row++) { 2927 ierr = MatSetValue(mat,rowIndices[offsets[f]+row],pInd[newOffsets[f]+row],1.,INSERT_VALUES);CHKERRQ(ierr); 2928 } 2929 } 2930 } 2931 else { 2932 PetscInt numRows = gDof, row; 2933 for (row = 0; row < numRows; row++) { 2934 ierr = MatSetValue(mat,rowIndices[row],pInd[row],1.,INSERT_VALUES);CHKERRQ(ierr); 2935 } 2936 } 2937 } 2938 else { /* interpolate from all */ 2939 if (numFields) { 2940 PetscInt f; 2941 for (f = 0; f < numFields; f++) { 2942 PetscInt numRows = offsets[f+1] - offsets[f]; 2943 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 2944 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],refPointFieldMats[childId - pRefStart][f],INSERT_VALUES);CHKERRQ(ierr); 2945 } 2946 } 2947 else { 2948 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,refPointFieldMats[childId - pRefStart][0],INSERT_VALUES);CHKERRQ(ierr); 2949 } 2950 } 2951 } 2952 else { /* interpolate from all */ 2953 PetscInt pMatOff; 2954 PetscScalar *pMat; 2955 2956 ierr = PetscSectionGetOffset(leafMatricesSec,p,&pMatOff);CHKERRQ(ierr); 2957 pMat = &leafMatrices[pMatOff]; 2958 if (childId < 0) { /* copy the incoming matrix */ 2959 if (numFields) { 2960 PetscInt f, count; 2961 for (f = 0, count = 0; f < numFields; f++) { 2962 PetscInt numRows = offsets[f+1]-offsets[f]; 2963 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2964 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2965 PetscScalar *inMat = &pMat[count]; 2966 2967 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],inMat,INSERT_VALUES);CHKERRQ(ierr); 2968 count += numCols * numInRows; 2969 } 2970 } 2971 else { 2972 ierr = MatSetValues(mat,gDof,rowIndices,numColIndices,pInd,pMat,INSERT_VALUES);CHKERRQ(ierr); 2973 } 2974 } 2975 else { /* multiply the incoming matrix by the child interpolation */ 2976 if (numFields) { 2977 PetscInt f, count; 2978 for (f = 0, count = 0; f < numFields; f++) { 2979 PetscInt numRows = offsets[f+1]-offsets[f]; 2980 PetscInt numCols = newOffsets[f+1]-newOffsets[f]; 2981 PetscInt numInRows = rowOffsets[f+1]-rowOffsets[f]; 2982 PetscScalar *inMat = &pMat[count]; 2983 PetscInt i, j, k; 2984 if (refPointFieldN[childId - pRefStart][f] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 2985 for (i = 0; i < numRows; i++) { 2986 for (j = 0; j < numCols; j++) { 2987 PetscScalar val = 0.; 2988 for (k = 0; k < numInRows; k++) { 2989 val += refPointFieldMats[childId - pRefStart][f][i * numInRows + k] * inMat[k * numCols + j]; 2990 } 2991 pointWork[i * numCols + j] = val; 2992 } 2993 } 2994 ierr = MatSetValues(mat,numRows,&rowIndices[offsets[f]],numCols,&pInd[newOffsets[f]],pointWork,INSERT_VALUES);CHKERRQ(ierr); 2995 count += numCols * numInRows; 2996 } 2997 } 2998 else { /* every dof gets a full row */ 2999 PetscInt numRows = gDof; 3000 PetscInt numCols = numColIndices; 3001 PetscInt numInRows = matSize / numColIndices; 3002 PetscInt i, j, k; 3003 if (refPointFieldN[childId - pRefStart][0] != numInRows) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point constraint matrix multiply dimension mismatch"); 3004 for (i = 0; i < numRows; i++) { 3005 for (j = 0; j < numCols; j++) { 3006 PetscScalar val = 0.; 3007 for (k = 0; k < numInRows; k++) { 3008 val += refPointFieldMats[childId - pRefStart][0][i * numInRows + k] * pMat[k * numCols + j]; 3009 } 3010 pointWork[i * numCols + j] = val; 3011 } 3012 } 3013 ierr = MatSetValues(mat,numRows,rowIndices,numCols,pInd,pointWork,INSERT_VALUES);CHKERRQ(ierr); 3014 } 3015 } 3016 } 3017 } 3018 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 3019 ierr = DMRestoreWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 3020 ierr = PetscFree(pointWork);CHKERRQ(ierr); 3021 } 3022 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3023 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3024 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr); 3025 ierr = PetscSectionDestroy(&leafMatricesSec);CHKERRQ(ierr); 3026 ierr = PetscFree2(leafIndices,leafMatrices);CHKERRQ(ierr); 3027 ierr = PetscFree2(perms,flips);CHKERRQ(ierr); 3028 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr); 3029 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 3030 PetscFunctionReturn(0); 3031 } 3032 3033 #undef __FUNCT__ 3034 #define __FUNCT__ "DMPlexComputeInjectorReferenceTree" 3035 /* 3036 * Assuming a nodal basis (w.r.t. the dual basis) basis: 3037 * 3038 * for each coarse dof \phi^c_i: 3039 * for each quadrature point (w_l,x_l) in the dual basis definition of \phi^c_i: 3040 * for each fine dof \phi^f_j; 3041 * a_{i,j} = 0; 3042 * for each fine dof \phi^f_k: 3043 * a_{i,j} += interp_{i,k} * \phi^f_k(x_l) * \phi^f_j(x_l) * w_l 3044 * [^^^ this is = \phi^c_i ^^^] 3045 */ 3046 PetscErrorCode DMPlexComputeInjectorReferenceTree(DM refTree, Mat *inj) 3047 { 3048 PetscDS ds; 3049 PetscSection section, cSection; 3050 DMLabel canonical, depth; 3051 Mat cMat, mat; 3052 PetscInt *nnz; 3053 PetscInt f, dim, numFields, numSecFields, p, pStart, pEnd, cStart, cEnd; 3054 PetscInt m, n; 3055 PetscScalar *pointScalar; 3056 PetscReal *v0, *v0parent, *vtmp, *J, *Jparent, *invJ, *pointRef, detJ, detJparent; 3057 PetscErrorCode ierr; 3058 3059 PetscFunctionBegin; 3060 ierr = DMGetDefaultSection(refTree,§ion);CHKERRQ(ierr); 3061 ierr = DMGetDimension(refTree, &dim);CHKERRQ(ierr); 3062 ierr = PetscMalloc6(dim,&v0,dim,&v0parent,dim,&vtmp,dim*dim,&J,dim*dim,&Jparent,dim*dim,&invJ);CHKERRQ(ierr); 3063 ierr = PetscMalloc2(dim,&pointScalar,dim,&pointRef);CHKERRQ(ierr); 3064 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3065 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3066 ierr = PetscSectionGetNumFields(section,&numSecFields);CHKERRQ(ierr); 3067 ierr = DMGetLabel(refTree,"canonical",&canonical);CHKERRQ(ierr); 3068 ierr = DMGetLabel(refTree,"depth",&depth);CHKERRQ(ierr); 3069 ierr = DMGetDefaultConstraints(refTree,&cSection,&cMat);CHKERRQ(ierr); 3070 ierr = DMPlexGetChart(refTree, &pStart, &pEnd);CHKERRQ(ierr); 3071 ierr = DMPlexGetHeightStratum(refTree, 0, &cStart, &cEnd);CHKERRQ(ierr); 3072 ierr = MatGetSize(cMat,&n,&m);CHKERRQ(ierr); /* the injector has transpose sizes from the constraint matrix */ 3073 /* Step 1: compute non-zero pattern. A proper subset of constraint matrix non-zero */ 3074 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 3075 for (p = pStart; p < pEnd; p++) { /* a point will have non-zeros if it is canonical, it has dofs, and its children have dofs */ 3076 const PetscInt *children; 3077 PetscInt numChildren; 3078 PetscInt i, numChildDof, numSelfDof; 3079 3080 if (canonical) { 3081 PetscInt pCanonical; 3082 ierr = DMLabelGetValue(canonical,p,&pCanonical);CHKERRQ(ierr); 3083 if (p != pCanonical) continue; 3084 } 3085 ierr = DMPlexGetTreeChildren(refTree,p,&numChildren,&children);CHKERRQ(ierr); 3086 if (!numChildren) continue; 3087 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3088 PetscInt child = children[i]; 3089 PetscInt dof; 3090 3091 ierr = PetscSectionGetDof(section,child,&dof);CHKERRQ(ierr); 3092 numChildDof += dof; 3093 } 3094 ierr = PetscSectionGetDof(section,p,&numSelfDof);CHKERRQ(ierr); 3095 if (!numChildDof || !numSelfDof) continue; 3096 for (f = 0; f < numFields; f++) { 3097 PetscInt selfOff; 3098 3099 if (numSecFields) { /* count the dofs for just this field */ 3100 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3101 PetscInt child = children[i]; 3102 PetscInt dof; 3103 3104 ierr = PetscSectionGetFieldDof(section,child,f,&dof);CHKERRQ(ierr); 3105 numChildDof += dof; 3106 } 3107 ierr = PetscSectionGetFieldDof(section,p,f,&numSelfDof);CHKERRQ(ierr); 3108 ierr = PetscSectionGetFieldOffset(section,p,f,&selfOff);CHKERRQ(ierr); 3109 } 3110 else { 3111 ierr = PetscSectionGetOffset(section,p,&selfOff);CHKERRQ(ierr); 3112 } 3113 for (i = 0; i < numSelfDof; i++) { 3114 nnz[selfOff + i] = numChildDof; 3115 } 3116 } 3117 } 3118 ierr = MatCreateAIJ(PETSC_COMM_SELF,m,n,m,n,-1,nnz,-1,NULL,&mat);CHKERRQ(ierr); 3119 ierr = PetscFree(nnz);CHKERRQ(ierr); 3120 /* Setp 2: compute entries */ 3121 for (p = pStart; p < pEnd; p++) { 3122 const PetscInt *children; 3123 PetscInt numChildren; 3124 PetscInt i, numChildDof, numSelfDof; 3125 3126 /* same conditions about when entries occur */ 3127 if (canonical) { 3128 PetscInt pCanonical; 3129 ierr = DMLabelGetValue(canonical,p,&pCanonical);CHKERRQ(ierr); 3130 if (p != pCanonical) continue; 3131 } 3132 ierr = DMPlexGetTreeChildren(refTree,p,&numChildren,&children);CHKERRQ(ierr); 3133 if (!numChildren) continue; 3134 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3135 PetscInt child = children[i]; 3136 PetscInt dof; 3137 3138 ierr = PetscSectionGetDof(section,child,&dof);CHKERRQ(ierr); 3139 numChildDof += dof; 3140 } 3141 ierr = PetscSectionGetDof(section,p,&numSelfDof);CHKERRQ(ierr); 3142 if (!numChildDof || !numSelfDof) continue; 3143 3144 for (f = 0; f < numFields; f++) { 3145 PetscInt selfOff, fComp, numSelfShapes, numChildShapes, parentCell; 3146 PetscInt cellShapeOff; 3147 PetscObject disc; 3148 PetscDualSpace dsp; 3149 PetscClassId classId; 3150 PetscScalar *pointMat; 3151 PetscInt *matRows, *matCols; 3152 PetscInt pO = PETSC_MIN_INT; 3153 const PetscInt *depthNumDof; 3154 3155 if (numSecFields) { 3156 for (i = 0, numChildDof = 0; i < numChildren; i++) { 3157 PetscInt child = children[i]; 3158 PetscInt dof; 3159 3160 ierr = PetscSectionGetFieldDof(section,child,f,&dof);CHKERRQ(ierr); 3161 numChildDof += dof; 3162 } 3163 ierr = PetscSectionGetFieldDof(section,p,f,&numSelfDof);CHKERRQ(ierr); 3164 ierr = PetscSectionGetFieldOffset(section,p,f,&selfOff);CHKERRQ(ierr); 3165 ierr = PetscSectionGetFieldComponents(section,f,&fComp);CHKERRQ(ierr); 3166 } 3167 else { 3168 ierr = PetscSectionGetOffset(section,p,&selfOff);CHKERRQ(ierr); 3169 fComp = 1; 3170 } 3171 numSelfShapes = numSelfDof / fComp; 3172 numChildShapes = numChildDof / fComp; 3173 3174 /* find a cell whose closure contains p */ 3175 if (p >= cStart && p < cEnd) { 3176 parentCell = p; 3177 } 3178 else { 3179 PetscInt *star = NULL; 3180 PetscInt numStar; 3181 3182 parentCell = -1; 3183 ierr = DMPlexGetTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3184 for (i = numStar - 1; i >= 0; i--) { 3185 PetscInt c = star[2 * i]; 3186 3187 if (c >= cStart && c < cEnd) { 3188 parentCell = c; 3189 break; 3190 } 3191 } 3192 ierr = DMPlexRestoreTransitiveClosure(refTree,p,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3193 } 3194 /* determine the offset of p's shape functions withing parentCell's shape functions */ 3195 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 3196 ierr = PetscObjectGetClassId(disc,&classId);CHKERRQ(ierr); 3197 if (classId == PETSCFE_CLASSID) { 3198 ierr = PetscFEGetDualSpace((PetscFE)disc,&dsp);CHKERRQ(ierr); 3199 } 3200 else if (classId == PETSCFV_CLASSID) { 3201 ierr = PetscFVGetDualSpace((PetscFV)disc,&dsp);CHKERRQ(ierr); 3202 } 3203 else { 3204 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported discretization object");CHKERRQ(ierr); 3205 } 3206 ierr = PetscDualSpaceGetNumDof(dsp,&depthNumDof);CHKERRQ(ierr); 3207 { 3208 PetscInt *closure = NULL; 3209 PetscInt numClosure; 3210 3211 ierr = DMPlexGetTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3212 for (i = 0, cellShapeOff = 0; i < numClosure; i++) { 3213 PetscInt point = closure[2 * i], pointDepth; 3214 3215 pO = closure[2 * i + 1]; 3216 if (point == p) break; 3217 ierr = DMLabelGetValue(depth,point,&pointDepth);CHKERRQ(ierr); 3218 cellShapeOff += depthNumDof[pointDepth]; 3219 } 3220 ierr = DMPlexRestoreTransitiveClosure(refTree,parentCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3221 } 3222 3223 ierr = DMGetWorkArray(refTree, numSelfShapes * numChildShapes, PETSC_SCALAR,&pointMat);CHKERRQ(ierr); 3224 ierr = DMGetWorkArray(refTree, numSelfShapes + numChildShapes, PETSC_INT,&matRows);CHKERRQ(ierr); 3225 matCols = matRows + numSelfShapes; 3226 for (i = 0; i < numSelfShapes; i++) { 3227 matRows[i] = selfOff + i * fComp; 3228 } 3229 { 3230 PetscInt colOff = 0; 3231 3232 for (i = 0; i < numChildren; i++) { 3233 PetscInt child = children[i]; 3234 PetscInt dof, off, j; 3235 3236 if (numSecFields) { 3237 ierr = PetscSectionGetFieldDof(cSection,child,f,&dof);CHKERRQ(ierr); 3238 ierr = PetscSectionGetFieldOffset(cSection,child,f,&off);CHKERRQ(ierr); 3239 } 3240 else { 3241 ierr = PetscSectionGetDof(cSection,child,&dof);CHKERRQ(ierr); 3242 ierr = PetscSectionGetOffset(cSection,child,&off);CHKERRQ(ierr); 3243 } 3244 3245 for (j = 0; j < dof / fComp; j++) { 3246 matCols[colOff++] = off + j * fComp; 3247 } 3248 } 3249 } 3250 if (classId == PETSCFE_CLASSID) { 3251 PetscFE fe = (PetscFE) disc; 3252 PetscInt fSize; 3253 3254 ierr = PetscFEGetDualSpace(fe,&dsp);CHKERRQ(ierr); 3255 ierr = PetscDualSpaceGetDimension(dsp,&fSize);CHKERRQ(ierr); 3256 for (i = 0; i < numSelfShapes; i++) { /* for every shape function */ 3257 PetscQuadrature q; 3258 PetscInt dim, numPoints, j, k; 3259 const PetscReal *points; 3260 const PetscReal *weights; 3261 PetscInt *closure = NULL; 3262 PetscInt numClosure; 3263 PetscInt parentCellShapeDof = cellShapeOff + (pO < 0 ? (numSelfShapes - 1 - i) : i); 3264 PetscReal *Bparent; 3265 3266 ierr = PetscDualSpaceGetFunctional(dsp,parentCellShapeDof,&q);CHKERRQ(ierr); 3267 ierr = PetscQuadratureGetData(q,&dim,&numPoints,&points,&weights);CHKERRQ(ierr); 3268 ierr = PetscFEGetTabulation(fe,numPoints,points,&Bparent,NULL,NULL);CHKERRQ(ierr); /* I'm expecting a nodal basis: weights[:]' * Bparent[:,cellShapeDof] = 1. */ 3269 for (k = 0; k < numChildShapes; k++) { 3270 pointMat[numChildShapes * i + k] = 0.; 3271 } 3272 for (j = 0; j < numPoints; j++) { 3273 PetscInt childCell = -1; 3274 PetscReal parentValAtPoint; 3275 const PetscReal *pointReal = &points[dim * j]; 3276 const PetscScalar *point; 3277 PetscReal *Bchild; 3278 PetscInt childCellShapeOff, pointMatOff; 3279 #if defined(PETSC_USE_COMPLEX) 3280 PetscInt d; 3281 3282 for (d = 0; d < dim; d++) { 3283 pointScalar[d] = points[dim * j + d]; 3284 } 3285 point = pointScalar; 3286 #else 3287 point = pointReal; 3288 #endif 3289 3290 parentValAtPoint = Bparent[(fSize * j + parentCellShapeDof) * fComp]; 3291 3292 for (k = 0; k < numChildren; k++) { /* locate the point in a child's star cell*/ 3293 PetscInt child = children[k]; 3294 PetscInt *star = NULL; 3295 PetscInt numStar, s; 3296 3297 ierr = DMPlexGetTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3298 for (s = numStar - 1; s >= 0; s--) { 3299 PetscInt c = star[2 * s]; 3300 3301 if (c < cStart || c >= cEnd) continue; 3302 ierr = DMPlexLocatePoint_Internal(refTree,dim,point,c,&childCell);CHKERRQ(ierr); 3303 if (childCell >= 0) break; 3304 } 3305 ierr = DMPlexRestoreTransitiveClosure(refTree,child,PETSC_FALSE,&numStar,&star);CHKERRQ(ierr); 3306 if (childCell >= 0) break; 3307 } 3308 if (childCell < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Could not locate quadrature point");CHKERRQ(ierr); 3309 ierr = DMPlexComputeCellGeometryFEM(refTree, childCell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 3310 ierr = DMPlexComputeCellGeometryFEM(refTree, parentCell, NULL, v0parent, Jparent, NULL, &detJparent);CHKERRQ(ierr); 3311 CoordinatesRefToReal(dim, dim, v0parent, Jparent, pointReal, vtmp); 3312 CoordinatesRealToRef(dim, dim, v0, invJ, vtmp, pointRef); 3313 3314 ierr = PetscFEGetTabulation(fe,1,pointRef,&Bchild,NULL,NULL);CHKERRQ(ierr); 3315 ierr = DMPlexGetTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3316 for (k = 0, pointMatOff = 0; k < numChildren; k++) { /* point is located in cell => child dofs support at point are in closure of cell */ 3317 PetscInt child = children[k], childDepth, childDof, childO = PETSC_MIN_INT; 3318 PetscInt l; 3319 3320 ierr = DMLabelGetValue(depth,child,&childDepth);CHKERRQ(ierr); 3321 childDof = depthNumDof[childDepth]; 3322 for (l = 0, childCellShapeOff = 0; l < numClosure; l++) { 3323 PetscInt point = closure[2 * l]; 3324 PetscInt pointDepth; 3325 3326 childO = closure[2 * l + 1]; 3327 if (point == child) break; 3328 ierr = DMLabelGetValue(depth,point,&pointDepth);CHKERRQ(ierr); 3329 childCellShapeOff += depthNumDof[pointDepth]; 3330 } 3331 if (l == numClosure) { 3332 pointMatOff += childDof; 3333 continue; /* child is not in the closure of the cell: has nothing to contribute to this point */ 3334 } 3335 for (l = 0; l < childDof; l++) { 3336 PetscInt childCellDof = childCellShapeOff + (childO ? (childDof - 1 - l) : l); 3337 PetscReal childValAtPoint = Bchild[childCellDof * fComp]; 3338 3339 pointMat[i * numChildShapes + pointMatOff + l] += weights[j] * parentValAtPoint * childValAtPoint; 3340 } 3341 pointMatOff += childDof; 3342 } 3343 ierr = DMPlexRestoreTransitiveClosure(refTree,childCell,PETSC_TRUE,&numClosure,&closure);CHKERRQ(ierr); 3344 ierr = PetscFERestoreTabulation(fe,1,pointRef,&Bchild,NULL,NULL);CHKERRQ(ierr); 3345 } 3346 ierr = PetscFERestoreTabulation(fe,numPoints,points,&Bparent,NULL,NULL);CHKERRQ(ierr); 3347 } 3348 } 3349 else { /* just the volume-weighted averages of the children */ 3350 PetscInt childShapeOff; 3351 PetscReal parentVol; 3352 3353 ierr = DMPlexComputeCellGeometryFVM(refTree, p, &parentVol, NULL, NULL);CHKERRQ(ierr); 3354 for (i = 0, childShapeOff = 0; i < numChildren; i++) { 3355 PetscInt child = children[i]; 3356 PetscReal childVol; 3357 3358 if (child < cStart || child >= cEnd) continue; 3359 ierr = DMPlexComputeCellGeometryFVM(refTree, child, &childVol, NULL, NULL);CHKERRQ(ierr); 3360 pointMat[childShapeOff] = childVol / parentVol; 3361 childShapeOff++; 3362 } 3363 } 3364 /* Insert pointMat into mat */ 3365 for (i = 0; i < fComp; i++) { 3366 PetscInt j; 3367 ierr = MatSetValues(mat,numSelfShapes,matRows,numChildShapes,matCols,pointMat,INSERT_VALUES);CHKERRQ(ierr); 3368 3369 for (j = 0; j < numSelfShapes; j++) { 3370 matRows[j]++; 3371 } 3372 for (j = 0; j < numChildShapes; j++) { 3373 matCols[j]++; 3374 } 3375 } 3376 ierr = DMRestoreWorkArray(refTree, numSelfShapes + numChildShapes, PETSC_INT,&matRows);CHKERRQ(ierr); 3377 ierr = DMRestoreWorkArray(refTree, numSelfShapes * numChildShapes, PETSC_SCALAR,&pointMat);CHKERRQ(ierr); 3378 } 3379 } 3380 ierr = PetscFree6(v0,v0parent,vtmp,J,Jparent,invJ);CHKERRQ(ierr); 3381 ierr = PetscFree2(pointScalar,pointRef);CHKERRQ(ierr); 3382 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3383 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3384 *inj = mat; 3385 PetscFunctionReturn(0); 3386 } 3387 3388 #undef __FUNCT__ 3389 #define __FUNCT__ "DMPlexReferenceTreeGetChildrenMatrices_Injection" 3390 static PetscErrorCode DMPlexReferenceTreeGetChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3391 { 3392 PetscDS ds; 3393 PetscInt numFields, f, pRefStart, pRefEnd, p, *rows, *cols, maxDof; 3394 PetscScalar ***refPointFieldMats; 3395 PetscSection refConSec, refSection; 3396 PetscErrorCode ierr; 3397 3398 PetscFunctionBegin; 3399 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3400 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3401 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 3402 ierr = DMGetDefaultSection(refTree,&refSection);CHKERRQ(ierr); 3403 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3404 ierr = PetscMalloc1(pRefEnd-pRefStart,&refPointFieldMats);CHKERRQ(ierr); 3405 ierr = PetscSectionGetMaxDof(refConSec,&maxDof);CHKERRQ(ierr); 3406 ierr = PetscMalloc1(maxDof,&rows);CHKERRQ(ierr); 3407 ierr = PetscMalloc1(maxDof*maxDof,&cols);CHKERRQ(ierr); 3408 for (p = pRefStart; p < pRefEnd; p++) { 3409 PetscInt parent, pDof, parentDof; 3410 3411 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 3412 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 3413 ierr = PetscSectionGetDof(refSection,parent,&parentDof);CHKERRQ(ierr); 3414 if (!pDof || !parentDof || parent == p) continue; 3415 3416 ierr = PetscMalloc1(numFields,&refPointFieldMats[p-pRefStart]);CHKERRQ(ierr); 3417 for (f = 0; f < numFields; f++) { 3418 PetscInt cDof, cOff, numCols, r, fComp; 3419 PetscObject disc; 3420 PetscClassId id; 3421 PetscFE fe = NULL; 3422 PetscFV fv = NULL; 3423 3424 ierr = PetscDSGetDiscretization(ds,f,&disc);CHKERRQ(ierr); 3425 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 3426 if (id == PETSCFE_CLASSID) { 3427 fe = (PetscFE) disc; 3428 ierr = PetscFEGetNumComponents(fe,&fComp);CHKERRQ(ierr); 3429 } 3430 else if (id == PETSCFV_CLASSID) { 3431 fv = (PetscFV) disc; 3432 ierr = PetscFVGetNumComponents(fv,&fComp);CHKERRQ(ierr); 3433 } 3434 else SETERRQ1(PetscObjectComm(disc),PETSC_ERR_ARG_UNKNOWN_TYPE, "PetscDS discretization id %d not recognized.", id); 3435 3436 if (numFields > 1) { 3437 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 3438 ierr = PetscSectionGetFieldOffset(refConSec,p,f,&cOff);CHKERRQ(ierr); 3439 } 3440 else { 3441 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 3442 ierr = PetscSectionGetOffset(refConSec,p,&cOff);CHKERRQ(ierr); 3443 } 3444 3445 for (r = 0; r < cDof; r++) { 3446 rows[r] = cOff + r; 3447 } 3448 numCols = 0; 3449 { 3450 PetscInt aDof, aOff, j; 3451 3452 if (numFields > 1) { 3453 ierr = PetscSectionGetFieldDof(refSection,parent,f,&aDof);CHKERRQ(ierr); 3454 ierr = PetscSectionGetFieldOffset(refSection,parent,f,&aOff);CHKERRQ(ierr); 3455 } 3456 else { 3457 ierr = PetscSectionGetDof(refSection,parent,&aDof);CHKERRQ(ierr); 3458 ierr = PetscSectionGetOffset(refSection,parent,&aOff);CHKERRQ(ierr); 3459 } 3460 3461 for (j = 0; j < aDof; j++) { 3462 cols[numCols++] = aOff + j; 3463 } 3464 } 3465 ierr = PetscMalloc1(cDof*numCols,&refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 3466 /* transpose of constraint matrix */ 3467 ierr = MatGetValues(inj,numCols,cols,cDof,rows,refPointFieldMats[p-pRefStart][f]);CHKERRQ(ierr); 3468 } 3469 } 3470 *childrenMats = refPointFieldMats; 3471 ierr = PetscFree(rows);CHKERRQ(ierr); 3472 ierr = PetscFree(cols);CHKERRQ(ierr); 3473 PetscFunctionReturn(0); 3474 } 3475 3476 #undef __FUNCT__ 3477 #define __FUNCT__ "DMPlexReferenceTreeRestoreChildrenMatrices_Injection" 3478 static PetscErrorCode DMPlexReferenceTreeRestoreChildrenMatrices_Injection(DM refTree, Mat inj, PetscScalar ****childrenMats) 3479 { 3480 PetscDS ds; 3481 PetscScalar ***refPointFieldMats; 3482 PetscInt numFields, pRefStart, pRefEnd, p, f; 3483 PetscSection refConSec, refSection; 3484 PetscErrorCode ierr; 3485 3486 PetscFunctionBegin; 3487 refPointFieldMats = *childrenMats; 3488 *childrenMats = NULL; 3489 ierr = DMGetDS(refTree,&ds);CHKERRQ(ierr); 3490 ierr = DMGetDefaultSection(refTree,&refSection);CHKERRQ(ierr); 3491 ierr = PetscDSGetNumFields(ds,&numFields);CHKERRQ(ierr); 3492 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 3493 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3494 for (p = pRefStart; p < pRefEnd; p++) { 3495 PetscInt parent, pDof, parentDof; 3496 3497 ierr = DMPlexGetTreeParent(refTree,p,&parent,NULL);CHKERRQ(ierr); 3498 ierr = PetscSectionGetDof(refConSec,p,&pDof);CHKERRQ(ierr); 3499 ierr = PetscSectionGetDof(refSection,parent,&parentDof);CHKERRQ(ierr); 3500 if (!pDof || !parentDof || parent == p) continue; 3501 3502 for (f = 0; f < numFields; f++) { 3503 PetscInt cDof; 3504 3505 if (numFields > 1) { 3506 ierr = PetscSectionGetFieldDof(refConSec,p,f,&cDof);CHKERRQ(ierr); 3507 } 3508 else { 3509 ierr = PetscSectionGetDof(refConSec,p,&cDof);CHKERRQ(ierr); 3510 } 3511 3512 ierr = PetscFree(refPointFieldMats[p - pRefStart][f]);CHKERRQ(ierr); 3513 } 3514 ierr = PetscFree(refPointFieldMats[p - pRefStart]);CHKERRQ(ierr); 3515 } 3516 ierr = PetscFree(refPointFieldMats);CHKERRQ(ierr); 3517 PetscFunctionReturn(0); 3518 } 3519 3520 #undef __FUNCT__ 3521 #define __FUNCT__ "DMPlexReferenceTreeGetInjector" 3522 static PetscErrorCode DMPlexReferenceTreeGetInjector(DM refTree,Mat *injRef) 3523 { 3524 Mat cMatRef; 3525 PetscObject injRefObj; 3526 PetscErrorCode ierr; 3527 3528 PetscFunctionBegin; 3529 ierr = DMGetDefaultConstraints(refTree,NULL,&cMatRef);CHKERRQ(ierr); 3530 ierr = PetscObjectQuery((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",&injRefObj);CHKERRQ(ierr); 3531 *injRef = (Mat) injRefObj; 3532 if (!*injRef) { 3533 ierr = DMPlexComputeInjectorReferenceTree(refTree,injRef);CHKERRQ(ierr); 3534 ierr = PetscObjectCompose((PetscObject)cMatRef,"DMPlexComputeInjectorTree_refTree",(PetscObject)*injRef);CHKERRQ(ierr); 3535 /* there is now a reference in cMatRef, which should be the only one for symmetry with the above case */ 3536 ierr = PetscObjectDereference((PetscObject)*injRef);CHKERRQ(ierr); 3537 } 3538 PetscFunctionReturn(0); 3539 } 3540 3541 #undef __FUNCT__ 3542 #define __FUNCT__ "DMPlexTransferInjectorTree" 3543 static PetscErrorCode DMPlexTransferInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, const PetscInt *childIds, Vec fineVec, PetscInt numFields, PetscInt *offsets, PetscSection *rootMultiSec, PetscSection *multiLeafSec, PetscInt **gatheredIndices, PetscScalar **gatheredValues) 3544 { 3545 PetscInt pStartF, pEndF, pStartC, pEndC, p, maxDof, numMulti; 3546 PetscSection globalCoarse, globalFine; 3547 PetscSection localCoarse, localFine, leafIndicesSec; 3548 PetscSection multiRootSec, rootIndicesSec; 3549 PetscInt *leafInds, *rootInds = NULL; 3550 const PetscInt *rootDegrees; 3551 PetscScalar *leafVals = NULL, *rootVals = NULL; 3552 PetscSF coarseToFineEmbedded; 3553 PetscErrorCode ierr; 3554 3555 PetscFunctionBegin; 3556 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3557 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 3558 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 3559 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 3560 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafIndicesSec);CHKERRQ(ierr); 3561 ierr = PetscSectionSetChart(leafIndicesSec,pStartF, pEndF);CHKERRQ(ierr); 3562 ierr = PetscSectionGetMaxDof(localFine,&maxDof);CHKERRQ(ierr); 3563 { /* winnow fine points that don't have global dofs out of the sf */ 3564 PetscInt l, nleaves, dof, cdof, numPointsWithDofs, offset, *pointsWithDofs, numIndices; 3565 const PetscInt *leaves; 3566 3567 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 3568 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 3569 p = leaves ? leaves[l] : l; 3570 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 3571 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 3572 if ((dof - cdof) > 0) { 3573 numPointsWithDofs++; 3574 3575 ierr = PetscSectionGetDof(localFine,p,&dof);CHKERRQ(ierr); 3576 ierr = PetscSectionSetDof(leafIndicesSec,p,dof + 1);CHKERRQ(ierr); 3577 } 3578 } 3579 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 3580 ierr = PetscSectionSetUp(leafIndicesSec);CHKERRQ(ierr); 3581 ierr = PetscSectionGetStorageSize(leafIndicesSec,&numIndices);CHKERRQ(ierr); 3582 ierr = PetscMalloc1(gatheredIndices ? numIndices : (maxDof + 1),&leafInds);CHKERRQ(ierr); 3583 if (gatheredValues) {ierr = PetscMalloc1(numIndices,&leafVals);CHKERRQ(ierr);} 3584 for (l = 0, offset = 0; l < nleaves; l++) { 3585 p = leaves ? leaves[l] : l; 3586 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 3587 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 3588 if ((dof - cdof) > 0) { 3589 PetscInt off, gOff; 3590 PetscInt *pInd; 3591 PetscScalar *pVal = NULL; 3592 3593 pointsWithDofs[offset++] = l; 3594 3595 ierr = PetscSectionGetOffset(leafIndicesSec,p,&off);CHKERRQ(ierr); 3596 3597 pInd = gatheredIndices ? (&leafInds[off + 1]) : leafInds; 3598 if (gatheredValues) { 3599 PetscInt i; 3600 3601 pVal = &leafVals[off + 1]; 3602 for (i = 0; i < dof; i++) pVal[i] = 0.; 3603 } 3604 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 3605 3606 offsets[0] = 0; 3607 if (numFields) { 3608 PetscInt f; 3609 3610 for (f = 0; f < numFields; f++) { 3611 PetscInt fDof; 3612 ierr = PetscSectionGetFieldDof(localFine,p,f,&fDof);CHKERRQ(ierr); 3613 offsets[f + 1] = fDof + offsets[f]; 3614 } 3615 DMPlexGetIndicesPointFields_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,-1,pInd);CHKERRQ(ierr); 3616 } 3617 else { 3618 DMPlexGetIndicesPoint_Internal(localFine,p,gOff < 0 ? -(gOff + 1) : gOff,offsets,PETSC_FALSE,NULL,pInd);CHKERRQ(ierr); 3619 } 3620 if (gatheredValues) {ierr = VecGetValues(fineVec,dof,pInd,pVal);CHKERRQ(ierr);} 3621 } 3622 } 3623 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 3624 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 3625 } 3626 3627 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3628 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 3629 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 3630 3631 { /* there may be the case where an sf root has a parent: broadcast parents back to children */ 3632 MPI_Datatype threeInt; 3633 PetscMPIInt rank; 3634 PetscInt (*parentNodeAndIdCoarse)[3]; 3635 PetscInt (*parentNodeAndIdFine)[3]; 3636 PetscInt p, nleaves, nleavesToParents; 3637 PetscSF pointSF, sfToParents; 3638 const PetscInt *ilocal; 3639 const PetscSFNode *iremote; 3640 PetscSFNode *iremoteToParents; 3641 PetscInt *ilocalToParents; 3642 3643 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);CHKERRQ(ierr); 3644 ierr = MPI_Type_contiguous(3,MPIU_INT,&threeInt);CHKERRQ(ierr); 3645 ierr = MPI_Type_commit(&threeInt);CHKERRQ(ierr); 3646 ierr = PetscMalloc2(pEndC-pStartC,&parentNodeAndIdCoarse,pEndF-pStartF,&parentNodeAndIdFine);CHKERRQ(ierr); 3647 ierr = DMGetPointSF(coarse,&pointSF);CHKERRQ(ierr); 3648 ierr = PetscSFGetGraph(pointSF,NULL,&nleaves,&ilocal,&iremote);CHKERRQ(ierr); 3649 for (p = pStartC; p < pEndC; p++) { 3650 PetscInt parent, childId; 3651 ierr = DMPlexGetTreeParent(coarse,p,&parent,&childId);CHKERRQ(ierr); 3652 parentNodeAndIdCoarse[p - pStartC][0] = rank; 3653 parentNodeAndIdCoarse[p - pStartC][1] = parent - pStartC; 3654 parentNodeAndIdCoarse[p - pStartC][2] = (p == parent) ? -1 : childId; 3655 if (nleaves > 0) { 3656 PetscInt leaf = -1; 3657 3658 if (ilocal) { 3659 ierr = PetscFindInt(parent,nleaves,ilocal,&leaf);CHKERRQ(ierr); 3660 } 3661 else { 3662 leaf = p - pStartC; 3663 } 3664 if (leaf >= 0) { 3665 parentNodeAndIdCoarse[p - pStartC][0] = iremote[leaf].rank; 3666 parentNodeAndIdCoarse[p - pStartC][1] = iremote[leaf].index; 3667 } 3668 } 3669 } 3670 for (p = pStartF; p < pEndF; p++) { 3671 parentNodeAndIdFine[p - pStartF][0] = -1; 3672 parentNodeAndIdFine[p - pStartF][1] = -1; 3673 parentNodeAndIdFine[p - pStartF][2] = -1; 3674 } 3675 ierr = PetscSFBcastBegin(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3676 ierr = PetscSFBcastEnd(coarseToFineEmbedded,threeInt,parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3677 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3678 PetscInt dof; 3679 3680 ierr = PetscSectionGetDof(leafIndicesSec,p,&dof);CHKERRQ(ierr); 3681 if (dof) { 3682 PetscInt off; 3683 3684 ierr = PetscSectionGetOffset(leafIndicesSec,p,&off);CHKERRQ(ierr); 3685 if (gatheredIndices) { 3686 leafInds[off] = PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3687 } else if (gatheredValues) { 3688 leafVals[off] = (PetscScalar) PetscMax(childIds[p-pStartF],parentNodeAndIdFine[p-pStartF][2]); 3689 } 3690 } 3691 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3692 nleavesToParents++; 3693 } 3694 } 3695 ierr = PetscMalloc1(nleavesToParents,&ilocalToParents);CHKERRQ(ierr); 3696 ierr = PetscMalloc1(nleavesToParents,&iremoteToParents);CHKERRQ(ierr); 3697 for (p = pStartF, nleavesToParents = 0; p < pEndF; p++) { 3698 if (parentNodeAndIdFine[p-pStartF][0] >= 0) { 3699 ilocalToParents[nleavesToParents] = p - pStartF; 3700 iremoteToParents[nleavesToParents].rank = parentNodeAndIdFine[p-pStartF][0]; 3701 iremoteToParents[nleavesToParents].index = parentNodeAndIdFine[p-pStartF][1]; 3702 nleavesToParents++; 3703 } 3704 } 3705 ierr = PetscSFCreate(PetscObjectComm((PetscObject)coarse),&sfToParents);CHKERRQ(ierr); 3706 ierr = PetscSFSetGraph(sfToParents,pEndC-pStartC,nleavesToParents,ilocalToParents,PETSC_OWN_POINTER,iremoteToParents,PETSC_OWN_POINTER);CHKERRQ(ierr); 3707 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3708 3709 coarseToFineEmbedded = sfToParents; 3710 3711 ierr = PetscFree2(parentNodeAndIdCoarse,parentNodeAndIdFine);CHKERRQ(ierr); 3712 ierr = MPI_Type_free(&threeInt);CHKERRQ(ierr); 3713 } 3714 3715 { /* winnow out coarse points that don't have dofs */ 3716 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 3717 PetscSF sfDofsOnly; 3718 3719 for (p = pStartC, numPointsWithDofs = 0; p < pEndC; p++) { 3720 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3721 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3722 if ((dof - cdof) > 0) { 3723 numPointsWithDofs++; 3724 } 3725 } 3726 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 3727 for (p = pStartC, offset = 0; p < pEndC; p++) { 3728 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3729 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3730 if ((dof - cdof) > 0) { 3731 pointsWithDofs[offset++] = p - pStartC; 3732 } 3733 } 3734 ierr = PetscSFCreateEmbeddedSF(coarseToFineEmbedded, numPointsWithDofs, pointsWithDofs, &sfDofsOnly);CHKERRQ(ierr); 3735 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3736 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 3737 coarseToFineEmbedded = sfDofsOnly; 3738 } 3739 3740 /* communicate back to the coarse mesh which coarse points have children (that may require injection) */ 3741 ierr = PetscSFComputeDegreeBegin(coarseToFineEmbedded,&rootDegrees);CHKERRQ(ierr); 3742 ierr = PetscSFComputeDegreeEnd(coarseToFineEmbedded,&rootDegrees);CHKERRQ(ierr); 3743 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&multiRootSec);CHKERRQ(ierr); 3744 ierr = PetscSectionSetChart(multiRootSec,pStartC,pEndC);CHKERRQ(ierr); 3745 for (p = pStartC; p < pEndC; p++) { 3746 ierr = PetscSectionSetDof(multiRootSec,p,rootDegrees[p-pStartC]);CHKERRQ(ierr); 3747 } 3748 ierr = PetscSectionSetUp(multiRootSec);CHKERRQ(ierr); 3749 ierr = PetscSectionGetStorageSize(multiRootSec,&numMulti);CHKERRQ(ierr); 3750 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootIndicesSec);CHKERRQ(ierr); 3751 { /* distribute the leaf section */ 3752 PetscSF multi, multiInv, indicesSF; 3753 PetscInt *remoteOffsets, numRootIndices; 3754 3755 ierr = PetscSFGetMultiSF(coarseToFineEmbedded,&multi);CHKERRQ(ierr); 3756 ierr = PetscSFCreateInverseSF(multi,&multiInv);CHKERRQ(ierr); 3757 ierr = PetscSFDistributeSection(multiInv,leafIndicesSec,&remoteOffsets,rootIndicesSec);CHKERRQ(ierr); 3758 ierr = PetscSFCreateSectionSF(multiInv,leafIndicesSec,remoteOffsets,rootIndicesSec,&indicesSF);CHKERRQ(ierr); 3759 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 3760 ierr = PetscSFDestroy(&multiInv);CHKERRQ(ierr); 3761 ierr = PetscSectionGetStorageSize(rootIndicesSec,&numRootIndices);CHKERRQ(ierr); 3762 if (gatheredIndices) { 3763 ierr = PetscMalloc1(numRootIndices,&rootInds);CHKERRQ(ierr); 3764 ierr = PetscSFBcastBegin(indicesSF,MPIU_INT,leafInds,rootInds);CHKERRQ(ierr); 3765 ierr = PetscSFBcastEnd(indicesSF,MPIU_INT,leafInds,rootInds);CHKERRQ(ierr); 3766 } 3767 if (gatheredValues) { 3768 ierr = PetscMalloc1(numRootIndices,&rootVals);CHKERRQ(ierr); 3769 ierr = PetscSFBcastBegin(indicesSF,MPIU_SCALAR,leafVals,rootVals);CHKERRQ(ierr); 3770 ierr = PetscSFBcastEnd(indicesSF,MPIU_SCALAR,leafVals,rootVals);CHKERRQ(ierr); 3771 } 3772 ierr = PetscSFDestroy(&indicesSF);CHKERRQ(ierr); 3773 } 3774 ierr = PetscSectionDestroy(&leafIndicesSec);CHKERRQ(ierr); 3775 ierr = PetscFree(leafInds);CHKERRQ(ierr); 3776 ierr = PetscFree(leafVals);CHKERRQ(ierr); 3777 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 3778 *rootMultiSec = multiRootSec; 3779 *multiLeafSec = rootIndicesSec; 3780 if (gatheredIndices) *gatheredIndices = rootInds; 3781 if (gatheredValues) *gatheredValues = rootVals; 3782 PetscFunctionReturn(0); 3783 } 3784 3785 #undef __FUNCT__ 3786 #define __FUNCT__ "DMPlexComputeInjectorTree" 3787 PetscErrorCode DMPlexComputeInjectorTree(DM coarse, DM fine, PetscSF coarseToFine, PetscInt *childIds, Mat mat) 3788 { 3789 DM refTree; 3790 PetscSection multiRootSec, rootIndicesSec; 3791 PetscSection globalCoarse, globalFine; 3792 PetscSection localCoarse, localFine; 3793 PetscSection cSecRef; 3794 PetscInt *rootIndices, *parentIndices, pRefStart, pRefEnd; 3795 Mat injRef; 3796 PetscInt numFields, maxDof; 3797 PetscInt pStartC, pEndC, pStartF, pEndF, p; 3798 PetscInt *offsets, *offsetsCopy, *rowOffsets; 3799 PetscLayout rowMap, colMap; 3800 PetscInt rowStart, rowEnd, colStart, colEnd, *nnzD, *nnzO; 3801 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 3802 PetscErrorCode ierr; 3803 3804 PetscFunctionBegin; 3805 3806 /* get the templates for the fine-to-coarse injection from the reference tree */ 3807 ierr = DMPlexGetReferenceTree(coarse,&refTree);CHKERRQ(ierr); 3808 ierr = DMGetDefaultConstraints(refTree,&cSecRef,NULL);CHKERRQ(ierr); 3809 ierr = PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);CHKERRQ(ierr); 3810 ierr = DMPlexReferenceTreeGetInjector(refTree,&injRef);CHKERRQ(ierr); 3811 3812 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 3813 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 3814 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 3815 ierr = PetscSectionGetNumFields(localFine,&numFields);CHKERRQ(ierr); 3816 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 3817 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 3818 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 3819 ierr = PetscSectionGetMaxDof(localCoarse,&maxDof);CHKERRQ(ierr); 3820 { 3821 PetscInt maxFields = PetscMax(1,numFields) + 1; 3822 ierr = PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);CHKERRQ(ierr); 3823 } 3824 3825 ierr = DMPlexTransferInjectorTree(coarse,fine,coarseToFine,childIds,NULL,numFields,offsets,&multiRootSec,&rootIndicesSec,&rootIndices,NULL);CHKERRQ(ierr); 3826 3827 ierr = PetscMalloc1(maxDof,&parentIndices);CHKERRQ(ierr); 3828 3829 /* count indices */ 3830 ierr = MatGetLayouts(mat,&rowMap,&colMap);CHKERRQ(ierr); 3831 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 3832 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 3833 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 3834 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 3835 ierr = PetscCalloc2(rowEnd-rowStart,&nnzD,rowEnd-rowStart,&nnzO);CHKERRQ(ierr); 3836 for (p = pStartC; p < pEndC; p++) { 3837 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3838 3839 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3840 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3841 if ((dof - cdof) <= 0) continue; 3842 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 3843 3844 rowOffsets[0] = 0; 3845 offsetsCopy[0] = 0; 3846 if (numFields) { 3847 PetscInt f; 3848 3849 for (f = 0; f < numFields; f++) { 3850 PetscInt fDof; 3851 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 3852 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3853 } 3854 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 3855 } 3856 else { 3857 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 3858 rowOffsets[1] = offsetsCopy[0]; 3859 } 3860 3861 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 3862 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 3863 leafEnd = leafStart + numLeaves; 3864 for (l = leafStart; l < leafEnd; l++) { 3865 PetscInt numIndices, childId, offset; 3866 const PetscInt *childIndices; 3867 3868 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 3869 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 3870 childId = rootIndices[offset++]; 3871 childIndices = &rootIndices[offset]; 3872 numIndices--; 3873 3874 if (childId == -1) { /* equivalent points: scatter */ 3875 PetscInt i; 3876 3877 for (i = 0; i < numIndices; i++) { 3878 PetscInt colIndex = childIndices[i]; 3879 PetscInt rowIndex = parentIndices[i]; 3880 if (rowIndex < 0) continue; 3881 if (colIndex < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unconstrained fine and constrained coarse"); 3882 if (colIndex >= colStart && colIndex < colEnd) { 3883 nnzD[rowIndex - rowStart] = 1; 3884 } 3885 else { 3886 nnzO[rowIndex - rowStart] = 1; 3887 } 3888 } 3889 } 3890 else { 3891 PetscInt parentId, f, lim; 3892 3893 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 3894 3895 lim = PetscMax(1,numFields); 3896 offsets[0] = 0; 3897 if (numFields) { 3898 PetscInt f; 3899 3900 for (f = 0; f < numFields; f++) { 3901 PetscInt fDof; 3902 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 3903 3904 offsets[f + 1] = fDof + offsets[f]; 3905 } 3906 } 3907 else { 3908 PetscInt cDof; 3909 3910 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 3911 offsets[1] = cDof; 3912 } 3913 for (f = 0; f < lim; f++) { 3914 PetscInt parentStart = rowOffsets[f], parentEnd = rowOffsets[f + 1]; 3915 PetscInt childStart = offsets[f], childEnd = offsets[f + 1]; 3916 PetscInt i, numD = 0, numO = 0; 3917 3918 for (i = childStart; i < childEnd; i++) { 3919 PetscInt colIndex = childIndices[i]; 3920 3921 if (colIndex < 0) continue; 3922 if (colIndex >= colStart && colIndex < colEnd) { 3923 numD++; 3924 } 3925 else { 3926 numO++; 3927 } 3928 } 3929 for (i = parentStart; i < parentEnd; i++) { 3930 PetscInt rowIndex = parentIndices[i]; 3931 3932 if (rowIndex < 0) continue; 3933 nnzD[rowIndex - rowStart] += numD; 3934 nnzO[rowIndex - rowStart] += numO; 3935 } 3936 } 3937 } 3938 } 3939 } 3940 /* preallocate */ 3941 ierr = MatXAIJSetPreallocation(mat,1,nnzD,nnzO,NULL,NULL);CHKERRQ(ierr); 3942 ierr = PetscFree2(nnzD,nnzO);CHKERRQ(ierr); 3943 /* insert values */ 3944 ierr = DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 3945 for (p = pStartC; p < pEndC; p++) { 3946 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 3947 3948 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 3949 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 3950 if ((dof - cdof) <= 0) continue; 3951 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 3952 3953 rowOffsets[0] = 0; 3954 offsetsCopy[0] = 0; 3955 if (numFields) { 3956 PetscInt f; 3957 3958 for (f = 0; f < numFields; f++) { 3959 PetscInt fDof; 3960 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 3961 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 3962 } 3963 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 3964 } 3965 else { 3966 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 3967 rowOffsets[1] = offsetsCopy[0]; 3968 } 3969 3970 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 3971 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 3972 leafEnd = leafStart + numLeaves; 3973 for (l = leafStart; l < leafEnd; l++) { 3974 PetscInt numIndices, childId, offset; 3975 const PetscInt *childIndices; 3976 3977 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 3978 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 3979 childId = rootIndices[offset++]; 3980 childIndices = &rootIndices[offset]; 3981 numIndices--; 3982 3983 if (childId == -1) { /* equivalent points: scatter */ 3984 PetscInt i; 3985 3986 for (i = 0; i < numIndices; i++) { 3987 ierr = MatSetValue(mat,parentIndices[i],childIndices[i],1.,INSERT_VALUES);CHKERRQ(ierr); 3988 } 3989 } 3990 else { 3991 PetscInt parentId, f, lim; 3992 3993 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 3994 3995 lim = PetscMax(1,numFields); 3996 offsets[0] = 0; 3997 if (numFields) { 3998 PetscInt f; 3999 4000 for (f = 0; f < numFields; f++) { 4001 PetscInt fDof; 4002 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 4003 4004 offsets[f + 1] = fDof + offsets[f]; 4005 } 4006 } 4007 else { 4008 PetscInt cDof; 4009 4010 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 4011 offsets[1] = cDof; 4012 } 4013 for (f = 0; f < lim; f++) { 4014 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 4015 PetscInt *rowIndices = &parentIndices[rowOffsets[f]]; 4016 const PetscInt *colIndices = &childIndices[offsets[f]]; 4017 4018 ierr = MatSetValues(mat,rowOffsets[f+1]-rowOffsets[f],rowIndices,offsets[f+1]-offsets[f],colIndices,childMat,INSERT_VALUES);CHKERRQ(ierr); 4019 } 4020 } 4021 } 4022 } 4023 ierr = PetscSectionDestroy(&multiRootSec);CHKERRQ(ierr); 4024 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 4025 ierr = PetscFree(parentIndices);CHKERRQ(ierr); 4026 ierr = DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 4027 ierr = PetscFree(rootIndices);CHKERRQ(ierr); 4028 ierr = PetscFree3(offsets,offsetsCopy,rowOffsets);CHKERRQ(ierr); 4029 4030 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4031 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4032 PetscFunctionReturn(0); 4033 } 4034 4035 #undef __FUNCT__ 4036 #define __FUNCT__ "DMPlexTransferVecTree_Interpolate" 4037 static PetscErrorCode DMPlexTransferVecTree_Interpolate(DM coarse, Vec vecCoarseLocal, DM fine, Vec vecFine, PetscSF coarseToFine, PetscInt *cids, Vec grad, Vec cellGeom) 4038 { 4039 PetscDS ds; 4040 PetscSF coarseToFineEmbedded; 4041 PetscSection globalCoarse, globalFine; 4042 PetscSection localCoarse, localFine; 4043 PetscSection aSec, cSec; 4044 PetscSection rootValuesSec; 4045 PetscSection leafValuesSec; 4046 PetscScalar *rootValues, *leafValues; 4047 IS aIS; 4048 const PetscInt *anchors; 4049 Mat cMat; 4050 PetscInt numFields; 4051 PetscInt pStartC, pEndC, pStartF, pEndF, p, cellStart, cellEnd, cellEndInterior; 4052 PetscInt aStart, aEnd, cStart, cEnd; 4053 PetscInt *maxChildIds; 4054 PetscInt *offsets, *newOffsets, *offsetsCopy, *newOffsetsCopy, *rowOffsets, *numD, *numO; 4055 PetscFV fv = NULL; 4056 PetscInt dim, numFVcomps = -1, fvField = -1; 4057 DM cellDM = NULL, gradDM = NULL; 4058 const PetscScalar *cellGeomArray = NULL; 4059 const PetscScalar *gradArray = NULL; 4060 PetscErrorCode ierr; 4061 4062 PetscFunctionBegin; 4063 ierr = VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4064 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 4065 ierr = DMPlexGetHeightStratum(coarse,0,&cellStart,&cellEnd);CHKERRQ(ierr); 4066 ierr = DMPlexGetHybridBounds(coarse,&cellEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 4067 cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior; 4068 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 4069 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 4070 ierr = DMGetCoordinateDim(coarse,&dim);CHKERRQ(ierr); 4071 { /* winnow fine points that don't have global dofs out of the sf */ 4072 PetscInt nleaves, l; 4073 const PetscInt *leaves; 4074 PetscInt dof, cdof, numPointsWithDofs, offset, *pointsWithDofs; 4075 4076 ierr = PetscSFGetGraph(coarseToFine,NULL,&nleaves,&leaves,NULL);CHKERRQ(ierr); 4077 4078 for (l = 0, numPointsWithDofs = 0; l < nleaves; l++) { 4079 PetscInt p = leaves ? leaves[l] : l; 4080 4081 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 4082 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 4083 if ((dof - cdof) > 0) { 4084 numPointsWithDofs++; 4085 } 4086 } 4087 ierr = PetscMalloc1(numPointsWithDofs,&pointsWithDofs);CHKERRQ(ierr); 4088 for (l = 0, offset = 0; l < nleaves; l++) { 4089 PetscInt p = leaves ? leaves[l] : l; 4090 4091 ierr = PetscSectionGetDof(globalFine,p,&dof);CHKERRQ(ierr); 4092 ierr = PetscSectionGetConstraintDof(globalFine,p,&cdof);CHKERRQ(ierr); 4093 if ((dof - cdof) > 0) { 4094 pointsWithDofs[offset++] = l; 4095 } 4096 } 4097 ierr = PetscSFCreateEmbeddedLeafSF(coarseToFine, numPointsWithDofs, pointsWithDofs, &coarseToFineEmbedded);CHKERRQ(ierr); 4098 ierr = PetscFree(pointsWithDofs);CHKERRQ(ierr); 4099 } 4100 /* communicate back to the coarse mesh which coarse points have children (that may require interpolation) */ 4101 ierr = PetscMalloc1(pEndC-pStartC,&maxChildIds);CHKERRQ(ierr); 4102 for (p = pStartC; p < pEndC; p++) { 4103 maxChildIds[p - pStartC] = -2; 4104 } 4105 ierr = PetscSFReduceBegin(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 4106 ierr = PetscSFReduceEnd(coarseToFineEmbedded,MPIU_INT,cids,maxChildIds,MPIU_MAX);CHKERRQ(ierr); 4107 4108 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 4109 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 4110 4111 ierr = DMPlexGetAnchors(coarse,&aSec,&aIS);CHKERRQ(ierr); 4112 ierr = ISGetIndices(aIS,&anchors);CHKERRQ(ierr); 4113 ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr); 4114 4115 ierr = DMGetDefaultConstraints(coarse,&cSec,&cMat);CHKERRQ(ierr); 4116 ierr = PetscSectionGetChart(cSec,&cStart,&cEnd);CHKERRQ(ierr); 4117 4118 /* create sections that will send to children the indices and matrices they will need to construct the interpolator */ 4119 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)coarse),&rootValuesSec);CHKERRQ(ierr); 4120 ierr = PetscSectionSetChart(rootValuesSec,pStartC,pEndC);CHKERRQ(ierr); 4121 ierr = PetscSectionGetNumFields(localCoarse,&numFields);CHKERRQ(ierr); 4122 { 4123 PetscInt maxFields = PetscMax(1,numFields) + 1; 4124 ierr = PetscMalloc7(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&newOffsets,maxFields,&newOffsetsCopy,maxFields,&rowOffsets,maxFields,&numD,maxFields,&numO);CHKERRQ(ierr); 4125 } 4126 if (grad) { 4127 PetscInt i; 4128 4129 ierr = VecGetDM(cellGeom,&cellDM);CHKERRQ(ierr); 4130 ierr = VecGetArrayRead(cellGeom,&cellGeomArray);CHKERRQ(ierr); 4131 ierr = VecGetDM(grad,&gradDM);CHKERRQ(ierr); 4132 ierr = VecGetArrayRead(grad,&gradArray);CHKERRQ(ierr); 4133 for (i = 0; i < PetscMax(1,numFields); i++) { 4134 PetscObject obj; 4135 PetscClassId id; 4136 4137 ierr = DMGetField(coarse, i, &obj);CHKERRQ(ierr); 4138 ierr = PetscObjectGetClassId(obj,&id);CHKERRQ(ierr); 4139 if (id == PETSCFV_CLASSID) { 4140 fv = (PetscFV) obj; 4141 ierr = PetscFVGetNumComponents(fv,&numFVcomps);CHKERRQ(ierr); 4142 fvField = i; 4143 break; 4144 } 4145 } 4146 } 4147 4148 for (p = pStartC; p < pEndC; p++) { /* count the sizes of the indices and matrices */ 4149 PetscInt dof; 4150 PetscInt maxChildId = maxChildIds[p - pStartC]; 4151 PetscInt numValues = 0; 4152 4153 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 4154 if (dof < 0) { 4155 dof = -(dof + 1); 4156 } 4157 offsets[0] = 0; 4158 newOffsets[0] = 0; 4159 if (maxChildId >= 0) { /* this point has children (with dofs) that will need to be interpolated from the closure of p */ 4160 PetscInt *closure = NULL, closureSize, cl; 4161 4162 ierr = DMPlexGetTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 4163 for (cl = 0; cl < closureSize; cl++) { /* get the closure */ 4164 PetscInt c = closure[2 * cl], clDof; 4165 4166 ierr = PetscSectionGetDof(localCoarse,c,&clDof);CHKERRQ(ierr); 4167 numValues += clDof; 4168 } 4169 ierr = DMPlexRestoreTransitiveClosure(coarse,p,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 4170 } 4171 else if (maxChildId == -1) { 4172 ierr = PetscSectionGetDof(localCoarse,p,&numValues);CHKERRQ(ierr); 4173 } 4174 /* we will pack the column indices with the field offsets */ 4175 if (maxChildId >= 0 && grad && p >= cellStart && p < cellEnd) { 4176 /* also send the centroid, and the gradient */ 4177 numValues += dim * (1 + numFVcomps); 4178 } 4179 ierr = PetscSectionSetDof(rootValuesSec,p,numValues);CHKERRQ(ierr); 4180 } 4181 ierr = PetscSectionSetUp(rootValuesSec);CHKERRQ(ierr); 4182 { 4183 PetscInt numRootValues; 4184 const PetscScalar *coarseArray; 4185 4186 ierr = PetscSectionGetStorageSize(rootValuesSec,&numRootValues);CHKERRQ(ierr); 4187 ierr = PetscMalloc1(numRootValues,&rootValues);CHKERRQ(ierr); 4188 ierr = VecGetArrayRead(vecCoarseLocal,&coarseArray);CHKERRQ(ierr); 4189 for (p = pStartC; p < pEndC; p++) { 4190 PetscInt numValues; 4191 PetscInt pValOff; 4192 PetscScalar *pVal; 4193 PetscInt maxChildId = maxChildIds[p - pStartC]; 4194 4195 ierr = PetscSectionGetDof(rootValuesSec,p,&numValues);CHKERRQ(ierr); 4196 if (!numValues) { 4197 continue; 4198 } 4199 ierr = PetscSectionGetOffset(rootValuesSec,p,&pValOff);CHKERRQ(ierr); 4200 pVal = &(rootValues[pValOff]); 4201 if (maxChildId >= 0) { /* build an identity matrix, apply matrix constraints on the right */ 4202 PetscInt closureSize = numValues; 4203 ierr = DMPlexVecGetClosure(coarse,NULL,vecCoarseLocal,p,&closureSize,&pVal);CHKERRQ(ierr); 4204 if (grad && p >= cellStart && p < cellEnd) { 4205 PetscFVCellGeom *cg; 4206 PetscScalar *gradVals = NULL; 4207 PetscInt i; 4208 4209 pVal += (numValues - dim * (1 + numFVcomps)); 4210 4211 ierr = DMPlexPointLocalRead(cellDM,p,cellGeomArray,(void *) &cg);CHKERRQ(ierr); 4212 for (i = 0; i < dim; i++) pVal[i] = cg->centroid[i]; 4213 pVal += dim; 4214 ierr = DMPlexPointGlobalRead(gradDM,p,gradArray,(void *) &gradVals);CHKERRQ(ierr); 4215 for (i = 0; i < dim * numFVcomps; i++) pVal[i] = gradVals[i]; 4216 } 4217 } 4218 else if (maxChildId == -1) { 4219 PetscInt lDof, lOff, i; 4220 4221 ierr = PetscSectionGetDof(localCoarse,p,&lDof);CHKERRQ(ierr); 4222 ierr = PetscSectionGetOffset(localCoarse,p,&lOff);CHKERRQ(ierr); 4223 for (i = 0; i < lDof; i++) pVal[i] = coarseArray[lOff + i]; 4224 } 4225 } 4226 ierr = VecRestoreArrayRead(vecCoarseLocal,&coarseArray);CHKERRQ(ierr); 4227 ierr = PetscFree(maxChildIds);CHKERRQ(ierr); 4228 } 4229 { 4230 PetscSF valuesSF; 4231 PetscInt *remoteOffsetsValues, numLeafValues; 4232 4233 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)fine),&leafValuesSec);CHKERRQ(ierr); 4234 ierr = PetscSFDistributeSection(coarseToFineEmbedded,rootValuesSec,&remoteOffsetsValues,leafValuesSec);CHKERRQ(ierr); 4235 ierr = PetscSFCreateSectionSF(coarseToFineEmbedded,rootValuesSec,remoteOffsetsValues,leafValuesSec,&valuesSF);CHKERRQ(ierr); 4236 ierr = PetscSFDestroy(&coarseToFineEmbedded);CHKERRQ(ierr); 4237 ierr = PetscFree(remoteOffsetsValues);CHKERRQ(ierr); 4238 ierr = PetscSectionGetStorageSize(leafValuesSec,&numLeafValues);CHKERRQ(ierr); 4239 ierr = PetscMalloc1(numLeafValues,&leafValues);CHKERRQ(ierr); 4240 ierr = PetscSFBcastBegin(valuesSF,MPIU_SCALAR,rootValues,leafValues);CHKERRQ(ierr); 4241 ierr = PetscSFBcastEnd(valuesSF,MPIU_SCALAR,rootValues,leafValues);CHKERRQ(ierr); 4242 ierr = PetscSFDestroy(&valuesSF);CHKERRQ(ierr); 4243 ierr = PetscFree(rootValues);CHKERRQ(ierr); 4244 ierr = PetscSectionDestroy(&rootValuesSec);CHKERRQ(ierr); 4245 } 4246 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 4247 { 4248 PetscInt maxDof; 4249 PetscInt *rowIndices; 4250 DM refTree; 4251 PetscInt **refPointFieldN; 4252 PetscScalar ***refPointFieldMats; 4253 PetscSection refConSec, refAnSec; 4254 PetscInt pRefStart,pRefEnd,leafStart,leafEnd; 4255 PetscScalar *pointWork; 4256 4257 ierr = PetscSectionGetMaxDof(globalFine,&maxDof);CHKERRQ(ierr); 4258 ierr = DMGetWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 4259 ierr = DMGetWorkArray(fine,maxDof,PETSC_SCALAR,&pointWork);CHKERRQ(ierr); 4260 ierr = DMPlexGetReferenceTree(fine,&refTree);CHKERRQ(ierr); 4261 ierr = DMGetDS(fine,&ds);CHKERRQ(ierr); 4262 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 4263 ierr = DMPlexReferenceTreeGetChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 4264 ierr = DMGetDefaultConstraints(refTree,&refConSec,NULL);CHKERRQ(ierr); 4265 ierr = DMPlexGetAnchors(refTree,&refAnSec,NULL);CHKERRQ(ierr); 4266 ierr = PetscSectionGetChart(refConSec,&pRefStart,&pRefEnd);CHKERRQ(ierr); 4267 ierr = PetscSectionGetChart(leafValuesSec,&leafStart,&leafEnd);CHKERRQ(ierr); 4268 ierr = DMPlexGetHeightStratum(fine,0,&cellStart,&cellEnd);CHKERRQ(ierr); 4269 ierr = DMPlexGetHybridBounds(fine,&cellEndInterior,NULL,NULL,NULL);CHKERRQ(ierr); 4270 cellEnd = (cellEndInterior < 0) ? cellEnd : cellEndInterior; 4271 for (p = leafStart; p < leafEnd; p++) { 4272 PetscInt gDof, gcDof, gOff, lDof; 4273 PetscInt numValues, pValOff; 4274 PetscInt childId; 4275 const PetscScalar *pVal; 4276 const PetscScalar *fvGradData = NULL; 4277 4278 ierr = PetscSectionGetDof(globalFine,p,&gDof);CHKERRQ(ierr); 4279 ierr = PetscSectionGetDof(localFine,p,&lDof);CHKERRQ(ierr); 4280 ierr = PetscSectionGetConstraintDof(globalFine,p,&gcDof);CHKERRQ(ierr); 4281 if ((gDof - gcDof) <= 0) { 4282 continue; 4283 } 4284 ierr = PetscSectionGetOffset(globalFine,p,&gOff);CHKERRQ(ierr); 4285 ierr = PetscSectionGetDof(leafValuesSec,p,&numValues);CHKERRQ(ierr); 4286 if (!numValues) continue; 4287 ierr = PetscSectionGetOffset(leafValuesSec,p,&pValOff);CHKERRQ(ierr); 4288 pVal = &leafValues[pValOff]; 4289 offsets[0] = 0; 4290 offsetsCopy[0] = 0; 4291 newOffsets[0] = 0; 4292 newOffsetsCopy[0] = 0; 4293 childId = cids[p - pStartF]; 4294 if (numFields) { 4295 PetscInt f; 4296 for (f = 0; f < numFields; f++) { 4297 PetscInt rowDof; 4298 4299 ierr = PetscSectionGetFieldDof(localFine,p,f,&rowDof);CHKERRQ(ierr); 4300 offsets[f + 1] = offsets[f] + rowDof; 4301 offsetsCopy[f + 1] = offsets[f + 1]; 4302 /* TODO: closure indices */ 4303 newOffsets[f + 1] = newOffsets[f] + ((childId == -1) ? rowDof : refPointFieldN[childId - pRefStart][f]); 4304 } 4305 ierr = DMPlexGetIndicesPointFields_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,-1,rowIndices);CHKERRQ(ierr); 4306 } 4307 else { 4308 offsets[0] = 0; 4309 offsets[1] = lDof; 4310 newOffsets[0] = 0; 4311 newOffsets[1] = (childId == -1) ? lDof : refPointFieldN[childId - pRefStart][0]; 4312 ierr = DMPlexGetIndicesPoint_Internal(localFine,p,gOff,offsetsCopy,PETSC_FALSE,NULL,rowIndices);CHKERRQ(ierr); 4313 } 4314 if (childId == -1) { /* no child interpolation: one nnz per */ 4315 ierr = VecSetValues(vecFine,numValues,rowIndices,pVal,INSERT_VALUES);CHKERRQ(ierr); 4316 } else { 4317 PetscInt f; 4318 4319 if (grad && p >= cellStart && p < cellEnd) { 4320 numValues -= (dim * (1 + numFVcomps)); 4321 fvGradData = &pVal[numValues]; 4322 } 4323 for (f = 0; f < PetscMax(1,numFields); f++) { 4324 const PetscScalar *childMat = refPointFieldMats[childId - pRefStart][f]; 4325 PetscInt numRows = offsets[f+1] - offsets[f]; 4326 PetscInt numCols = newOffsets[f + 1] - newOffsets[f]; 4327 const PetscScalar *cVal = &pVal[newOffsets[f]]; 4328 PetscScalar *rVal = &pointWork[offsets[f]]; 4329 PetscInt i, j; 4330 4331 #if 0 4332 ierr = PetscInfo6(coarse,"field %D, numFields %D, childId %D, numRows %D, numCols %D, refPointFieldN %D\n",f,numFields,childId,numRows,numCols,refPointFieldN[childId - pRefStart][f]);CHKERRQ(ierr); 4333 #endif 4334 for (i = 0; i < numRows; i++) { 4335 PetscScalar val = 0.; 4336 for (j = 0; j < numCols; j++) { 4337 val += childMat[i * numCols + j] * cVal[j]; 4338 } 4339 rVal[i] = val; 4340 } 4341 if (f == fvField && p >= cellStart && p < cellEnd) { 4342 PetscReal centroid[3]; 4343 PetscScalar diff[3]; 4344 const PetscScalar *parentCentroid = &fvGradData[0]; 4345 const PetscScalar *gradient = &fvGradData[dim]; 4346 4347 ierr = DMPlexComputeCellGeometryFVM(fine,p,NULL,centroid,NULL);CHKERRQ(ierr); 4348 for (i = 0; i < dim; i++) { 4349 diff[i] = centroid[i] - parentCentroid[i]; 4350 } 4351 for (i = 0; i < numFVcomps; i++) { 4352 PetscScalar val = 0.; 4353 4354 for (j = 0; j < dim; j++) { 4355 val += gradient[dim * i + j] * diff[j]; 4356 } 4357 rVal[i] += val; 4358 } 4359 } 4360 ierr = VecSetValues(vecFine,numRows,&rowIndices[offsets[f]],rVal,INSERT_VALUES);CHKERRQ(ierr); 4361 } 4362 } 4363 } 4364 ierr = DMPlexReferenceTreeRestoreChildrenMatrices(refTree,&refPointFieldMats,&refPointFieldN);CHKERRQ(ierr); 4365 ierr = DMRestoreWorkArray(fine,maxDof,PETSC_INT,&rowIndices);CHKERRQ(ierr); 4366 ierr = DMRestoreWorkArray(fine,maxDof,PETSC_SCALAR,&pointWork);CHKERRQ(ierr); 4367 } 4368 ierr = PetscFree(leafValues);CHKERRQ(ierr); 4369 ierr = PetscSectionDestroy(&leafValuesSec);CHKERRQ(ierr); 4370 ierr = PetscFree7(offsets,offsetsCopy,newOffsets,newOffsetsCopy,rowOffsets,numD,numO);CHKERRQ(ierr); 4371 ierr = ISRestoreIndices(aIS,&anchors);CHKERRQ(ierr); 4372 PetscFunctionReturn(0); 4373 } 4374 4375 #undef __FUNCT__ 4376 #define __FUNCT__ "DMPlexTransferVecTree_Inject" 4377 static PetscErrorCode DMPlexTransferVecTree_Inject(DM fine, Vec vecFine, DM coarse, Vec vecCoarse, PetscSF coarseToFine, PetscInt *cids) 4378 { 4379 PetscDS ds; 4380 DM refTree; 4381 PetscSection multiRootSec, rootIndicesSec; 4382 PetscSection globalCoarse, globalFine; 4383 PetscSection localCoarse, localFine; 4384 PetscSection cSecRef; 4385 PetscInt *parentIndices, pRefStart, pRefEnd; 4386 PetscScalar *rootValues, *parentValues; 4387 Mat injRef; 4388 PetscInt numFields, maxDof; 4389 PetscInt pStartC, pEndC, pStartF, pEndF, p; 4390 PetscInt *offsets, *offsetsCopy, *rowOffsets; 4391 PetscLayout rowMap, colMap; 4392 PetscInt rowStart, rowEnd, colStart, colEnd; 4393 PetscScalar ***childrenMats=NULL ; /* gcc -O gives 'may be used uninitialized' warning'. Initializing to suppress this warning */ 4394 PetscErrorCode ierr; 4395 4396 PetscFunctionBegin; 4397 4398 /* get the templates for the fine-to-coarse injection from the reference tree */ 4399 ierr = VecSetOption(vecFine,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4400 ierr = VecSetOption(vecCoarse,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);CHKERRQ(ierr); 4401 ierr = DMPlexGetReferenceTree(coarse,&refTree);CHKERRQ(ierr); 4402 ierr = DMGetDS(coarse,&ds);CHKERRQ(ierr); 4403 ierr = DMSetDS(refTree,ds);CHKERRQ(ierr); 4404 ierr = DMGetDefaultConstraints(refTree,&cSecRef,NULL);CHKERRQ(ierr); 4405 ierr = PetscSectionGetChart(cSecRef,&pRefStart,&pRefEnd);CHKERRQ(ierr); 4406 ierr = DMPlexReferenceTreeGetInjector(refTree,&injRef);CHKERRQ(ierr); 4407 4408 ierr = DMPlexGetChart(fine,&pStartF,&pEndF);CHKERRQ(ierr); 4409 ierr = DMGetDefaultSection(fine,&localFine);CHKERRQ(ierr); 4410 ierr = DMGetDefaultGlobalSection(fine,&globalFine);CHKERRQ(ierr); 4411 ierr = PetscSectionGetNumFields(localFine,&numFields);CHKERRQ(ierr); 4412 ierr = DMPlexGetChart(coarse,&pStartC,&pEndC);CHKERRQ(ierr); 4413 ierr = DMGetDefaultSection(coarse,&localCoarse);CHKERRQ(ierr); 4414 ierr = DMGetDefaultGlobalSection(coarse,&globalCoarse);CHKERRQ(ierr); 4415 ierr = PetscSectionGetMaxDof(localCoarse,&maxDof);CHKERRQ(ierr); 4416 { 4417 PetscInt maxFields = PetscMax(1,numFields) + 1; 4418 ierr = PetscMalloc3(maxFields,&offsets,maxFields,&offsetsCopy,maxFields,&rowOffsets);CHKERRQ(ierr); 4419 } 4420 4421 ierr = DMPlexTransferInjectorTree(coarse,fine,coarseToFine,cids,vecFine,numFields,offsets,&multiRootSec,&rootIndicesSec,NULL,&rootValues);CHKERRQ(ierr); 4422 4423 ierr = PetscMalloc2(maxDof,&parentIndices,maxDof,&parentValues);CHKERRQ(ierr); 4424 4425 /* count indices */ 4426 ierr = VecGetLayout(vecFine,&colMap);CHKERRQ(ierr); 4427 ierr = VecGetLayout(vecCoarse,&rowMap);CHKERRQ(ierr); 4428 ierr = PetscLayoutSetUp(rowMap);CHKERRQ(ierr); 4429 ierr = PetscLayoutSetUp(colMap);CHKERRQ(ierr); 4430 ierr = PetscLayoutGetRange(rowMap,&rowStart,&rowEnd);CHKERRQ(ierr); 4431 ierr = PetscLayoutGetRange(colMap,&colStart,&colEnd);CHKERRQ(ierr); 4432 /* insert values */ 4433 ierr = DMPlexReferenceTreeGetChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 4434 for (p = pStartC; p < pEndC; p++) { 4435 PetscInt numLeaves, leafStart, leafEnd, l, dof, cdof, gOff; 4436 PetscBool contribute = PETSC_FALSE; 4437 4438 ierr = PetscSectionGetDof(globalCoarse,p,&dof);CHKERRQ(ierr); 4439 ierr = PetscSectionGetConstraintDof(globalCoarse,p,&cdof);CHKERRQ(ierr); 4440 if ((dof - cdof) <= 0) continue; 4441 ierr = PetscSectionGetDof(localCoarse,p,&dof);CHKERRQ(ierr); 4442 ierr = PetscSectionGetOffset(globalCoarse,p,&gOff);CHKERRQ(ierr); 4443 4444 rowOffsets[0] = 0; 4445 offsetsCopy[0] = 0; 4446 if (numFields) { 4447 PetscInt f; 4448 4449 for (f = 0; f < numFields; f++) { 4450 PetscInt fDof; 4451 ierr = PetscSectionGetFieldDof(localCoarse,p,f,&fDof);CHKERRQ(ierr); 4452 rowOffsets[f + 1] = offsetsCopy[f + 1] = fDof + rowOffsets[f]; 4453 } 4454 DMPlexGetIndicesPointFields_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,-1,parentIndices);CHKERRQ(ierr); 4455 } 4456 else { 4457 DMPlexGetIndicesPoint_Internal(localCoarse,p,gOff < 0 ? -(gOff + 1) : gOff,offsetsCopy,PETSC_FALSE,NULL,parentIndices);CHKERRQ(ierr); 4458 rowOffsets[1] = offsetsCopy[0]; 4459 } 4460 4461 ierr = PetscSectionGetDof(multiRootSec,p,&numLeaves);CHKERRQ(ierr); 4462 ierr = PetscSectionGetOffset(multiRootSec,p,&leafStart);CHKERRQ(ierr); 4463 leafEnd = leafStart + numLeaves; 4464 for (l = 0; l < dof; l++) parentValues[l] = 0.; 4465 for (l = leafStart; l < leafEnd; l++) { 4466 PetscInt numIndices, childId, offset; 4467 const PetscScalar *childValues; 4468 4469 ierr = PetscSectionGetDof(rootIndicesSec,l,&numIndices);CHKERRQ(ierr); 4470 ierr = PetscSectionGetOffset(rootIndicesSec,l,&offset);CHKERRQ(ierr); 4471 childId = (PetscInt) PetscRealPart(rootValues[offset++]); 4472 childValues = &rootValues[offset]; 4473 numIndices--; 4474 4475 if (childId == -2) { /* skip */ 4476 continue; 4477 } else if (childId == -1) { /* equivalent points: scatter */ 4478 PetscInt m; 4479 4480 contribute = PETSC_TRUE; 4481 for (m = 0; m < numIndices; m++) parentValues[m] = childValues[m]; 4482 } else { /* contributions from children: sum with injectors from reference tree */ 4483 PetscInt parentId, f, lim; 4484 4485 contribute = PETSC_TRUE; 4486 ierr = DMPlexGetTreeParent(refTree,childId,&parentId,NULL);CHKERRQ(ierr); 4487 4488 lim = PetscMax(1,numFields); 4489 offsets[0] = 0; 4490 if (numFields) { 4491 PetscInt f; 4492 4493 for (f = 0; f < numFields; f++) { 4494 PetscInt fDof; 4495 ierr = PetscSectionGetFieldDof(cSecRef,childId,f,&fDof);CHKERRQ(ierr); 4496 4497 offsets[f + 1] = fDof + offsets[f]; 4498 } 4499 } 4500 else { 4501 PetscInt cDof; 4502 4503 ierr = PetscSectionGetDof(cSecRef,childId,&cDof);CHKERRQ(ierr); 4504 offsets[1] = cDof; 4505 } 4506 for (f = 0; f < lim; f++) { 4507 PetscScalar *childMat = &childrenMats[childId - pRefStart][f][0]; 4508 PetscInt n = offsets[f+1]-offsets[f]; 4509 PetscInt i, j; 4510 const PetscScalar *colValues = &childValues[offsets[f]]; 4511 4512 for (i = rowOffsets[f]; i < rowOffsets[f + 1]; i++) { 4513 PetscScalar val = 0.; 4514 for (j = 0; j < n; j++) { 4515 val += childMat[n * i + j] * colValues[j]; 4516 } 4517 parentValues[i] += val; 4518 } 4519 } 4520 } 4521 } 4522 if (contribute) {ierr = VecSetValues(vecCoarse,dof,parentIndices,parentValues,INSERT_VALUES);CHKERRQ(ierr);} 4523 } 4524 ierr = PetscSectionDestroy(&multiRootSec);CHKERRQ(ierr); 4525 ierr = PetscSectionDestroy(&rootIndicesSec);CHKERRQ(ierr); 4526 ierr = PetscFree2(parentIndices,parentValues);CHKERRQ(ierr); 4527 ierr = DMPlexReferenceTreeRestoreChildrenMatrices_Injection(refTree,injRef,&childrenMats);CHKERRQ(ierr); 4528 ierr = PetscFree(rootValues);CHKERRQ(ierr); 4529 ierr = PetscFree3(offsets,offsetsCopy,rowOffsets);CHKERRQ(ierr); 4530 PetscFunctionReturn(0); 4531 } 4532 4533 #undef __FUNCT__ 4534 #define __FUNCT__ "DMPlexTransferVecTree" 4535 /*@ 4536 DMPlexTransferVecTree - transfer a vector between two meshes that differ from each other by refinement/coarsening 4537 that can be represented by a common reference tree used by both. This routine can be used for a combination of 4538 coarsening and refinement at the same time. 4539 4540 collective 4541 4542 Input Parameters: 4543 + dmIn - The DMPlex mesh for the input vector 4544 . vecIn - The input vector 4545 . sfRefine - A star forest indicating points in the mesh dmIn (roots in the star forest) that are parents to points in 4546 the mesh dmOut (leaves in the star forest), i.e. where dmOut is more refined than dmIn 4547 . sfCoarsen - A star forest indicating points in the mesh dmOut (roots in the star forest) that are parents to points in 4548 the mesh dmIn (leaves in the star forest), i.e. where dmOut is more coarsened than dmIn 4549 . cidsRefine - The childIds of the points in dmOut. These childIds relate back to the reference tree: childid[j] = k implies 4550 that mesh point j of dmOut was refined from a point in dmIn just as the mesh point k in the reference 4551 tree was refined from its parent. childid[j] = -1 indicates that the point j in dmOut is exactly 4552 equivalent to its root in dmIn, so no interpolation is necessary. childid[j] = -2 indicates that this 4553 point j in dmOut is not a leaf of sfRefine. 4554 . cidsCoarsen - The childIds of the points in dmIn. These childIds relate back to the reference tree: childid[j] = k implies 4555 that mesh point j of dmIn coarsens to a point in dmOut just as the mesh point k in the reference 4556 tree coarsens to its parent. childid[j] = -2 indicates that point j in dmOut is not a leaf in sfCoarsen. 4557 . useBCs - PETSC_TRUE indicates that boundary values should be inserted into vecIn before transfer. 4558 - time - Used if boundary values are time dependent. 4559 4560 Output Parameters: 4561 . vecOut - Using interpolation and injection operators calculated on the reference tree, the transfered 4562 projection of vecIn from dmIn to dmOut. Note that any field discretized with a PetscFV finite volume 4563 method that uses gradient reconstruction will use reconstructed gradients when interpolating from 4564 coarse points to fine points. 4565 4566 Level: developer 4567 4568 .seealso(): DMPlexSetReferenceTree(), DMPlexGetReferenceTree(), PetscFVGetComputeGradients() 4569 @*/ 4570 PetscErrorCode DMPlexTransferVecTree(DM dmIn, Vec vecIn, DM dmOut, Vec vecOut, PetscSF sfRefine, PetscSF sfCoarsen, PetscInt *cidsRefine, PetscInt *cidsCoarsen, PetscBool useBCs, PetscReal time) 4571 { 4572 PetscErrorCode ierr; 4573 4574 PetscFunctionBegin; 4575 ierr = VecSet(vecOut,0.0);CHKERRQ(ierr); 4576 if (sfRefine) { 4577 Vec vecInLocal; 4578 DM dmGrad = NULL; 4579 Vec faceGeom = NULL, cellGeom = NULL, grad = NULL; 4580 4581 ierr = DMGetLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr); 4582 ierr = VecSet(vecInLocal,0.0);CHKERRQ(ierr); 4583 { 4584 PetscInt numFields, i; 4585 4586 ierr = DMGetNumFields(dmIn, &numFields);CHKERRQ(ierr); 4587 for (i = 0; i < numFields; i++) { 4588 PetscObject obj; 4589 PetscClassId classid; 4590 4591 ierr = DMGetField(dmIn, i, &obj);CHKERRQ(ierr); 4592 ierr = PetscObjectGetClassId(obj, &classid);CHKERRQ(ierr); 4593 if (classid == PETSCFV_CLASSID) { 4594 ierr = DMPlexGetDataFVM(dmIn,(PetscFV)obj,&cellGeom,&faceGeom,&dmGrad);CHKERRQ(ierr); 4595 break; 4596 } 4597 } 4598 } 4599 if (useBCs) { 4600 ierr = DMPlexInsertBoundaryValues(dmIn,PETSC_TRUE,vecInLocal,time,faceGeom,cellGeom,NULL);CHKERRQ(ierr); 4601 } 4602 ierr = DMGlobalToLocalBegin(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr); 4603 ierr = DMGlobalToLocalEnd(dmIn,vecIn,INSERT_VALUES,vecInLocal);CHKERRQ(ierr); 4604 if (dmGrad) { 4605 ierr = DMGetGlobalVector(dmGrad,&grad);CHKERRQ(ierr); 4606 ierr = DMPlexReconstructGradientsFVM(dmIn,vecInLocal,grad);CHKERRQ(ierr); 4607 } 4608 ierr = DMPlexTransferVecTree_Interpolate(dmIn,vecInLocal,dmOut,vecOut,sfRefine,cidsRefine,grad,cellGeom);CHKERRQ(ierr); 4609 ierr = DMRestoreLocalVector(dmIn,&vecInLocal);CHKERRQ(ierr); 4610 if (dmGrad) { 4611 ierr = DMRestoreGlobalVector(dmGrad,&grad);CHKERRQ(ierr); 4612 } 4613 } 4614 if (sfCoarsen) { 4615 ierr = DMPlexTransferVecTree_Inject(dmIn,vecIn,dmOut,vecOut,sfCoarsen,cidsCoarsen);CHKERRQ(ierr); 4616 } 4617 ierr = VecAssemblyBegin(vecOut);CHKERRQ(ierr); 4618 ierr = VecAssemblyEnd(vecOut);CHKERRQ(ierr); 4619 PetscFunctionReturn(0); 4620 } 4621