<div dir="ltr">On Wed, Oct 30, 2013 at 11:13 AM, Bernhard Reinhardt <span dir="ltr"><<a href="mailto:b.reinhardt@physik.uni-muenchen.de" target="_blank">b.reinhardt@physik.uni-muenchen.de</a>></span> wrote:<br><div class="gmail_extra">
<div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">On 28.10.2013 16:43, Matthew Knepley wrote::<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
On Mon, Oct 28, 2013 at 6:17 AM, Bernhard Reinhardt<br>
<<a href="mailto:b.reinhardt@physik.uni-muenchen.de" target="_blank">b.reinhardt@physik.uni-<u></u>muenchen.de</a><br>
<mailto:<a href="mailto:b.reinhardt@physik.uni-muenchen.de" target="_blank">b.reinhardt@physik.<u></u>uni-muenchen.de</a>>> wrote:<br>
<br>
Hi there,<br>
<br>
I am trying to familiarize myself with unstructured grids. For this<br>
I started with the example given in the petsc-manual (mesh of 2<br>
triangles) and tried to modify it according to my needs.<br>
<br>
However, the example seems to be incomplete. At least I was not able<br>
to create a matrix without applying a DMPlexSetDimension.<br>
<br>
In the end I decided to set up an an example as simple as possible:<br>
1D transport of a scalar in a bar with periodic<br>
</blockquote>
<br>
aries.<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<br>
Lets say I have a bar consisting of 3 volumes and 3 faces.<br>
You could write this down as a vector of length 3 containing the<br>
scalar in the volumes and a 3x3 matrix containing the transport<br>
coefficients on the faces.<br>
<br>
I have already tried many different ways of setting this up with<br>
plex but always failed. Attached you find the version that makes<br>
most sense to me. It gives me vectors of length 3 and a 3x3 matrix.<br>
However, the matrix contains no writable entries and I think the<br>
returned closures are not correct.<br>
<br>
Is there a simple example with good documentation of the individual<br>
steps like in the manual which is complete and working?<br>
<br>
<br>
I have modified your example and it is attached. Here are a few points:<br>
<br>
1) You ask for the closure of point 4 in your Vec, and get 1.0. This is<br>
correct.<br>
<br>
The closure of 4 is {4, 1, 2}, the cell + its vertices. There are<br>
only values on cells,<br>
so you get the value on 4, which is 1.0.<br>
<br>
2) The default for preallocation is FEM-style topology, namely that<br>
unknowns are connected to other<br>
unknowns with overlapping support. The support of your unknowns is<br>
only the cell itself, so there<br>
are no connections by default. I put in a call to<br>
DMPlexSetPreallocationCenterDi<u></u>mension(), which<br>
is the dimension of intermediate points. In FEM these are cells,<br>
but in FVM they are vertices (or<br>
faces). I set your topology to FVM, and you should get the<br>
preallocation you want.<br>
</blockquote>
<br>
Dear Matt,<br>
<br>
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.<br>
<br>
However, I still don't get what DMPlexMatSetClosure does.</blockquote><div><br></div><div>The words used are precise. Closure means the transitive closure of a point in the graph. For example, the</div><div>closure of point 7 is</div>
<div><br></div><div> closure(7) = {7, 2, 3}</div><div><br></div><div>MatSet() takes a set of indices and a set of values and updates a Mat. Now MatSetClosure(), first calculates the indices</div><div><br></div><div> for q in closure(p):</div>
<div> PetscSectionGetOffset(s, q, &off)</div><div> PetscSectionGetDof(s, q, &dof)</div><div> indices += [ind, ind+dof)</div><div> MatSetValues(indices, values)</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
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).<br>
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.<br></blockquote><div><br></div><div>What you seem to be suggesting is FD, which you can do on a structured grid using DMDA. On an unstructured grid, what I am familiar</div>
<div>with is FV, which means looping over faces instead of cells, and doing something like MatSetStar() instead of MatSetClosure().</div><div><br></div><div>If you are using some unstructured method which sets whole rows of the Jacobian at once (like FD), then you would want something like</div>
<div>MatSetStarClosure(), and then you would need to add the flag for ghost cells in DMPlexDistribute(). I have never done this since all</div><div>methods I have implemented are more local.</div><div><br></div><div>It would help to understand what you really want to do.</div>
<div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
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.<br>
<br>
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<br>
<br>
- - - - - (without any faces)<br>
or<br>
<br>
|-||-||-||-|-| (with double faces)<br>
or<br>
_ _ _ _ _<br>
|-||-||-||-|-| (with double faces and an additional "internal" face)<br>
<br>
In the best case I get the foreseen allocation, but am stil am not able to set the right values in the matrix.<br>
<br>
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?<br>
<br>
Best regards!<span class="HOEnZb"><font color="#888888"><br>
<br>
Bernhard<br>
</font></span></blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>