<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Wed, Nov 15, 2017 at 3:11 AM, Matteo Semplice <span dir="ltr"><<a href="mailto:matteo.semplice@unito.it" target="_blank">matteo.semplice@unito.it</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi.<br>
<br>
I am struggling with indices into matrices associated to a DMPLex mesh. I can explain my problem to the following minimal example.<br>
<br>
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<br>
<br>
- loop over the faces of the DMPlex (height=1)<br>
<br>
- 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)<br>
<br>
- 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<br>
<br>
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...<br></blockquote><div><br></div><div>Everything always reduces to raw Section calls. For instance, for cell c</div><div><br></div><div>  PetscSectionGetDof(sec, c, &dof);</div><div>  PetscSectionGetOffset(sec, c, &off);</div><div><br></div><div>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).</div><div>The Closure stuff just calls this for every point in the closure.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
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?</blockquote><div><br></div><div>There are separate functions for reconstructing the gradient since it is so expensive (and can be inaccurate for some unstructured cases at boundaries).</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> But even if so, this wouldn't sovle the problem of where to put the elements in the matrix, right?</blockquote><div><br></div><div>No.</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"><span class="HOEnZb"><font color="#888888"><br>
Matteo</font></span></blockquote></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div>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><br></div><div><a href="http://www.caam.rice.edu/~mk51/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div>
</div></div>