[petsc-dev] PetscSection

[Karl Rupp]

On 11/10/2012 02:29 PM, Jed Brown wrote:
> On Sat, Nov 10, 2012 at 2:09 PM, Blaise A Bourdin <bourdin at lsu.edu
> <mailto:bourdin at lsu.edu>> wrote:
>
>     Nope... and please, don't break that...
>
>     Here is an example of a mesh with 4 element blocks. In an analysis
>     code, for each block , I would specify the element type, type of
>     physics, coefficients etc in a petscBag. That said, I may be doing
>     it wrong...
>
>     block 10 is made of tri
>     block 22 is made of quad
>     blocks 100 and 102 are made of bars
>
>
> I'm not suggesting any demotion of this ability. I'm trying to
> understand what is the advantage of having a stratum with unsorted mixed
> dimension. If the blocks were sorted by dimension, you would still
> access them using a label (which gives you an index set, perhaps
> contiguous). You would still be able to do all the normal operations
> with those points. If you want to have one loop over all elements of the
> various types, you would union or concatenate the index sets of each
> type. Am I missing something?

May I add something that is not specific to the current implementation 
in PETSc, yet still backs up the 'arrays of elements of the same type' 
approach:
  a) When thinking of accelerators, having, say, all tets 
stored/referenced contiguously in one array and all hexahedra stored 
contiguously in another (separate) array, makes it a lot easier to move 
all tets to a kernel working on tets, and to feed all hexes to a second 
kernel working on hexes.

  b) In a strongly typed setting (my experience here is in C++), one can 
do pretty nice things by keeping the two arrays separate rather than 
trying to unify them in some base class. For simplicity, suppose you 
have a mesh and two containers holding the tets and hexes. A FEM 
assembly could then look like this:

   for_each(mesh_cells.begin(), mesh_cells.end(), assembly_functor());

where assembly_functor statically dispatches into the cell type:
   struct assembly_functor{
     void operator()(tetrahedron const & t) { ... }
     void operator()(hexahedron const & t) { ... }
   }
and for_each is implemented such that it internally expands to a 
for_each over each of the two containers for tets and hexes. Thus, 
instead of dispatching for each element over and over again, the 
dispatch is done once for the whole array.

Just my 2 cents,
Karli