[petsc-users] DMPlex distribution with custom adjacency

[Lawrence Mitchell]

> On 5 Jun 2017, at 23:01, Matthew Knepley <knepley at gmail.com> wrote:
> 
> To do a FEM integral on a facet f I need:
> 
> i) to evaluate coefficients at quadrature points (on the facet)
> ii) to evaluate basis functions at quadrature points (on the facet)
> 
> for (i), I need all the dofs in closure(support(f)).
> 
> So this is a jump term, since I need both sides. Is it nonlinear? If its linear, things are easy
> and just compute from both sides and add, but of course that does not work for nonlinear things.

I can't guarantee what downstream users will write.  Most of the time it will probably be nonlinear in some way.  Hence wanting to compute on the facet (getting both sides of the contribution at once).  Rather than spinning over cells, computing the one-sided facet contribution and adding in.


> 
> for (ii), I just need the definition of the finite element.
> 
> So now, my model for how I want to global assembly of a facet integral is:
> 
> loop over all facets:
>    gather from global coefficient to local data
>    evaluate coefficient at quad points
>    perform integral
>    local to global
> 
> In parallel, I just make a partition of the facets (so that each facet is integrated exactly once).
> 
> OK, so what data do I need in parallel?
> 
> Exactly the dofs that correspond to closure(support(facet)) for all owned facets in my partition.
> 
> So I was hoping to be able to grow a distributed cell partition by exactly that means: add in those remote cells which are in the support of owned facets (I'm happy if this is symmetrically grown, although I think I can do it with one-sided growth).
> 
> So that's my rationale for wanting this "strange" adjacency.  I can get enough adjacency by using the current FEM adjacency and filtering which entities I iterate over, but it seems a bit wasteful.
> 
> I see that you have edges you do not necessarily want, but do they mess up your loop? It seems like you will not encounter them looping over facets.
> This is exactly what happens to me in FV, where I just ignore them.

Remember in 2D the edges are the facets.  I haven't checked what happens in 3D (but I expect it will be similar), because I can't draw 3D meshes.

> If you do really want to prune them, then I guess overriding the DMPlexGetAdjacency() as you propose is probably the best way. I would
> be willing to put it in. Please send me a reminder email since this week is pretty heinous for me.

Sure.  I think this is quite a cute usage because I can make the ghost region "one-sided" quite easily by only growing the adjacency through the facets on the ranks that I need.  So the halo exchange between two processes is not symmetric.  The code I sketched that did this seems to work properly, once the adjacency computation is right.

Lawrence