Subject: [Fenics] systematic names option for elements based on FEEC

From: Douglas Arnold <[email-masked]>
Date: Wed, 29 Sep 2010 20:31:51 -0500

FEniCS/UFL uses historical element names, like Raviart-Thomas,
Nedelec 1st kind H(div), etc.  These have to be learned, although
this is fortunately ameliorated somewhat by having some
synonyms, like N1curl.  It can also cause confusion sometimes,
as the names are not always used consistently in the literature
(and I think that "Discontinuous Lagrange" is a contradiction
in terms).

However there is a systematic, consistent way to refer to
the majority of elements available, so I am making the
suggesting that this systematic form be added as
an additional option.  The systematic way is what is
used in the Finite Element Exterior Calculus, namely the
two families, which are in denoted in LaTeX by
$P_r\Lambda^k$ and $P_r^-\Lambda^k$, defined for
polynomial degree r=1,2,... and form degree k between
0 and n on simplices of dimension n.  Thus I am suggesting
adding something like

   V = FunctionSpace(mesh, "P Lambda", r, k)
and
   V = FunctionSpace(mesh, "P- Lambda", r, k)

The translation to current names is:

   {"P Lambda", r, 0} and {"P- Lambda", r, 0} coincide and
      are equal to {"Lagrange", r}
   {"P Lambda", r, n} is the same as {"Discontinuous Lagrange", r}
   {"P- Lambda", r, n} is the same as {"Discontinuous Lagrange", r-1}
   {"P- Lambda", r, 1} is {"Nedelec 1st kind H(curl)", r}
   {"P Lambda", r, 1} is {"Nedelec 2nd kind H(curl)", r}
   {"P- Lambda", r, n-1} is {"Nedelec 1st kind H(div)", r}
   {"P Lambda", r, n-1} is {"Nedelec 2nd kind H(div)", r}

I want to stress that there is nothing arbitrary about this.
The whole family of P Lambda and P- Lambda spaces can be
defined at once for all r, k, n.  This even works for n=1
and n>3 (although the latter are not yet implemented in
FEniCS).

There is a small issue in 2D where n-1 = 1, so the last
two cases conflict with the preceding two.  This is the
choice between the rotated Raviart-Thomas or Brezzi-Douglas-Marini
elements or the usual ones in 2D.  I would suggest always
thinking of 1-forms as the images of gradients, and hence
1-forms always correspond to H(curl), also in 2D.

What do you all think?

   -- Doug

P.S. On Marie's suggestion, I am posting this both to
the UFL and the FEniCS mailing lists.

---

From: Anders Logg <[email-masked]>
Date: Fri, 1 Oct 2010 15:54:14 +0200

>   {"P Lambda", r, n-1} is {"Nedelec 2nd kind H(div)", r}

Sounds good to me.

> and n>3 (although the latter are not yet implemented in

(They are, but only in the experimental Exterior package.)

> the UFL and the FEniCS mailing lists.

I have added a blueprint for this:

https://blueprints.launchpad.net/ufl/+spec/exterior-calculus-names/

And assigned it to the appropriate person ;-)

--
Anders

---

From: Marie Rognes <[email-masked]>
Date: Fri, 22 Oct 2010 15:12:47 +0200

> 1-forms always correspond to H(curl), also in 2D.

This has now been implemented in UFL using the notation

   V = FiniteElement("P Lambda", cell, r, k)
   V = FiniteElement("P- Lambda", cell, r, k)

where cell can be one of [interval, triangle, tetrahedron], r >= 1 [cf
footnote 1] and k in [0, n] where n is the topological dimension of the
cell.

While at it, I also added an "exterior_derivative" operator, mapping to
grad, curl, div, id based on the finite element space. An example of how
this can be used in DOLFIN is included below. (Hopefully, the original
requester will recognize the problem.)

------------------------------------------------------------------------------------------------------------------------
from dolfin import *
from ufl import exterior_derivative as d

meshes = [UnitInterval(2), UnitSquare(2, 2), UnitCube(2, 2, 2)]
families = ["P Lambda", "P- Lambda"]

r = 1
for (n, mesh) in enumerate(meshes):
  for family in families:
      for k in range(1, n+1):

          S = FunctionSpace(mesh, "P- Lambda", r+1, k-1)
          V = FunctionSpace(mesh, "P Lambda", r, k)
          W = S * V

          (s, u) = TrialFunctions(W)
          (t, v) = TestFunctions(W)

          a  = (- s, t) + (d(t), u) + (d(s), v) + (d(u), d(v))
          A = assemble(a)
          info(A)

------------------------------------------------------------------------------------------------------------------------

Footnotes:

[1] "P Lambda", "tetrahedron", r, 1 (aka Nedelec 2nd kind H(curl)
elements) is only supported by FIAT for r == 1.

--
Marie

---

From: Anders Logg <[email-masked]>
Date: Mon, 25 Oct 2010 09:52:09 +0200

> requester will recognize the problem.)

Very nice! :-)

--
Anders