Subject: [FFC-dev] General / Special quadrature rules in FFC

From: Peter Brune <[email-masked]>
Date: Mon, 6 Jul 2009 05:02:17 -0500

Some recent discussions and problems that have cropped up have led me to
need/want to add some specialized quadrature rules to FFC.  These include
"optimal" combinatorially-derived rules for given order polynomials for the
n-simplex, as well as rules that can be explicitly guaranteed to have
certain properties (positivity, all-internal-nodes, etc)

Does anyone else need anything like this?  I already have implementations of
a few nice rules for higher-dimensional applications, but have started to
adapt them for this purpose.  Ideally we would be able to specify a rule in
the ffc command.  I'm looking at where to put this.

There are other potential benefits to this move.  If we were to switch to
rules we have definite information about, we could use, say, group symmetry
properties, combined with symmetry properties of the unknowns, for further
optimization of the form compilation by quadrature.  We would already win,
especially in 3D, by getting away from the squashed rules and towards rules
that are more optimal.

- Peter

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 06 Jul 2009 12:25:04 +0200

> the ffc command.  I'm looking at where to put this.

So on the command line we should be able to do:

ffc --quadrature-rule some_string
or
ffc -q some_string

and then in the file ffc/fem/quadrature.py in the function
make_quadrature() there should be an extra argument

make_quadrature(shape, n, quad_rule)

with default None? Then you can just put your code and some switches in this
file.
I can add the command line option and propagate it to the make_quadrature()
function, at least for quadrature representation, if nobody objects.

> especially in 3D, by getting away from the squashed rules and towards rules

This would be nice.

Kristian

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 06 Jul 2009 13:03:55 +0200

> function, at least for quadrature representation, if nobody objects.

I added this to FFC, just add your rules to the make_quadrature() function.

Kristian

> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev
>

---

From: Robert Kirby <[email-masked]>
Date: Mon, 6 Jul 2009 08:13:14 -0500

Such rules may best be put into FIAT, and then propogated into ffc.

Rob

On Mon, Jul 6, 2009 at 6:03 AM, Kristian Oelgaard
<[email-masked]>wrote:

> > > Does anyone else need anything like this?  I already have
> implementations
> > of

---

From: Peter Brune <[email-masked]>
Date: Mon, 6 Jul 2009 08:28:43 -0500

It seems that this would just be another level through which one has to pass
an option.  I suppose that this would allow for EVERYTHING FENICS to use the
new rules immediately.

On that note, I have the Möller-Grundmann rules coded up in such a way that
FFC can use it, and have done some brief tests.  It seems to generate more
points than the squashed rule in 2D, but significantly less points in 3D (as
I could have guessed).  The FFC "number of points in each dimension" way of
specifying rules sort of doesn't apply here, either, so I basically have to
just go for the equivalent 2n + 1 order rule to match things up.  There's
also the issue of exploiting the symmetry.

- Peter

On Mon, Jul 6, 2009 at 8:13 AM, Robert Kirby <[email-masked]>wrote:

> On Mon, Jul 6, 2009 at 6:03 AM, Kristian Oelgaard <[email-masked]
> > wrote:
>

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 06 Jul 2009 15:34:29 +0200

> Such rules may best be put into FIAT, and then propogated into ffc.

That would be great, the infrastructure is in place on the FFC side, for
quadrature representation, it is just a matter of deciding on a naming
convention and then have FIAT accepting an additional string argument when
generating the rules.

I also know that Anders would like to move a lot of other things (like the
finite element classes) from the ffc/fem module to FIAT in order to shrink the
code base, but I'll let him comment on the details.

Kristian

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 06 Jul 2009 15:40:42 +0200

> an option.  I suppose that this would allow for EVERYTHING FENICS to use the

Yes, but I think it makes sense to let FIAT compute the quadrature scheme when
we rely on it for tabulating the basis at those points later on anyway.

> also the issue of exploiting the symmetry.

The 'number of points in each dimension' way should be changed, and I think that
by changing FIAT such that it takes a (shape, order, rule) set of arguments when
creating the rule will do the trick.

Kristian

---

From: Peter Brune <[email-masked]>
Date: Mon, 6 Jul 2009 08:45:05 -0500

This already seems to exist in FIAT through the make_quadrature_by_degree
function.  This could just be swapped out in fem/quadrature.py

- Peter

> Kristian
>
> > - Peter
> >
> > On Mon, Jul 6, 2009 at 8:13 AM, Robert Kirby <[email-masked]
> >wrote:
> >
> > >> > > Some recent discussions and problems that have cropped up have led
> me
> > >> to
> > >> > > "optimal" combinatorially-derived rules for given order
> polynomials
> > >> for the
> > >> > > a few nice rules for higher-dimensional applications, but have
> started
> > >> to
> > >> > > adapt them for this purpose.  Ideally we would be able to specify
> a
> > >> rule in
> > >> > with default None? Then you can just put your code and some switches
> in
> > >> this
> > >> > file.
> > >> > > There are other potential benefits to this move.  If we were to
> switch
> > >> to
> > >> > > optimization of the form compilation by quadrature.  We would
> already
> > >> win,
> > >> > > especially in 3D, by getting away from the squashed rules and
> towards
> > >> rules
> > >> > > that are more optimal.
> > >> >
> > >> > This would be nice.
> > >> >
> > >> > Kristian
> > >> >
> > >> > > - Peter
> > >> > >
> > >> >
> > >> >

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 06 Jul 2009 16:11:47 +0200

> function.  This could just be swapped out in fem/quadrature.py

Yes, it looks like it. However, we still need to put the 'rule' argument in
there. I'll wait for this before changing things.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Mon, 6 Jul 2009 16:35:45 +0200

> code base, but I'll let him comment on the details.

Almost. I would like to be able to use the FIAT FiniteElement class
directly from FFC without wrapping it. For that to work, some (small)
changes will be needed on the FIAT side (and some on the FFC side).

One thing that we currently handle in FFC which FIAT can't handle is
the MixedElement class. This needs to tabulate basis functions for a
list of elements with zeros inserted (taking the tensor product). Rob
has not wanted to put this into FIAT before, but maybe he can be
convinced... :-)

In general, it would be good if we could reduce the number of
FiniteElement classes in FEniCS. There are currently 7 (seven)!

  FIAT
  FFC
  UFL
  UFC
  DOLFIN
  SyFi
  (Exterior)
  
-- 
Anders

---

From: Anders Logg <[email-masked]>
Date: Mon, 6 Jul 2009 16:37:40 +0200

> has not wanted to put this into FIAT before, but maybe he can be

And while we're on the subject, should it be

  TH = P2 + P1 (what we have now)

or

  TH = P2 * P1

?

Someone mentioned it to me recently, I think it was Kristian.

-- 
Anders

---

From: Robert Kirby <[email-masked]>
Date: Tue, 7 Jul 2009 09:40:58 -0500

Hi, I'm about to be on vacation with the family for a while, but let me say
there are two issues here:1.) It is relatively easy to add rules into FIAT
(when one is not otherwise occupied) as a list of weights
and points.  There would be a small amount of work to add an option to use
them in ffc.  I originally
did the current rules because they were only used for creating integral
moments in FIAT.  I definitely see
the need for other rules, and would be glad to add them or have someone else
do so.

2.) Optimally exploiting quadrature rules is still a research topic.  When
ffc better supports matrix actions,
techniques like sum factorization on the smashed rules give big speedups.
 There may be ways of doing
things with group structure in symmetric rules.

Regarding mixed elements, I don't really want to be convinced to tabulate a
bunch of zeros.  The semantics of
the tensor product of bases can be achieved without blowing these things up
the memory usage and having to
strip out the zeros later.

On Mon, Jul 6, 2009 at 9:37 AM, Anders Logg <[email-masked]> wrote:

> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAkpSDDQACgkQTuwUCDsYZdGd/QCghMVQO77GvzGsVOl2LTUF2ZsQ
> BLQAn2IHvOem0ZU66wTl//sZT/ude4K2
> =nbji
> -----END PGP SIGNATURE-----
>
>

---

From: Peter Brune <[email-masked]>
Date: Tue, 7 Jul 2009 09:51:11 -0500

On Tue, Jul 7, 2009 at 9:40 AM, Robert Kirby <[email-masked]>wrote:

> the need for other rules, and would be glad to add them or have someone
> else do so.
>

I already put one into (my local copy) of FFC... I'll move it downstream and
send you the patch when appropriate.  Ideally we would have a whole forest
of potential rules with different properties.

> things with group structure in symmetric rules.

Also -- looking at this. :)  You should be able to look at some barycentric
representation of both the quadrature and basis and do something nice in
terms of eliminating permutation-equivalent things.  Python makes all this
EXCEEDINGLY easy to test!

- Peter