[petsc-users] Re-ordering in DMPlexCreateFromCellListParallelPetsc

[Nicolas Barral]

Thanks for your answers Matt.

On 22/03/2021 13:45, Matthew Knepley wrote:
> On Mon, Mar 22, 2021 at 6:22 AM Nicolas Barral 
> <nicolas.barral at math.u-bordeaux.fr 
> <mailto:nicolas.barral at math.u-bordeaux.fr>> wrote:
> 
>     On 21/03/2021 21:29, Matthew Knepley wrote:
>      > On Sat, Mar 20, 2021 at 10:07 AM Nicolas Barral
>      > <nicolas.barral at math.u-bordeaux.fr
>     <mailto:nicolas.barral at math.u-bordeaux.fr>
>      > <mailto:nicolas.barral at math.u-bordeaux.fr
>     <mailto:nicolas.barral at math.u-bordeaux.fr>>> wrote:
>      >
>      >     Hi all,
>      >
>      >     I'm building a plex from elements arrays using
>      >     DMPlexCreateFromCellListParallelPetsc. Once the plex is built, I
>      >     need to
>      >     set up boundary labels. I have an array of faces containing a
>     series of
>      >     3 vertex local indices. To rebuild boundary labels, I need to
>     loop over
>      >     the array and get the join of 3 consecutive points to find the
>      >     corresponding face point in the DAG.
>      >
>      >
>      > This is very common. We should have a built-in thing that does this.
>      >
>     Ordering apart, it's not very complicated once you figure out "join" is
>     the right operation to use. We need more doc on graph layout and
>     operations in DMPlex, I think I'm going to make pictures when I'm done
>     with that code because I waste too much time every time. Is there a
>     starting point you like ?
> 
> 
> I think the discussion in
> 
> @article{LangeMitchellKnepleyGorman2015,
>    title     = {Efficient mesh management in {Firedrake} using 
> {PETSc-DMPlex}},
>    author    = {Michael Lange and Lawrence Mitchell and Matthew G. 
> Knepley and Gerard J. Gorman},
>    journal   = {SIAM Journal on Scientific Computing},
>    volume    = {38},
>    number    = {5},
>    pages     = {S143--S155},
>    eprint    = {http://arxiv.org/abs/1506.07749},
>    doi       = {10.1137/15M1026092},
>    year      = {2016}
> }
> 
> is pretty good.
> 
It is, I'd just like to have something more complete (meet, join, 
height, depth...) and with more 2D & 3D pictures. It's all information 
available somewhere, but I would find it convenient to be all at the 
same place. Are sources of the paper available somewhere ?


>      >     Problem, vertices get reordered by
>      >     DMPlexCreateFromCellListParallelPetsc
>      >     so that locally owned vertices are before remote ones, so
>     local indices
>      >     are changed and the indices in the face array are not good
>     anymore.
>      >
>      >
>      > This is not exactly what happens. I will talk through the
>     algorithm so
>      > that maybe we can find a good
>      > interface. I can probably write the code quickly:
>      >
>      > 1) We take in cells[numCells, numCorners], which is a list of all
>     the
>      > vertices in each cell
>      >
>      >       The vertex numbers do not have to be a contiguous set. You can
>      > have any integers you want.
>      >
>      > 2) We create a sorted list of the unique vertex numbers on each
>     process.
>      > The new local vertex numbers
>      >       are the locations in this list.
>      >
> 
>     Ok It took me re-writing this email a couple times but I think I
>     understand. I was too focused on local/global indices. But if I get
>     this
>     right, you still make an assumption: that the numVertices*dim
>     coordinates passed in vertexCoords are the coordinates of the
>     numvertices first vertices in the sorted list. Is that true ?
> 
> 
> No. It can be arbitrary. That is why we make the vertexSF, so we can map 
> those coordinates back to the right processes.
> 
So how do you know which coordinates correspond to which vertex since no 
map is explicitly provided ?



>      > Here is my proposed interface. We preserve this list of unique
>     vertices,
>      > just as we preserve the vertexSF.
>      > Then after DMPlexCreateFromCellListParallelPetsc(), can
>      > DMPlexInterpolate(), you could call
>      >
>      >    DMPlexBuildFaceLabelsFromCellList(dm, numFaces, faces, labelName,
>      > labelValues)
>      >
>      > Would that work for you? I think I could do that in a couple of
>     hours.
>      >
> 
>     So that function would be very helpful. But is it as simple as you're
>     thinking ? The sorted list gives local index -> unique identifier, but
> 
>     what we need is the other way round, isn't it ?
> 
> 
> You get the other way using search.
> 
Do you mean a linear search for each vertex ?

Thanks,

-- 
Nicolas

>    Thanks,
> 
>       Matt
> 
>     Once I understand better, I can have a first draft in my plexadapt code
>     and we pull it out later.
> 
>     Thanks
> 
>     -- 
>     Nicolas
> 
> 
>      >    Thanks,
>      >
>      >       Matt
>      >
>      >     Is there a way to track this renumbering ? For owned
>     vertices, I can
>      >     find the local index from the global one (so do old local
>     index ->
>      >     global index -> new local index). For the remote ones, I'm
>     not sure. I
>      >     can hash global indices, but is there a more idiomatic way ?
>      >
>      >     Thanks,
>      >
>      >     --
>      >     Nicolas
>      >
>      >
>      >
>      > --
>      > 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
>      >
>      > https://www.cse.buffalo.edu/~knepley/
>     <http://www.cse.buffalo.edu/~knepley/>
> 
> 
> 
> -- 
> 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
> 
> https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>