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