On Sat, Jun 12, 2021 at 5:25 AM Lawrence Mitchell <[email protected]> wrote:
> > > On 11 Jun 2021, at 13:19, Matthew Knepley <[email protected]> wrote: > > > > Before doing what Lawrence suggests, I would like to talk this through. > I don't think we should need communication. So, > > > > Goal: Section with dofs on each face that has support 2 (separates 2 > cells) > > > > Okay, suppose we loop over our local faces and get the support size of > each. Boundary faces are disqualified, which is good, > > but if a face is on the boundary with another process, it will also have > support size 1. Thus, let us distribute with overlap = 1. > > Now if we loop over that same face set, we will get the "right answer" > for the support size. > > > > However, using overlap = 1 put in a bunch of new faces. We do not care > about the ones on the process boundary. They will be > > handled by the other process. We do care about faces between two ghost > cells, since they will be a false positive. Luckily, these > > are labeled by the "ghost" label. So I think this is our algorithm: > > Ah, I missed that the mesh is already overlapped. > > > 1) DMClone(dm, &fdm); > > > > 2) DMGetLocalSection(fdm, &ls); > > > > 3) Mark each face with 2 neighbors that is not a ghost > > > > DMGetLabel(dm, "ghost", &ghostLabel); > > This label is only created if one calls DMPlexConstructGhostCells, Shoot, that is the wrong workflow. > I don't understand this code really, but it seems like if I call > DMPlexDistribute(dm, overlap=1, &newdm) and then DMPlexConstructGhostCells, > the "depth 2" faces/etc... will be marked with the ghost label, but the > "depth 1" faces will not be. > Yes, I only marked those because those are the only ones you might calculate fluxes on. > That said, I agree with Matt's algorithm, and if the ghost label isn't > available or set up right you can replace > > DMLabelGetValue(ghostLabel, f, &ghost) > > by checking if the point f is in the ilocal array of the DM's pointSF, I > think. > It is a little more complicated. The check for "ghost" is that you have 2 cells in the support, but they are both not owned, meaning both are in the SF. Thanks, Matt > > DMGetHeightStratum(fdm, 1, &fStart, &fEnd); > > for (f = fStartl f < fEnd; ++f) { > > PetscInt supportSize, ghost; > > > > DMGetSupportSize(fdm, f, &supportSize); > > DMLabelGetValue(ghostLabel, f, &ghost); > > if (ghost < 0 && supportSize == 2) PetscSectionSetDof(ls, f, 1); > > } > > > > 4) The global section will be created automatically when we ask for it > > > > DMGetGlobalSection(fdm, &s); > > > > Let me know if this does not work. > > > Lawrence -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
