ffc team mailing list archive

[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.

-- 
Anders

Attachment: signature.asc
Description: Digital signature