[petsc-users] flux vector

[Lawrence Mitchell]

> On 11 Jun 2021, at 11:29, Mark Adams <mfadams at lbl.gov> wrote:
> 
> This is a Matt question, but You can set a "Boundary" label (https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexMarkBoundaryFaces.html) and that will label boundaries.
> Next you want the dual of that ...

Note that if you want this to be the _global_ boundary you must call it before mesh distribution.

You need to reduce the support size over the point SF, where you would make a contribution to the count only if the cell is in the owned part of your local cells.

So (pseudocode):

sf = DMGetPointSF(dm);

// since you're only going to reduce over facet points you could create an embedded sf here, but meh.

nroots, nleaves, ilocal, iremote = SFGetGraph(sf);

cstart, cend = DMPlexGetChart(dm);

count = PetscCalloc(nroots); // nroots should be == cend - cstart 

fstart, fend = DMPlexGetDepthStratum(dm, 1); // codim 1 entities

for f = fstart; f < fend; f++ {
  supportsize = DMPlexGetSupportSize(dm, f);
  support = DMPlexGetSupport(dm, f);
  for c in support {
     if (c is owned) { // c not in ilocal
         count[f]++;
     }
  }
}

gcount = PetscMalloc(nroots);
PetscSFReduceBegin(sf, count, gcount, MPI_SUM);
PetscSFReduceEnd(sf, count, gcount, MPI_SUM);
PetscSFBcastBegin(sf, gcount, count, MPI_REPLACE);
PetscSFBcastEnd(sf, gcount, count, MPI_REPLACE);
interior_faces = DMLabelCreate("interior_faces")
for f = fstart; f < fend; f++ {
   if count[f] == 2 {
       DMLabelSetValue(interior_faces, f, 1);
   }
}


(Obviously) untested.

Lawrence