dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

A first version of the new UFC-based assembly is now implemented.
A bunch of features are still missing, but we can now assemble the
bilinear form for Poisson's equation and get correct results... :-)

Here are some notes if anyone wants to dig in. For example, it might
be a good idea to start thinking of merging in the parallel assembler
at this point.

 - The new assembler is implemented by the class Assembler in
   src/kernel/fem/.

 - A small test program is in src/sandbox/assembly/.

 - Assembler reuses data so one can save some time by doing

     Assembler assembler;
     assemble.assemble(A, form, mesh);
     assemble.assemble(A, form, mesh);

 - There is a simple function call assemble() for simple
   assembly without needing to create an Assembler:

     assemble(A, form, mesh);

 - What syntax do we want? Here are two options:

     assemble(A, form, mesh);
     assemble(form, A, mesh);

 - The new code is considerably shorter than before, currently 117
   lines in Assembler.cpp compared to 1177 in FEM.cpp. This is a very
   unfair comparison since many features are missing and there is a
   lot of extra testing code in FEM.cpp, but it feels much cleaner
   than the current assembler.

 - Local data structures like the local tensor, mapped degrees of
   freedom, local dimensions etc are stored in a data structure called
   UFC. We had something similar before with test_nodes and
   trial_nodes in the class BilinearForm. This way we can remove all
   the data allocation new/delete and initialization from Assembler.cpp.

 - There are two simple classes UFCMesh and UFCCell that provide a
   layer between DOLFIN::Mesh, DOLFIN::Cell and ufc::mesh, ufc::cell.
   The UFC data structures are not meant to be used for storing a mesh
   or a cell, but are just place holders/views for data that we need
   to send to the UFC functions. During assembly, we need to fill in
   the values on each cell and this can be done by UFCCell::update().
   A UFCCell is stored in the data structure UFC, so we need to do

       for (CellIterator cell...)
          ufc.cell.update(*cell);

   and we can then access the UFCCell by ufc.cell.

 - I have also added a new class GenericTensor which specifies the
   interface for a tensor of arbitrary rank. This is used to keep the
   number of functions and cases down. There is now only one
   assemble() function, namely assemble a form of arity r into a
   tensor of rank r. Previously, we had different classes for
   different forms (BilinearForm, LinearForm, Functional) and
   different arguments for the data structures (GenericMatrix,
   GenericVector, real). Now, it's just assembly of

       ufc::form --> GenericTensor

   I have added a simple sub class of GenericTensor called
   AssemblyMatrix. This is just something used for testing but it
   might be useful (it's in src/sandbox/assembly). It represents a
   sparse matrix as

       std::vector<std::map<uint, real> > A;

   Some benchmarks we did earlier indicated that this might be even
   faster than assembly into a PETSc or uBLAS matrix.

   The GenericTensor interface has only the functions needed to do
   assembly, currently

     /// Initialize zero tensor of given rank and dimensions
     virtual void init(uint rank, uint* dims) = 0;

     /// Return size of given dimension
     virtual uint size(const uint dim) const = 0;

     /// Add entries to tensor
     virtual void add(real* block, uint* size, uint** entries) = 0;

   Actually, size() is not used so it's currently only init() and
   add(). We might want to add some more functions as we refine the
   assembly.

   I'm thinking about letting GenericMatrix and GenericVector be sub
   classes of GenericTensor.

Discussion?

/Anders