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