ffc team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

This change breaks installation of FFC:

ffc/ext/time_elements_interface.cpp:23:10: fatal error:
'numpy/arrayobject.h' file not found
#include <numpy/arrayobject.h>

I assume that to commit doesn't check for the numpy path.

It's also a lot of C++ code for FFC, which is a Python library.

Garth


---------- Forwarded message ----------
From:  <noreply@xxxxxxxxxxxxx>
Date: 5 April 2013 03:24
Subject: [Branch ~ffc-core/ffc/main] Rev 1866: Merge work of Benjamin
Kehlet on including Radau and Lobatto
To: Garth Wells <gnw20@xxxxxxxxx>


Merge authors:
  Marie Rognes (meg-simula)
------------------------------------------------------------
revno: 1866 [merge]
committer: Marie E. Rognes <meg@xxxxxxxxx>
branch nick: trunk
timestamp: Thu 2013-04-04 21:21:30 +0200
message:
  Merge work of Benjamin Kehlet on including Radau and Lobatto
  elements. (These are elements defined on the interval via Lobatto and
  Radau quadrature points as degrees of freedom.)
added:
  ffc/ext/
  ffc/ext/Legendre.cpp
  ffc/ext/Legendre.h
  ffc/ext/LobattoQuadrature.cpp
  ffc/ext/LobattoQuadrature.h
  ffc/ext/RadauQuadrature.cpp
  ffc/ext/RadauQuadrature.h
  ffc/ext/time_elements.cpp
  ffc/ext/time_elements.h
  ffc/ext/time_elements_interface.cpp
  ffc/timeelements.py
  test/unit/elements/
  test/unit/elements/test.py
modified:
  ffc/fiatinterface.py
  setup.py


--
lp:ffc
https://code.launchpad.net/~ffc-core/ffc/main

Your team FFC Core Team is subscribed to branch lp:ffc.
To unsubscribe from this branch go to
https://code.launchpad.net/~ffc-core/ffc/main/+edit-subscription

=== added directory 'ffc/ext'
=== added file 'ffc/ext/Legendre.cpp'
--- ffc/ext/Legendre.cpp        1970-01-01 00:00:00 +0000
+++ ffc/ext/Legendre.cpp        2013-04-04 18:23:04 +0000
@@ -0,0 +1,117 @@
+// Copyright (C) 2003-2008 Anders Logg
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// Modified by Benjamin Kehlet 2011-2012
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+
+#include "Legendre.h"
+#include <cmath>
+
+//-----------------------------------------------------------------------------
+Legendre::Legendre(unsigned int n) : n(n), cache_x(0.0), cache(n + 1)
+{
+  cache[0] = 1.0; //constant value
+
+  // eval to initialize cache
+  eval(n, -1.0);
+}
+//-----------------------------------------------------------------------------
+double Legendre::operator() (double x)
+{
+  return eval(n, x);
+}
+//-----------------------------------------------------------------------------
+double Legendre::ddx(double x)
+{
+  return ddx(n, x);
+}
+//-----------------------------------------------------------------------------
+double Legendre::d2dx(double x)
+{
+  return d2dx(n, x);
+}
+//-----------------------------------------------------------------------------
+double Legendre::eval(unsigned int nn, double x)
+{
+  //recursive formula, BETA page 254
+  //return ( (2.0*nn-1.0)*x*eval(nn-1, x) - (nn-1.0)*eval(nn-2, x) ) / nn;
+
+  //The special cases
+  if (n == 0)
+    return 1.0;
+  else if (n == 1)
+    return x;
+
+  //check cache
+  if (x != cache_x)
+  {
+    cache[1] = x;
+    for (unsigned int i = 2; i <= n; ++i)
+    {
+      double ii = i;
+      cache[i] = ( (2.0*ii-1.0)*x*cache[i-1] - (ii-1.0)*cache[i-2] ) / ii;
+    }
+    cache_x = x;
+  }
+
+  return cache[nn];
+}
+//-----------------------------------------------------------------------------
+double Legendre::ddx(unsigned int n, double x)
+{
+  // Special cases
+  if (n == 0)
+    return 0.0;
+  else if (n == 1)
+    return 1.0;
+
+  // Avoid division by zero
+  if (std::abs(x - 1.0) < EPSILON)
+    x -= 2.0*EPSILON;
+
+  if (std::abs(x + 1.0) < EPSILON)
+    x += 2.0*EPSILON;
+
+  // Formula, BETA page 254
+  const double nn = n;
+  return nn * (x*eval(n, x) - eval(n-1, x)) / (x*x - 1.0);
+}
+//-----------------------------------------------------------------------------
+double Legendre::d2dx(unsigned int, double x)
+{
+  // Special case n = 0
+  if (n == 0)
+    return 0.0;
+
+  // Special case n = 1
+  if (n == 1)
+    return 0.0;
+
+  // Avoid division by zero
+  if (std::abs(x - 1.0) < EPSILON)
+    x -= 2.0*EPSILON;
+  if (std::abs(x + 1.0) < EPSILON)
+    x += 2.0*EPSILON;
+
+  // Formula, BETA page 254
+  const double nn = double(n);
+  return (2.0*x*ddx(n, x) - nn*(nn+1)*eval(n, x)) / (1.0-x*x);
+}
+//-----------------------------------------------------------------------------

=== added file 'ffc/ext/Legendre.h'
--- ffc/ext/Legendre.h  1970-01-01 00:00:00 +0000
+++ ffc/ext/Legendre.h  2013-04-04 18:23:04 +0000
@@ -0,0 +1,65 @@
+// Copyright (C) 2003-2009 Anders Logg
+//
+// This file is part of FFC.
+//
+// DOLFIN is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// DOLFIN is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with DOLFIN.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+#ifndef __LEGENDRE_H
+#define __LEGENDRE_H
+
+/// Legendre polynomial of given degree n on the interval [-1,1].
+///
+///   P0(x) = 1
+///   P1(x) = x
+///   P2(x) = (3x^2 - 1) / 2
+///   ...
+///
+/// The function values and derivatives are computed using
+/// three-term recurrence formulas.
+
+#include <vector>
+#define EPSILON 10e-15
+
+class Legendre
+{
+ public:
+
+  Legendre(unsigned int n);
+
+  /// Evaluation at given point
+  double operator() (double x);
+
+  /// Evaluation of derivative at given point
+  double ddx(double x);
+
+  /// Evaluation of second derivative at given point
+  double d2dx(double x);
+
+  /// Evaluation of arbitrary order, nn <= n (useful ie in RadauQuadrature)
+  double eval(unsigned int nn, double x);
+
+  double ddx(unsigned int n, double x);
+  double d2dx(unsigned int n, double x);
+
+ private:
+
+  const unsigned int n;
+  double cache_x;
+  std::vector<double> cache;
+
+};
+#endif

=== added file 'ffc/ext/LobattoQuadrature.cpp'
--- ffc/ext/LobattoQuadrature.cpp       1970-01-01 00:00:00 +0000
+++ ffc/ext/LobattoQuadrature.cpp       2013-04-04 18:23:04 +0000
@@ -0,0 +1,93 @@
+// Copyright (C) 2003-2006 Anders Logg
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+#include "LobattoQuadrature.h"
+#include "Legendre.h"
+
+#include <cmath>
+
+//-----------------------------------------------------------------------------
+LobattoQuadrature::LobattoQuadrature(unsigned int n)
+  : points(n)
+{
+  // FIXME: Do proper arguement checking
+  // if (n < 2)
+  //   error("Lobatto quadrature requires at least 2 points.");
+
+  // init();
+
+  // if (!check(2*n - 3))
+  //   error("Lobatto quadrature not ok, check failed.");
+
+
+
+  // Compute the Lobatto quadrature points in [-1,1] as the endpoints
+  // and the zeroes of the derivatives of the Legendre polynomials
+  // using Newton's method
+
+  //const unsigned int n = points.size();
+
+  // Special case n = 1 (should not be used)
+  if (n == 1)
+  {
+    points[0] = 0.0;
+    return;
+  }
+
+  // Special case n = 2
+  if (n == 2)
+  {
+    points[0] = -1.0;
+    points[1] = 1.0;
+    return;
+  }
+
+  Legendre p(n - 1);
+  double x, dx;
+
+  // Set the first and last nodal points which are 0 and 1
+  points[0] = -1.0;
+  points[n - 1] = 1.0;
+
+  // Compute the rest of the nodes by Newton's method
+  for (unsigned int i = 1; i <= ((n-1)/2); i++)
+  {
+
+    // Initial guess
+    x = cos(3.1415926*double(i)/double(n - 1));
+
+    // Newton's method
+    do
+    {
+      dx = -p.ddx(x)/p.d2dx(x);
+      x  = x + dx;
+    } while (std::abs(dx) > EPSILON);
+
+    // Save the value using the symmetry of the points
+    points[i] = -x;
+    points[n - 1 - i] = x;
+
+  }
+
+  // Fix the middle node
+  if ((n % 2) != 0)
+    points[n/2] = 0.0;
+}
+

=== added file 'ffc/ext/LobattoQuadrature.h'
--- ffc/ext/LobattoQuadrature.h 1970-01-01 00:00:00 +0000
+++ ffc/ext/LobattoQuadrature.h 2013-04-04 18:23:04 +0000
@@ -0,0 +1,47 @@
+// Copyright (C) 2003-2009 Anders Logg
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+#ifndef __LOBATTO_QUADRATURE_H
+#define __LOBATTO_QUADRATURE_H
+
+/// Lobatto (Gauss-Lobatto) quadrature on the interval [-1,1].
+/// The n quadrature points are given by the end-points -1 and 1,
+/// and the zeros of P{n-1}'(x), where P{n-1}(x) is the (n-1):th
+/// Legendre polynomial.
+///
+/// The quadrature points are computed using Newton's method, and
+/// the quadrature weights are computed by solving a linear system
+/// determined by the condition that Lobatto quadrature with n points
+/// should be exact for polynomials of degree 2n-3.
+
+#include <vector>
+
+class LobattoQuadrature
+{
+ public:
+
+  /// Create Lobatto quadrature with n points
+  LobattoQuadrature(unsigned int n);
+  //~LobattoQuadrature();
+
+  std::vector<double> points;
+};
+
+#endif

=== added file 'ffc/ext/RadauQuadrature.cpp'
--- ffc/ext/RadauQuadrature.cpp 1970-01-01 00:00:00 +0000
+++ ffc/ext/RadauQuadrature.cpp 2013-04-04 18:23:04 +0000
@@ -0,0 +1,84 @@
+// Copyright (C) 2003-2006 Anders Logg
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+#include "RadauQuadrature.h"
+#include "Legendre.h"
+
+#include <cmath>
+
+//-----------------------------------------------------------------------------
+RadauQuadrature::RadauQuadrature(unsigned int n)
+  : points(n+1)
+{
+  // Compute the Radau quadrature points in [-1,1] as -1 and the zeros
+  // of ( Pn-1(x) + Pn(x) ) / (1+x) where Pn is the n:th Legendre
+  // polynomial. Computation is a little different than for Gauss and
+  // Lobatto quadrature, since we don't know of any good initial
+  // approximation for the Newton iterations.
+
+  // Special case n = 1
+  if (n == 1)
+  {
+    points[0] = -1.0;
+    return;
+  }
+
+  Legendre p(n);
+  double x, dx, step, sign;
+
+  // Set size of stepping for seeking starting points
+  step = 1.0/(double(n - 1)*15.0);
+
+  // Set the first nodal point which is -1
+  points[0] = -1.0;
+
+  // Start at -1 + step
+  x = -1.0 + step;
+
+  // Set the sign at -1 + epsilon
+  sign = ((p.eval(n - 1, x) + p(x)) > 0 ? 1.0 : -1.0);
+
+  // Compute the rest of the nodes by Newton's method
+  for (unsigned int i = 1; i < n; i++)
+  {
+
+    // Step to a sign change
+    while ((p.eval(n - 1, x) + p(x))*sign > 0.0)
+      x += step;
+
+    // Newton's method
+    do
+    {
+      dx = -(p.eval(n-1, x) + p(x))/(p.ddx(n - 1, x) + p.ddx(x));
+      x  = x + dx;
+    } while (std::abs(dx) > EPSILON);
+
+    // Set the node value
+    points[i] = x;
+
+    // Fix step so that it's not too large
+    if (step > (points[i] - points[i-1])/10.0)
+      step = (points[i] - points[i-1])/10.0;
+
+    // Step forward
+    sign = -sign;
+    x += step;
+  }
+}

=== added file 'ffc/ext/RadauQuadrature.h'
--- ffc/ext/RadauQuadrature.h   1970-01-01 00:00:00 +0000
+++ ffc/ext/RadauQuadrature.h   2013-04-04 18:23:04 +0000
@@ -0,0 +1,49 @@
+// Copyright (C) 2003-2009 Anders Logg
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+#ifndef __RADAU_QUADRATURE_H
+#define __RADAU_QUADRATURE_H
+
+/// Radau (Gauss-Radau) quadrature on the interval [-1,1].
+/// The n quadrature points are given by the zeros of
+///
+///     ( Pn-1(x) + Pn(x) ) / (1+x)
+///
+/// where Pn is the n:th Legendre polynomial.
+///
+/// The quadrature points are computed using Newton's method, and
+/// the quadrature weights are computed by solving a linear system
+/// determined by the condition that Radau quadrature with n points
+/// should be exact for polynomials of degree 2n-2.
+
+#include <vector>
+
+class RadauQuadrature
+{
+ public:
+
+  /// Create Radau quadrature with n points
+  RadauQuadrature(unsigned int n);
+
+  std::vector<double> points;
+
+};
+
+#endif

=== added file 'ffc/ext/time_elements.cpp'
--- ffc/ext/time_elements.cpp   1970-01-01 00:00:00 +0000
+++ ffc/ext/time_elements.cpp   2013-04-04 18:23:04 +0000
@@ -0,0 +1,52 @@
+// Copyright (C) 2012 Benjamin Kehlet
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+#include <iostream>
+#include "LobattoQuadrature.h"
+#include "RadauQuadrature.h"
+#include "Legendre.h"
+
+void compute_lobatto_points(double* points, const unsigned int degree)
+{
+  // Compute the nodal basis
+  LobattoQuadrature lobatto(degree + 1);
+  for (unsigned int i = 0; i < degree +1; i++)
+    points[i] = (lobatto.points[i] + 1.0) / 2.0;
+}
+
+void compute_radau_points  (double* points, const unsigned int degree)
+{
+  RadauQuadrature radau(degree+1);
+
+  for (unsigned int i = 0; i < degree+1; i++)
+    points[degree-i] = (-radau.points[i] + 1.0) / 2.0;
+}
+
+void compute_legendre_coeffs(double* coeffs, const double *points,
const unsigned int num_points, const unsigned int degree)
+{
+  for (unsigned int i = 0; i < degree; i++)
+  {
+    Legendre p(i);
+    for (unsigned int j = 0; j < num_points; j++)
+    {
+      coeffs[i*num_points + j] = p(points[j]*2.0 -1.0);
+    }
+  }
+}

=== added file 'ffc/ext/time_elements.h'
--- ffc/ext/time_elements.h     1970-01-01 00:00:00 +0000
+++ ffc/ext/time_elements.h     2013-04-04 18:23:04 +0000
@@ -0,0 +1,29 @@
+// Copyright (C) 2012 Benjamin Kehlet
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+
+void compute_lobatto_points(double* points, const unsigned int degree);
+
+void compute_radau_points  (double* points, const unsigned int degree);
+
+void compute_legendre_coeffs(double* coeffs,
+                            const double *points,
+                            const unsigned int num_points,
+                            const unsigned int degree);

=== added file 'ffc/ext/time_elements_interface.cpp'
--- ffc/ext/time_elements_interface.cpp 1970-01-01 00:00:00 +0000
+++ ffc/ext/time_elements_interface.cpp 2013-04-04 18:23:04 +0000
@@ -0,0 +1,124 @@
+// Copyright (C) 2012 Benjamin Kehlet
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+
+#include <Python.h>
+#include <numpy/arrayobject.h>
+
+#include "time_elements.h"
+
+static PyObject *compute_lobatto_interface(PyObject *dummy, PyObject *args)
+{
+  int degree;
+
+  /* parse argument tuple */
+  if (!PyArg_ParseTuple(args, "i", &degree))
+  {
+    return NULL; /* PyArg_ParseTuple has raised an exception */
+  }
+
+  npy_intp n = degree+1;
+
+  PyArrayObject *py_array_points =  (PyArrayObject*)
PyArray_SimpleNew(1, &n, NPY_DOUBLE);
+
+  double *points = (double*)  PyArray_DATA(py_array_points);
+
+  compute_lobatto_points(points, degree);
+
+  // return values
+  return (PyObject*) py_array_points;
+}
+
+static PyObject *compute_radau_interface(PyObject *dummy, PyObject *args)
+{
+  int degree;
+
+  /* parse argument tuple */
+  if (!PyArg_ParseTuple(args, "i", &degree))
+  {
+    return NULL; /* PyArg_ParseTuple has raised an exception */
+  }
+
+  npy_intp n = degree+1;
+
+  PyArrayObject *py_array_points =  (PyArrayObject*)
PyArray_SimpleNew(1, &n, NPY_DOUBLE);
+
+  double *points = (double*)  PyArray_DATA(py_array_points);
+
+  compute_radau_points(points, degree);
+
+  // return values
+  return (PyObject*) py_array_points;
+}
+
+static PyObject *compute_legendre_coeffs_interface(PyObject *dummy,
PyObject *args)
+{
+  PyArrayObject *points_array;
+
+  /* parse argument tuple */
+  if (!PyArg_ParseTuple(args, "O!",
+                       &PyArray_Type, &points_array))
+  {
+    return NULL; /* PyArg_ParseTuple has raised an exception */
+  }
+
+  const npy_intp num_points = PyArray_DIMS(points_array)[0];
+  npy_intp dims[2] = { num_points, num_points };
+  PyArrayObject *py_array_coeffs =  (PyArrayObject*)
PyArray_SimpleNew(2, dims, NPY_DOUBLE);
+
+  double *coeffs = (double*)  PyArray_DATA(py_array_coeffs);
+
+  compute_legendre_coeffs(coeffs, (double*)
PyArray_DATA(points_array), num_points, num_points);
+
+  // return values
+  return (PyObject*) py_array_coeffs;
+}
+
+static char compute_lobatto_doc[] = \
+  "Doc string for compute_lobatto_points";
+static char compute_radau_doc[] = \
+  "Doc string for compute_radau_points";
+static char compute_legendre_coeffs_doc[] = \
+  "Doc string for compute_legendre_coeffs";
+
+
+static PyMethodDef time_elements_ext_methods[] = {
+  {"compute_lobatto_points",
+   compute_lobatto_interface,
+   METH_VARARGS,
+   compute_lobatto_doc},
+  {"compute_radau_points",
+   compute_radau_interface,
+   METH_VARARGS,
+   compute_radau_doc},
+  {"compute_legendre_coeffs",
+   compute_legendre_coeffs_interface,
+   METH_VARARGS,
+   compute_legendre_coeffs_doc},
+
+  {NULL, NULL}
+};
+
+PyMODINIT_FUNC
+inittime_elements_ext(void)
+{
+   (void)Py_InitModule("time_elements_ext", time_elements_ext_methods);
+   import_array();
+}

=== modified file 'ffc/fiatinterface.py'
--- ffc/fiatinterface.py        2013-03-13 10:58:55 +0000
+++ ffc/fiatinterface.py        2013-04-04 19:07:05 +0000
@@ -16,7 +16,7 @@
 # along with FFC. If not, see <http://www.gnu.org/licenses/>.
 #
 # Modified by Garth N. Wells, 2009.
-# Modified by Marie Rognes, 2009-2010.
+# Modified by Marie Rognes, 2009-2013.
 # Modified by Martin Alnaes, 2013
 #
 # First added:  2009-03-06
@@ -32,6 +32,8 @@
 # FFC modules
 from ffc.log import debug, error, ffc_assert
 from ffc.quadratureelement import QuadratureElement as FFCQuadratureElement
+from ffc.timeelements import LobattoElement as FFCLobattoElement
+from ffc.timeelements import RadauElement as FFCRadauElement

 from ffc.mixedelement import MixedElement
 from ffc.restrictedelement import RestrictedElement
@@ -58,8 +60,10 @@
                       "Crouzeix-Raviart",
                       "Discontinuous Lagrange",
                       "Lagrange",
+                      "Lobatto",
                       "Nedelec 1st kind H(curl)",
                       "Nedelec 2nd kind H(curl)",
+                      "Radau",
                       "Raviart-Thomas",
                       "Real",
                       "Bubble",
@@ -139,6 +143,12 @@
         constant = _create_fiat_element(dg0_element)
         return SpaceOfReals(constant)

+    # Handle the specialized time elements
+    if family == "Lobatto" :
+        return FFCLobattoElement(ufl_element.degree())
+    if family == "Radau" :
+        return FFCRadauElement(ufl_element.degree())
+
     # FIXME: AL: Should this really be here?
     # Handle QuadratureElement
     if family == "Quadrature":

=== added file 'ffc/timeelements.py'
--- ffc/timeelements.py 1970-01-01 00:00:00 +0000
+++ ffc/timeelements.py 2013-04-04 18:23:04 +0000
@@ -0,0 +1,128 @@
+# Copyright (C) 2012 Benjamin Kehlet
+#
+# This file is part of FFC.
+#
+# FFC is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# FFC is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with FFC. If not, see <http://www.gnu.org/licenses/>.
+#
+# Modified by Marie E. Rognes, 2012
+#
+# First added:  2012-08-15
+# Last changed: 2012-09-07
+
+
+from FIAT import finite_element, polynomial_set, dual_set,
functional, reference_element
+import time_elements_ext as ext
+import numpy
+
+class TimeElementDualSet(dual_set.DualSet):
+  """. """
+  def __init__(self, family, degree):
+    assert(family == "Lobatto" or family == "Radau"), \
+        "Unknown time element '%s'" % family
+    if family == "Lobatto" :
+        assert (degree > 0), "Lobatto not defined for degree < 1!"
+    else :
+        assert(degree >= 0), "Degree must be >= 0"
+
+    # ids is a map from mesh entity (numbers) to dof numbers
+    ids = {}
+
+    # dofs is a list of functionals
+    dofs = []
+
+    # Only defined in 1D (on an inteval)
+    cell = reference_element.UFCInterval()
+
+    self.coords = (ext.compute_lobatto_points(degree) if family == "Lobatto"
+                   else ext.compute_radau_points(degree))
+    points = [(c,) for c in self.coords]
+
+    # Create dofs from points
+    dofs = [functional.PointEvaluation(cell, point)
+            for point in points]
+
+    # Create ids
+    if family == "Lobatto":
+      ids[0] = {0: [0], 1: [len(points)-1]}
+      ids[1] = {0: range(1, len(points)-1)}
+    elif family == "Radau":
+      ids[0] = {0: [], 1: []}
+      ids[1] = {0: range(len(points))} # Treat all Radau points as internal
+    else:
+      error("Undefined family: %s" % family)
+
+    # Initialize dual set
+    dual_set.DualSet.__init__(self, dofs, cell, ids)
+
+class TimeElement(finite_element.FiniteElement):
+  """."""
+  def __init__(self, family, degree):
+    "Create time element with given (polynomial degree)."
+
+    # Only defined in 1D (on an inteval)
+    cell = reference_element.UFCInterval()
+
+    # Initialize polynomial space of degree 'degree'
+    polynomial_space = polynomial_set.ONPolynomialSet(cell, degree)
+
+    # Create dual (degrees of freedom)
+    dual = TimeElementDualSet(family, degree)
+
+    # Initialize super class
+    finite_element.FiniteElement.__init__(self,
+                                          polynomial_space,
+                                          dual,
+                                          degree
+                                          )
+
+  def compute_quadrature_weights(self) :
+    """Compute the quadrature weights by solving a linear system of equations
+    for exact integration of polynomials. We compute the integrals over
+    [-1,1] of the Legendre polynomials of degree <= n - 1; These integrals
+    are all zero, except for the integral of P0 which is 2.
+
+    This requires that the n-point quadrature rule is exact at least for
+    polynomials of degree n-1."""
+
+    n = len(self.dual.coords)
+
+    # Special case n = 0
+    if n == 0 :
+      weights[0] = 2.0;
+      return weights
+
+    # Initialize linear system
+    A = ext.compute_legendre_coeffs(self.dual.coords)
+
+    b = numpy.zeros(n)
+    b[0] = 2.0
+
+    weights = numpy.linalg.solve(A, b)
+
+    # Weights are computed on interval [-1, 1]. Scale to reference interval
+    return weights/2.0
+
+
+class LobattoElement(TimeElement):
+  """."""
+  def __init__(self, degree):
+      "Create Lobatto element with given (polynomial degree)."
+      TimeElement.__init__(self, "Lobatto", degree)
+
+class RadauElement(TimeElement):
+  """."""
+  def __init__(self, degree):
+      "Create Radau element with given (polynomial degree)."
+      TimeElement.__init__(self, "Radau", degree)
+

=== modified file 'setup.py'
--- setup.py    2013-03-25 09:20:08 +0000
+++ setup.py    2013-04-04 18:23:04 +0000
@@ -1,7 +1,7 @@
 #!/usr/bin/env python

 import sys, platform
-from distutils.core import setup
+from distutils.core import setup, Extension
 from os import chdir
 from os.path import join, split

@@ -20,6 +20,13 @@
         batch_files.append(batch_file)
     scripts.extend(batch_files)

+ext = Extension("ffc.time_elements_ext",
+                ["ffc/ext/time_elements_interface.cpp",
+                 "ffc/ext/time_elements.cpp",
+                 "ffc/ext/LobattoQuadrature.cpp",
+                 "ffc/ext/RadauQuadrature.cpp",
+                 "ffc/ext/Legendre.cpp"])
+
 setup(name = "FFC",
       version = "1.2.0",
       description = "The FEniCS Form Compiler",
@@ -32,5 +39,6 @@
                   "ffc.dolfin"],
       package_dir={"ffc": "ffc"},
       scripts = scripts,
+      ext_modules = [ext],
       data_files = [(join("share", "man", "man1"),
                      [join("doc", "man", "man1", "ffc.1.gz")])])

=== added directory 'test/unit/elements'
=== added file 'test/unit/elements/test.py'
--- test/unit/elements/test.py  1970-01-01 00:00:00 +0000
+++ test/unit/elements/test.py  2013-04-04 19:07:05 +0000
@@ -0,0 +1,43 @@
+"""Unit tests for FFC finite elements"""
+
+# Copyright (C) 2013 Marie E. Rognes
+#
+# This file is part of FFC.
+#
+# FFC is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# FFC is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with FFC. If not, see <http://www.gnu.org/licenses/>.
+
+import unittest
+
+from ufl import interval
+from ufl import FiniteElement
+
+from ffc import compile_element
+
+class TestCompileElements(unittest.TestCase):
+
+    def testRadau(self):
+        "Test that Radau elements compile."
+        for degree in range(3):
+            element = FiniteElement("Radau", interval, degree)
+            compile_element(element)
+
+    def testLobatto(self):
+        "Test that Lobatto elements compile."
+        for degree in range(1, 4):
+            element = FiniteElement("Lobatto", interval, degree)
+            compile_element(element)
+
+
+if __name__ == "__main__":
+    unittest.main()

=== added directory 'ffc/ext'
=== added file 'ffc/ext/Legendre.cpp'
--- ffc/ext/Legendre.cpp	1970-01-01 00:00:00 +0000
+++ ffc/ext/Legendre.cpp	2013-04-04 18:23:04 +0000
@@ -0,0 +1,117 @@
+// Copyright (C) 2003-2008 Anders Logg
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// Modified by Benjamin Kehlet 2011-2012
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+
+#include "Legendre.h"
+#include <cmath>
+
+//-----------------------------------------------------------------------------
+Legendre::Legendre(unsigned int n) : n(n), cache_x(0.0), cache(n + 1)
+{
+  cache[0] = 1.0; //constant value
+
+  // eval to initialize cache
+  eval(n, -1.0);
+}
+//-----------------------------------------------------------------------------
+double Legendre::operator() (double x)
+{
+  return eval(n, x);
+}
+//-----------------------------------------------------------------------------
+double Legendre::ddx(double x)
+{
+  return ddx(n, x);
+}
+//-----------------------------------------------------------------------------
+double Legendre::d2dx(double x)
+{
+  return d2dx(n, x);
+}
+//-----------------------------------------------------------------------------
+double Legendre::eval(unsigned int nn, double x)
+{
+  //recursive formula, BETA page 254
+  //return ( (2.0*nn-1.0)*x*eval(nn-1, x) - (nn-1.0)*eval(nn-2, x) ) / nn;
+
+  //The special cases
+  if (n == 0)
+    return 1.0;
+  else if (n == 1)
+    return x;
+
+  //check cache
+  if (x != cache_x)
+  {
+    cache[1] = x;
+    for (unsigned int i = 2; i <= n; ++i)
+    {
+      double ii = i;
+      cache[i] = ( (2.0*ii-1.0)*x*cache[i-1] - (ii-1.0)*cache[i-2] ) / ii;
+    }
+    cache_x = x;
+  }
+
+  return cache[nn];
+}
+//-----------------------------------------------------------------------------
+double Legendre::ddx(unsigned int n, double x)
+{
+  // Special cases
+  if (n == 0)
+    return 0.0;
+  else if (n == 1)
+    return 1.0;
+
+  // Avoid division by zero
+  if (std::abs(x - 1.0) < EPSILON)
+    x -= 2.0*EPSILON;
+
+  if (std::abs(x + 1.0) < EPSILON)
+    x += 2.0*EPSILON;
+
+  // Formula, BETA page 254
+  const double nn = n;
+  return nn * (x*eval(n, x) - eval(n-1, x)) / (x*x - 1.0);
+}
+//-----------------------------------------------------------------------------
+double Legendre::d2dx(unsigned int, double x)
+{
+  // Special case n = 0
+  if (n == 0)
+    return 0.0;
+
+  // Special case n = 1
+  if (n == 1)
+    return 0.0;
+
+  // Avoid division by zero
+  if (std::abs(x - 1.0) < EPSILON)
+    x -= 2.0*EPSILON;
+  if (std::abs(x + 1.0) < EPSILON)
+    x += 2.0*EPSILON;
+
+  // Formula, BETA page 254
+  const double nn = double(n);
+  return (2.0*x*ddx(n, x) - nn*(nn+1)*eval(n, x)) / (1.0-x*x);
+}
+//-----------------------------------------------------------------------------

=== added file 'ffc/ext/Legendre.h'
--- ffc/ext/Legendre.h	1970-01-01 00:00:00 +0000
+++ ffc/ext/Legendre.h	2013-04-04 18:23:04 +0000
@@ -0,0 +1,65 @@
+// Copyright (C) 2003-2009 Anders Logg
+//
+// This file is part of FFC.
+//
+// DOLFIN is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// DOLFIN is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with DOLFIN.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+#ifndef __LEGENDRE_H
+#define __LEGENDRE_H
+
+/// Legendre polynomial of given degree n on the interval [-1,1].
+///
+///   P0(x) = 1
+///   P1(x) = x
+///   P2(x) = (3x^2 - 1) / 2
+///   ...
+///
+/// The function values and derivatives are computed using
+/// three-term recurrence formulas.
+
+#include <vector>
+#define EPSILON 10e-15
+
+class Legendre
+{
+ public:
+
+  Legendre(unsigned int n);
+
+  /// Evaluation at given point
+  double operator() (double x);
+
+  /// Evaluation of derivative at given point
+  double ddx(double x);
+
+  /// Evaluation of second derivative at given point
+  double d2dx(double x);
+
+  /// Evaluation of arbitrary order, nn <= n (useful ie in RadauQuadrature)
+  double eval(unsigned int nn, double x);
+
+  double ddx(unsigned int n, double x);
+  double d2dx(unsigned int n, double x);
+ 
+ private:
+
+  const unsigned int n;
+  double cache_x;
+  std::vector<double> cache;
+
+};
+#endif

=== added file 'ffc/ext/LobattoQuadrature.cpp'
--- ffc/ext/LobattoQuadrature.cpp	1970-01-01 00:00:00 +0000
+++ ffc/ext/LobattoQuadrature.cpp	2013-04-04 18:23:04 +0000
@@ -0,0 +1,93 @@
+// Copyright (C) 2003-2006 Anders Logg
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+#include "LobattoQuadrature.h"
+#include "Legendre.h"
+
+#include <cmath>
+
+//-----------------------------------------------------------------------------
+LobattoQuadrature::LobattoQuadrature(unsigned int n)
+  : points(n)
+{
+  // FIXME: Do proper arguement checking
+  // if (n < 2)
+  //   error("Lobatto quadrature requires at least 2 points.");
+
+  // init();
+
+  // if (!check(2*n - 3))
+  //   error("Lobatto quadrature not ok, check failed.");
+
+
+
+  // Compute the Lobatto quadrature points in [-1,1] as the endpoints
+  // and the zeroes of the derivatives of the Legendre polynomials
+  // using Newton's method
+
+  //const unsigned int n = points.size();
+
+  // Special case n = 1 (should not be used)
+  if (n == 1)
+  {
+    points[0] = 0.0;
+    return;
+  }
+
+  // Special case n = 2
+  if (n == 2)
+  {
+    points[0] = -1.0;
+    points[1] = 1.0;
+    return;
+  }
+
+  Legendre p(n - 1);
+  double x, dx;
+
+  // Set the first and last nodal points which are 0 and 1
+  points[0] = -1.0;
+  points[n - 1] = 1.0;
+
+  // Compute the rest of the nodes by Newton's method
+  for (unsigned int i = 1; i <= ((n-1)/2); i++)
+  {
+
+    // Initial guess
+    x = cos(3.1415926*double(i)/double(n - 1));
+
+    // Newton's method
+    do
+    {
+      dx = -p.ddx(x)/p.d2dx(x);
+      x  = x + dx;
+    } while (std::abs(dx) > EPSILON);
+
+    // Save the value using the symmetry of the points
+    points[i] = -x;
+    points[n - 1 - i] = x;
+
+  }
+
+  // Fix the middle node
+  if ((n % 2) != 0)
+    points[n/2] = 0.0;
+}
+

=== added file 'ffc/ext/LobattoQuadrature.h'
--- ffc/ext/LobattoQuadrature.h	1970-01-01 00:00:00 +0000
+++ ffc/ext/LobattoQuadrature.h	2013-04-04 18:23:04 +0000
@@ -0,0 +1,47 @@
+// Copyright (C) 2003-2009 Anders Logg
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+#ifndef __LOBATTO_QUADRATURE_H
+#define __LOBATTO_QUADRATURE_H
+
+/// Lobatto (Gauss-Lobatto) quadrature on the interval [-1,1].
+/// The n quadrature points are given by the end-points -1 and 1,
+/// and the zeros of P{n-1}'(x), where P{n-1}(x) is the (n-1):th
+/// Legendre polynomial.
+///
+/// The quadrature points are computed using Newton's method, and
+/// the quadrature weights are computed by solving a linear system
+/// determined by the condition that Lobatto quadrature with n points
+/// should be exact for polynomials of degree 2n-3.
+
+#include <vector>
+
+class LobattoQuadrature
+{
+ public:
+
+  /// Create Lobatto quadrature with n points
+  LobattoQuadrature(unsigned int n);
+  //~LobattoQuadrature();
+
+  std::vector<double> points;
+};
+
+#endif

=== added file 'ffc/ext/RadauQuadrature.cpp'
--- ffc/ext/RadauQuadrature.cpp	1970-01-01 00:00:00 +0000
+++ ffc/ext/RadauQuadrature.cpp	2013-04-04 18:23:04 +0000
@@ -0,0 +1,84 @@
+// Copyright (C) 2003-2006 Anders Logg
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+#include "RadauQuadrature.h"
+#include "Legendre.h"
+
+#include <cmath>
+
+//-----------------------------------------------------------------------------
+RadauQuadrature::RadauQuadrature(unsigned int n)
+  : points(n+1)
+{
+  // Compute the Radau quadrature points in [-1,1] as -1 and the zeros
+  // of ( Pn-1(x) + Pn(x) ) / (1+x) where Pn is the n:th Legendre
+  // polynomial. Computation is a little different than for Gauss and
+  // Lobatto quadrature, since we don't know of any good initial
+  // approximation for the Newton iterations.
+
+  // Special case n = 1
+  if (n == 1)
+  {
+    points[0] = -1.0;
+    return;
+  }
+
+  Legendre p(n);
+  double x, dx, step, sign;
+
+  // Set size of stepping for seeking starting points
+  step = 1.0/(double(n - 1)*15.0);
+
+  // Set the first nodal point which is -1
+  points[0] = -1.0;
+
+  // Start at -1 + step
+  x = -1.0 + step;
+
+  // Set the sign at -1 + epsilon
+  sign = ((p.eval(n - 1, x) + p(x)) > 0 ? 1.0 : -1.0);
+
+  // Compute the rest of the nodes by Newton's method
+  for (unsigned int i = 1; i < n; i++)
+  {
+
+    // Step to a sign change
+    while ((p.eval(n - 1, x) + p(x))*sign > 0.0)
+      x += step;
+
+    // Newton's method
+    do
+    {
+      dx = -(p.eval(n-1, x) + p(x))/(p.ddx(n - 1, x) + p.ddx(x));
+      x  = x + dx;
+    } while (std::abs(dx) > EPSILON);
+
+    // Set the node value
+    points[i] = x;
+
+    // Fix step so that it's not too large
+    if (step > (points[i] - points[i-1])/10.0)
+      step = (points[i] - points[i-1])/10.0;
+
+    // Step forward
+    sign = -sign;
+    x += step;
+  }
+}

=== added file 'ffc/ext/RadauQuadrature.h'
--- ffc/ext/RadauQuadrature.h	1970-01-01 00:00:00 +0000
+++ ffc/ext/RadauQuadrature.h	2013-04-04 18:23:04 +0000
@@ -0,0 +1,49 @@
+// Copyright (C) 2003-2009 Anders Logg
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+#ifndef __RADAU_QUADRATURE_H
+#define __RADAU_QUADRATURE_H
+
+/// Radau (Gauss-Radau) quadrature on the interval [-1,1].
+/// The n quadrature points are given by the zeros of
+///
+///     ( Pn-1(x) + Pn(x) ) / (1+x)
+///
+/// where Pn is the n:th Legendre polynomial.
+///
+/// The quadrature points are computed using Newton's method, and
+/// the quadrature weights are computed by solving a linear system
+/// determined by the condition that Radau quadrature with n points
+/// should be exact for polynomials of degree 2n-2.
+
+#include <vector>
+
+class RadauQuadrature
+{
+ public:
+
+  /// Create Radau quadrature with n points
+  RadauQuadrature(unsigned int n);
+
+  std::vector<double> points;
+
+};
+
+#endif

=== added file 'ffc/ext/time_elements.cpp'
--- ffc/ext/time_elements.cpp	1970-01-01 00:00:00 +0000
+++ ffc/ext/time_elements.cpp	2013-04-04 18:23:04 +0000
@@ -0,0 +1,52 @@
+// Copyright (C) 2012 Benjamin Kehlet
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+#include <iostream>
+#include "LobattoQuadrature.h"
+#include "RadauQuadrature.h"
+#include "Legendre.h"
+
+void compute_lobatto_points(double* points, const unsigned int degree)
+{
+  // Compute the nodal basis
+  LobattoQuadrature lobatto(degree + 1);
+  for (unsigned int i = 0; i < degree +1; i++)
+    points[i] = (lobatto.points[i] + 1.0) / 2.0;
+}
+
+void compute_radau_points  (double* points, const unsigned int degree)
+{
+  RadauQuadrature radau(degree+1);
+
+  for (unsigned int i = 0; i < degree+1; i++)
+    points[degree-i] = (-radau.points[i] + 1.0) / 2.0;
+}
+
+void compute_legendre_coeffs(double* coeffs, const double *points, const unsigned int num_points, const unsigned int degree)
+{
+  for (unsigned int i = 0; i < degree; i++)
+  {
+    Legendre p(i);
+    for (unsigned int j = 0; j < num_points; j++)
+    {
+      coeffs[i*num_points + j] = p(points[j]*2.0 -1.0);
+    }
+  }
+}

=== added file 'ffc/ext/time_elements.h'
--- ffc/ext/time_elements.h	1970-01-01 00:00:00 +0000
+++ ffc/ext/time_elements.h	2013-04-04 18:23:04 +0000
@@ -0,0 +1,29 @@
+// Copyright (C) 2012 Benjamin Kehlet
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+
+void compute_lobatto_points(double* points, const unsigned int degree);
+
+void compute_radau_points  (double* points, const unsigned int degree);
+
+void compute_legendre_coeffs(double* coeffs, 
+			     const double *points, 
+			     const unsigned int num_points, 
+			     const unsigned int degree);

=== added file 'ffc/ext/time_elements_interface.cpp'
--- ffc/ext/time_elements_interface.cpp	1970-01-01 00:00:00 +0000
+++ ffc/ext/time_elements_interface.cpp	2013-04-04 18:23:04 +0000
@@ -0,0 +1,124 @@
+// Copyright (C) 2012 Benjamin Kehlet
+//
+// This file is part of FFC.
+//
+// FFC is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// FFC is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
+//
+// First added:  2012-08-20
+// Last changed: 2012-09-05
+
+
+#include <Python.h>
+#include <numpy/arrayobject.h>
+
+#include "time_elements.h"
+
+static PyObject *compute_lobatto_interface(PyObject *dummy, PyObject *args)
+{
+  int degree;
+
+  /* parse argument tuple */
+  if (!PyArg_ParseTuple(args, "i", &degree)) 
+  {
+    return NULL; /* PyArg_ParseTuple has raised an exception */
+  } 
+
+  npy_intp n = degree+1;
+
+  PyArrayObject *py_array_points =  (PyArrayObject*) PyArray_SimpleNew(1, &n, NPY_DOUBLE);
+
+  double *points = (double*)  PyArray_DATA(py_array_points);
+
+  compute_lobatto_points(points, degree);
+
+  // return values
+  return (PyObject*) py_array_points;
+}
+
+static PyObject *compute_radau_interface(PyObject *dummy, PyObject *args)
+{
+  int degree;
+
+  /* parse argument tuple */
+  if (!PyArg_ParseTuple(args, "i", &degree)) 
+  {
+    return NULL; /* PyArg_ParseTuple has raised an exception */
+  } 
+
+  npy_intp n = degree+1;
+
+  PyArrayObject *py_array_points =  (PyArrayObject*) PyArray_SimpleNew(1, &n, NPY_DOUBLE);
+
+  double *points = (double*)  PyArray_DATA(py_array_points);
+
+  compute_radau_points(points, degree);
+
+  // return values
+  return (PyObject*) py_array_points;
+}
+
+static PyObject *compute_legendre_coeffs_interface(PyObject *dummy, PyObject *args)
+{
+  PyArrayObject *points_array;
+
+  /* parse argument tuple */
+  if (!PyArg_ParseTuple(args, "O!", 
+			&PyArray_Type, &points_array))
+  {
+    return NULL; /* PyArg_ParseTuple has raised an exception */
+  } 
+
+  const npy_intp num_points = PyArray_DIMS(points_array)[0];
+  npy_intp dims[2] = { num_points, num_points };
+  PyArrayObject *py_array_coeffs =  (PyArrayObject*) PyArray_SimpleNew(2, dims, NPY_DOUBLE);
+
+  double *coeffs = (double*)  PyArray_DATA(py_array_coeffs);
+
+  compute_legendre_coeffs(coeffs, (double*) PyArray_DATA(points_array), num_points, num_points);
+
+  // return values
+  return (PyObject*) py_array_coeffs;
+}
+
+static char compute_lobatto_doc[] = \
+  "Doc string for compute_lobatto_points";
+static char compute_radau_doc[] = \
+  "Doc string for compute_radau_points";
+static char compute_legendre_coeffs_doc[] = \
+  "Doc string for compute_legendre_coeffs";
+
+
+static PyMethodDef time_elements_ext_methods[] = {
+  {"compute_lobatto_points",
+   compute_lobatto_interface, 
+   METH_VARARGS,  
+   compute_lobatto_doc},
+  {"compute_radau_points",
+   compute_radau_interface,
+   METH_VARARGS,
+   compute_radau_doc},
+  {"compute_legendre_coeffs",
+   compute_legendre_coeffs_interface,
+   METH_VARARGS,
+   compute_legendre_coeffs_doc},
+
+  {NULL, NULL}    
+};
+
+PyMODINIT_FUNC
+inittime_elements_ext(void)
+{
+   (void)Py_InitModule("time_elements_ext", time_elements_ext_methods);
+   import_array();
+}

=== modified file 'ffc/fiatinterface.py'
--- ffc/fiatinterface.py	2013-03-13 10:58:55 +0000
+++ ffc/fiatinterface.py	2013-04-04 19:07:05 +0000
@@ -16,7 +16,7 @@
 # along with FFC. If not, see <http://www.gnu.org/licenses/>.
 #
 # Modified by Garth N. Wells, 2009.
-# Modified by Marie Rognes, 2009-2010.
+# Modified by Marie Rognes, 2009-2013.
 # Modified by Martin Alnaes, 2013
 #
 # First added:  2009-03-06
@@ -32,6 +32,8 @@
 # FFC modules
 from ffc.log import debug, error, ffc_assert
 from ffc.quadratureelement import QuadratureElement as FFCQuadratureElement
+from ffc.timeelements import LobattoElement as FFCLobattoElement
+from ffc.timeelements import RadauElement as FFCRadauElement
 
 from ffc.mixedelement import MixedElement
 from ffc.restrictedelement import RestrictedElement
@@ -58,8 +60,10 @@
                       "Crouzeix-Raviart",
                       "Discontinuous Lagrange",
                       "Lagrange",
+                      "Lobatto",
                       "Nedelec 1st kind H(curl)",
                       "Nedelec 2nd kind H(curl)",
+                      "Radau",
                       "Raviart-Thomas",
                       "Real",
                       "Bubble",
@@ -139,6 +143,12 @@
         constant = _create_fiat_element(dg0_element)
         return SpaceOfReals(constant)
 
+    # Handle the specialized time elements
+    if family == "Lobatto" :
+        return FFCLobattoElement(ufl_element.degree())
+    if family == "Radau" :
+        return FFCRadauElement(ufl_element.degree())
+
     # FIXME: AL: Should this really be here?
     # Handle QuadratureElement
     if family == "Quadrature":

=== added file 'ffc/timeelements.py'
--- ffc/timeelements.py	1970-01-01 00:00:00 +0000
+++ ffc/timeelements.py	2013-04-04 18:23:04 +0000
@@ -0,0 +1,128 @@
+# Copyright (C) 2012 Benjamin Kehlet
+#
+# This file is part of FFC.
+#
+# FFC is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# FFC is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with FFC. If not, see <http://www.gnu.org/licenses/>.
+#
+# Modified by Marie E. Rognes, 2012
+#
+# First added:  2012-08-15
+# Last changed: 2012-09-07
+
+
+from FIAT import finite_element, polynomial_set, dual_set, functional, reference_element
+import time_elements_ext as ext
+import numpy
+
+class TimeElementDualSet(dual_set.DualSet):
+  """. """
+  def __init__(self, family, degree):
+    assert(family == "Lobatto" or family == "Radau"), \
+        "Unknown time element '%s'" % family
+    if family == "Lobatto" :
+        assert (degree > 0), "Lobatto not defined for degree < 1!"
+    else :
+        assert(degree >= 0), "Degree must be >= 0"
+
+    # ids is a map from mesh entity (numbers) to dof numbers
+    ids = {}
+
+    # dofs is a list of functionals
+    dofs = []
+
+    # Only defined in 1D (on an inteval)
+    cell = reference_element.UFCInterval()
+
+    self.coords = (ext.compute_lobatto_points(degree) if family == "Lobatto"
+                   else ext.compute_radau_points(degree))
+    points = [(c,) for c in self.coords]
+
+    # Create dofs from points
+    dofs = [functional.PointEvaluation(cell, point)
+            for point in points]
+
+    # Create ids
+    if family == "Lobatto":
+      ids[0] = {0: [0], 1: [len(points)-1]}
+      ids[1] = {0: range(1, len(points)-1)}
+    elif family == "Radau":
+      ids[0] = {0: [], 1: []}
+      ids[1] = {0: range(len(points))} # Treat all Radau points as internal
+    else:
+      error("Undefined family: %s" % family)
+
+    # Initialize dual set
+    dual_set.DualSet.__init__(self, dofs, cell, ids)
+
+class TimeElement(finite_element.FiniteElement):
+  """."""
+  def __init__(self, family, degree):
+    "Create time element with given (polynomial degree)."
+
+    # Only defined in 1D (on an inteval)
+    cell = reference_element.UFCInterval()
+
+    # Initialize polynomial space of degree 'degree'
+    polynomial_space = polynomial_set.ONPolynomialSet(cell, degree)
+
+    # Create dual (degrees of freedom)
+    dual = TimeElementDualSet(family, degree)
+
+    # Initialize super class
+    finite_element.FiniteElement.__init__(self,
+                                          polynomial_space,
+                                          dual,
+                                          degree
+                                          )
+
+  def compute_quadrature_weights(self) :
+    """Compute the quadrature weights by solving a linear system of equations
+    for exact integration of polynomials. We compute the integrals over
+    [-1,1] of the Legendre polynomials of degree <= n - 1; These integrals
+    are all zero, except for the integral of P0 which is 2.
+
+    This requires that the n-point quadrature rule is exact at least for
+    polynomials of degree n-1."""
+
+    n = len(self.dual.coords)
+
+    # Special case n = 0
+    if n == 0 :
+      weights[0] = 2.0;
+      return weights
+
+    # Initialize linear system
+    A = ext.compute_legendre_coeffs(self.dual.coords)
+
+    b = numpy.zeros(n)
+    b[0] = 2.0
+
+    weights = numpy.linalg.solve(A, b)
+
+    # Weights are computed on interval [-1, 1]. Scale to reference interval
+    return weights/2.0
+
+
+class LobattoElement(TimeElement):
+  """."""
+  def __init__(self, degree):
+      "Create Lobatto element with given (polynomial degree)."
+      TimeElement.__init__(self, "Lobatto", degree)
+
+class RadauElement(TimeElement):
+  """."""
+  def __init__(self, degree):
+      "Create Radau element with given (polynomial degree)."
+      TimeElement.__init__(self, "Radau", degree)
+

=== modified file 'setup.py'
--- setup.py	2013-03-25 09:20:08 +0000
+++ setup.py	2013-04-04 18:23:04 +0000
@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 
 import sys, platform
-from distutils.core import setup
+from distutils.core import setup, Extension
 from os import chdir
 from os.path import join, split
 
@@ -20,6 +20,13 @@
         batch_files.append(batch_file)
     scripts.extend(batch_files)
 
+ext = Extension("ffc.time_elements_ext",
+                ["ffc/ext/time_elements_interface.cpp",
+                 "ffc/ext/time_elements.cpp",
+                 "ffc/ext/LobattoQuadrature.cpp",
+                 "ffc/ext/RadauQuadrature.cpp",
+                 "ffc/ext/Legendre.cpp"])
+
 setup(name = "FFC",
       version = "1.2.0",
       description = "The FEniCS Form Compiler",
@@ -32,5 +39,6 @@
                   "ffc.dolfin"],
       package_dir={"ffc": "ffc"},
       scripts = scripts,
+      ext_modules = [ext],
       data_files = [(join("share", "man", "man1"),
                      [join("doc", "man", "man1", "ffc.1.gz")])])

=== added directory 'test/unit/elements'
=== added file 'test/unit/elements/test.py'
--- test/unit/elements/test.py	1970-01-01 00:00:00 +0000
+++ test/unit/elements/test.py	2013-04-04 19:07:05 +0000
@@ -0,0 +1,43 @@
+"""Unit tests for FFC finite elements"""
+
+# Copyright (C) 2013 Marie E. Rognes
+#
+# This file is part of FFC.
+#
+# FFC is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# FFC is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with FFC. If not, see <http://www.gnu.org/licenses/>.
+
+import unittest
+
+from ufl import interval
+from ufl import FiniteElement
+
+from ffc import compile_element
+
+class TestCompileElements(unittest.TestCase):
+
+    def testRadau(self):
+        "Test that Radau elements compile."
+        for degree in range(3):
+            element = FiniteElement("Radau", interval, degree)
+            compile_element(element)
+
+    def testLobatto(self):
+        "Test that Lobatto elements compile."
+        for degree in range(1, 4):
+            element = FiniteElement("Lobatto", interval, degree)
+            compile_element(element)
+
+
+if __name__ == "__main__":
+    unittest.main()