<div dir="ltr"><div dir="ltr">On Fri, Jun 11, 2021 at 7:38 AM Lawrence Mitchell <<a href="mailto:wence@gmx.li">wence@gmx.li</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br>
<br>
> On 11 Jun 2021, at 11:29, Mark Adams <<a href="mailto:mfadams@lbl.gov" target="_blank">mfadams@lbl.gov</a>> wrote:<br>
> <br>
> This is a Matt question, but You can set a "Boundary" label (<a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexMarkBoundaryFaces.html" rel="noreferrer" target="_blank">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMPLEX/DMPlexMarkBoundaryFaces.html</a>) and that will label boundaries.<br>
> Next you want the dual of that ...<br></blockquote><div><br></div><div>Before doing what Lawrence suggests, I would like to talk this through. I don't think we should need communication. So,</div><div><br></div><div> Goal: Section with dofs on each face that has support 2 (separates 2 cells)</div><div><br></div><div>Okay, suppose we loop over our local faces and get the support size of each. Boundary faces are disqualified, which is good,</div><div>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.</div><div>Now if we loop over that same face set, we will get the "right answer" for the support size.</div><div><br></div><div>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</div><div>handled by the other process. We do care about faces between two ghost cells, since they will be a false positive. Luckily, these</div><div>are labeled by the "ghost" label. So I think this is our algorithm:</div><div><br></div><div>1) DMClone(dm, &fdm);</div><div><br></div><div>2) DMGetLocalSection(fdm, &ls);</div><div><br></div><div>3) Mark each face with 2 neighbors that is not a ghost</div><div><br></div><div> DMGetLabel(dm, "ghost", &ghostLabel);</div><div> DMGetHeightStratum(fdm, 1, &fStart, &fEnd);</div><div> for (f = fStartl f < fEnd; ++f) {</div><div> PetscInt supportSize, ghost;<br></div><div><br></div><div> DMGetSupportSize(fdm, f, &supportSize);</div><div> DMLabelGetValue(ghostLabel, f, &ghost);</div><div> if (ghost < 0 && supportSize == 2) PetscSectionSetDof(ls, f, 1);</div><div> }</div><div><br></div><div>4) The global section will be created automatically when we ask for it</div><div><br></div><div> DMGetGlobalSection(fdm, &s);</div><div><br></div><div>Let me know if this does not work.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Note that if you want this to be the _global_ boundary you must call it before mesh distribution.<br>
<br>
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.<br>
<br>
So (pseudocode):<br>
<br>
sf = DMGetPointSF(dm);<br>
<br>
// since you're only going to reduce over facet points you could create an embedded sf here, but meh.<br>
<br>
nroots, nleaves, ilocal, iremote = SFGetGraph(sf);<br>
<br>
cstart, cend = DMPlexGetChart(dm);<br>
<br>
count = PetscCalloc(nroots); // nroots should be == cend - cstart <br>
<br>
fstart, fend = DMPlexGetDepthStratum(dm, 1); // codim 1 entities<br>
<br>
for f = fstart; f < fend; f++ {<br>
supportsize = DMPlexGetSupportSize(dm, f);<br>
support = DMPlexGetSupport(dm, f);<br>
for c in support {<br>
if (c is owned) { // c not in ilocal<br>
count[f]++;<br>
}<br>
}<br>
}<br>
<br>
gcount = PetscMalloc(nroots);<br>
PetscSFReduceBegin(sf, count, gcount, MPI_SUM);<br>
PetscSFReduceEnd(sf, count, gcount, MPI_SUM);<br>
PetscSFBcastBegin(sf, gcount, count, MPI_REPLACE);<br>
PetscSFBcastEnd(sf, gcount, count, MPI_REPLACE);<br>
interior_faces = DMLabelCreate("interior_faces")<br>
for f = fstart; f < fend; f++ {<br>
if count[f] == 2 {<br>
DMLabelSetValue(interior_faces, f, 1);<br>
}<br>
}<br>
<br>
<br>
(Obviously) untested.<br>
<br>
Lawrence<br>
<br>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>