Re: [petsc-users] indices into Vec/Mat associated to a DMPlex

[Matteo Semplice]

On 15/11/2017 11:39, Matthew Knepley wrote:
On Wed, Nov 15, 2017 at 3:11 AM, Matteo Semplice 
<matteo.sempl...@unito.it <mailto:matteo.sempl...@unito.it>> wrote:

    Hi.

    I am struggling with indices into matrices associated to a DMPLex
    mesh. I can explain my problem to the following minimal example.

    Let's say I want to assemble the matrix to solve an equation (say
    Laplace) with data attached to cells and the finite volume method.
    In principle I

    - loop over the faces of the DMPlex (height=1)

    - for each face, find neighbouring cells (say points i and j in
    the DMPlex) and compute the contributions coming from that face
    (in the example it would be something like (u_j-u_i)*<x_j-x_i,n>
    with x=centroids, n=scaled normal to face and <,> the inner product)

    - insert the contributions +/-<x_j-x_i,n> in rows/columns n(i) and
    n(j) of the matrix where n(k) is the index into Vec/Mat of the
    unknown associated to the k-th cell

    My problem is how to find out n(k). I assume that the Section
    should be able to tell me, but I cannot find the correct function
    to call. I see that FEM methods can use DMPlexMatSetClosure but
    here we'd need a DMPlexMatSetCone...


Everything always reduces to raw Section calls. For instance, for cell c

  PetscSectionGetDof(sec, c, &dof);
  PetscSectionGetOffset(sec, c, &off);

give the number of degrees of freedom on the cell, and the offset into 
storage (a local vector for the local section, and global vector for 
the global section).
The Closure stuff just calls this for every point in the closure.

All right, so the offset is also the (global) index and that is to be 
used in calls to MatSetValues.

Thanks a lot!

    Matteo