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