dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Apr 11, 2006 at 08:58:38AM +0200, Garth N. Wells wrote:
> 
> 
> 
> Quoting 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.
> > 
> 
> I saw that there is a lot of duplicate code in the mesh iterators ripe for some
> simplifications. I didn't look into it deeply because I didn't see the point if
> the whole mesh data structure will be changed shortly in favour of Sieve. I
> think that we should stick with a simple solution in terms of cleaning up
> FEM.cpp to get the boundary conditions sorted out, which I think the
> BoundaryIterator template offers (feel free to change the name), and hold off
> making any changes in the mesh files until we can get in idea of when Sieve will
> be ready. 
> 
> Garth

Agree. Until we know how Sieve looks and how we'll use it, we shouldn't
do any major work on the mesh classes.

I expect that we won't need to template over dimension with Sieve
since it will provide a dimension-independent interface, so then the
templates may go away, but it seems a very practical solution right
now to remove a lot of duplication.

/Anders


> 
> > /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/