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