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