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. 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()
signature.asc
Description: OpenPGP digital signature
