On Wed, Nov 15, 2017 at 6:09 AM, Matteo Semplice <matteo.sempl...@unito.it> wrote:
> On 15/11/2017 11:39, Matthew Knepley wrote: > > On Wed, Nov 15, 2017 at 3:11 AM, Matteo Semplice <matteo.sempl...@unito.it > > 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)*<x_j-x_i,n> with x=centroids, n=scaled >> normal to face and <,> the inner product) >> >> - 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 >> >> 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, &dof); > PetscSectionGetOffset(sec, c, &off); > > 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/ <http://www.caam.rice.edu/~mk51/>