dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

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