[petsc-users] DMPlex distribution with FVM adjacency

[Lawrence Mitchell]

On 25/05/17 16:25, Matthew Knepley wrote:
> On Thu, May 25, 2017 at 9:27 AM, Lawrence Mitchell
> <lawrence.mitchell at imperial.ac.uk
> <mailto:lawrence.mitchell at imperial.ac.uk>> 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()

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 473 bytes
Desc: OpenPGP digital signature
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170525/11c00049/attachment.pgp>