<div class="gmail_extra">On Fri, Nov 9, 2012 at 5:12 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div id=":4g2"> You are "down defining" IS so you can say it doesn't support doing these things. I am proposing "up defining" IS to be able to do these things.<br>
<br>
      I consider IS to be _THE_  PETSc tool for "indexing  stuff" (and want it to be the only tool, having two tools is redundant and confusing). When IS doesn't have a capability for some "kind of indexing" we wish to do then I think the right thing to do is generalize IS as needed, not introduce a different class.<br>
</div></blockquote><div><br></div><div>Matt and I have been debating a similar point independently. I think we concluded that we could make do with index sets, DMImpl-based vector access, and and sub-DMs (for splitting out a piece of a larger problem, e.g., to solve an auxiliary problem along a fault).</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div id=":4g2">
<br>
      Consider the normal scenario for finite elements:  (simple case)<br>
<br>
          Loop over elements:<br>
               Get the vector entries Ue from U for the element<br>
               Compute with Ue producing Fe<br>
               Put the Fe into the output vector F.<br>
<br>
        this is an example of indexing.  Plain old ISGeneral doesn't directly provide a way to do this,  Having a concept of a collection of groups of indices (of varying sizes for different element types) makes this kind of finite element scenario straight forward.  Your PetscSection seems like yet a different type of useful indexing.<br>

<br>
    What I propose is we "start at the beginning" using all the ideas we've gained over the years and design a new IS (which we can call PetscIS (or something else)) using good abstractions of "indexing stuff". Once we have this good design<br>

IS will quickly be removed and replaced with the new IS construct. I am open to any idea on what the model (set of abstractions, whatever ..) we use for the new IS.<br></div></blockquote><div><br></div><div>The key feature is accessing a range by block number. In the case of Section, the linear scalar indices (from concatenating the scalar entries for each block) is not used directly.</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div id=":4g2">
<br>
    Though you use fancy abstract words with PetscSection like Section and Atlas you actually have only one concrete specific implementation that is merely block size and ranges. Is that really all we want/need for a very general and useful concept of efficient indexing?  Perhaps it is, but let's explore it. For example in pseudo code like above show me how you would do the normal scenario for finite elements using your PetscSection?</div>
</blockquote></div><br></div><div class="gmail_extra">Matt does not explicitly store the vector indices with support on the element. Instead he has the "closure" operation which, given an cell number, gathers all the lower-dimensional topological "points" (faces, edges, and vertices) using a breadth-first traversal. (This is something that might make sense to cache explicitly for each element. It uses more memory, but it makes the access pattern simpler so might pay off. Most FEM codes do that.)<br>
<br>Each of those "points" in the closure has zero or more dofs in the Vec, as indicated by the "default" Section.</div><div class="gmail_extra"><br></div><div class="gmail_extra">for each element e:</div>
<div class="gmail_extra">    PetscScalar *cu,*cf;</div><div class="gmail_extra">    // need workspace for cf</div><div class="gmail_extra">    DMPlexVecGetClosure(dm,U,e,&csize,&cu); // I'm just going to pretend that he has abandoned the colliding name</div>
<div class="gmail_extra">    ComputeElementResidual(csize,cu,cf); // geometry and function space conveniently tucked under the rug</div><div class="gmail_extra">    DMPlexVecRestoreClosure(dm,U,e,&csize,&cu);</div>
<div class="gmail_extra">    DMPlexVecSetClosure(dm,F,e,csize,cf,ADD_VALUES);</div><div class="gmail_extra"><br></div><div class="gmail_extra"><br></div><div class="gmail_extra"><br></div><div class="gmail_extra">Note that there is no PetscSection above because its use is tucked inside the convenience routines which internally rely on DMComplexGetTransitiveClosure() and then figure out what part of the Vec to access using PetscSectionGetOffset() and PetscSectionGetDof().</div>