dolfin team mailing list archive

["Garth N. Wells" <g.n.wells@xxxxxxxxxx>]

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


> /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/
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
>