[petsc-users] flux vector

[Lawrence Mitchell]

> On 11 Jun 2021, at 13:19, Matthew Knepley <knepley at gmail.com> 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, 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.

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.


>   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