dolfin team mailing list archive

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

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

> 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