<div class="gmail_quote">On Fri, Jun 17, 2011 at 23:00, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
It is possible that list of elements is not correct for periodic boundary conditions; I've never used it. I'm forwarding this to petsc-dev so Matt can check it and see if the elements make sense and what indices they routine (are the local, global? I don't know) if they are local then maybe you don't need to the MatSetValuesStencil() and just call MatSetValuesLocal().</blockquote>
</div><br><div>The function returns local indices, so you would use MatSetValuesLocal(). I put a link to these in the documentation.</div><div><br></div><div>But the implementation is not correct for periodic domains. It should not be a great deal of work to implement. You would just need to teach that function about the "wrap-around" elements. But this is actually not trivial for a mapped grid. In no case can you just use the coordinates that were set on the DMDA, you have to extrapolate. But it's unclear how to extrapolate for a non-uniform mesh in order to make the periodic problem consistent. I don't think this is easy to rectify without changing the API for managing coordinates, but you can certainly make it work for regular Cartesian grids.</div>
<div><br></div><div>By the way, I think the whole premise of DMDAGetElements() is flawed. Shape functions are not provided so the user has to implement those themselves and manage consistency. The actual logic inside DMDAGetElements() is really small compared to defining quadratures and basis functions, so it's no particular effort for the user to manage it on their own. Additionally, managing it on their own will perform somewhat better (no intermediate indexing structures) and gives more flexibility in boundary conditions.</div>
<div><br></div><div>It might not be a bad idea to offer an easy-to-use, first order only, structured grid FEM API, but this ain't it.</div>