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

[Matthew Knepley]

On Wed, Nov 15, 2017 at 6:09 AM, Matteo Semplice 
wrote:

> On 15/11/2017 11:39, Matthew Knepley wrote:
>
> On Wed, Nov 15, 2017 at 3:11 AM, Matteo Semplice  > 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)* with x=centroids, n=scaled
>> normal to face and <,> the inner product)
>>
>> - insert the contributions +/- 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, );
>   PetscSectionGetOffset(sec, c, );
>
> 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.
>

Not quite. For all owned dofs, it is (in the global section). For all
unowned dofs, it is -(off+1), so you have
to convert it if you want to set off-process values.

  Thanks,

Matt


> Thanks a lot!
>
> Matteo
>



-- 
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/ 

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 
> 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)*
with x=centroids, n=scaled normal to face and <,> the inner product)

- insert the contributions +/- 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, );
  PetscSectionGetOffset(sec, c, );

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

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

[Matthew Knepley]

On Wed, Nov 15, 2017 at 3:11 AM, Matteo Semplice 
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)* with x=centroids, n=scaled
> normal to face and <,> the inner product)
>
> - insert the contributions +/- 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, );
  PetscSectionGetOffset(sec, c, );

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   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
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/