[petsc-users] DMPlex distribution with FVM adjacency

[Lawrence Mitchell]

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?

For example the following petsc4py code breaks when run on 3 processes:

$ mpiexec -n 3 python bork.py
[1] DMPlexGetOrdering() line 133 in
/data/lmitche1/src/deps/petsc/src/dm/impls/plex/plexreorder.c
[1] DMPlexCreateOrderingClosure_Static() line 41 in
/data/lmitche1/src/deps/petsc/src/dm/impls/plex/plexreorder.c
[1] Petsc has generated inconsistent data
[1] Number of depth 2 faces 34 does not match permuted nubmer 29
: error code 77
[2] DMPlexGetOrdering() line 133 in
/data/lmitche1/src/deps/petsc/src/dm/impls/plex/plexreorder.c
[2] DMPlexCreateOrderingClosure_Static() line 41 in
/data/lmitche1/src/deps/petsc/src/dm/impls/plex/plexreorder.c
[2] Petsc has generated inconsistent data
[2] Number of depth 2 faces 33 does not match permuted nubmer 28
: error code 77
[0] DMPlexGetOrdering() line 133 in
/data/lmitche1/src/deps/petsc/src/dm/impls/plex/plexreorder.c
[0] DMPlexCreateOrderingClosure_Static() line 41 in
/data/lmitche1/src/deps/petsc/src/dm/impls/plex/plexreorder.c
[0] Petsc has generated inconsistent data
[0] Number of depth 2 faces 33 does not match permuted nubmer 31

$ cat > bork.py<<\EOF
from petsc4py import PETSc
import numpy as np
Lx = Ly = 1
nx = ny = 4

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, comm=comm)
else:
    dm = PETSc.DMPlex().createFromCellList(2, np.zeros((0, 4),
dtype=PETSc.IntType),
                                           np.zeros((0, 2),
dtype=PETSc.RealType),
                                           comm=comm)

dm.setAdjacencyUseClosure(False)
dm.setAdjacencyUseCone(True)

dm.distribute(overlap=1)

dm.getOrdering(PETSc.Mat.OrderingType.RCM)

dm.view()
EOF

Am I doing something wrong?  Is this not expected to work?

Cheers,

Lawrence

-------------- 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/7fb54dbd/attachment.pgp>