dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Mon, Apr 10, 2006 at 09:04:15PM +0200, Garth N. Wells wrote:
> On Tue, 2006-04-04 at 12:16 -0500, Anders Logg wrote:
> > On Tue, Apr 04, 2006 at 09:08:20AM +0200, Johan Hoffman wrote:
> > 
> > > >> >  1. General boundary conditions including Neumann
> > >
> > > [...]
> > > >
> > > > It seems that we might want to start with (1) and then do 2-3-4 when
> > > > Sieve is released.
> > > >
> > > > I expect implementing (1) on both the FFC and DOLFIN sides would take
> > > > me a week to finish if I didn't have anything else on my hands. Who
> > > > wants to help?
> > > 
> > > I think this is a first priority as well. I can help, I have some time
> > > next week, and possibly some this week.
> > 
> > Very good!
> > 
> > Here is a quick sketch of what we need to do:
> > 
> > 1. Generate correct code for boundary integrals in FFC.
> > 
> > FFC generates the following two functions:
> > 
> >   void eval(real block[], const AffineMap& map) const;
> >   void eval(real block[], const AffineMap& map, unsigned int boundary) const;
> > 
> > I'm not sure "eval" is the best name, but "computeElementTensor" is a
> > bit long...
> > 
> > The second of the eval functions is for the boundary, so this is the
> > function we need to call during the iteration over the boundary.
> > 
> > The extra argument should tell you which piece of the boundary you are
> > currently integrating over (one of three edges for the triangle and
> > one of four faces for the tetrahedron).
> > 
> > 2. Iterate over the boundary in DOLFIN during assembly
> > 
> > While we fix this, we might want to think of a way to reduce the
> > amount of code in FEM.cpp which has grown pretty big. There is a lot
> > of repetition of code that we should try to remove. One of the
> > problems is that we have different versions for 2D and 3D where the
> > only difference is that in one case we iterate over edges and in the
> > other case over faces. With a modified mesh interface, we should be
> > able to do both cases (and also 1D) with the same code. Sieve has a
> > dimension-independent interface so we should look at that.
> > 
> 
> I've added a template class BoundaryIterator<class U, class V> to help
> with this. For 2D problems
> 
>   BoundaryIterator<EdgeIterator,Edge> iterator(boundary);
> 
> and for 3D problems
> 
>   BoundaryIterator<FaceIterator,Face> iterator(boundary);
> 
> I've tested it, but haven't checked in the changes to FEM.cpp yet. Will
> do soon. It will allow a big chunk of the code in FEM.cpp to be cut.
> 
> Garth

This is a step in the right direction, but I have a few concerns:

1. The other iterators for the mesh classes are named FooIterator
where Foo is the name of a single entity in the set of entities that
you iterate over: a CellIterator iterates over Cells.

We have an iterator called MeshIterator which iterates over the Meshes
in a MeshHierarchy, so a BoundaryIterator would indicate that you
iterate over Boundaries.

I don't know what to call edges and faces with a common name. The only
name I can think of is boundary entities, so we would have
BoundaryEntityIterator which is a bit long...

2. There is a lot of code duplication in all the mesh iterators:
VertexIterator, EdgeIterator, FaceIterator, CellIterator,
MeshIterator. All of this could probably be templated and then we can
typedef all the individual iterators for convenience.

I suggest we create a template MeshEntityIterator and use it to
replace all the different mesh iterators. If we did this,
then we could probably use it for iterating over edges of faces on the
boundary as well.

/Anders


> 
> > It would be great if you could start experimenting. I have added a new
> > test problem to FFC that we can use as a test (NeumannProblem.form):
> > 
> >     element = FiniteElement("Lagrange", "triangle", 1)
> > 
> >     v = TestFunction(element)
> >     U = TrialFunction(element)
> >     f = Function(element)
> >     g = Function(element)
> > 
> >     a = dot(grad(v), grad(U))*dx
> >     L = v*f*dx + v*g*ds
> > 
> > We could solve the following problem on the unit square:
> > 
> >    - div grad u(x, y) = f(x, y)
> > 
> >    f(x, y) = x*sin(y)
> > 
> >    u(x, y)     = 0         for x = 0 or y = 0
> >    du/dn(x, y) = sin(y)    for x = 1
> >    du/dn(x, y) = x*cos(1)  for y = 1
> > 
> > which has the exact solution u(x, y) = x*sin(y) (but please
> > double-check).
> > 
> > I expect I will need to fix the implementation on the FFC side, but
> > it would be great if I could get some help on the DOLFIN side.
> > 
> > /Anders
> > 
> > _______________________________________________
> > DOLFIN-dev mailing list
> > DOLFIN-dev@xxxxxxxxxx
> > http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
> 
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
> 

-- 
Anders Logg
Research Assistant Professor
Toyota Technological Institute at Chicago
http://www.tti-c.org/logg/