<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, May 25, 2017 at 11:27 AM, Lawrence Mitchell <span dir="ltr"><<a href="mailto:lawrence.mitchell@imperial.ac.uk" target="_blank">lawrence.mitchell@imperial.ac.uk</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class="">On 25/05/17 16:25, Matthew Knepley wrote:<br>
> On Thu, May 25, 2017 at 9:27 AM, Lawrence Mitchell<br>
> <<a href="mailto:lawrence.mitchell@imperial.ac.uk">lawrence.mitchell@imperial.<wbr>ac.uk</a><br>
</span><div><div class="h5">> <mailto:<a href="mailto:lawrence.mitchell@imperial.ac.uk">lawrence.mitchell@<wbr>imperial.ac.uk</a>>> wrote:<br>
><br>
> Dear petsc-users,<br>
><br>
> I am trying to distribute a triangle mesh with a cell halo defined by<br>
> FVM adjacency (i.e. if I have a facet in the initial (0-overlap)<br>
> distribution, I want the cell on the other side of it).<br>
><br>
> Reading the documentation, I think I do:<br>
><br>
> DMPlexSetAdjacencyUseCone(<wbr>PETSC_TRUE)<br>
> DMPlexSetAdjacencyUseClosure(<wbr>PETSC_FALSE)<br>
><br>
> and then<br>
> DMPlexDistribute(..., ovelap=1)<br>
><br>
> If I do this for a simple mesh and then try and do anything on it, I<br>
> run into all sorts of problems because I have a plex where I have some<br>
> facets, but not even one cell in the support of the facet. Is this to<br>
> be expected?<br>
><br>
><br>
> Hmm. I don't think so. You should have at least one cell in the<br>
> support of every facet.<br>
> TS ex11 works exactly this way.<br>
><br>
> When using that adjacency, the overlap cells you get will not have<br>
> anything but the<br>
> facet connecting them to that partition. Although, if you have<br>
> adjacent cells in that overlap layer,<br>
> you can get ghost faces between those.<br>
><br>
> With the code below, do you get an interpolated mesh when you create<br>
> it there. That call in C<br>
> has another argument<br>
><br>
> <a href="http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/DMPLEX/DMPlexCreateFromCellList.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-master/docs/manualpages/<wbr>DMPLEX/<wbr>DMPlexCreateFromCellList.html</a><br>
<br>
</div></div>The mesh is interpolated.<br>
<br>
<br>
OK, so let's see if I can understand what the different adjacency<br>
relations are:<br>
<br>
usecone=False, useclosure=False:<br>
<br>
adj(p) => cone(p) + cone(support(p))<br>
<br>
usecone=True, useclosure=False:<br>
<br>
adj(p) => support(p) + support(cone(p))<br>
<br>
usecone=False, useclosure=True<br>
<br>
adj(p) => closure(star(p))<br>
<br>
usecone=True, useclosure=True<br>
<br>
adj(p) => star(closure(p))<br>
<br>
So let's imagine I have a facet f, the adjacent points are the<br>
support(cone(f)) so the support of the vertices in 2D, so those are<br>
some new facets.<br></blockquote><div><br></div><div>If you want that, is there a reason you cannot use the FEM style FALSE+TRUE?</div><div>If you already want the closure, usually the star is not really adding anything new.</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
So now, following <a href="https://arxiv.org/pdf/1506.06194.pdf" rel="noreferrer" target="_blank">https://arxiv.org/pdf/1506.<wbr>06194.pdf</a>, I need to<br>
complete this new mesh, so I ask for the closure of these new facets.<br>
But that might mean I won't ask for cells, right? So I think I would<br>
end up with some facets that don't have any support. And empirically<br>
I observe that:<br>
<br>
e.g. the code attached:<br>
<br>
$ mpiexec -n 3 python bar.py<br>
[0] 7 [0]<br>
[0] 8 [0]<br>
[0] 9 [0 1]<br>
[0] 10 [1]<br>
[0] 11 [1]<br>
[0] 12 []<br>
[1] 10 [0 2]<br>
[1] 11 [0 1]<br>
[1] 12 [0]<br>
[1] 13 [1]<br>
[1] 14 [2]<br>
[1] 15 [2]<br>
[1] 16 [1 3]<br>
[1] 17 [3]<br>
[1] 18 [3]<br>
[2] 7 [0 1]<br>
[2] 8 [0]<br>
[2] 9 [0]<br>
[2] 10 [1]<br>
[2] 11 []<br>
[2] 12 [1]<br>
<br>
<br>
What I would like (although I'm not sure if this is supported right<br>
now), is the overlap to contain closure(support(facet)) for all shared<br>
facets. I think that's equivalent to closure(support(p)) \forall p.<br>
<br>
That way on any shared facets, I have both cells and their closure.<br>
<br>
Is that easy to do?<br>
<br>
Lawrence<br>
<br>
import sys, petsc4py<br>
petsc4py.init(sys.argv)<br>
<span class="">from petsc4py import PETSc<br>
import numpy as np<br>
Lx = Ly = 1<br>
</span>nx = 1<br>
ny = 2<br>
<span class=""><br>
xcoords = np.linspace(0.0, Lx, nx + 1, dtype=PETSc.RealType)<br>
ycoords = np.linspace(0.0, Ly, ny + 1, dtype=PETSc.RealType)<br>
coords = np.asarray(np.meshgrid(<wbr>xcoords, ycoords)).swapaxes(0,<br>
2).reshape(-1, 2)<br>
<br>
# cell vertices<br>
i, j = np.meshgrid(np.arange(nx, dtype=PETSc.IntType), np.arange(ny,<br>
dtype=PETSc.IntType))<br>
cells = [i*(ny+1) + j, i*(ny+1) + j+1, (i+1)*(ny+1) + j+1,<br>
(i+1)*(ny+1) + j]<br>
cells = np.asarray(cells, dtype=PETSc.IntType).swapaxes(<wbr>0,<br>
2).reshape(-1, 4)<br>
idx = [0, 1, 3, 1, 2, 3]<br>
cells = cells[:, idx].reshape(-1, 3)<br>
<br>
comm = PETSc.COMM_WORLD<br>
if comm.rank == 0:<br>
dm = PETSc.DMPlex().<wbr>createFromCellList(2, cells, coords,<br>
</span>interpolate=True, comm=comm)<br>
<span class="">else:<br>
dm = PETSc.DMPlex().<wbr>createFromCellList(2, np.zeros((0,<br>
</span>cells.shape[1]), dtype=PETSc.IntType),<br>
np.zeros((0, 2),<br>
dtype=PETSc.RealType),<br>
interpolate=True,<br>
comm=comm)<br>
<br>
dm.setAdjacencyUseClosure(<wbr>False)<br>
dm.setAdjacencyUseCone(True)<br>
<br>
dm.distribute(overlap=1)<br>
sf = dm.getPointSF()<br>
<br>
for p in range(*dm.getDepthStratum(dm.<wbr>getDepth()-1)):<br>
PETSc.Sys.syncPrint("[%d] %d %s" % (comm.rank, p, dm.getSupport(p)))<br>
<br>
PETSc.Sys.syncFlush()<br>
<br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature"><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.caam.rice.edu/~mk51/" target="_blank">http://www.caam.rice.edu/~mk51/</a><br></div></div></div>
</div></div>