Subject: [Fenics] Patch for changed Python interface of Expressions

From: Johan Hake <[email-masked]>
Date: Mon, 6 Sep 2010 08:44:59 -0700

I also draged in some changes done in the C++ source code of la/eigenvalue.

Johan
# Bazaar merge directive format 2 (Bazaar 0.90)
# revision_id: [email-masked]
# target_branch: bzr+ssh://bazaar.launchpad.net/~fenics-core/fenics-\
#   doc/main/
# testament_sha1: 624bff4ed92bd2ceeef59bda34f7782bf777a0e4
# timestamp: 2010-09-06 08:42:44 -0700
# base_revision_id: [email-masked]\
#   ciqdciv4nc1e145w
# 
# Begin patch
=== modified file 'source/demos/la/eigenvalue/cpp/StiffnessMatrix.h'
--- source/demos/la/eigenvalue/cpp/StiffnessMatrix.h	2010-09-01 15:13:40 +0000
+++ source/demos/la/eigenvalue/cpp/StiffnessMatrix.h	2010-09-06 15:40:49 +0000
@@ -1,5 +1,5 @@
-// This code conforms with the UFC specification version 1.4
-// and was automatically generated by FFC version 0.9.2+.
+// This code conforms with the UFC specification version 1.4.2
+// and was automatically generated by FFC version 0.9.4+.
 //
 // This code was generated with the option '-l dolfin' and
 // contains DOLFIN-specific wrappers that depend on DOLFIN.
@@ -9,12 +9,13 @@
 //   cache_dir:                      ''
 //   convert_exceptions_to_warnings: False
 //   cpp_optimize:                   False
+//   cpp_optimize_flags:             '-O2'
 //   epsilon:                        1e-14
 //   form_postfix:                   True
 //   format:                         'dolfin'
-//   log_level:                      10
+//   log_level:                      20
 //   log_prefix:                     ''
-//   optimize:                       True
+//   optimize:                       False
 //   output_dir:                     '.'
 //   precision:                      15
 //   quadrature_degree:              'auto'
@@ -51,19 +52,19 @@
   /// Return a string identifying the finite element
   virtual const char* signature() const
   {
-    return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
+    return "FiniteElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1)";
   }
 
   /// Return the cell shape
   virtual ufc::shape cell_shape() const
   {
-    return ufc::triangle;
+    return ufc::tetrahedron;
   }
 
   /// Return the dimension of the finite element function space
   virtual unsigned int space_dimension() const
   {
-    return 3;
+    return 4;
   }
 
   /// Return the rank of the value space
@@ -90,21 +91,40 @@
     // Compute Jacobian of affine map from reference cell
     const double J_00 = x[1][0] - x[0][0];
     const double J_01 = x[2][0] - x[0][0];
+    const double J_02 = x[3][0] - x[0][0];
     const double J_10 = x[1][1] - x[0][1];
     const double J_11 = x[2][1] - x[0][1];
+    const double J_12 = x[3][1] - x[0][1];
+    const double J_20 = x[1][2] - x[0][2];
+    const double J_21 = x[2][2] - x[0][2];
+    const double J_22 = x[3][2] - x[0][2];
+    
+    // Compute sub determinants
+    const double d_00 = J_11*J_22 - J_12*J_21;
+    const double d_01 = J_12*J_20 - J_10*J_22;
+    const double d_02 = J_10*J_21 - J_11*J_20;
+    const double d_10 = J_02*J_21 - J_01*J_22;
+    const double d_11 = J_00*J_22 - J_02*J_20;
+    const double d_12 = J_01*J_20 - J_00*J_21;
+    const double d_20 = J_01*J_12 - J_02*J_11;
+    const double d_21 = J_02*J_10 - J_00*J_12;
+    const double d_22 = J_00*J_11 - J_01*J_10;
     
     // Compute determinant of Jacobian
-    double detJ = J_00*J_11 - J_01*J_10;
+    double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
     
     // Compute inverse of Jacobian
     
     // Compute constants
-    const double C0 = x[1][0] + x[2][0];
-    const double C1 = x[1][1] + x[2][1];
+    const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0];
+    const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1];
+    const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2];
     
     // Get coordinates and map to the reference (FIAT) element
-    double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
-    double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
+    double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ;
+    double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ;
+    double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ;
+    
     
     // Reset values.
     *values = 0.000000000000000;
@@ -114,37 +134,49 @@
       {
         
       // Array of basisvalues.
-      double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
+      double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000};
       
       // Declare helper variables.
       unsigned int rr = 0;
       unsigned int ss = 0;
-      double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
+      double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
       
       // Compute basisvalues.
       basisvalues[0] = 1.000000000000000;
       basisvalues[1] = tmp0;
       for (unsigned int r = 0; r < 1; r++)
       {
-        rr = (r + 1)*(r + 1 + 1)/2 + 1;
-        ss = r*(r + 1)/2;
-        basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
+        rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
+        ss = r*(r + 1)*(r + 2)/6;
+        basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
+      }// end loop over 'r'
+      for (unsigned int r = 0; r < 1; r++)
+      {
+        for (unsigned int s = 0; s < 1 - r; s++)
+        {
+          rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
+          ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
+          basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
+        }// end loop over 's'
       }// end loop over 'r'
       for (unsigned int r = 0; r < 2; r++)
       {
         for (unsigned int s = 0; s < 2 - r; s++)
         {
-          rr = (r + s)*(r + s + 1)/2 + s;
-          basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
+          for (unsigned int t = 0; t < 2 - r - s; t++)
+          {
+            rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
+            basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
+          }// end loop over 't'
         }// end loop over 's'
       }// end loop over 'r'
       
       // Table(s) of coefficients.
-      static const double coefficients0[3] = \
-      {0.471404520791032, -0.288675134594813, -0.166666666666667};
+      static const double coefficients0[4] = \
+      {0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993};
       
       // Compute value(s).
-      for (unsigned int r = 0; r < 3; r++)
+      for (unsigned int r = 0; r < 4; r++)
       {
         *values += coefficients0[r]*basisvalues[r];
       }// end loop over 'r'
@@ -154,37 +186,49 @@
       {
         
       // Array of basisvalues.
-      double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
+      double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000};
       
       // Declare helper variables.
       unsigned int rr = 0;
       unsigned int ss = 0;
-      double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
+      double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
       
       // Compute basisvalues.
       basisvalues[0] = 1.000000000000000;
       basisvalues[1] = tmp0;
       for (unsigned int r = 0; r < 1; r++)
       {
-        rr = (r + 1)*(r + 1 + 1)/2 + 1;
-        ss = r*(r + 1)/2;
-        basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
+        rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
+        ss = r*(r + 1)*(r + 2)/6;
+        basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
+      }// end loop over 'r'
+      for (unsigned int r = 0; r < 1; r++)
+      {
+        for (unsigned int s = 0; s < 1 - r; s++)
+        {
+          rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
+          ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
+          basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
+        }// end loop over 's'
       }// end loop over 'r'
       for (unsigned int r = 0; r < 2; r++)
       {
         for (unsigned int s = 0; s < 2 - r; s++)
         {
-          rr = (r + s)*(r + s + 1)/2 + s;
-          basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
+          for (unsigned int t = 0; t < 2 - r - s; t++)
+          {
+            rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
+            basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
+          }// end loop over 't'
         }// end loop over 's'
       }// end loop over 'r'
       
       // Table(s) of coefficients.
-      static const double coefficients0[3] = \
-      {0.471404520791032, 0.288675134594813, -0.166666666666667};
+      static const double coefficients0[4] = \
+      {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993};
       
       // Compute value(s).
-      for (unsigned int r = 0; r < 3; r++)
+      for (unsigned int r = 0; r < 4; r++)
       {
         *values += coefficients0[r]*basisvalues[r];
       }// end loop over 'r'
@@ -194,37 +238,101 @@
       {
         
       // Array of basisvalues.
-      double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
-      
-      // Declare helper variables.
-      unsigned int rr = 0;
-      unsigned int ss = 0;
-      double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
-      
-      // Compute basisvalues.
-      basisvalues[0] = 1.000000000000000;
-      basisvalues[1] = tmp0;
-      for (unsigned int r = 0; r < 1; r++)
-      {
-        rr = (r + 1)*(r + 1 + 1)/2 + 1;
-        ss = r*(r + 1)/2;
-        basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
-      }// end loop over 'r'
-      for (unsigned int r = 0; r < 2; r++)
-      {
-        for (unsigned int s = 0; s < 2 - r; s++)
-        {
-          rr = (r + s)*(r + s + 1)/2 + s;
-          basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
-        }// end loop over 's'
-      }// end loop over 'r'
-      
-      // Table(s) of coefficients.
-      static const double coefficients0[3] = \
-      {0.471404520791032, 0.000000000000000, 0.333333333333333};
-      
-      // Compute value(s).
-      for (unsigned int r = 0; r < 3; r++)
+      double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000};
+      
+      // Declare helper variables.
+      unsigned int rr = 0;
+      unsigned int ss = 0;
+      double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
+      
+      // Compute basisvalues.
+      basisvalues[0] = 1.000000000000000;
+      basisvalues[1] = tmp0;
+      for (unsigned int r = 0; r < 1; r++)
+      {
+        rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
+        ss = r*(r + 1)*(r + 2)/6;
+        basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
+      }// end loop over 'r'
+      for (unsigned int r = 0; r < 1; r++)
+      {
+        for (unsigned int s = 0; s < 1 - r; s++)
+        {
+          rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
+          ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
+          basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
+        }// end loop over 's'
+      }// end loop over 'r'
+      for (unsigned int r = 0; r < 2; r++)
+      {
+        for (unsigned int s = 0; s < 2 - r; s++)
+        {
+          for (unsigned int t = 0; t < 2 - r - s; t++)
+          {
+            rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
+            basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
+          }// end loop over 't'
+        }// end loop over 's'
+      }// end loop over 'r'
+      
+      // Table(s) of coefficients.
+      static const double coefficients0[4] = \
+      {0.288675134594813, 0.000000000000000, 0.210818510677892, -0.074535599249993};
+      
+      // Compute value(s).
+      for (unsigned int r = 0; r < 4; r++)
+      {
+        *values += coefficients0[r]*basisvalues[r];
+      }// end loop over 'r'
+        break;
+      }
+    case 3:
+      {
+        
+      // Array of basisvalues.
+      double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000};
+      
+      // Declare helper variables.
+      unsigned int rr = 0;
+      unsigned int ss = 0;
+      double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
+      
+      // Compute basisvalues.
+      basisvalues[0] = 1.000000000000000;
+      basisvalues[1] = tmp0;
+      for (unsigned int r = 0; r < 1; r++)
+      {
+        rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
+        ss = r*(r + 1)*(r + 2)/6;
+        basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
+      }// end loop over 'r'
+      for (unsigned int r = 0; r < 1; r++)
+      {
+        for (unsigned int s = 0; s < 1 - r; s++)
+        {
+          rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
+          ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
+          basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
+        }// end loop over 's'
+      }// end loop over 'r'
+      for (unsigned int r = 0; r < 2; r++)
+      {
+        for (unsigned int s = 0; s < 2 - r; s++)
+        {
+          for (unsigned int t = 0; t < 2 - r - s; t++)
+          {
+            rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
+            basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
+          }// end loop over 't'
+        }// end loop over 's'
+      }// end loop over 'r'
+      
+      // Table(s) of coefficients.
+      static const double coefficients0[4] = \
+      {0.288675134594813, 0.000000000000000, 0.000000000000000, 0.223606797749979};
+      
+      // Compute value(s).
+      for (unsigned int r = 0; r < 4; r++)
       {
         *values += coefficients0[r]*basisvalues[r];
       }// end loop over 'r'
@@ -243,7 +351,7 @@
     double dof_values = 0.000000000000000;
     
     // Loop dofs and call evaluate_basis.
-    for (unsigned int r = 0; r < 3; r++)
+    for (unsigned int r = 0; r < 4; r++)
     {
       evaluate_basis(r, &dof_values, coordinates, c);
       values[r] = dof_values;
@@ -263,31 +371,55 @@
     // Compute Jacobian of affine map from reference cell
     const double J_00 = x[1][0] - x[0][0];
     const double J_01 = x[2][0] - x[0][0];
+    const double J_02 = x[3][0] - x[0][0];
     const double J_10 = x[1][1] - x[0][1];
     const double J_11 = x[2][1] - x[0][1];
+    const double J_12 = x[3][1] - x[0][1];
+    const double J_20 = x[1][2] - x[0][2];
+    const double J_21 = x[2][2] - x[0][2];
+    const double J_22 = x[3][2] - x[0][2];
+    
+    // Compute sub determinants
+    const double d_00 = J_11*J_22 - J_12*J_21;
+    const double d_01 = J_12*J_20 - J_10*J_22;
+    const double d_02 = J_10*J_21 - J_11*J_20;
+    const double d_10 = J_02*J_21 - J_01*J_22;
+    const double d_11 = J_00*J_22 - J_02*J_20;
+    const double d_12 = J_01*J_20 - J_00*J_21;
+    const double d_20 = J_01*J_12 - J_02*J_11;
+    const double d_21 = J_02*J_10 - J_00*J_12;
+    const double d_22 = J_00*J_11 - J_01*J_10;
     
     // Compute determinant of Jacobian
-    double detJ = J_00*J_11 - J_01*J_10;
+    double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
     
     // Compute inverse of Jacobian
-    const double K_00 =  J_11 / detJ;
-    const double K_01 = -J_01 / detJ;
-    const double K_10 = -J_10 / detJ;
-    const double K_11 =  J_00 / detJ;
+    const double K_00 = d_00 / detJ;
+    const double K_01 = d_10 / detJ;
+    const double K_02 = d_20 / detJ;
+    const double K_10 = d_01 / detJ;
+    const double K_11 = d_11 / detJ;
+    const double K_12 = d_21 / detJ;
+    const double K_20 = d_02 / detJ;
+    const double K_21 = d_12 / detJ;
+    const double K_22 = d_22 / detJ;
     
     // Compute constants
-    const double C0 = x[1][0] + x[2][0];
-    const double C1 = x[1][1] + x[2][1];
+    const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0];
+    const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1];
+    const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2];
     
     // Get coordinates and map to the reference (FIAT) element
-    double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
-    double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
+    double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ;
+    double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ;
+    double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ;
+    
     
     // Compute number of derivatives.
     unsigned int num_derivatives = 1;
     for (unsigned int r = 0; r < n; r++)
     {
-      num_derivatives *= 2;
+      num_derivatives *= 3;
     }// end loop over 'r'
     
     // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
@@ -306,7 +438,7 @@
       {
         for (unsigned int col = n-1; col+1 > 0; col--)
         {
-          if (combinations[row][col] + 1 > 1)
+          if (combinations[row][col] + 1 > 2)
             combinations[row][col] = 0;
           else
           {
@@ -318,7 +450,7 @@
     }
     
     // Compute inverse of Jacobian
-    const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
+    const double Jinv[3][3] = {{K_00, K_01, K_02}, {K_10, K_11, K_12}, {K_20, K_21, K_22}};
     
     // Declare transformation matrix
     // Declare pointer to two dimensional array and initialise
@@ -353,45 +485,65 @@
       {
         
       // Array of basisvalues.
-      double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
+      double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000};
       
       // Declare helper variables.
       unsigned int rr = 0;
       unsigned int ss = 0;
-      double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
+      double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
       
       // Compute basisvalues.
       basisvalues[0] = 1.000000000000000;
       basisvalues[1] = tmp0;
       for (unsigned int r = 0; r < 1; r++)
       {
-        rr = (r + 1)*(r + 1 + 1)/2 + 1;
-        ss = r*(r + 1)/2;
-        basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
+        rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
+        ss = r*(r + 1)*(r + 2)/6;
+        basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
+      }// end loop over 'r'
+      for (unsigned int r = 0; r < 1; r++)
+      {
+        for (unsigned int s = 0; s < 1 - r; s++)
+        {
+          rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
+          ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
+          basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
+        }// end loop over 's'
       }// end loop over 'r'
       for (unsigned int r = 0; r < 2; r++)
       {
         for (unsigned int s = 0; s < 2 - r; s++)
         {
-          rr = (r + s)*(r + s + 1)/2 + s;
-          basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
+          for (unsigned int t = 0; t < 2 - r - s; t++)
+          {
+            rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
+            basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
+          }// end loop over 't'
         }// end loop over 's'
       }// end loop over 'r'
       
       // Table(s) of coefficients.
-      static const double coefficients0[3] = \
-      {0.471404520791032, -0.288675134594813, -0.166666666666667};
+      static const double coefficients0[4] = \
+      {0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993};
       
       // Tables of derivatives of the polynomial base (transpose).
-      static const double dmats0[3][3] = \
-      {{0.000000000000000, 0.000000000000000, 0.000000000000000},
-      {4.898979485566356, 0.000000000000000, 0.000000000000000},
-      {0.000000000000000, 0.000000000000000, 0.000000000000000}};
-      
-      static const double dmats1[3][3] = \
-      {{0.000000000000000, 0.000000000000000, 0.000000000000000},
-      {2.449489742783178, 0.000000000000000, 0.000000000000000},
-      {4.242640687119285, 0.000000000000000, 0.000000000000000}};
+      static const double dmats0[4][4] = \
+      {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
+      
+      static const double dmats1[4][4] = \
+      {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
+      
+      static const double dmats2[4][4] = \
+      {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
       
       // Compute reference derivatives.
       // Declare pointer to array of derivatives on FIAT element.
@@ -402,24 +554,26 @@
       }// end loop over 'r'
       
       // Declare derivative matrix (of polynomial basis).
-      double dmats[3][3] = \
-      {{1.000000000000000, 0.000000000000000, 0.000000000000000},
-      {0.000000000000000, 1.000000000000000, 0.000000000000000},
-      {0.000000000000000, 0.000000000000000, 1.000000000000000}};
+      double dmats[4][4] = \
+      {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}};
       
       // Declare (auxiliary) derivative matrix (of polynomial basis).
-      double dmats_old[3][3] = \
-      {{1.000000000000000, 0.000000000000000, 0.000000000000000},
-      {0.000000000000000, 1.000000000000000, 0.000000000000000},
-      {0.000000000000000, 0.000000000000000, 1.000000000000000}};
+      double dmats_old[4][4] = \
+      {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}};
       
       // Loop possible derivatives.
       for (unsigned int r = 0; r < num_derivatives; r++)
       {
         // Resetting dmats values to compute next derivative.
-        for (unsigned int t = 0; t < 3; t++)
+        for (unsigned int t = 0; t < 4; t++)
         {
-          for (unsigned int u = 0; u < 3; u++)
+          for (unsigned int u = 0; u < 4; u++)
           {
             dmats[t][u] = 0.000000000000000;
             if (t == u)
@@ -434,9 +588,9 @@
         for (unsigned int s = 0; s < n; s++)
         {
           // Updating dmats_old with new values and resetting dmats.
-          for (unsigned int t = 0; t < 3; t++)
+          for (unsigned int t = 0; t < 4; t++)
           {
-            for (unsigned int u = 0; u < 3; u++)
+            for (unsigned int u = 0; u < 4; u++)
             {
               dmats_old[t][u] = dmats[t][u];
               dmats[t][u] = 0.000000000000000;
@@ -446,11 +600,11 @@
           // Update dmats using an inner product.
           if (combinations[r][s] == 0)
           {
-          for (unsigned int t = 0; t < 3; t++)
+          for (unsigned int t = 0; t < 4; t++)
           {
-            for (unsigned int u = 0; u < 3; u++)
+            for (unsigned int u = 0; u < 4; u++)
             {
-              for (unsigned int tu = 0; tu < 3; tu++)
+              for (unsigned int tu = 0; tu < 4; tu++)
               {
                 dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
               }// end loop over 'tu'
@@ -460,11 +614,11 @@
           
           if (combinations[r][s] == 1)
           {
-          for (unsigned int t = 0; t < 3; t++)
+          for (unsigned int t = 0; t < 4; t++)
           {
-            for (unsigned int u = 0; u < 3; u++)
+            for (unsigned int u = 0; u < 4; u++)
             {
-              for (unsigned int tu = 0; tu < 3; tu++)
+              for (unsigned int tu = 0; tu < 4; tu++)
               {
                 dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
               }// end loop over 'tu'
@@ -472,10 +626,24 @@
           }// end loop over 't'
           }
           
+          if (combinations[r][s] == 2)
+          {
+          for (unsigned int t = 0; t < 4; t++)
+          {
+            for (unsigned int u = 0; u < 4; u++)
+            {
+              for (unsigned int tu = 0; tu < 4; tu++)
+              {
+                dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
+              }// end loop over 'tu'
+            }// end loop over 'u'
+          }// end loop over 't'
+          }
+          
         }// end loop over 's'
-        for (unsigned int s = 0; s < 3; s++)
+        for (unsigned int s = 0; s < 4; s++)
         {
-          for (unsigned int t = 0; t < 3; t++)
+          for (unsigned int t = 0; t < 4; t++)
           {
             derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
           }// end loop over 't'
@@ -511,45 +679,65 @@
       {
         
       // Array of basisvalues.
-      double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
+      double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000};
       
       // Declare helper variables.
       unsigned int rr = 0;
       unsigned int ss = 0;
-      double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
+      double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
       
       // Compute basisvalues.
       basisvalues[0] = 1.000000000000000;
       basisvalues[1] = tmp0;
       for (unsigned int r = 0; r < 1; r++)
       {
-        rr = (r + 1)*(r + 1 + 1)/2 + 1;
-        ss = r*(r + 1)/2;
-        basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
+        rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
+        ss = r*(r + 1)*(r + 2)/6;
+        basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
+      }// end loop over 'r'
+      for (unsigned int r = 0; r < 1; r++)
+      {
+        for (unsigned int s = 0; s < 1 - r; s++)
+        {
+          rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
+          ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
+          basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
+        }// end loop over 's'
       }// end loop over 'r'
       for (unsigned int r = 0; r < 2; r++)
       {
         for (unsigned int s = 0; s < 2 - r; s++)
         {
-          rr = (r + s)*(r + s + 1)/2 + s;
-          basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
+          for (unsigned int t = 0; t < 2 - r - s; t++)
+          {
+            rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
+            basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
+          }// end loop over 't'
         }// end loop over 's'
       }// end loop over 'r'
       
       // Table(s) of coefficients.
-      static const double coefficients0[3] = \
-      {0.471404520791032, 0.288675134594813, -0.166666666666667};
+      static const double coefficients0[4] = \
+      {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993};
       
       // Tables of derivatives of the polynomial base (transpose).
-      static const double dmats0[3][3] = \
-      {{0.000000000000000, 0.000000000000000, 0.000000000000000},
-      {4.898979485566356, 0.000000000000000, 0.000000000000000},
-      {0.000000000000000, 0.000000000000000, 0.000000000000000}};
-      
-      static const double dmats1[3][3] = \
-      {{0.000000000000000, 0.000000000000000, 0.000000000000000},
-      {2.449489742783178, 0.000000000000000, 0.000000000000000},
-      {4.242640687119285, 0.000000000000000, 0.000000000000000}};
+      static const double dmats0[4][4] = \
+      {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
+      
+      static const double dmats1[4][4] = \
+      {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
+      
+      static const double dmats2[4][4] = \
+      {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
       
       // Compute reference derivatives.
       // Declare pointer to array of derivatives on FIAT element.
@@ -560,24 +748,26 @@
       }// end loop over 'r'
       
       // Declare derivative matrix (of polynomial basis).
-      double dmats[3][3] = \
-      {{1.000000000000000, 0.000000000000000, 0.000000000000000},
-      {0.000000000000000, 1.000000000000000, 0.000000000000000},
-      {0.000000000000000, 0.000000000000000, 1.000000000000000}};
+      double dmats[4][4] = \
+      {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}};
       
       // Declare (auxiliary) derivative matrix (of polynomial basis).
-      double dmats_old[3][3] = \
-      {{1.000000000000000, 0.000000000000000, 0.000000000000000},
-      {0.000000000000000, 1.000000000000000, 0.000000000000000},
-      {0.000000000000000, 0.000000000000000, 1.000000000000000}};
+      double dmats_old[4][4] = \
+      {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}};
       
       // Loop possible derivatives.
       for (unsigned int r = 0; r < num_derivatives; r++)
       {
         // Resetting dmats values to compute next derivative.
-        for (unsigned int t = 0; t < 3; t++)
+        for (unsigned int t = 0; t < 4; t++)
         {
-          for (unsigned int u = 0; u < 3; u++)
+          for (unsigned int u = 0; u < 4; u++)
           {
             dmats[t][u] = 0.000000000000000;
             if (t == u)
@@ -592,9 +782,9 @@
         for (unsigned int s = 0; s < n; s++)
         {
           // Updating dmats_old with new values and resetting dmats.
-          for (unsigned int t = 0; t < 3; t++)
+          for (unsigned int t = 0; t < 4; t++)
           {
-            for (unsigned int u = 0; u < 3; u++)
+            for (unsigned int u = 0; u < 4; u++)
             {
               dmats_old[t][u] = dmats[t][u];
               dmats[t][u] = 0.000000000000000;
@@ -604,11 +794,11 @@
           // Update dmats using an inner product.
           if (combinations[r][s] == 0)
           {
-          for (unsigned int t = 0; t < 3; t++)
+          for (unsigned int t = 0; t < 4; t++)
           {
-            for (unsigned int u = 0; u < 3; u++)
+            for (unsigned int u = 0; u < 4; u++)
             {
-              for (unsigned int tu = 0; tu < 3; tu++)
+              for (unsigned int tu = 0; tu < 4; tu++)
               {
                 dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
               }// end loop over 'tu'
@@ -618,11 +808,11 @@
           
           if (combinations[r][s] == 1)
           {
-          for (unsigned int t = 0; t < 3; t++)
+          for (unsigned int t = 0; t < 4; t++)
           {
-            for (unsigned int u = 0; u < 3; u++)
+            for (unsigned int u = 0; u < 4; u++)
             {
-              for (unsigned int tu = 0; tu < 3; tu++)
+              for (unsigned int tu = 0; tu < 4; tu++)
               {
                 dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
               }// end loop over 'tu'
@@ -630,10 +820,24 @@
           }// end loop over 't'
           }
           
+          if (combinations[r][s] == 2)
+          {
+          for (unsigned int t = 0; t < 4; t++)
+          {
+            for (unsigned int u = 0; u < 4; u++)
+            {
+              for (unsigned int tu = 0; tu < 4; tu++)
+              {
+                dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
+              }// end loop over 'tu'
+            }// end loop over 'u'
+          }// end loop over 't'
+          }
+          
         }// end loop over 's'
-        for (unsigned int s = 0; s < 3; s++)
+        for (unsigned int s = 0; s < 4; s++)
         {
-          for (unsigned int t = 0; t < 3; t++)
+          for (unsigned int t = 0; t < 4; t++)
           {
             derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
           }// end loop over 't'
@@ -669,129 +873,359 @@
       {
         
       // Array of basisvalues.
-      double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
-      
-      // Declare helper variables.
-      unsigned int rr = 0;
-      unsigned int ss = 0;
-      double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
-      
-      // Compute basisvalues.
-      basisvalues[0] = 1.000000000000000;
-      basisvalues[1] = tmp0;
-      for (unsigned int r = 0; r < 1; r++)
-      {
-        rr = (r + 1)*(r + 1 + 1)/2 + 1;
-        ss = r*(r + 1)/2;
-        basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
-      }// end loop over 'r'
-      for (unsigned int r = 0; r < 2; r++)
-      {
-        for (unsigned int s = 0; s < 2 - r; s++)
-        {
-          rr = (r + s)*(r + s + 1)/2 + s;
-          basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
-        }// end loop over 's'
-      }// end loop over 'r'
-      
-      // Table(s) of coefficients.
-      static const double coefficients0[3] = \
-      {0.471404520791032, 0.000000000000000, 0.333333333333333};
-      
-      // Tables of derivatives of the polynomial base (transpose).
-      static const double dmats0[3][3] = \
-      {{0.000000000000000, 0.000000000000000, 0.000000000000000},
-      {4.898979485566356, 0.000000000000000, 0.000000000000000},
-      {0.000000000000000, 0.000000000000000, 0.000000000000000}};
-      
-      static const double dmats1[3][3] = \
-      {{0.000000000000000, 0.000000000000000, 0.000000000000000},
-      {2.449489742783178, 0.000000000000000, 0.000000000000000},
-      {4.242640687119285, 0.000000000000000, 0.000000000000000}};
-      
-      // Compute reference derivatives.
-      // Declare pointer to array of derivatives on FIAT element.
-      double *derivatives = new double[num_derivatives];
-      for (unsigned int r = 0; r < num_derivatives; r++)
-      {
-        derivatives[r] = 0.000000000000000;
-      }// end loop over 'r'
-      
-      // Declare derivative matrix (of polynomial basis).
-      double dmats[3][3] = \
-      {{1.000000000000000, 0.000000000000000, 0.000000000000000},
-      {0.000000000000000, 1.000000000000000, 0.000000000000000},
-      {0.000000000000000, 0.000000000000000, 1.000000000000000}};
-      
-      // Declare (auxiliary) derivative matrix (of polynomial basis).
-      double dmats_old[3][3] = \
-      {{1.000000000000000, 0.000000000000000, 0.000000000000000},
-      {0.000000000000000, 1.000000000000000, 0.000000000000000},
-      {0.000000000000000, 0.000000000000000, 1.000000000000000}};
-      
-      // Loop possible derivatives.
-      for (unsigned int r = 0; r < num_derivatives; r++)
-      {
-        // Resetting dmats values to compute next derivative.
-        for (unsigned int t = 0; t < 3; t++)
-        {
-          for (unsigned int u = 0; u < 3; u++)
-          {
-            dmats[t][u] = 0.000000000000000;
-            if (t == u)
-            {
-            dmats[t][u] = 1.000000000000000;
-            }
-            
-          }// end loop over 'u'
-        }// end loop over 't'
-        
-        // Looping derivative order to generate dmats.
-        for (unsigned int s = 0; s < n; s++)
-        {
-          // Updating dmats_old with new values and resetting dmats.
-          for (unsigned int t = 0; t < 3; t++)
-          {
-            for (unsigned int u = 0; u < 3; u++)
-            {
-              dmats_old[t][u] = dmats[t][u];
-              dmats[t][u] = 0.000000000000000;
-            }// end loop over 'u'
-          }// end loop over 't'
-          
-          // Update dmats using an inner product.
-          if (combinations[r][s] == 0)
-          {
-          for (unsigned int t = 0; t < 3; t++)
-          {
-            for (unsigned int u = 0; u < 3; u++)
-            {
-              for (unsigned int tu = 0; tu < 3; tu++)
-              {
-                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
-              }// end loop over 'tu'
-            }// end loop over 'u'
-          }// end loop over 't'
-          }
-          
-          if (combinations[r][s] == 1)
-          {
-          for (unsigned int t = 0; t < 3; t++)
-          {
-            for (unsigned int u = 0; u < 3; u++)
-            {
-              for (unsigned int tu = 0; tu < 3; tu++)
-              {
-                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
-              }// end loop over 'tu'
-            }// end loop over 'u'
-          }// end loop over 't'
-          }
-          
-        }// end loop over 's'
-        for (unsigned int s = 0; s < 3; s++)
-        {
-          for (unsigned int t = 0; t < 3; t++)
+      double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000};
+      
+      // Declare helper variables.
+      unsigned int rr = 0;
+      unsigned int ss = 0;
+      double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
+      
+      // Compute basisvalues.
+      basisvalues[0] = 1.000000000000000;
+      basisvalues[1] = tmp0;
+      for (unsigned int r = 0; r < 1; r++)
+      {
+        rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
+        ss = r*(r + 1)*(r + 2)/6;
+        basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
+      }// end loop over 'r'
+      for (unsigned int r = 0; r < 1; r++)
+      {
+        for (unsigned int s = 0; s < 1 - r; s++)
+        {
+          rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
+          ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
+          basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
+        }// end loop over 's'
+      }// end loop over 'r'
+      for (unsigned int r = 0; r < 2; r++)
+      {
+        for (unsigned int s = 0; s < 2 - r; s++)
+        {
+          for (unsigned int t = 0; t < 2 - r - s; t++)
+          {
+            rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
+            basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
+          }// end loop over 't'
+        }// end loop over 's'
+      }// end loop over 'r'
+      
+      // Table(s) of coefficients.
+      static const double coefficients0[4] = \
+      {0.288675134594813, 0.000000000000000, 0.210818510677892, -0.074535599249993};
+      
+      // Tables of derivatives of the polynomial base (transpose).
+      static const double dmats0[4][4] = \
+      {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
+      
+      static const double dmats1[4][4] = \
+      {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
+      
+      static const double dmats2[4][4] = \
+      {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
+      
+      // Compute reference derivatives.
+      // Declare pointer to array of derivatives on FIAT element.
+      double *derivatives = new double[num_derivatives];
+      for (unsigned int r = 0; r < num_derivatives; r++)
+      {
+        derivatives[r] = 0.000000000000000;
+      }// end loop over 'r'
+      
+      // Declare derivative matrix (of polynomial basis).
+      double dmats[4][4] = \
+      {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}};
+      
+      // Declare (auxiliary) derivative matrix (of polynomial basis).
+      double dmats_old[4][4] = \
+      {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}};
+      
+      // Loop possible derivatives.
+      for (unsigned int r = 0; r < num_derivatives; r++)
+      {
+        // Resetting dmats values to compute next derivative.
+        for (unsigned int t = 0; t < 4; t++)
+        {
+          for (unsigned int u = 0; u < 4; u++)
+          {
+            dmats[t][u] = 0.000000000000000;
+            if (t == u)
+            {
+            dmats[t][u] = 1.000000000000000;
+            }
+            
+          }// end loop over 'u'
+        }// end loop over 't'
+        
+        // Looping derivative order to generate dmats.
+        for (unsigned int s = 0; s < n; s++)
+        {
+          // Updating dmats_old with new values and resetting dmats.
+          for (unsigned int t = 0; t < 4; t++)
+          {
+            for (unsigned int u = 0; u < 4; u++)
+            {
+              dmats_old[t][u] = dmats[t][u];
+              dmats[t][u] = 0.000000000000000;
+            }// end loop over 'u'
+          }// end loop over 't'
+          
+          // Update dmats using an inner product.
+          if (combinations[r][s] == 0)
+          {
+          for (unsigned int t = 0; t < 4; t++)
+          {
+            for (unsigned int u = 0; u < 4; u++)
+            {
+              for (unsigned int tu = 0; tu < 4; tu++)
+              {
+                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
+              }// end loop over 'tu'
+            }// end loop over 'u'
+          }// end loop over 't'
+          }
+          
+          if (combinations[r][s] == 1)
+          {
+          for (unsigned int t = 0; t < 4; t++)
+          {
+            for (unsigned int u = 0; u < 4; u++)
+            {
+              for (unsigned int tu = 0; tu < 4; tu++)
+              {
+                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
+              }// end loop over 'tu'
+            }// end loop over 'u'
+          }// end loop over 't'
+          }
+          
+          if (combinations[r][s] == 2)
+          {
+          for (unsigned int t = 0; t < 4; t++)
+          {
+            for (unsigned int u = 0; u < 4; u++)
+            {
+              for (unsigned int tu = 0; tu < 4; tu++)
+              {
+                dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
+              }// end loop over 'tu'
+            }// end loop over 'u'
+          }// end loop over 't'
+          }
+          
+        }// end loop over 's'
+        for (unsigned int s = 0; s < 4; s++)
+        {
+          for (unsigned int t = 0; t < 4; t++)
+          {
+            derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
+          }// end loop over 't'
+        }// end loop over 's'
+      }// end loop over 'r'
+      
+      // Transform derivatives back to physical element
+      for (unsigned int r = 0; r < num_derivatives; r++)
+      {
+        for (unsigned int s = 0; s < num_derivatives; s++)
+        {
+          values[r] += transform[r][s]*derivatives[s];
+        }// end loop over 's'
+      }// end loop over 'r'
+      
+      // Delete pointer to array of derivatives on FIAT element
+      delete [] derivatives;
+      
+      // Delete pointer to array of combinations of derivatives and transform
+      for (unsigned int r = 0; r < num_derivatives; r++)
+      {
+        delete [] combinations[r];
+      }// end loop over 'r'
+      delete [] combinations;
+      for (unsigned int r = 0; r < num_derivatives; r++)
+      {
+        delete [] transform[r];
+      }// end loop over 'r'
+      delete [] transform;
+        break;
+      }
+    case 3:
+      {
+        
+      // Array of basisvalues.
+      double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000};
+      
+      // Declare helper variables.
+      unsigned int rr = 0;
+      unsigned int ss = 0;
+      double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
+      
+      // Compute basisvalues.
+      basisvalues[0] = 1.000000000000000;
+      basisvalues[1] = tmp0;
+      for (unsigned int r = 0; r < 1; r++)
+      {
+        rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
+        ss = r*(r + 1)*(r + 2)/6;
+        basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
+      }// end loop over 'r'
+      for (unsigned int r = 0; r < 1; r++)
+      {
+        for (unsigned int s = 0; s < 1 - r; s++)
+        {
+          rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
+          ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
+          basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
+        }// end loop over 's'
+      }// end loop over 'r'
+      for (unsigned int r = 0; r < 2; r++)
+      {
+        for (unsigned int s = 0; s < 2 - r; s++)
+        {
+          for (unsigned int t = 0; t < 2 - r - s; t++)
+          {
+            rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
+            basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
+          }// end loop over 't'
+        }// end loop over 's'
+      }// end loop over 'r'
+      
+      // Table(s) of coefficients.
+      static const double coefficients0[4] = \
+      {0.288675134594813, 0.000000000000000, 0.000000000000000, 0.223606797749979};
+      
+      // Tables of derivatives of the polynomial base (transpose).
+      static const double dmats0[4][4] = \
+      {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
+      
+      static const double dmats1[4][4] = \
+      {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
+      
+      static const double dmats2[4][4] = \
+      {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
+      
+      // Compute reference derivatives.
+      // Declare pointer to array of derivatives on FIAT element.
+      double *derivatives = new double[num_derivatives];
+      for (unsigned int r = 0; r < num_derivatives; r++)
+      {
+        derivatives[r] = 0.000000000000000;
+      }// end loop over 'r'
+      
+      // Declare derivative matrix (of polynomial basis).
+      double dmats[4][4] = \
+      {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}};
+      
+      // Declare (auxiliary) derivative matrix (of polynomial basis).
+      double dmats_old[4][4] = \
+      {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000},
+      {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}};
+      
+      // Loop possible derivatives.
+      for (unsigned int r = 0; r < num_derivatives; r++)
+      {
+        // Resetting dmats values to compute next derivative.
+        for (unsigned int t = 0; t < 4; t++)
+        {
+          for (unsigned int u = 0; u < 4; u++)
+          {
+            dmats[t][u] = 0.000000000000000;
+            if (t == u)
+            {
+            dmats[t][u] = 1.000000000000000;
+            }
+            
+          }// end loop over 'u'
+        }// end loop over 't'
+        
+        // Looping derivative order to generate dmats.
+        for (unsigned int s = 0; s < n; s++)
+        {
+          // Updating dmats_old with new values and resetting dmats.
+          for (unsigned int t = 0; t < 4; t++)
+          {
+            for (unsigned int u = 0; u < 4; u++)
+            {
+              dmats_old[t][u] = dmats[t][u];
+              dmats[t][u] = 0.000000000000000;
+            }// end loop over 'u'
+          }// end loop over 't'
+          
+          // Update dmats using an inner product.
+          if (combinations[r][s] == 0)
+          {
+          for (unsigned int t = 0; t < 4; t++)
+          {
+            for (unsigned int u = 0; u < 4; u++)
+            {
+              for (unsigned int tu = 0; tu < 4; tu++)
+              {
+                dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
+              }// end loop over 'tu'
+            }// end loop over 'u'
+          }// end loop over 't'
+          }
+          
+          if (combinations[r][s] == 1)
+          {
+          for (unsigned int t = 0; t < 4; t++)
+          {
+            for (unsigned int u = 0; u < 4; u++)
+            {
+              for (unsigned int tu = 0; tu < 4; tu++)
+              {
+                dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
+              }// end loop over 'tu'
+            }// end loop over 'u'
+          }// end loop over 't'
+          }
+          
+          if (combinations[r][s] == 2)
+          {
+          for (unsigned int t = 0; t < 4; t++)
+          {
+            for (unsigned int u = 0; u < 4; u++)
+            {
+              for (unsigned int tu = 0; tu < 4; tu++)
+              {
+                dmats[t][u] += dmats2[t][tu]*dmats_old[tu][u];
+              }// end loop over 'tu'
+            }// end loop over 'u'
+          }// end loop over 't'
+          }
+          
+        }// end loop over 's'
+        for (unsigned int s = 0; s < 4; s++)
+        {
+          for (unsigned int t = 0; t < 4; t++)
           {
             derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
           }// end loop over 't'
@@ -837,7 +1271,7 @@
     unsigned int num_derivatives = 1;
     for (unsigned int r = 0; r < n; r++)
     {
-      num_derivatives *= 2;
+      num_derivatives *= 3;
     }// end loop over 'r'
     
     // Helper variable to hold values of a single dof.
@@ -848,7 +1282,7 @@
     }// end loop over 'r'
     
     // Loop dofs and call evaluate_basis_derivatives.
-    for (unsigned int r = 0; r < 3; r++)
+    for (unsigned int r = 0; r < 4; r++)
     {
       evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
       for (unsigned int s = 0; s < num_derivatives; s++)
@@ -870,7 +1304,7 @@
     double vals[1];
     
     // Declare variable for physical coordinates.
-    double y[2];
+    double y[3];
     const double * const * x = c.coordinates;
     switch (i)
     {
@@ -878,6 +1312,7 @@
       {
         y[0] = x[0][0];
       y[1] = x[0][1];
+      y[2] = x[0][2];
       f.evaluate(vals, y, c);
       return vals[0];
         break;
@@ -886,6 +1321,7 @@
       {
         y[0] = x[1][0];
       y[1] = x[1][1];
+      y[2] = x[1][2];
       f.evaluate(vals, y, c);
       return vals[0];
         break;
@@ -894,6 +1330,16 @@
       {
         y[0] = x[2][0];
       y[1] = x[2][1];
+      y[2] = x[2][2];
+      f.evaluate(vals, y, c);
+      return vals[0];
+        break;
+      }
+    case 3:
+      {
+        y[0] = x[3][0];
+      y[1] = x[3][1];
+      y[2] = x[3][2];
       f.evaluate(vals, y, c);
       return vals[0];
         break;
@@ -912,20 +1358,28 @@
     double vals[1];
     
     // Declare variable for physical coordinates.
-    double y[2];
+    double y[3];
     const double * const * x = c.coordinates;
     y[0] = x[0][0];
     y[1] = x[0][1];
+    y[2] = x[0][2];
     f.evaluate(vals, y, c);
     values[0] = vals[0];
     y[0] = x[1][0];
     y[1] = x[1][1];
+    y[2] = x[1][2];
     f.evaluate(vals, y, c);
     values[1] = vals[0];
     y[0] = x[2][0];
     y[1] = x[2][1];
+    y[2] = x[2][2];
     f.evaluate(vals, y, c);
     values[2] = vals[0];
+    y[0] = x[3][0];
+    y[1] = x[3][1];
+    y[2] = x[3][2];
+    f.evaluate(vals, y, c);
+    values[3] = vals[0];
   }
 
   /// Interpolate vertex values from dof values
@@ -937,6 +1391,7 @@
     vertex_values[0] = dof_values[0];
     vertex_values[1] = dof_values[1];
     vertex_values[2] = dof_values[2];
+    vertex_values[3] = dof_values[3];
   }
 
   /// Return the number of sub elements (for a mixed element)
@@ -978,7 +1433,7 @@
   /// Return a string identifying the dof map
   virtual const char* signature() const
   {
-    return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
+    return "FFC dofmap for FiniteElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1)";
   }
 
   /// Return true iff mesh entities of topological dimension d are needed
@@ -1001,6 +1456,11 @@
         return false;
         break;
       }
+    case 3:
+      {
+        return false;
+        break;
+      }
     }
     
     return false;
@@ -1035,25 +1495,25 @@
   /// Return the dimension of the local finite element function space for a cell
   virtual unsigned int local_dimension(const ufc::cell& c) const
   {
-    return 3;
+    return 4;
   }
 
   /// Return the maximum dimension of the local finite element function space
   virtual unsigned int max_local_dimension() const
   {
-    return 3;
+    return 4;
   }
 
   // Return the geometric dimension of the coordinates this dof map provides
   virtual unsigned int geometric_dimension() const
   {
-    return 2;
+    return 3;
   }
 
   /// Return the number of dofs on each cell facet
   virtual unsigned int num_facet_dofs() const
   {
-    return 2;
+    return 3;
   }
 
   /// Return the number of dofs associated with each cell entity of dimension d
@@ -1076,6 +1536,11 @@
         return 0;
         break;
       }
+    case 3:
+      {
+        return 0;
+        break;
+      }
     }
     
     return 0;
@@ -1089,6 +1554,7 @@
     dofs[0] = c.entity_indices[0][0];
     dofs[1] = c.entity_indices[0][1];
     dofs[2] = c.entity_indices[0][2];
+    dofs[3] = c.entity_indices[0][3];
   }
 
   /// Tabulate the local-to-local mapping from facet dofs to cell dofs
@@ -1101,18 +1567,28 @@
       {
         dofs[0] = 1;
       dofs[1] = 2;
+      dofs[2] = 3;
         break;
       }
     case 1:
       {
         dofs[0] = 0;
       dofs[1] = 2;
+      dofs[2] = 3;
         break;
       }
     case 2:
       {
         dofs[0] = 0;
       dofs[1] = 1;
+      dofs[2] = 3;
+        break;
+      }
+    case 3:
+      {
+        dofs[0] = 0;
+      dofs[1] = 1;
+      dofs[2] = 2;
         break;
       }
     }
@@ -1123,18 +1599,18 @@
   virtual void tabulate_entity_dofs(unsigned int* dofs,
                                     unsigned int d, unsigned int i) const
   {
-    if (d > 2)
+    if (d > 3)
     {
-    throw std::runtime_error("d is larger than dimension (2)");
+    throw std::runtime_error("d is larger than dimension (3)");
     }
     
     switch (d)
     {
     case 0:
       {
-        if (i > 2)
+        if (i > 3)
       {
-      throw std::runtime_error("i is larger than number of entities (2)");
+      throw std::runtime_error("i is larger than number of entities (3)");
       }
       
       switch (i)
@@ -1154,6 +1630,11 @@
           dofs[0] = 2;
           break;
         }
+      case 3:
+        {
+          dofs[0] = 3;
+          break;
+        }
       }
       
         break;
@@ -1168,6 +1649,11 @@
         
         break;
       }
+    case 3:
+      {
+        
+        break;
+      }
     }
     
   }
@@ -1180,10 +1666,16 @@
     
     coordinates[0][0] = x[0][0];
     coordinates[0][1] = x[0][1];
+    coordinates[0][2] = x[0][2];
     coordinates[1][0] = x[1][0];
     coordinates[1][1] = x[1][1];
+    coordinates[1][2] = x[1][2];
     coordinates[2][0] = x[2][0];
     coordinates[2][1] = x[2][1];
+    coordinates[2][2] = x[2][2];
+    coordinates[3][0] = x[3][0];
+    coordinates[3][1] = x[3][1];
+    coordinates[3][2] = x[3][2];
   }
 
   /// Return the number of sub dof maps (for a mixed element)
@@ -1225,10 +1717,10 @@
                                const double * const * w,
                                const ufc::cell& c) const
   {
-    // Number of operations (multiply-add pairs) for Jacobian data:      11
-    // Number of operations (multiply-add pairs) for geometry tensor:    8
-    // Number of operations (multiply-add pairs) for tensor contraction: 11
-    // Total number of operations (multiply-add pairs):                  30
+    // Number of operations (multiply-add pairs) for Jacobian data:      32
+    // Number of operations (multiply-add pairs) for geometry tensor:    27
+    // Number of operations (multiply-add pairs) for tensor contraction: 28
+    // Total number of operations (multiply-add pairs):                  87
     
     // Extract vertex coordinates
     const double * const * x = c.coordinates;
@@ -1236,37 +1728,70 @@
     // Compute Jacobian of affine map from reference cell
     const double J_00 = x[1][0] - x[0][0];
     const double J_01 = x[2][0] - x[0][0];
+    const double J_02 = x[3][0] - x[0][0];
     const double J_10 = x[1][1] - x[0][1];
     const double J_11 = x[2][1] - x[0][1];
+    const double J_12 = x[3][1] - x[0][1];
+    const double J_20 = x[1][2] - x[0][2];
+    const double J_21 = x[2][2] - x[0][2];
+    const double J_22 = x[3][2] - x[0][2];
+    
+    // Compute sub determinants
+    const double d_00 = J_11*J_22 - J_12*J_21;
+    const double d_01 = J_12*J_20 - J_10*J_22;
+    const double d_02 = J_10*J_21 - J_11*J_20;
+    const double d_10 = J_02*J_21 - J_01*J_22;
+    const double d_11 = J_00*J_22 - J_02*J_20;
+    const double d_12 = J_01*J_20 - J_00*J_21;
+    const double d_20 = J_01*J_12 - J_02*J_11;
+    const double d_21 = J_02*J_10 - J_00*J_12;
+    const double d_22 = J_00*J_11 - J_01*J_10;
     
     // Compute determinant of Jacobian
-    double detJ = J_00*J_11 - J_01*J_10;
+    double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
     
     // Compute inverse of Jacobian
-    const double K_00 =  J_11 / detJ;
-    const double K_01 = -J_01 / detJ;
-    const double K_10 = -J_10 / detJ;
-    const double K_11 =  J_00 / detJ;
+    const double K_00 = d_00 / detJ;
+    const double K_01 = d_10 / detJ;
+    const double K_02 = d_20 / detJ;
+    const double K_10 = d_01 / detJ;
+    const double K_11 = d_11 / detJ;
+    const double K_12 = d_21 / detJ;
+    const double K_20 = d_02 / detJ;
+    const double K_21 = d_12 / detJ;
+    const double K_22 = d_22 / detJ;
     
     // Set scale factor
     const double det = std::abs(detJ);
     
     // Compute geometry tensor
-    const double G0_0_0 = det*(K_00*K_00 + K_01*K_01);
-    const double G0_0_1 = det*(K_00*K_10 + K_01*K_11);
-    const double G0_1_0 = det*(K_10*K_00 + K_11*K_01);
-    const double G0_1_1 = det*(K_10*K_10 + K_11*K_11);
+    const double G0_0_0 = det*(K_00*K_00 + K_01*K_01 + K_02*K_02);
+    const double G0_0_1 = det*(K_00*K_10 + K_01*K_11 + K_02*K_12);
+    const double G0_0_2 = det*(K_00*K_20 + K_01*K_21 + K_02*K_22);
+    const double G0_1_0 = det*(K_10*K_00 + K_11*K_01 + K_12*K_02);
+    const double G0_1_1 = det*(K_10*K_10 + K_11*K_11 + K_12*K_12);
+    const double G0_1_2 = det*(K_10*K_20 + K_11*K_21 + K_12*K_22);
+    const double G0_2_0 = det*(K_20*K_00 + K_21*K_01 + K_22*K_02);
+    const double G0_2_1 = det*(K_20*K_10 + K_21*K_11 + K_22*K_12);
+    const double G0_2_2 = det*(K_20*K_20 + K_21*K_21 + K_22*K_22);
     
     // Compute element tensor
-    A[0] = 0.500000000000000*G0_0_0 + 0.500000000000000*G0_0_1 + 0.500000000000000*G0_1_0 + 0.500000000000000*G0_1_1;
-    A[1] = -0.500000000000000*G0_0_0 - 0.500000000000000*G0_1_0;
-    A[2] = -0.500000000000000*G0_0_1 - 0.500000000000000*G0_1_1;
-    A[3] = -0.500000000000000*G0_0_0 - 0.500000000000000*G0_0_1;
-    A[4] = 0.500000000000000*G0_0_0;
-    A[5] = 0.500000000000000*G0_0_1;
-    A[6] = -0.500000000000000*G0_1_0 - 0.500000000000000*G0_1_1;
-    A[7] = 0.500000000000000*G0_1_0;
-    A[8] = 0.500000000000000*G0_1_1;
+    A[0] = 0.166666666666666*G0_0_0 + 0.166666666666666*G0_0_1 + 0.166666666666666*G0_0_2 + 0.166666666666666*G0_1_0 + 0.166666666666666*G0_1_1 + 0.166666666666666*G0_1_2 + 0.166666666666666*G0_2_0 + 0.166666666666666*G0_2_1 + 0.166666666666666*G0_2_2;
+    A[1] = -0.166666666666666*G0_0_0 - 0.166666666666666*G0_1_0 - 0.166666666666666*G0_2_0;
+    A[2] = -0.166666666666666*G0_0_1 - 0.166666666666666*G0_1_1 - 0.166666666666666*G0_2_1;
+    A[3] = -0.166666666666666*G0_0_2 - 0.166666666666666*G0_1_2 - 0.166666666666666*G0_2_2;
+    A[4] = -0.166666666666666*G0_0_0 - 0.166666666666666*G0_0_1 - 0.166666666666666*G0_0_2;
+    A[5] = 0.166666666666666*G0_0_0;
+    A[6] = 0.166666666666666*G0_0_1;
+    A[7] = 0.166666666666666*G0_0_2;
+    A[8] = -0.166666666666666*G0_1_0 - 0.166666666666666*G0_1_1 - 0.166666666666666*G0_1_2;
+    A[9] = 0.166666666666666*G0_1_0;
+    A[10] = 0.166666666666667*G0_1_1;
+    A[11] = 0.166666666666666*G0_1_2;
+    A[12] = -0.166666666666666*G0_2_0 - 0.166666666666666*G0_2_1 - 0.166666666666666*G0_2_2;
+    A[13] = 0.166666666666666*G0_2_0;
+    A[14] = 0.166666666666666*G0_2_1;
+    A[15] = 0.166666666666666*G0_2_2;
   }
 
 };
@@ -1305,7 +1830,7 @@
   /// Return a string identifying the form
   virtual const char* signature() const
   {
-    return "Form([Integral(IndexSum(Product(Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 1), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2})), Measure('cell', 0, None))])";
+    return "Form([Integral(IndexSum(Product(Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1), 0), MultiIndex((Index(0),), {Index(0): 3})), MultiIndex((Index(0),), {Index(0): 3})), MultiIndex((Index(1),), {Index(1): 3})), Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1), 1), MultiIndex((Index(2),), {Index(2): 3})), MultiIndex((Index(2),), {Index(2): 3})), MultiIndex((Index(1),), {Index(1): 3}))), MultiIndex((Index(1),), {Index(1): 3})), Measure('cell', 0, None))])";
   }
 
   /// Return the rank of the global tensor (r)

=== modified file 'source/demos/la/eigenvalue/cpp/StiffnessMatrix.ufl'
--- source/demos/la/eigenvalue/cpp/StiffnessMatrix.ufl	2010-09-01 15:13:40 +0000
+++ source/demos/la/eigenvalue/cpp/StiffnessMatrix.ufl	2010-09-06 15:40:49 +0000
@@ -2,11 +2,11 @@
 # Licensed under the GNU LGPL Version 2.1
 #
 # First added:  2005-06-05
-# Last changed: 2006-03-28
+# Last changed: 2010-09-05
 #
 # The bilinear form for a stiffness matrix (Poisson).
 
-element = FiniteElement("Lagrange", "triangle", 1)
+element = FiniteElement("Lagrange", "tetrahedron", 1)
 
 v = TestFunction(element)
 u = TrialFunction(element)

=== modified file 'source/demos/la/eigenvalue/cpp/main.cpp'
--- source/demos/la/eigenvalue/cpp/main.cpp	2010-09-01 15:17:46 +0000
+++ source/demos/la/eigenvalue/cpp/main.cpp	2010-09-06 15:40:49 +0000
@@ -2,9 +2,10 @@
 // Licensed under the GNU LGPL Version 2.1.
 //
 // Modified by Anders Logg, 2008.
+// Modified by Marie E. Rognes, 2010.
 //
 // First added:  2007-03-08
-// Last changed: 2010-01-02
+// Last changed: 2010-09-05
 //
 // This simple program illustrates the use of the SLEPc eigenvalue solver.
 
@@ -21,34 +22,36 @@
   parameters["linear_algebra_backend"] = "PETSc";
 
   // Create mesh
-  UnitSquare mesh(64, 64);
+  Mesh mesh("box_with_dent.xml.gz");
 
-  // Build stiftness matrix
+  // Build stiffness matrix
   Matrix A;
   StiffnessMatrix::FunctionSpace V(mesh);
   StiffnessMatrix::BilinearForm a(V, V);
   assemble(A, a);
 
-  cout << A << endl;
-
   // Get PETSc matrix
   PETScMatrix& AA(A.down_cast<PETScMatrix>());
 
-  // Compute the first n eigenvalues
-  unsigned int n = 10;
+  // Create eigensolver
   SLEPcEigenSolver esolver;
-  esolver.parameters["spectrum"] = "largest magnitude";
-  esolver.solve(AA, n);
-
-  cout << "Solver converted in " << esolver.get_iteration_number() << " iterations" << endl;
-
-  // Display eigenvalues
-  for (unsigned int i = 0; i < n; i++)
-  {
-    double lr, lc;
-    esolver.get_eigenvalue(lr, lc, i);
-    cout << "Eigenvalue " << i << ": " << lr << endl;
-  }
+
+  // Compute all eigenvalues of A x = \lambda x
+  esolver.solve(AA);
+
+  // Extract largest (first) eigenpair
+  double r, c;
+  PETScVector rx(A.size(1));
+  PETScVector cx(A.size(1));
+  esolver.get_eigenpair(r, c, rx, cx, 0);
+
+  std::cout << "Largest eigenvalue: " << r << std::endl;
+
+  // Initialize function with eigenvector
+  Function u(V, rx);
+
+  // Plot eigenfunction
+  plot(u);
 
   #else

=== modified file 'source/demos/pde/biharmonic/python/demo.py'
--- source/demos/pde/biharmonic/python/demo.py	2010-09-01 20:32:22 +0000
+++ source/demos/pde/biharmonic/python/demo.py	2010-09-06 15:40:49 +0000
@@ -52,7 +52,7 @@
 h = CellSize(mesh)
 h_avg = (h('+') + h('-'))/2.0
 n = FacetNormal(mesh)
-f = Source(V)
+f = Source()
 
 # Penalty parameter
 alpha = Constant(8.0)

=== modified file 'source/demos/pde/biharmonic/python/documentation.rst'
--- source/demos/pde/biharmonic/python/documentation.rst	2010-09-01 22:59:53 +0000
+++ source/demos/pde/biharmonic/python/documentation.rst	2010-09-06 15:40:49 +0000
@@ -85,7 +85,7 @@
     h = CellSize(mesh)
     h_avg = (h('+') + h('-'))/2.0
     n = FacetNormal(mesh)
-    f = Source(V)
+    f = Source()
 
     # Penalty parameter
     alpha = Constant(8.0)

=== modified file 'source/demos/pde/cahn-hilliard/python/demo.py'
--- source/demos/pde/cahn-hilliard/python/demo.py	2010-09-01 20:32:22 +0000
+++ source/demos/pde/cahn-hilliard/python/demo.py	2010-09-06 15:40:49 +0000
@@ -18,8 +18,8 @@
     def eval(self, values, x):
         values[0] = 0.63 + 0.02*(0.5 - random.random())
         values[1] = 0.0
-    def dim(self):
-        return 2
+    def value_shape(self):
+        return (2,)
 
 # Class for interfacing with the Newton solver
 class CahnHilliardEquation(NonlinearProblem):

=== modified file 'source/demos/pde/cahn-hilliard/python/documentation.rst'
--- source/demos/pde/cahn-hilliard/python/documentation.rst	2010-09-01 20:10:55 +0000
+++ source/demos/pde/cahn-hilliard/python/documentation.rst	2010-09-06 15:40:49 +0000
@@ -34,8 +34,8 @@
         def eval(self, values, x):
             values[0] = 0.63 + 0.02*(0.5 - random.random())
             values[1] = 0.0
-        def dim(self):
-            return 2
+        def value_shape(self):
+            return (2,)
 
 It is a subclass of ``Expression``. In the constructor (``__init__``),
 the random number generator is seeded. If the program is run in
@@ -43,8 +43,8 @@
 number to ensure a different sequence of numbers on each process.  The
 function ``eval`` returns values for a function of dimension two.  For
 the first component of the function, a randomized value is returned.
-The function ``dim`` declares that the ``Expression`` is of dimension
-two.
+The method ``value_shape`` declares that the ``Expression`` is vector  
+valued with dimension two.
 
 .. index:: NonlinearProblem

=== modified file 'source/demos/pde/mixed-poisson/python/demo.py'
--- source/demos/pde/mixed-poisson/python/demo.py	2010-09-01 15:17:46 +0000
+++ source/demos/pde/mixed-poisson/python/demo.py	2010-09-06 15:40:49 +0000
@@ -53,10 +53,9 @@
         g = sin(5*data.x()[0])
         values[0] = g*data.normal()[0]
         values[1] = g*data.normal()[1]
-    def rank(self):
-        return 1
-    def dim(self):
-        return 2
+    def value_shape(self):
+        return (2,)
+
 G = BoundarySource()
 
 # Define essential boundary

=== modified file 'source/demos/pde/mixed-poisson/python/documentation.rst'
--- source/demos/pde/mixed-poisson/python/documentation.rst	2010-08-31 21:08:44 +0000
+++ source/demos/pde/mixed-poisson/python/documentation.rst	2010-09-06 15:40:49 +0000
@@ -104,7 +104,7 @@
 ``Expression`` class. Overloading the ``eval_data`` method (instead of
 the usual ``eval``) allows us to extract more geometry information
 such as the facet normals. Since this is a vector-valued expression,
-the methods ``rank`` and ``dim`` must also be specified.
+we need to overload the ``value_shape`` method.
 
 .. index:: Expression
 
@@ -116,10 +116,9 @@
             g = sin(5*data.x()[0])
             values[0] = g*data.normal()[0]
             values[1] = g*data.normal()[1]
-        def rank(self):
-            return 1
-        def dim(self):
-            return 2
+        def value_shape(self):
+            return (2,)
+
     G = BoundarySource()
 
 Specifying the relevant part of the boundary can be done as for the

=== modified file 'source/demos/undocumented/dielectric/python/demo.py'
--- source/demos/undocumented/dielectric/python/demo.py	2010-09-01 15:13:40 +0000
+++ source/demos/undocumented/dielectric/python/demo.py	2010-09-06 15:40:49 +0000
@@ -12,7 +12,7 @@
 #     u(x,y)      = V  for y = 0
 
 __author__ = "Kristen Kaasbjerg ([email-masked])"
-__date__ = "2008-02-14 -- 2009-10-07"
+__date__ = "2008-02-14 -- 2010-09-05"
 __copyright__ = ""
 __license__  = "GNU LGPL Version 2.1"
 
@@ -104,12 +104,12 @@
 v = TestFunction(V2)
 u = TrialFunction(V2)
 f = Constant(0.0)
-b = Coefficient(V0)
+b = Coefficient()
 a = b*dot(grad(v), grad(u))*dx
 L = v*f*dx
 
 # Define boundary condition
-u0 = DirichletFunction(V2)
+u0 = DirichletFunction()
 bc = DirichletBC(V2, u0, DirichletBoundary())
 
 # Solve PDE and plot solution
@@ -123,7 +123,7 @@
 # because it is an interpolation of the exact solution
 # in the finite element space!
 Pk = FunctionSpace(mesh, "CG", 5)
-exact = Exact(Pk)
+exact = Exact()
 
 e = u - exact
 L2_norm = e*e*dx

=== modified file 'source/demos/undocumented/elasticity/python/demo.py'
--- source/demos/undocumented/elasticity/python/demo.py	2010-09-01 22:59:53 +0000
+++ source/demos/undocumented/elasticity/python/demo.py	2010-09-06 15:40:49 +0000
@@ -3,7 +3,7 @@
 for a gear clamped at two of its ends and twisted 30 degrees."""
 
 __author__ = "Kristian B. Oelgaard ([email-masked])"
-__date__ = "2007-11-14 -- 2009-10-07"
+__date__ = "2007-11-14 -- 2010-09-05"
 __copyright__ = "Copyright (C) 2007 Kristian B. Oelgaard"
 __license__  = "GNU LGPL Version 2.1"
 
@@ -45,8 +45,8 @@
         values[1] = y - x[1]
         values[2] = z - x[2]
 
-    def dim(self):
-        return 3
+    def value_shape(self):
+        return (3,)
 
 # Sub domain for rotation at right end
 def right(x, on_boundary):

=== modified file 'source/demos/undocumented/elastodynamics/python/demo.py'
--- source/demos/undocumented/elastodynamics/python/demo.py	2010-09-01 22:59:53 +0000
+++ source/demos/undocumented/elastodynamics/python/demo.py	2010-09-06 15:40:49 +0000
@@ -41,23 +41,22 @@
         self.t   = t
         self.dt  = dt
         self.old = old
+
     def eval(self, values, x):
 
         # 'Shift' time for n-1 values
-        t_tmp = float(self.t)
+        t_tmp = self.t
         if self.old and t > 0.0:
             t_tmp -= self.dt
 
         cutoff_t = 10.0*1.0/32.0;
-        if t_tmp < cutoff_t:
-            values[0] = 1.0*t_tmp/cutoff_t
-            values[1] = 0.0
-        else:
-            values[0] = 1.0;
-            values[1] = 0.0;
-
-    def dim(self):
-        return 2
+        weight = t_tmp/cutoff_t if t_tmp < cutoff_t else 1.0
+
+        values[0] = 1.0*weight
+        values[1] = 0.0
+
+    def value_shape(self):
+        return (2,)
 
 # Sub domain for clamp at left end
 def left(x, on_boundary):
@@ -93,7 +92,7 @@
 beta    = 0.36
 gamma   = 0.7
 dt      = 1.0/32.0
-t       = Time(0.0)
+t       = 0.0
 T       = 10*dt
 
 # Some useful factors
@@ -144,11 +143,13 @@
 problem = VariationalProblem(a, L, bcs=bc, exterior_facet_domains=boundary_subdomains)
 
 vtk_file = File("elasticity.pvd")
-while t < T:
+while t <= T:
 
     t += dt
     print "Time: ", t
 
+    p.t = t
+    p0.t = t
     # Solve and update functions
     u = problem.solve()
     update(u, u0, v0, a0, beta, gamma, dt)

=== modified file 'source/demos/undocumented/periodic/python/demo.py'
--- source/demos/undocumented/periodic/python/demo.py	2010-09-01 15:13:40 +0000
+++ source/demos/undocumented/periodic/python/demo.py	2010-09-06 15:40:49 +0000
@@ -8,7 +8,7 @@
 # Original implementation: ../cpp/main.cpp by Anders Logg
 #
 __author__ = "Kristian B. Oelgaard ([email-masked])"
-__date__ = "2007-11-15 -- 2009-10-07"
+__date__ = "2007-11-15 -- 2010-09-05"
 __copyright__ = "Copyright (C) 2007 Kristian B. Oelgaard"
 __license__  = "GNU LGPL Version 2.1"
 
@@ -55,7 +55,7 @@
 # Define variational problem
 v = TestFunction(V)
 u = TrialFunction(V)
-f = Source(V)
+f = Source()
 
 a = dot(grad(v), grad(u))*dx
 L = v*f*dx

=== modified file 'source/demos/undocumented/sym-dirichlet-bc/python/demo.py'
--- source/demos/undocumented/sym-dirichlet-bc/python/demo.py	2010-09-01 15:13:40 +0000
+++ source/demos/undocumented/sym-dirichlet-bc/python/demo.py	2010-09-06 15:40:49 +0000
@@ -4,7 +4,7 @@
 # Modified by Kristian Oelgaard, 2008
 
 __author__ = "Kent-Andre Mardal ([email-masked])"
-__date__ = "2008-08-13 -- 2009-10-07"
+__date__ = "2008-08-13 -- 2010-09-05"
 __copyright__ = "Copyright (C) 2008 Kent-Andre Mardal"
 __license__  = "GNU LGPL Version 2.1"
 
@@ -37,8 +37,8 @@
 # Define variational problem
 v = TestFunction(V)
 u = TrialFunction(V)
-f = Source(V)
-g = Flux(V)
+f = Source()
+g = Flux()
 
 a = inner(grad(v), grad(u))*dx
 L = v*f*dx + v*g*ds

=== modified file 'source/programmers-reference/cpp/mesh/index.rst'
--- source/programmers-reference/cpp/mesh/index.rst	2010-08-30 13:52:07 +0000
+++ source/programmers-reference/cpp/mesh/index.rst	2010-09-06 15:40:49 +0000
@@ -34,7 +34,6 @@
     MeshPrimitive
     MeshTopology
     Point
-    Point
     PrimitiveIntersector
     PrimitiveTraits
     Rectangle

# Begin bundle

---

From: Anders Logg <[email-masked]>
Date: Mon, 6 Sep 2010 17:48:30 +0200

Patch applied

--
Anders

> IyBCYXphYXIgcmV2aXNpb24gYnVuZGxlIHY0CiMKQlpoOTFBWSZTWR7dD5kANcN/gF19GABY////
> f//ffv////pgKG4+oB94WB2jXzw1rnvSyqA33veFA9seQHbmus8Tru2vc9Sj733eCSS7BnsCtS0y
> UoXXd3F5Smlbu50Nr600KyW7jvb1EemvdtJUQmtOtUABVBQVRttbuEkiNVP1NpoyGp6NJtJ6I0ZG
> U8qeakGQaPSAAABieoJQEASE1MEIynlMhp6jTTQ9QAAAAAAAaaAhSTyn6U9QDQAAAAAAAAAAACTS
> SFMSaNAETQ0yaep6gGjR6g09QAGgAyaBoESU0CaTKYyVN4mjSYyaTap5NBGMp6h6INqDQ0ZANAVJ
> EBACAEaBNA0m1GE001TJ5T1NNqeUw1Gg9TNIyYSQgWfDlbKc/QGiFrN5CDabZIHV43EaLt6KLLmc
> oxXSs8jzYaTYfczju/HsfBwXfrTuZaf8cinY6ZufHmHDVOx4QWa5+tk9HE5CYxOHqTsdJpDhM1xt
> rVtWKB2J8X0vL5btZ4fJcy6+IfL184dc2bNfHPKSeDSqVqueNDdMZaFrrB/nKEw6Wb9PC7bGwUiz
> PHSGPK2Mh4xgogqoep5h4KkUXMsc6uzcslimW6llKXMYyFFFlz+eE0wW2sSnRZtGNjGLOy3nO+Sq
> iUor6iH8RQlEYQ3pSKK/2WU6KXRa9pz4YwbWLcvdZWwuxKMxh18uWPBpUGabIo81neUosyYQrrKN
> opvULKWYu1TDAZCtYk1F4Uj5GxcYxklIs2rDMlzMUuxKGAuLC/CqnrtgtXs88e4MqHpkFobjpjSV
> CVKtyd6r8DjtfaS7fty8N9JoApfDOAIJIakKXLqthRXZ1KAl3785vYlUb3aLsej5MqkFnLUuPRQg
> YpezDlbVsSzDJRUQiZl4tDLE3vO8RWk3HeWacMFjAWUlqITe9LUs9VvCqyOqYwMSEYZyXbEVNKKA
> WEAQYgqD/SCIo8UBFRwB3wERkQDxqEItQip+iWRhUpS/2nVPt9eoRePsFfJXzv8kaJvtyyPIxBoG
> LNqEmT/AXeQ7Q69vba1rW0GglDShcULiwobxQsLChYZjMXGAoYDEf3hpDCVKlSp9Yf3ihkNpLDsF
> BQ0ETaMQm0biVIc+aOmiXCBhMJhLLLJEDcDUE5ec73IJuGGRsNQ6huC4OUe1X7RcBZIRSAKSLCCI
> CJFgCkUWQWEFgHUdURJAgkQCDPJ6bbbbaLW21qtfW0TQFYCgEnwwgTp1vg7vzl1y9vSr29rC2eda
> nlxEeodNl+XFSwFShCHsZQBIgdvLpx1dEh3VTO3MMTZNIaRJAkNOiZFFEkyZmYV0ybYrUvD5YCu8
> kxN3e6irIovhRi0lqxVpCzTML2c4GGsj0bKuJw73e+FxFKDVtTfIDmCbRE1jZptvkjTxltuyrTAJ
> E+vYbAEj5B5A7HzkNj44bdOGcZVVoekHRaqOrKZVN8g6BrNTXM5gX3thmN9Iy7HJVoyXXKkLCnUD
> qcxOiJgL0rm/KCNNziu6zRXZ6RNMzLOpHsbp0u9Dkwbtmu5O0gB1wm3HPi9s4fUizLUaFUK6YKLK
> yzFVZSqvdhLMxdZMreREJN1gqUW8sEqrsrWxar3u1wrPjBeIl3GvecRcxFybYhgwoXSlL4uxiRYr
> dQIth55DyJ2ixZFFhP6A/0uwkrDxvbs/yP/D2vcyE0HskjcKsXgVIag3Asj92Fj7+F3P1t7c3NNM
> aesMgXyRAsNhvJbbbbcnQh4OcwDeM3jKbyzCyxggli7y8mryu88Ywgwg+MPDZBlpX4eFu5J5MELY
> 43uQwYwY3tNDBj1MabmNNNzBykLN6WWhZaWWhZaWX1cQ55zygMo86fbtMSrw1opqKUJShUXiSuFO
> 92zKvxKf0TE5qGgaNktV+tnkXGK3Y1NfTCHg7fRenI3iMRiM7xSylLKUspS0sstLOn87z63U5OP6
> 6jE5bIWnWJ65t/zwbFEQ3EEQQQSlhf1dAAnf+Q8ed3i+Xn8VV+ARy2OGkk9uyf6+nDyW8Lxnt9C9
> 1oST1+HRlvxCL5nSsvQ+Ekv3MVmA29mnBzKtKnGMbx1HDV1BHTPkcZJfLJMzNm9APYPwiTl1xJgO
> oYRJrkSm9vDM19TUI0kMJCwLSGg13eWRw3+TIb+Es2ibQX0BcS4KGzHf1cUks6duu3q1bhsNDRTL
> wxt5suBNhKJ20RYi11yLar2zah1dwwKCq48NCQu/37uRTSQiG4aGDBX8RMOAcgbeHHud1v7tne4g
> /v28bXckbSM/yDnwc5rOjTm9SqqqqEl1NryU0ikeKIm28d+y3Z20vmmZt7x/599yB4gQhtMyfsdF
> QUkBZAv3uXy/Cqr7JIsiy1lo+5SKQ79huP+m6NpJKKIz8dixloWH3SVrVsCjQfxF10sCjWLIlChS
> 5V1p4MmUmls4fbsb7bIkwChQKFAoKGFCwoKDkEQ3QCSnzTkeIfZMTv3MCW1ZUhfNaRHisxLxG9SH
> ll+oZpmngnqcO+Te+nJvT9mBUaVkMxk5dQwFxFChgKChWRVNtRbnJdabGSaLhvFBgI+lMbzLCq+b
> wv9b9eNWrtq49tiv6xxkNg5jrFw2MGiyyxZZ7B+Q64ZBoKEmgochQ8kZyHzUG2bJtlpaWlNMYNMY
> NMYxjTTTTTTTrAMyyk4KQ/LmQJAD2ACTrDudqqKdlKqirWXhLCyoVdSFIfTlI8SSdfwOD6T0uvv/
> byOozKFihuLFiZlClBFKQSlGzSMy4nGxfLrqqqwtY9TAXfXXyEiVBfC8z0sALcDAhLTtE18Wubg5
> naNi8TG+YtGuNyjE62HePDObWKwYjnDgji4jIx0ZPiwRijWaRD9Wxa/KrWub1FuMh7nu9m0uZI88
> nU5c/FBgdbxdfFuvGPGd/vw7OCblPNKSlJ+ZS7+WW6uiFqFChQoUKDiLlxQShSEUKEiiilSQUqlK
> KQsOwtCN9E6zpliUMAvYt2zuPCSHbOMg9vofJj6p98j+o7hL6zwWWPi014ezKbp3/lIk94lBH1z9
> 66j+1hUBhdR5yrNriw2tX7Au8wAnvgCYqqJmEdONJOR6TPodhLGAUVrxmtg4jEAheCaWLlttuQsp
> MKU6ySbpICn0XzsnUlTTN02QhiQ0qLMT22QNnmmyE6kwhaNoNRS6Km0D04OLokKNbH1mq80w7nP0
> OCP3QeDc5FJYkYtttt7SBkyZM7ocRjjjlZOUlWf2iKQKFERQpIJ7pBEYosUgsgqwoHMCwO/x6OQi
> o8SN11Ilwio7gHYICMdHXiUwgKnGU4psiVBDkKQR3CBoGO7ddxgWFUC6gh1uckwnB74Utqi1MnAw
> j6MGLbqXK8hANC7KtRwkewZHBUYGUF4MnqeFcQBtDmC27wg0G2libDJ0Er8BWsm9lBzYcBShsLX7
> SpM9lZjiqBWoMCIFZwNkknCEHcm4IiLT5gD/quEVTBPtViSSVFSTscTk/D9/Bi7nVTmspv6ThfHY
> 7DogU8WPFq3ZWlF6im6IaY7LrooZTGWvvq5ONgE9YBAOe29ttu0LKTQGjDB4dXjDl0Y6iw1TM24v
> BLKVcWgzxDYs6OjoMMZRU85jHDZrIuUavcBhsIoBGABGAYGBkUhSINttttuBTDmYbAT5QFkkoFCh
> 0U04KTfFclviAO+tUNk5gvr+czIOxqkki37zdMubjBPAqyPLjTRRoiYFKRpmRjVQUD0owgjI6Cpr
> dk3o5klESULjQ4iF06QYUlGCEcUFB0UooL1akpy9NLdeWJ3HlanEcJyg+o/KdNOqhebFcZDe4r91
> jpIxjC6wvLdMFqJwoGR/CsfRQCw11NnRGB+1vc79um+9uF6dna4QJVKy3oLOGFpz0ts0WOeNp5ql
> leQqimBYhadYwZF5t2WvIdi921XgrYpxmZhP+ih0Xcl0O3LQE9Ha5U4AiLC/F5hjnxbyR1e5xb9z
> yLpaZRHlJPanfBhUSZSTNzZOflxYNris2HA7pZhZi2PO3OdJD3u+Tbhn1DaV30a7sLs7WKqGGViy
> ow6PiFMpimIhGuGBmzOSNdnbO9MzYqohitaaamtOeEsmEBKNG5ufEWC7a58H+tExHxlAJDLLOqqq
> qEKKSxIxbbS3odYbQ882QBO9GcGiEOY1H/euUjOoOTwOfruykqBOPU3LLGo6TW6B9dKwOiHa8rbo
> SIncRJWqdE7vObPRkyJo85XAQBrvW8c2aQLKBBsokyJyOskh6RgZKbAVrojRGMBsrRQ94gx6Iud1
> U4WKorg7gp5GNT2nQTg9AxcvVqo/fIXDyASyS0KyKqC9C3X4uLVqti4xchpA8BqRBWOuyXXR3AXX
> wx4OSzTr92aGfcN6XnkGB7h7kjCiJlJNdex4P4PJZzfLO7orF07zzeRyXZJbWuvtnj2oq6IyDUhB
> UqqvwwzyeNoz5XjA65TfLZ6RNFAURrpskD4ALmkLDWVZoDRyWgBjMp0Ot9y8Rg9HRzfbwdWboX3O
> vm1MvM5LTWlqK08ehbuam4+LSTFrYNxrdeYSsMq2Cr7sWSgjdA6mx0upXFj1YlchBRKRKurdgrLZ
> V5PWEgB26FwrZKSUW4cetwD0SoIEEi4KFXxkPBZtY73NxODVzYMmiD2Tmp8H4PlOtdg47O9YN7t0
> AuJAuIQ46WpC3FjMy7OY2sRJL0WOyz2mHs9r0tjyuHammR62vZ1wYQ8bCAa34cDQGUAUZNqNliy+
> JBsDkMZM4I00pbBAxHaWO4hK1bXr5hrHqSgcC8al+xTQoUE0GyhyAwHgekhN3wu8S/AaMhEEgvlN
> IjwQyi8DHcCykEnXwDFMUyBxxPI5iAcFYO7HVB0p3Mm5nRyzYJ1uZg4G8NIg9GFCTsaOSnFy3sHR
> i5U8JuZPBzdPq+ueydtbXDqnUhGYhBoQoCcii6iqB4msmAaRMZ7QXzCml7tM2hr2awOs1Q9YF48w
> CUSxbIXAp2ud9a1i/DwPdwEkB2uCjMactSanuMOXOEsBDiM54xKkBlMCAZCSYOZCTYnjzg1VSiC3
> DgOIDYYgNCNqWe1yDuHKGKsC95gkwBQGlUdWVMgNJYhHXJADUTo1SwmDRUPLjYo5s0UHJyfAcGCx
> YPMfJPhOwLFTKLy3uF4lEAGS62ksGdZlQIL8YAANJk3s2zXDqdSsKpiqRhq0mJwVve4h4hNr3LhA
> Duq3FS7jriLqXqpuQWemQvrqrUxdzEW3117XBfbVy0QqWTuITPbrjODGRwyHBoECfRqzaAgOaZEL
> RDucmnDduZcxxOK+Cz9LfxrRWzMVys28HBWTcpdg7erZJQdg8Ojo0aq6mcwcihs1wUG83q72bR4h
> U+YaQ6p7iIHpSaz4moPAhYKNpZDpVtpmJKvazjTd6Y2gPh6WZrCrSybEZeFjSCGgbL20yzdSChjs
> 7JPAvrkVrdHWsXSpBoG0KaReC1Je0AeQuJVM4O0vUGJ0AqF9Hk9GeOvDhqa465C+rsb+nW6OHLjr
> wMm/iLJxcE543y0yc2rZTV3GurC2Ba6s2s6+uDw4NHTs4KFkOwUMHQ4O85Ea5SdGZVae556zzasQ
> PjMUtao0w04xa9LUapaa3shMVUsBh9FS+TB+NGz7m3G+QrONysVX22y7nsuUZuajKo4KSczdeYGa
> LkCaFT0iAWStrEwSNxWaqotDBlgPOLysoFzsPFD9+K89tgi5kllzDKYcrsDtY7MM3Aw6miDjuz8f
> GzrbXWzpwe18jv3OrmZckNKAqrroats3iYlbTD2nFHmZsojCAeC5pCCBogYQBih2lJW+SuInyDFq
> kHxoNwLyEguD5FfPIUftKiNC7ZTlssGQuB3FQ5nIx7q0DgqWiPfSpeI4HOAckg4M7gtX0ZMYEwTY
> 1F3e7+bwZOpwU0aN3q6Ox8sz7HTkt0QbQ4j3YL8oSGenwXb+c3HbZe5ycFA3g7dvIW4eOiQ7ceqX
> vrFSdphYMLQVpmtkFhtQ0trbW0Ed+hqCkWcdV0U4HqaNvIYmZErUCCASZUDreaIhYpYdH4wMzMzD
> Mw4YAVAkAhISWMaGmzYJfCBqE3wH2hllAwgX0YXBZ8smSzAV1Duy/Rwx3jXQvRcYbRcYe8bUsi/f
> WyJOacE0GwtYwkMHR5FlDHAFxLCxLEoiiP4+uCSlHqVkgA/WT+UQEKIBAj7iIJfBU/MMhEBJMAZL
> osAnX3/GViqiCsRFT2QtAoDG92SXFQDYh6/uySzzdxvD0KB6VL/b87976wCixrDerBJGT6rvuUJP
> O7jAnuCb/7GF/jyvV8j4cD4WvMa+X15UFiwUFigoKKDvPTJueoYxESx8eKjQiYsQgwYMQiEGAxgL
> Eid6ym03yYAxBiZMyYAxBiZMyYAxDl8ANk9oST1M4Q94QKbH4vNNeQU3uqwz+lZvFIUKRPgb0qNb
> SW/eqVg9CCkjomBSMLSXQgagYyhSMhqhdCxVGSGoGM1ZhouoKAoqp54+sZEYj4HvMC/tqqUyjftj
> T0shw+0wFgTUDKQK8ww5mEnMTYCgIjki/AB4MTEyKKKLewCWGw2GkpKU9jN8P0mehUKnt+P2+2uM
> 0qJYnAPIImACtBaBWUMpeWmQw5D3R9tg19NmwaZ1NirVKqyNspzhZLbbdjM4lzjghucQ4HeDMEwQ
> yDMEwYyzLcwy0ttywHFvMIfRmhnmI5gfUGOrVRRVVW9/yxIQ+cGmMCwRKEoWCWMpQlDWjaAW0pbb
> tNTaWFEpLGLGMsKJSUElhSxlKSxjrRtAKq3Jom0KWUssbKSxJYUspQoSwpYyjKRzTNisk8D2Gh/d
> 5E3YlSofdS1VIp7KQ+3kH7KywVvAiv5yPiqekLap4G43eBTMAMjA6TmOd2rrO1wz7WTzblM7tTPl
> Ifu6zexU3qatx2C3h8q0mMGG9Oe9azeq7iN4OwIoKDEbsZZKIP9j0WTF7ER0TcP9POQdyZejS568
> Pifp9NpEhc/SocgbfErpr0mhNdBTSQi3vcuAzwIGhnScoFHPef0xM3w83IUQZuvfvrynJh8Fy/MT
> Hk6EWxW21GRkTkIyjlSqiAJ+BAXz9sy6A3nAvnYW9hgMpmkIhjYVkhrOQpo6jR7nrYLMphIXfIpw
> avQ9/fE/gKH3/zH+I/kMAw2uDY7F25wc6dHJZxcmTKMV3N2u19WRPEhPmUg0KUnwcJ+v+rtk7urx
> ne9SEjcR17MDdLPQrDue+bHg72ayl3cMMVlRnKwmK7pnLhISjDSz/8VeGSk9a8EtL5S+Yu5zfF5z
> FdzfN4s7HaWp4M3e9fg9v7OfNwZt8w8JAC7pe87C1lKvV9fkVCnyodFPMgMAijEfINQZBgZWw6GG
> HmTTR3O0KzDjHHv8KzDjjYkk3aBHMp6NmN1XdtndNPF7/T8jPxZrbmr3qeS7hwOzL3qXg1er1PW5
> QU4Po9rTYbTtB4FXcHIL4HmTmIkIkKaBYwXLCUD9x5wYnMRzaMJUeKlDYPM7qSbpDtNbMNdA9tVV
> VVXbt8h14Bh4ihKB9XqNYakgeeQeiRPSruDni2iYLA0/1G161MBZgPmvxnKrOP19I6qYQSyWFpBa
> USwqhZZ1vRKT18zCyD0bn4PyfT2Ke5UP7rJpFE+SD1JREoCoB2RR+LI1lyYBATUzIiZlCPhC/tqO
> KSQVdBUD6DDF/ZdCrhgU1mRo9aavXzZOcgzknwc53LOHHfz+ExXWfSeUkKyLWQMIgIa03ZE9257n
> 3Ya5iQuDG3aXyUccWTtaOTx3YWknpknSY1DxotaO5RgxvJO0gkAykxiJxxVmzAeZozlMBAkKDES+
> 0RzWL+Cem2um0BJZIegsLSssIisVwykW2Wq1oYRZDmNSS5wp0ouajgwXTuXOR46Gi9jjNIYWKEWG
> pKgOgp/ceBp6kawncuEfP6T308VeP7FXHuppTXvXkl02ff2TNJjBsObCCPONWDabXymS304ZTE4u
> 3JO5QNM56lpn8hxu+WlKV2rRzudzCkTQMDkTIPmSWCAxiBijAlVBTP/Zo+Hc7m2WwWOtflDfUiXt
> JaXbkroZoYSDuZXQ0MITP02XXkk03xMW5rAE9jOK4KYh5XybSB09dhUe5orO2HOaBxKMdRgwfqfn
> 0b2rJ1+DXB6deT5rMVOe1fFsDi5dHtXmw5QCysAxVn2K86u1X5tj1Z5L3+6odAge9IHvAKCwRkOg
> R9MWIkOsPnwXVqRtHUvfDza3R7X/p2Wkqa+4/BlJ65Cmx95Z1B13HOZrR61Q+QKAVu4XpNmsdRq7
> DfJaymV6uD8X0BUk19m3Ny+zWajCPsppvTM/H9PYz17mbVJp0jQu732TyQR8vXnk7/QtJ3jzjTXd
> 2hvLAknj7cCld53+tVigWtyK4niQ6+rqOXWAYLjXGvYQ2gK86odnAr5oK/qWbGqjORtmk7ZiJ2k+
> KHj8td3p2PNuT9scAivUHT7fJ16stTJAkSSIECQA5XjPlg7u3fMz9KHBYD2x8yVt8K3T0T2iLiKk
> YPk0WqSVbtQepm+x+dg9k2vtWlO4vg+1B/Y7iR+XqQ/JI+pgnKnlaLWfbJDq6HvmybxDtDNvTyaj
> u71Nf4blO+6G3EA7EVaVkgsVhFdWj7khYffpd8R2mJtbk8XqLHDf07e6sAE/j0agwSTTJTlLlFPH
> 5ek3hjmFDwGcaFCGTzEUKuhYH4lQ83bOi9w4U5w54DZGyMYlpJ2Owp1jFSk02443Z2++ybZ2cNMn
> x+hn6G+EvJOG8MZD6DdrI+J2OD8XBuAu9qwZyxMtmTBNSRWH0BCbRkuiqj51nUqCpJhJDqnF9LZ2
> FwixyhgIfYiRjLWzvAVzrSHaYN/OuA5YfjUZXRrDudWGEbqfjI5pg1g0jiJJeSy4JaZCSVIM7GK8
> YXxl5JPybh+5TGJMr64E+VIwJGDGAtheDT26Eqh3xZg9DxLIGFQ9RU04qhsH/H8PnTX4h9BdLyQ0
> ol/DsG7YsIvxvNxx2sTiSXKQMyBdED/3o4DVaCoE+46BH7y4mnhPBCN5ey3V7et0CC9yqkzyLNhD
> EyAcD0Rdt70UGw9mPOBB/ab9sEb3bXPObXvshG5P68IVR3IK0jDabPJR8aklQXfyMw2js6G+Pztb
> P3IPvexh5lHapNrltXdsnl3Rtn2sXVsmDg2pNo+UT7pHCYTfGW2XOA+DljnVgigk+cI4igoUKChS
> Q5ElRJeMI7KWUKo/fJFjPEZiohVYx5T6vsmeJNg+soKie1ghSb/ZUW9tRtkzTb4nJNUSpFtR9O5o
> xQWd73T9n7VVLRVfo9X5Zhof4KNnw8L3EdUiuc+gpLSuiGH/Dm+RNlu1q8PfHWP3DtHQxE6OZzdC
> nvLUPiLuvtWNZ+aQw7NVlpj+YpkPaKEoKEoUJQtaQ+QlSFlhUwztu0pjfGkaGHv94ds2m0TEsveR
> M/g+eyNQgQE9cF95+Y679kAzA2lDxRJoUhUVRVUk1KMqqFQlqgoO4ZBS5C5Au8x0N+W8kA3mJdxD
> r5MQ1pGlds7yre/HIST+yR7lu5Xehk5EHoayTlB7epkY6cXLBctTiPWLnsGU7grtUnFXFqCW5RK9
> 5GsT+v41W9C8rTUDyYcq2vf83TATSgd1QPZno09DR8VAnVXi00EjbQMD3RFo7BWNKzx3CdkKBYYz
> EnqPzHwmAGmQLgcoLVK7BW/mmOOvWiYjCaR1RjFLKks17YtGWqjzzmAnh0PDq0ST27EtDR7CfV57
> eeklmlRFa8/uUcm97HM0eFVVKiqDYZLCPzvKdJI7ZmbUM5211pkehskYYVLFFFCilRE5fQvB0xMx
> O4UTyC6HzL49CWIolEZDQFOkakjJIDAPB2KG41nqD9IeVNgBoKRQydU8hB9EXVIKkiPvb/mn8P2S
> nzqfB7FoPDzI+pnwnSOP0ScyL1hFCBAB63Pr6yuxDqX+VJte3OlP1exMHukHPe2LqXpP+yvb/CzC
> nFPjRVH0fVVh91IjE3iOcUnKi+tb7TDoEUrBsySTJD3e/3PNg0MNUSrWZA61fUL1PSuXlfvCNFXq
> SC9MGopgwNQntgvD0STMwbMe0rEBCcOChX0wwMzsseBvgW+9dn5+5Ld4puWDAu8Xoerwr5FVSlSF
> UZ7xyCm891IxLXpU/mRsMaEaipwgOK6uopSU7YfOKkVPyUhSPlGB/iTwFhFggDS1CgSQwiYAYHjX
> +oK/WXqmABuCMk5AJwQ34X+9tlth9u0wEzKByiBH/i7kinChID26HzI=

> _______________________________________________
> Mailing list: https://launchpad.net/~fenics
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~fenics
> More help   : https://help.launchpad.net/ListHelp

--
Anders