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

Matthew Knepley knepley at gmail.com
Wed Nov 15 04:39:09 CST 2017

```On Wed, Nov 15, 2017 at 3:11 AM, Matteo Semplice <matteo.semplice at 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.

> On a side note, I noticed that the grad field in PetscFVFaceGeom is not
> computed by DMPlexComputeGeometryFVM. Is it meant or could it be used to
> store the precomputed  <x_j-x_i,n> values?

There are separate functions for reconstructing the gradient since it is so
expensive (and can be inaccurate for some unstructured cases at boundaries).

> But even if so, this wouldn't sovle the problem of where to put the
> elements in the matrix, right?

No.

Matt

>
> Matteo

--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their