Subject: [FFC-dev] Constant() and quadrature representation

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 27 May 2008 11:29:16 +0200

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()?

What do You have in mind for UFL?

Kristian

---

From: Anders Logg <[email-masked]>
Date: Sat, 7 Jun 2008 00:25:50 +0200

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

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Sat, 7 Jun 2008 11:54:06 +0200

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

Also note that I have done nothing to support DG in UFL, which
features are needed for that?

-- 
Martin

---

From: Anders Logg <[email-masked]>
Date: Sat, 7 Jun 2008 12:21:03 +0200

> are unique. If dofmap->global_dimension() == value size, then this dofmap can be

Yes, we should add this. It's a fairly simple fix and we can hide it
in the UFC class.

> 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

---

From: Kristian Oelgaard <[email-masked]>
Date: Sat, 07 Jun 2008 12:28:28 +0200

> 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

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Sun, 8 Jun 2008 15:36:32 +0200

> facet. In FFC, one may do

Ok, that should be possible to add with few lines of code.

> start looking more closely at basis functions.

Ok. I'm happy except for a couple of FIXMEs in there, mostly
related to the value rank and shape of mixed elements.

--
Martin