ufl team mailing list archive

[Kristian Oelgaard <k.b.oelgaard@xxxxxxxxxx>]

Quoting Anders Logg <logg@xxxxxxxxx>:

> On Sat, Jun 07, 2008 at 11:54:06AM +0200, Martin Sandve Alnæs wrote:
> > 2008/6/7 Anders Logg <logg@xxxxxxxxx>:
> > > On Tue, May 27, 2008 at 11:29:16AM +0200, Kristian Oelgaard wrote:
> > >>
> > >>
> > >> Hi,
> > >>
> > >> I have some problems with the built-in function Constant() when
> generating code
> > >> using the quadrature option. The problem is that it is restricted by
> default to
> > >> ('+/-'). Could someone remind me why this is the case?
> > >>
> > >> The reason this causes problems is that the map for
> restriction.CONSTANT
> > >> (ufcformat.py) is "0". So if one does:
> > >>
> > >> el = FiniteElement("Lagrange", "triangle", 1)
> > >> DG = FiniteElement("DG", "triangle", 0)
> > >>
> > >> v = TestFunction(el)
> > >> u = TrialFunction(el)
> > >> C = Constant("triangle")
> > >> Cu = C*u
> > >>
> > >> # 1
> > >> a = dot(grad(v), grad(Cu))*dx
> > >>
> > >> Transforms of the kind Jinv0_00 will appear in the header file. These
> transforms
> > >> are of course not defined for a cell integral which will result in a
> compile
> > >> error in DOLFIN. Tensor representation works because terms involving
> Jinv0_00
> > >> fall away.
> > >>
> > >> A second form that will fail (at FFC compile time) is:
> > >>
> > >> # 2
> > >> a = dot(jump(v), jump(Cu))*dS
> > >>
> > >> Because 'C' is already restricted.
> > >>
> > >> However
> > >> # 3
> > >> a = C*dot(jump(v), jump(u))*dS
> > >>
> > >> will work.
> > >>
> > >> If 'C' is defined as
> > >> C = Function(DG)
> > >> things work because it is no longer restricted and the map for None is
> "", so we
> > >> only have transforms of the kind Jinv_00.
> > >>
> > >> Now, why would someone take derivatives of a constant (like in #1) since
> the
> > >> following form is equivalent (and works)?
> > >>
> > >> # 4
> > >> a = C*dot(grad(v), grad(u))*dx
> > >>
> > >> Well, for very complex forms involving a lot of functions and constants
> it is
> > >> just easier to let FFC worry about which terms vanishes and which don't.
> An
> > >> optimisation of quadrature would naturally be to remove terms where the
> > >> basisfunctions are zero, this will take care of the problem as terms
> involving
> > >> Jinv0_00 disappears.
> > >>
> > >> But since form #2 also fails I think that it would make more sense not
> to
> > >> restrict Constant() by default and let the user handle restrictions.
> > >>
> > >> The way things are the user has to define forms like #3 and #4 or
> define
> > >> constants like:
> > >>
> > >> C = Function(DG)
> > >>
> > >> but then, why do we support Constant()?
> > >
> > > I don't know why we had the +/-. I removed it and all unit tests seem
> > > to run,
> > >
> > >> What do You have in mind for UFL?
> > >
> > > So far, the same as in FFC. Constant is just a shortcut for creating a
> > > DG0 element. I think Martin mentioned not long ago that the overhead
> > > for using DG0 elements as opposed to true constants is negligible. So
> > > we can use DG0 elements which simplifies both the UFC interface and
> > > the implementation in UFL.
> > >
> > 
> > I understand Constant(polygon) as a shortcut for
> > Function(FiniteElement("DG", polygon, 0)).
> > 
> > The overhead is neglible yet O(n) (a few operations per element).
> > 
> > A related optimization trick for the assembler is to enumerate which
> > of the dofmaps
> > are unique before iterating over the mesh, then call tabulate_dofs
> > only for those which
> > are unique. If dofmap->global_dimension() == value size, then this dofmap
> can be
> > omitted as well.
> 
> Yes, we should add this. It's a fairly simple fix and we can hide it
> in the UFC class.
> 
> > Also note that I have done nothing to support DG in UFL, which
> > features are needed for that?
> 
> We need two things:
> 
> 1. The ability to restrict things to either the + or the - side of a
> facet. In FFC, one may do
> 
>   v('+') - v('-')
> 
> 2. Operators that expand into restrictions:
> 
>   jump(v), avg(v)
> 
> or
> 
>   jump(v, n), avg(v, n)
> 
> There are a bunch of different jump operators defined differently for
> scalars and vectors, so it's important we get these right.
> 
> I can start looking at this. If we're both happy with the UFL finite
> element classes (still need to look at your latest fixes), I will
> start looking more closely at basis functions.

.. and we need to be able to integrate on interior facets (dS), as well as the
'macro element' to handle the dof_map related to restrictions.

Kristian

> -- 
> Anders
>