dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Apr 05, 2007 at 10:47:47AM +0200, Anders Logg wrote:
> On Thu, Apr 05, 2007 at 10:45:35AM +0200, Garth N. Wells wrote:
> > 
> > 
> > Anders Logg wrote:
> > > So, the new UFC-based assembly seems to be mostly complete and I think
> > > we've reached a point where we need to break everything in order to
> > > finish up.
> > > 
> > > This should be the last major rewrite we do for a while. We've done
> > > quite a few, but not needing to do another one soon is the whole point
> > > of the UFC interface.
> > > 
> > > Where should we do this? We can do it in dolfin-dev (dolfin), or we
> > > can set up a separate branch. Doing it in dolfin-dev would be the most
> > > practical, since we will need to touch the code in quite a few places
> > > to get things done.
> > > 
> > > Meanwhile, we have dolfin-stable for those who need something that
> > > works. (We should write something in the wiki about this!)
> > > 
> > > Here's a summary of what we have and what we need to do:
> > > 
> > > - The new UFC-based assembly has been implemented in the class
> > > Assembler and is mostly finished.
> > > 
> > > - There is a small test program in src/sandbox/assembly.
> > > 
> > > - We've been able to remove quite a lot of code. As mentioned before,
> > > Assembler.cpp is (still) considerably shorter than FEM.cpp.
> > > 
> > > - A whole lot of code can be removed from the Function classes. For
> > > example, here is what the implementation of UserFunction looks like:
> > > 
> > >   void UserFunction::interpolate(real* coefficients,
> > >                                  const ufc::cell& cell,
> > >                                  const ufc::finite_element& finite_element)
> > >   {
> > >     for (uint i = 0; i < finite_element.space_dimension(); i++)
> > >       coefficients[i] = finite_element.evaluate_dof(i, *this, cell);
> > >   }
> > > 
> > > And here is what the implementation of DiscreteFunction looks like:
> > > 
> > >   void NewDiscreteFunction::interpolate(real* coefficients,
> > >                                         const ufc::cell& cell,
> > >                                         const ufc::finite_element& finite_element)
> > >   {
> > >     dof_map->tabulate_dofs(dofs, cell);
> > >     x.get(coefficients, dof_map->local_dimension(), dofs);
> > >   }
> > > 
> > > - The class NewDiscreteFunction will replace DiscreteFunction and I
> > > would like us to try hard to avoid trying to make too much of it. In
> > > particular, it owns all its data (except for the mesh which is a
> > > reference) so we don't need to mess with smart pointers or tracking
> > > who owns the vector etc.
> > > 
> > 
> > Any problem adding a function to DiscreteFunction which returns the 
> > DofMap? For I/O, I need to get the mesh (which I can get from DofMap).
> > 
> > Garth
> 
> I might have an objection to this. I need to run now but will get
> back tonight and say more.
> 
> /Anders

For I/O, one would need access to the mesh, the vector, the dof map
and the finite element (to use finite_element::interpolate_vertex_values).

This would require putting in 4 special functions in the Function
interface that are only valid for one type of representation. It would
also require a lot of knowledge from the user (like the guy
implementing MFile.cpp or VTKFile.cpp) to combine these four objects
into something that gives the vertex values.

Instead, I'd like something like

    /// Interpolate function at vertices of mesh
    void interpolate(real* values, Mesh& mesh);

in Function. So the I/O would just rely on the Function class dumping
the values into an array.

There are some gaps here though, like where does the mesh come in. If
we want to do

    file << u;

then u must know about the mesh. For this to work, I think we should
require that the constructor of Function always must have a Mesh.

So the interpolate() function in DiscreteFunction would make use of
the vector and the dof map and the finite_element, while the
interpolate() function in UserFunction and ConstantFunction would
just evaluate at the vertices of the mesh.

And here's another gap: how do UserFunction and ConstantFunction know
the dimension of the value type (scalar or vector)?

/Anders