ffc team mailing list archive

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

2008/6/7 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('-')

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


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

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

--
Martin