ffc team mailing list archive

["Martin Sandve Alnæs" <martinal@xxxxxxxxx>]

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

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