[petsc-users] DMPlex example with manual grid layout

[Bernhard Reinhardt]

On 28.10.2013 16:43, Matthew Knepley wrote::
> On Mon, Oct 28, 2013 at 6:17 AM, Bernhard Reinhardt
> <b.reinhardt at physik.uni-muenchen.de
> <mailto:b.reinhardt at physik.uni-muenchen.de>> wrote:
>
>     Hi there,
>
>     I am trying to familiarize myself with unstructured grids. For this
>     I started with the example given in the petsc-manual (mesh of 2
>     triangles) and tried to modify it according to my needs.
>
>     However, the example seems to be incomplete. At least I was not able
>     to create a matrix without applying a DMPlexSetDimension.
>
>     In the end I decided to set up an an example as simple as possible:
>     1D transport of a scalar in a bar with periodic

aries.
>
>     Lets say I have a bar consisting of 3 volumes and 3 faces.
>     You could write this down as a vector of length 3 containing the
>     scalar in the volumes and a 3x3 matrix containing the transport
>     coefficients on the faces.
>
>     I have already tried many different ways of setting this up with
>     plex but always failed. Attached you find the version that makes
>     most sense to me. It gives me vectors of length 3 and a 3x3 matrix.
>     However, the matrix contains no writable entries and I think the
>     returned closures are not correct.
>
>     Is there a simple example with good documentation of the individual
>     steps like in the manual which is complete and working?
>
>
> I have modified your example and it is attached. Here are a few points:
>
> 1) You ask for the closure of point 4 in your Vec, and get 1.0. This is
> correct.
>
>    The closure of 4 is {4, 1, 2}, the cell + its vertices. There are
> only values on cells,
>    so you get the value on 4, which is 1.0.
>
> 2) The default for preallocation is FEM-style topology, namely that
> unknowns are connected to other
>       unknowns with overlapping support. The support of your unknowns is
> only the cell itself, so there
>       are no connections by default. I put in a call to
> DMPlexSetPreallocationCenterDimension(), which
>       is the dimension of intermediate points. In FEM these are cells,
> but in FVM they are vertices (or
>       faces). I set your topology to FVM, and you should get the
> preallocation you want.

Dear Matt,

thank you for your advice. Now I get the right allocation for my 
exemplary problem (although I have to admit I´m not entirely convinced 
for the right reasons ;-). I expanded the bar from 3 to 5 volumes, so 
that the allocation becomes clearer. See attached file.

However, I still don't get what DMPlexMatSetClosure does. What I would 
like to do, is set the operator into the matrix which computes the 
change of the value at point p. So let's say in my example the Laplacian 
operator (1, -2, 1).
E.g. Delta_p7 = A_21 p6 + A_22 p7 + A_23 p8; with A_21=1, A_22 = -2, 
A_23 = 1. However, in my example MatClosure(p7) consits only of A_22.

Probably my sieve-layout is not appropriate. I recognize that the 
Laplacian consists of 3 elements, but my volumes have only a cone size 
of 2. Nevertheless I hoped to be able to set at least A_21 and A_23.

I tried all combinations of dof=1 or dof=0 at the faces or volumes and 
PreallocationCenterDimension 0 or 1 without elucidation. I also tried 
different mesh layouts without much succes. That is

- - - - - (without any faces)
or

|-||-||-||-|-| (with double faces)
or
  _  _  _  _ _
|-||-||-||-|-| (with double faces and an additional "internal" face)

In the best case I get the foreseen allocation, but am stil am not able 
to set the right values in the matrix.

This leads me to the question: Is the sieve-layout appropriate 
(|-|-|-|-|-)? If yes, what is the envisaged way to set the operator for 
point p in the matrix?

Best regards!

Bernhard
-------------- next part --------------
A non-text attachment was scrubbed...
Name: main.f90
Type: text/x-fortran
Size: 5714 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131030/5a2c3a83/attachment.bin>