dolfin team mailing list archive

["Johan Hoffman" <jhoffman@xxxxxxxxxx>]

> 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

Apart from the issue of duplicating code with or without templates, what
do we need to figure out to get the assembly of boundary integrals
correct?

Do we have a clear way to obtain the FE basis for the boundary of a cell
given the basis of the cell? Do we use the corresponding basis for for d-1
dimensions?

/Johan





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