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