~chaffra/ffc/main-old : contents of test/regression/references/r

~chaffra/ffc/main-old : (revision 1683)
// This code conforms with the UFC specification version 2.1.0+
// and was automatically generated by FFC version 1.1.0+.
// 
// This code was generated with the following parameters:
// 
//   cache_dir:                      ''
//   convert_exceptions_to_warnings: True
//   cpp_optimize:                   False
//   cpp_optimize_flags:             '-O2'
//   epsilon:                        1e-14
//   error_control:                  False
//   form_postfix:                   True
//   format:                         'ufc'
//   log_level:                      20
//   log_prefix:                     ''
//   optimize:                       False
//   output_dir:                     '.'
//   precision:                      '8'
//   quadrature_degree:              'auto'
//   quadrature_rule:                'auto'
//   representation:                 'quadrature'
//   split:                          False
//   swig_binary:                    'swig'
//   swig_path:                      ''

#ifndef __HEAT_H
#define __HEAT_H

#include <cmath>
#include <stdexcept>
#include <fstream>
#include <ufc.h>

/// This class defines the interface for a finite element.

class heat_finite_element_0: public ufc::finite_element
{
public:

  /// Constructor
  heat_finite_element_0() : ufc::finite_element()
  {
    // Do nothing
  }

  /// Destructor
  virtual ~heat_finite_element_0()
  {
    // Do nothing
  }

  /// Return a string identifying the finite element
  virtual const char* signature() const
  {
    return "FiniteElement('Real', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 0, None)";
  }

  /// Return the cell shape
  virtual ufc::shape cell_shape() const
  {
    return ufc::triangle;
  }

  /// Return the topological dimension of the cell shape
  virtual std::size_t topological_dimension() const
  {
    return 2;
  }

  /// Return the geometric dimension of the cell shape
  virtual std::size_t geometric_dimension() const
  {
    return 2;
  }

  /// Return the dimension of the finite element function space
  virtual std::size_t space_dimension() const
  {
    return 1;
  }

  /// Return the rank of the value space
  virtual std::size_t value_rank() const
  {
    return 0;
  }

  /// Return the dimension of the value space for axis i
  virtual std::size_t value_dimension(std::size_t i) const
  {
    return 1;
  }

  /// Evaluate basis function i at given point in cell
  virtual void evaluate_basis(std::size_t i,
                              double* values,
                              const double* coordinates,
                              const ufc::cell& c) const
  {
    // Extract vertex coordinates
    
    // Compute Jacobian of affine map from reference cell
    
    // Compute determinant of Jacobian
    
    // Compute inverse of Jacobian
    
    // Compute constants
    
    // Get coordinates and map to the reference (FIAT) element
    
    // Reset values.
    *values = 0.0;
    
    // Array of basisvalues.
    double basisvalues[1] = {0.0};
    
    // Declare helper variables.
    
    // Compute basisvalues.
    basisvalues[0] = 1.0;
    
    // Table(s) of coefficients.
    static const double coefficients0[1] = \
    {1.0};
    
    // Compute value(s).
    for (unsigned int r = 0; r < 1; r++)
    {
      *values += coefficients0[r]*basisvalues[r];
    }// end loop over 'r'
  }

  /// Evaluate all basis functions at given point in cell
  virtual void evaluate_basis_all(double* values,
                                  const double* coordinates,
                                  const ufc::cell& c) const
  {
    // Element is constant, calling evaluate_basis.
    evaluate_basis(0, values, coordinates, c);
  }

  /// Evaluate order n derivatives of basis function i at given point in cell
  virtual void evaluate_basis_derivatives(std::size_t i,
                                          std::size_t n,
                                          double* values,
                                          const double* coordinates,
                                          const ufc::cell& c) const
  {
    // Extract vertex coordinates
    const double * const * x = c.coordinates;
    
    // 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_10 = x[1][1] - x[0][1];
    const double J_11 = x[2][1] - x[0][1];
    
    // Compute determinant of Jacobian
    const double detJ = J_00*J_11 - J_01*J_10;
    
    // 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;
    
    // Compute constants
    
    // Get coordinates and map to the reference (FIAT) element
    
    // Compute number of derivatives.
    unsigned int num_derivatives = 1;
    for (unsigned int r = 0; r < n; r++)
    {
      num_derivatives *= 2;
    }// end loop over 'r'
    
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
    unsigned int **combinations = new unsigned int *[num_derivatives];
    for (unsigned int row = 0; row < num_derivatives; row++)
    {
      combinations[row] = new unsigned int [n];
      for (unsigned int col = 0; col < n; col++)
        combinations[row][col] = 0;
    }
    
    // Generate combinations of derivatives
    for (unsigned int row = 1; row < num_derivatives; row++)
    {
      for (unsigned int num = 0; num < row; num++)
      {
        for (unsigned int col = n-1; col+1 > 0; col--)
        {
          if (combinations[row][col] + 1 > 1)
            combinations[row][col] = 0;
          else
          {
            combinations[row][col] += 1;
            break;
          }
        }
      }
    }
    
    // Compute inverse of Jacobian
    const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
    
    // Declare transformation matrix
    // Declare pointer to two dimensional array and initialise
    double **transform = new double *[num_derivatives];
    
    for (unsigned int j = 0; j < num_derivatives; j++)
    {
      transform[j] = new double [num_derivatives];
      for (unsigned int k = 0; k < num_derivatives; k++)
        transform[j][k] = 1;
    }
    
    // Construct transformation matrix
    for (unsigned int row = 0; row < num_derivatives; row++)
    {
      for (unsigned int col = 0; col < num_derivatives; col++)
      {
        for (unsigned int k = 0; k < n; k++)
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
      }
    }
    
    // Reset values. Assuming that values is always an array.
    for (unsigned int r = 0; r < num_derivatives; r++)
    {
      values[r] = 0.0;
    }// end loop over 'r'
    
    
    // Array of basisvalues.
    double basisvalues[1] = {0.0};
    
    // Declare helper variables.
    
    // Compute basisvalues.
    basisvalues[0] = 1.0;
    
    // Table(s) of coefficients.
    static const double coefficients0[1] = \
    {1.0};
    
    // Tables of derivatives of the polynomial base (transpose).
    static const double dmats0[1][1] = \
    {{0.0}};
    
    static const double dmats1[1][1] = \
    {{0.0}};
    
    // 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.0;
    }// end loop over 'r'
    
    // Declare derivative matrix (of polynomial basis).
    double dmats[1][1] = \
    {{1.0}};
    
    // Declare (auxiliary) derivative matrix (of polynomial basis).
    double dmats_old[1][1] = \
    {{1.0}};
    
    // 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 < 1; t++)
      {
        for (unsigned int u = 0; u < 1; u++)
        {
          dmats[t][u] = 0.0;
          if (t == u)
          {
          dmats[t][u] = 1.0;
          }
          
        }// 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 < 1; t++)
        {
          for (unsigned int u = 0; u < 1; u++)
          {
            dmats_old[t][u] = dmats[t][u];
            dmats[t][u] = 0.0;
          }// 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 < 1; t++)
        {
          for (unsigned int u = 0; u < 1; u++)
          {
            for (unsigned int tu = 0; tu < 1; 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 < 1; t++)
        {
          for (unsigned int u = 0; u < 1; u++)
          {
            for (unsigned int tu = 0; tu < 1; 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 < 1; s++)
      {
        for (unsigned int t = 0; t < 1; 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;
  }

  /// Evaluate order n derivatives of all basis functions at given point in cell
  virtual void evaluate_basis_derivatives_all(std::size_t n,
                                              double* values,
                                              const double* coordinates,
                                              const ufc::cell& c) const
  {
    // Element is constant, calling evaluate_basis_derivatives.
    evaluate_basis_derivatives(0, n, values, coordinates, c);
  }

  /// Evaluate linear functional for dof i on the function f
  virtual double evaluate_dof(std::size_t i,
                              const ufc::function& f,
                              const ufc::cell& c) const
  {
    // Declare variables for result of evaluation.
    double vals[1];
    
    // Declare variable for physical coordinates.
    double y[2];
    const double * const * x = c.coordinates;
    switch (i)
    {
    case 0:
      {
        y[0] = 0.33333333*x[0][0] + 0.33333333*x[1][0] + 0.33333333*x[2][0];
      y[1] = 0.33333333*x[0][1] + 0.33333333*x[1][1] + 0.33333333*x[2][1];
      f.evaluate(vals, y, c);
      return vals[0];
        break;
      }
    }
    
    return 0.0;
  }

  /// Evaluate linear functionals for all dofs on the function f
  virtual void evaluate_dofs(double* values,
                             const ufc::function& f,
                             const ufc::cell& c) const
  {
    // Declare variables for result of evaluation.
    double vals[1];
    
    // Declare variable for physical coordinates.
    double y[2];
    const double * const * x = c.coordinates;
    y[0] = 0.33333333*x[0][0] + 0.33333333*x[1][0] + 0.33333333*x[2][0];
    y[1] = 0.33333333*x[0][1] + 0.33333333*x[1][1] + 0.33333333*x[2][1];
    f.evaluate(vals, y, c);
    values[0] = vals[0];
  }

  /// Interpolate vertex values from dof values
  virtual void interpolate_vertex_values(double* vertex_values,
                                         const double* dof_values,
                                         const ufc::cell& c) const
  {
    // Evaluate function and change variables
    vertex_values[0] = dof_values[0];
    vertex_values[1] = dof_values[0];
    vertex_values[2] = dof_values[0];
  }

  /// Map coordinate xhat from reference cell to coordinate x in cell
  virtual void map_from_reference_cell(double* x,
                                       const double* xhat,
                                       const ufc::cell& c) const
  {
    std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
  }

  /// Map from coordinate x in cell to coordinate xhat in reference cell
  virtual void map_to_reference_cell(double* xhat,
                                     const double* x,
                                     const ufc::cell& c) const
  {
    std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
  }

  /// Return the number of sub elements (for a mixed element)
  virtual std::size_t num_sub_elements() const
  {
    return 0;
  }

  /// Create a new finite element for sub element i (for a mixed element)
  virtual ufc::finite_element* create_sub_element(std::size_t i) const
  {
    return 0;
  }

  /// Create a new class instance
  virtual ufc::finite_element* create() const
  {
    return new heat_finite_element_0();
  }

};

/// This class defines the interface for a finite element.

class heat_finite_element_1: public ufc::finite_element
{
public:

  /// Constructor
  heat_finite_element_1() : ufc::finite_element()
  {
    // Do nothing
  }

  /// Destructor
  virtual ~heat_finite_element_1()
  {
    // Do nothing
  }

  /// Return a string identifying the finite element
  virtual const char* signature() const
  {
    return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
  }

  /// Return the cell shape
  virtual ufc::shape cell_shape() const
  {
    return ufc::triangle;
  }

  /// Return the topological dimension of the cell shape
  virtual std::size_t topological_dimension() const
  {
    return 2;
  }

  /// Return the geometric dimension of the cell shape
  virtual std::size_t geometric_dimension() const
  {
    return 2;
  }

  /// Return the dimension of the finite element function space
  virtual std::size_t space_dimension() const
  {
    return 3;
  }

  /// Return the rank of the value space
  virtual std::size_t value_rank() const
  {
    return 0;
  }

  /// Return the dimension of the value space for axis i
  virtual std::size_t value_dimension(std::size_t i) const
  {
    return 1;
  }

  /// Evaluate basis function i at given point in cell
  virtual void evaluate_basis(std::size_t i,
                              double* values,
                              const double* coordinates,
                              const ufc::cell& c) const
  {
    // Extract vertex coordinates
    const double * const * x = c.coordinates;
    
    // 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_10 = x[1][1] - x[0][1];
    const double J_11 = x[2][1] - x[0][1];
    
    // Compute determinant of Jacobian
    const double detJ = J_00*J_11 - J_01*J_10;
    
    // Compute inverse of Jacobian
    
    // Compute constants
    const double C0 = x[1][0] + x[2][0];
    const double C1 = x[1][1] + x[2][1];
    
    // 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;
    
    // Reset values.
    *values = 0.0;
    switch (i)
    {
    case 0:
      {
        
      // Array of basisvalues.
      double basisvalues[3] = {0.0, 0.0, 0.0};
      
      // Declare helper variables.
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
      
      // Compute basisvalues.
      basisvalues[0] = 1.0;
      basisvalues[1] = tmp0;
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
      basisvalues[0] *= std::sqrt(0.5);
      basisvalues[2] *= std::sqrt(1.0);
      basisvalues[1] *= std::sqrt(3.0);
      
      // Table(s) of coefficients.
      static const double coefficients0[3] = \
      {0.47140452, -0.28867513, -0.16666667};
      
      // Compute value(s).
      for (unsigned int r = 0; r < 3; r++)
      {
        *values += coefficients0[r]*basisvalues[r];
      }// end loop over 'r'
        break;
      }
    case 1:
      {
        
      // Array of basisvalues.
      double basisvalues[3] = {0.0, 0.0, 0.0};
      
      // Declare helper variables.
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
      
      // Compute basisvalues.
      basisvalues[0] = 1.0;
      basisvalues[1] = tmp0;
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
      basisvalues[0] *= std::sqrt(0.5);
      basisvalues[2] *= std::sqrt(1.0);
      basisvalues[1] *= std::sqrt(3.0);
      
      // Table(s) of coefficients.
      static const double coefficients0[3] = \
      {0.47140452, 0.28867513, -0.16666667};
      
      // Compute value(s).
      for (unsigned int r = 0; r < 3; r++)
      {
        *values += coefficients0[r]*basisvalues[r];
      }// end loop over 'r'
        break;
      }
    case 2:
      {
        
      // Array of basisvalues.
      double basisvalues[3] = {0.0, 0.0, 0.0};
      
      // Declare helper variables.
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
      
      // Compute basisvalues.
      basisvalues[0] = 1.0;
      basisvalues[1] = tmp0;
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
      basisvalues[0] *= std::sqrt(0.5);
      basisvalues[2] *= std::sqrt(1.0);
      basisvalues[1] *= std::sqrt(3.0);
      
      // Table(s) of coefficients.
      static const double coefficients0[3] = \
      {0.47140452, 0.0, 0.33333333};
      
      // Compute value(s).
      for (unsigned int r = 0; r < 3; r++)
      {
        *values += coefficients0[r]*basisvalues[r];
      }// end loop over 'r'
        break;
      }
    }
    
  }

  /// Evaluate all basis functions at given point in cell
  virtual void evaluate_basis_all(double* values,
                                  const double* coordinates,
                                  const ufc::cell& c) const
  {
    // Helper variable to hold values of a single dof.
    double dof_values = 0.0;
    
    // Loop dofs and call evaluate_basis.
    for (unsigned int r = 0; r < 3; r++)
    {
      evaluate_basis(r, &dof_values, coordinates, c);
      values[r] = dof_values;
    }// end loop over 'r'
  }

  /// Evaluate order n derivatives of basis function i at given point in cell
  virtual void evaluate_basis_derivatives(std::size_t i,
                                          std::size_t n,
                                          double* values,
                                          const double* coordinates,
                                          const ufc::cell& c) const
  {
    // Extract vertex coordinates
    const double * const * x = c.coordinates;
    
    // 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_10 = x[1][1] - x[0][1];
    const double J_11 = x[2][1] - x[0][1];
    
    // Compute determinant of Jacobian
    const double detJ = J_00*J_11 - J_01*J_10;
    
    // 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;
    
    // Compute constants
    const double C0 = x[1][0] + x[2][0];
    const double C1 = x[1][1] + x[2][1];
    
    // 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;
    
    // Compute number of derivatives.
    unsigned int num_derivatives = 1;
    for (unsigned int r = 0; r < n; r++)
    {
      num_derivatives *= 2;
    }// end loop over 'r'
    
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
    unsigned int **combinations = new unsigned int *[num_derivatives];
    for (unsigned int row = 0; row < num_derivatives; row++)
    {
      combinations[row] = new unsigned int [n];
      for (unsigned int col = 0; col < n; col++)
        combinations[row][col] = 0;
    }
    
    // Generate combinations of derivatives
    for (unsigned int row = 1; row < num_derivatives; row++)
    {
      for (unsigned int num = 0; num < row; num++)
      {
        for (unsigned int col = n-1; col+1 > 0; col--)
        {
          if (combinations[row][col] + 1 > 1)
            combinations[row][col] = 0;
          else
          {
            combinations[row][col] += 1;
            break;
          }
        }
      }
    }
    
    // Compute inverse of Jacobian
    const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
    
    // Declare transformation matrix
    // Declare pointer to two dimensional array and initialise
    double **transform = new double *[num_derivatives];
    
    for (unsigned int j = 0; j < num_derivatives; j++)
    {
      transform[j] = new double [num_derivatives];
      for (unsigned int k = 0; k < num_derivatives; k++)
        transform[j][k] = 1;
    }
    
    // Construct transformation matrix
    for (unsigned int row = 0; row < num_derivatives; row++)
    {
      for (unsigned int col = 0; col < num_derivatives; col++)
      {
        for (unsigned int k = 0; k < n; k++)
          transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
      }
    }
    
    // Reset values. Assuming that values is always an array.
    for (unsigned int r = 0; r < num_derivatives; r++)
    {
      values[r] = 0.0;
    }// end loop over 'r'
    
    switch (i)
    {
    case 0:
      {
        
      // Array of basisvalues.
      double basisvalues[3] = {0.0, 0.0, 0.0};
      
      // Declare helper variables.
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
      
      // Compute basisvalues.
      basisvalues[0] = 1.0;
      basisvalues[1] = tmp0;
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
      basisvalues[0] *= std::sqrt(0.5);
      basisvalues[2] *= std::sqrt(1.0);
      basisvalues[1] *= std::sqrt(3.0);
      
      // Table(s) of coefficients.
      static const double coefficients0[3] = \
      {0.47140452, -0.28867513, -0.16666667};
      
      // Tables of derivatives of the polynomial base (transpose).
      static const double dmats0[3][3] = \
      {{0.0, 0.0, 0.0},
      {4.8989795, 0.0, 0.0},
      {0.0, 0.0, 0.0}};
      
      static const double dmats1[3][3] = \
      {{0.0, 0.0, 0.0},
      {2.4494897, 0.0, 0.0},
      {4.2426407, 0.0, 0.0}};
      
      // 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.0;
      }// end loop over 'r'
      
      // Declare derivative matrix (of polynomial basis).
      double dmats[3][3] = \
      {{1.0, 0.0, 0.0},
      {0.0, 1.0, 0.0},
      {0.0, 0.0, 1.0}};
      
      // Declare (auxiliary) derivative matrix (of polynomial basis).
      double dmats_old[3][3] = \
      {{1.0, 0.0, 0.0},
      {0.0, 1.0, 0.0},
      {0.0, 0.0, 1.0}};
      
      // 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.0;
            if (t == u)
            {
            dmats[t][u] = 1.0;
            }
            
          }// 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.0;
            }// 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++)
          {
            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 1:
      {
        
      // Array of basisvalues.
      double basisvalues[3] = {0.0, 0.0, 0.0};
      
      // Declare helper variables.
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
      
      // Compute basisvalues.
      basisvalues[0] = 1.0;
      basisvalues[1] = tmp0;
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
      basisvalues[0] *= std::sqrt(0.5);
      basisvalues[2] *= std::sqrt(1.0);
      basisvalues[1] *= std::sqrt(3.0);
      
      // Table(s) of coefficients.
      static const double coefficients0[3] = \
      {0.47140452, 0.28867513, -0.16666667};
      
      // Tables of derivatives of the polynomial base (transpose).
      static const double dmats0[3][3] = \
      {{0.0, 0.0, 0.0},
      {4.8989795, 0.0, 0.0},
      {0.0, 0.0, 0.0}};
      
      static const double dmats1[3][3] = \
      {{0.0, 0.0, 0.0},
      {2.4494897, 0.0, 0.0},
      {4.2426407, 0.0, 0.0}};
      
      // 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.0;
      }// end loop over 'r'
      
      // Declare derivative matrix (of polynomial basis).
      double dmats[3][3] = \
      {{1.0, 0.0, 0.0},
      {0.0, 1.0, 0.0},
      {0.0, 0.0, 1.0}};
      
      // Declare (auxiliary) derivative matrix (of polynomial basis).
      double dmats_old[3][3] = \
      {{1.0, 0.0, 0.0},
      {0.0, 1.0, 0.0},
      {0.0, 0.0, 1.0}};
      
      // 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.0;
            if (t == u)
            {
            dmats[t][u] = 1.0;
            }
            
          }// 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.0;
            }// 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++)
          {
            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 2:
      {
        
      // Array of basisvalues.
      double basisvalues[3] = {0.0, 0.0, 0.0};
      
      // Declare helper variables.
      double tmp0 = (1.0 + Y + 2.0*X)/2.0;
      
      // Compute basisvalues.
      basisvalues[0] = 1.0;
      basisvalues[1] = tmp0;
      basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
      basisvalues[0] *= std::sqrt(0.5);
      basisvalues[2] *= std::sqrt(1.0);
      basisvalues[1] *= std::sqrt(3.0);
      
      // Table(s) of coefficients.
      static const double coefficients0[3] = \
      {0.47140452, 0.0, 0.33333333};
      
      // Tables of derivatives of the polynomial base (transpose).
      static const double dmats0[3][3] = \
      {{0.0, 0.0, 0.0},
      {4.8989795, 0.0, 0.0},
      {0.0, 0.0, 0.0}};
      
      static const double dmats1[3][3] = \
      {{0.0, 0.0, 0.0},
      {2.4494897, 0.0, 0.0},
      {4.2426407, 0.0, 0.0}};
      
      // 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.0;
      }// end loop over 'r'
      
      // Declare derivative matrix (of polynomial basis).
      double dmats[3][3] = \
      {{1.0, 0.0, 0.0},
      {0.0, 1.0, 0.0},
      {0.0, 0.0, 1.0}};
      
      // Declare (auxiliary) derivative matrix (of polynomial basis).
      double dmats_old[3][3] = \
      {{1.0, 0.0, 0.0},
      {0.0, 1.0, 0.0},
      {0.0, 0.0, 1.0}};
      
      // 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.0;
            if (t == u)
            {
            dmats[t][u] = 1.0;
            }
            
          }// 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.0;
            }// 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++)
          {
            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;
      }
    }
    
  }

  /// Evaluate order n derivatives of all basis functions at given point in cell
  virtual void evaluate_basis_derivatives_all(std::size_t n,
                                              double* values,
                                              const double* coordinates,
                                              const ufc::cell& c) const
  {
    // Compute number of derivatives.
    unsigned int num_derivatives = 1;
    for (unsigned int r = 0; r < n; r++)
    {
      num_derivatives *= 2;
    }// end loop over 'r'
    
    // Helper variable to hold values of a single dof.
    double *dof_values = new double[num_derivatives];
    for (unsigned int r = 0; r < num_derivatives; r++)
    {
      dof_values[r] = 0.0;
    }// end loop over 'r'
    
    // Loop dofs and call evaluate_basis_derivatives.
    for (unsigned int r = 0; r < 3; r++)
    {
      evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
      for (unsigned int s = 0; s < num_derivatives; s++)
      {
        values[r*num_derivatives + s] = dof_values[s];
      }// end loop over 's'
    }// end loop over 'r'
    
    // Delete pointer.
    delete [] dof_values;
  }

  /// Evaluate linear functional for dof i on the function f
  virtual double evaluate_dof(std::size_t i,
                              const ufc::function& f,
                              const ufc::cell& c) const
  {
    // Declare variables for result of evaluation.
    double vals[1];
    
    // Declare variable for physical coordinates.
    double y[2];
    const double * const * x = c.coordinates;
    switch (i)
    {
    case 0:
      {
        y[0] = x[0][0];
      y[1] = x[0][1];
      f.evaluate(vals, y, c);
      return vals[0];
        break;
      }
    case 1:
      {
        y[0] = x[1][0];
      y[1] = x[1][1];
      f.evaluate(vals, y, c);
      return vals[0];
        break;
      }
    case 2:
      {
        y[0] = x[2][0];
      y[1] = x[2][1];
      f.evaluate(vals, y, c);
      return vals[0];
        break;
      }
    }
    
    return 0.0;
  }

  /// Evaluate linear functionals for all dofs on the function f
  virtual void evaluate_dofs(double* values,
                             const ufc::function& f,
                             const ufc::cell& c) const
  {
    // Declare variables for result of evaluation.
    double vals[1];
    
    // Declare variable for physical coordinates.
    double y[2];
    const double * const * x = c.coordinates;
    y[0] = x[0][0];
    y[1] = x[0][1];
    f.evaluate(vals, y, c);
    values[0] = vals[0];
    y[0] = x[1][0];
    y[1] = x[1][1];
    f.evaluate(vals, y, c);
    values[1] = vals[0];
    y[0] = x[2][0];
    y[1] = x[2][1];
    f.evaluate(vals, y, c);
    values[2] = vals[0];
  }

  /// Interpolate vertex values from dof values
  virtual void interpolate_vertex_values(double* vertex_values,
                                         const double* dof_values,
                                         const ufc::cell& c) const
  {
    // Evaluate function and change variables
    vertex_values[0] = dof_values[0];
    vertex_values[1] = dof_values[1];
    vertex_values[2] = dof_values[2];
  }

  /// Map coordinate xhat from reference cell to coordinate x in cell
  virtual void map_from_reference_cell(double* x,
                                       const double* xhat,
                                       const ufc::cell& c) const
  {
    std::cerr << "*** FFC warning: " << "map_from_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
  }

  /// Map from coordinate x in cell to coordinate xhat in reference cell
  virtual void map_to_reference_cell(double* xhat,
                                     const double* x,
                                     const ufc::cell& c) const
  {
    std::cerr << "*** FFC warning: " << "map_to_reference_cell not yet implemented (introduced in UFC 2.0)." << std::endl;
  }

  /// Return the number of sub elements (for a mixed element)
  virtual std::size_t num_sub_elements() const
  {
    return 0;
  }

  /// Create a new finite element for sub element i (for a mixed element)
  virtual ufc::finite_element* create_sub_element(std::size_t i) const
  {
    return 0;
  }

  /// Create a new class instance
  virtual ufc::finite_element* create() const
  {
    return new heat_finite_element_1();
  }

};

/// This class defines the interface for a local-to-global mapping of
/// degrees of freedom (dofs).

class heat_dofmap_0: public ufc::dofmap
{
public:

  /// Constructor
  heat_dofmap_0() : ufc::dofmap()
  {
    // Do nothing
  }

  /// Destructor
  virtual ~heat_dofmap_0()
  {
    // Do nothing
  }

  /// Return a string identifying the dofmap
  virtual const char* signature() const
  {
    return "FFC dofmap for FiniteElement('Real', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 0, None)";
  }

  /// Return true iff mesh entities of topological dimension d are needed
  virtual bool needs_mesh_entities(std::size_t d) const
  {
    switch (d)
    {
    case 0:
      {
        return false;
        break;
      }
    case 1:
      {
        return false;
        break;
      }
    case 2:
      {
        return false;
        break;
      }
    }
    
    return false;
  }

  /// Return the topological dimension of the associated cell shape
  virtual std::size_t topological_dimension() const
  {
    return 2;
  }

  /// Return the geometric dimension of the associated cell shape
  virtual std::size_t geometric_dimension() const
  {
    return 2;
  }

  /// Return the dimension of the global finite element function space
  virtual std::size_t global_dimension(const std::vector<std::size_t>&
                                       num_global_entities) const
  {
    return 1;
  }

  /// Return the dimension of the local finite element function space for a cell
  virtual std::size_t local_dimension(const ufc::cell& c) const
  {
    return 1;
  }

  /// Return the maximum dimension of the local finite element function space
  virtual std::size_t max_local_dimension() const
  {
    return 1;
  }

  /// Return the number of dofs on each cell facet
  virtual std::size_t num_facet_dofs() const
  {
    return 0;
  }

  /// Return the number of dofs associated with each cell entity of dimension d
  virtual std::size_t num_entity_dofs(std::size_t d) const
  {
    switch (d)
    {
    case 0:
      {
        return 0;
        break;
      }
    case 1:
      {
        return 0;
        break;
      }
    case 2:
      {
        return 1;
        break;
      }
    }
    
    return 0;
  }

  /// Tabulate the local-to-global mapping of dofs on a cell
  virtual void tabulate_dofs(std::size_t* dofs,
                             const std::vector<std::size_t>& num_global_entities,
                             const ufc::cell& c) const
  {
    dofs[0] = 0;
  }

  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
  virtual void tabulate_facet_dofs(std::size_t* dofs,
                                   std::size_t facet) const
  {
    switch (facet)
    {
    case 0:
      {
        
        break;
      }
    case 1:
      {
        
        break;
      }
    case 2:
      {
        
        break;
      }
    }
    
  }

  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
  virtual void tabulate_entity_dofs(std::size_t* dofs,
                                    std::size_t d, std::size_t i) const
  {
    if (d > 2)
    {
    std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
    }
    
    switch (d)
    {
    case 0:
      {
        
        break;
      }
    case 1:
      {
        
        break;
      }
    case 2:
      {
        if (i > 0)
      {
      std::cerr << "*** FFC warning: " << "i is larger than number of entities (0)" << std::endl;
      }
      
      dofs[0] = 0;
        break;
      }
    }
    
  }

  /// Tabulate the coordinates of all dofs on a cell
  virtual void tabulate_coordinates(double** coordinates,
                                    const ufc::cell& c) const
  {
    const double * const * x = c.coordinates;
    
    coordinates[0][0] = 0.33333333*x[0][0] + 0.33333333*x[1][0] + 0.33333333*x[2][0];
    coordinates[0][1] = 0.33333333*x[0][1] + 0.33333333*x[1][1] + 0.33333333*x[2][1];
  }

  /// Return the number of sub dofmaps (for a mixed element)
  virtual std::size_t num_sub_dofmaps() const
  {
    return 0;
  }

  /// Create a new dofmap for sub dofmap i (for a mixed element)
  virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
  {
    return 0;
  }

  /// Create a new class instance
  virtual ufc::dofmap* create() const
  {
    return new heat_dofmap_0();
  }

};

/// This class defines the interface for a local-to-global mapping of
/// degrees of freedom (dofs).

class heat_dofmap_1: public ufc::dofmap
{
public:

  /// Constructor
  heat_dofmap_1() : ufc::dofmap()
  {
    // Do nothing
  }

  /// Destructor
  virtual ~heat_dofmap_1()
  {
    // Do nothing
  }

  /// Return a string identifying the dofmap
  virtual const char* signature() const
  {
    return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
  }

  /// Return true iff mesh entities of topological dimension d are needed
  virtual bool needs_mesh_entities(std::size_t d) const
  {
    switch (d)
    {
    case 0:
      {
        return true;
        break;
      }
    case 1:
      {
        return false;
        break;
      }
    case 2:
      {
        return false;
        break;
      }
    }
    
    return false;
  }

  /// Return the topological dimension of the associated cell shape
  virtual std::size_t topological_dimension() const
  {
    return 2;
  }

  /// Return the geometric dimension of the associated cell shape
  virtual std::size_t geometric_dimension() const
  {
    return 2;
  }

  /// Return the dimension of the global finite element function space
  virtual std::size_t global_dimension(const std::vector<std::size_t>&
                                       num_global_entities) const
  {
    return num_global_entities[0];
  }

  /// Return the dimension of the local finite element function space for a cell
  virtual std::size_t local_dimension(const ufc::cell& c) const
  {
    return 3;
  }

  /// Return the maximum dimension of the local finite element function space
  virtual std::size_t max_local_dimension() const
  {
    return 3;
  }

  /// Return the number of dofs on each cell facet
  virtual std::size_t num_facet_dofs() const
  {
    return 2;
  }

  /// Return the number of dofs associated with each cell entity of dimension d
  virtual std::size_t num_entity_dofs(std::size_t d) const
  {
    switch (d)
    {
    case 0:
      {
        return 1;
        break;
      }
    case 1:
      {
        return 0;
        break;
      }
    case 2:
      {
        return 0;
        break;
      }
    }
    
    return 0;
  }

  /// Tabulate the local-to-global mapping of dofs on a cell
  virtual void tabulate_dofs(std::size_t* dofs,
                             const std::vector<std::size_t>& num_global_entities,
                             const ufc::cell& c) const
  {
    dofs[0] = c.entity_indices[0][0];
    dofs[1] = c.entity_indices[0][1];
    dofs[2] = c.entity_indices[0][2];
  }

  /// Tabulate the local-to-local mapping from facet dofs to cell dofs
  virtual void tabulate_facet_dofs(std::size_t* dofs,
                                   std::size_t facet) const
  {
    switch (facet)
    {
    case 0:
      {
        dofs[0] = 1;
      dofs[1] = 2;
        break;
      }
    case 1:
      {
        dofs[0] = 0;
      dofs[1] = 2;
        break;
      }
    case 2:
      {
        dofs[0] = 0;
      dofs[1] = 1;
        break;
      }
    }
    
  }

  /// Tabulate the local-to-local mapping of dofs on entity (d, i)
  virtual void tabulate_entity_dofs(std::size_t* dofs,
                                    std::size_t d, std::size_t i) const
  {
    if (d > 2)
    {
    std::cerr << "*** FFC warning: " << "d is larger than dimension (2)" << std::endl;
    }
    
    switch (d)
    {
    case 0:
      {
        if (i > 2)
      {
      std::cerr << "*** FFC warning: " << "i is larger than number of entities (2)" << std::endl;
      }
      
      switch (i)
      {
      case 0:
        {
          dofs[0] = 0;
          break;
        }
      case 1:
        {
          dofs[0] = 1;
          break;
        }
      case 2:
        {
          dofs[0] = 2;
          break;
        }
      }
      
        break;
      }
    case 1:
      {
        
        break;
      }
    case 2:
      {
        
        break;
      }
    }
    
  }

  /// Tabulate the coordinates of all dofs on a cell
  virtual void tabulate_coordinates(double** coordinates,
                                    const ufc::cell& c) const
  {
    const double * const * x = c.coordinates;
    
    coordinates[0][0] = x[0][0];
    coordinates[0][1] = x[0][1];
    coordinates[1][0] = x[1][0];
    coordinates[1][1] = x[1][1];
    coordinates[2][0] = x[2][0];
    coordinates[2][1] = x[2][1];
  }

  /// Return the number of sub dofmaps (for a mixed element)
  virtual std::size_t num_sub_dofmaps() const
  {
    return 0;
  }

  /// Create a new dofmap for sub dofmap i (for a mixed element)
  virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
  {
    return 0;
  }

  /// Create a new class instance
  virtual ufc::dofmap* create() const
  {
    return new heat_dofmap_1();
  }

};

/// This class defines the interface for the tabulation of the cell
/// tensor corresponding to the local contribution to a form from
/// the integral over a cell.

class heat_cell_integral_0_0: public ufc::cell_integral
{
public:

  /// Constructor
  heat_cell_integral_0_0() : ufc::cell_integral()
  {
    // Do nothing
  }

  /// Destructor
  virtual ~heat_cell_integral_0_0()
  {
    // Do nothing
  }

  /// Tabulate the tensor for the contribution from a local cell
  virtual void tabulate_tensor(double* A,
                               const double * const * w,
                               const ufc::cell& c) const
  {
    // Extract vertex coordinates
    const double * const * x = c.coordinates;
    
    // 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_10 = x[1][1] - x[0][1];
    const double J_11 = x[2][1] - x[0][1];
    
    // Compute determinant of Jacobian
    const double detJ = J_00*J_11 - J_01*J_10;
    
    // 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;
    
    // Set scale factor
    const double det = std::abs(detJ);
    
    // Cell volume.
    
    // Compute circumradius of triangle in 2D.
    
    
    // Facet Area.
    
    // Array of quadrature weights.
    static const double W3[3] = {0.16666667, 0.16666667, 0.16666667};
    // Quadrature points on the UFC reference element: (0.16666667, 0.16666667), (0.16666667, 0.66666667), (0.66666667, 0.16666667)
    
    // Value of basis functions at quadrature points.
    static const double FE0[3][3] = \
    {{0.66666667, 0.16666667, 0.16666667},
    {0.16666667, 0.16666667, 0.66666667},
    {0.16666667, 0.66666667, 0.16666667}};
    
    static const double FE0_D01[3][3] = \
    {{-1.0, 0.0, 1.0},
    {-1.0, 0.0, 1.0},
    {-1.0, 0.0, 1.0}};
    
    static const double FE0_D10[3][3] = \
    {{-1.0, 1.0, 0.0},
    {-1.0, 1.0, 0.0},
    {-1.0, 1.0, 0.0}};
    
    // Reset values in the element tensor.
    for (unsigned int r = 0; r < 9; r++)
    {
      A[r] = 0.0;
    }// end loop over 'r'
    
    // Compute element tensor using UFL quadrature representation
    // Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
    
    // Loop quadrature points for integral.
    // Number of operations to compute element tensor for following IP loop = 612
    for (unsigned int ip = 0; ip < 3; ip++)
    {
      
      // Coefficient declarations.
      double F0 = 0.0;
      
      // Total number of operations to compute function values = 6
      for (unsigned int r = 0; r < 3; r++)
      {
        F0 += FE0[ip][r]*w[0][r];
      }// end loop over 'r'
      
      // Number of operations for primary indices: 198
      for (unsigned int j = 0; j < 3; j++)
      {
        for (unsigned int k = 0; k < 3; k++)
        {
          // Number of operations to compute entry: 22
          A[j*3 + k] += (FE0[ip][j]*FE0[ip][k] + ((((K_00*FE0_D10[ip][j] + K_10*FE0_D01[ip][j]))*((K_00*FE0_D10[ip][k] + K_10*FE0_D01[ip][k])) + ((K_01*FE0_D10[ip][j] + K_11*FE0_D01[ip][j]))*((K_01*FE0_D10[ip][k] + K_11*FE0_D01[ip][k]))))*F0*w[1][0])*W3[ip]*det;
        }// end loop over 'k'
      }// end loop over 'j'
    }// end loop over 'ip'
  }

  /// Tabulate the tensor for the contribution from a local cell
  /// using the specified reference cell quadrature points/weights
  virtual void tabulate_tensor(double* A,
                               const double * const * w,
                               const ufc::cell& c,
                               std::size_t num_quadrature_points,
                               const double * const * quadrature_points,
                               const double* quadrature_weights) const
  {
    std::cerr << "*** FFC warning: " << "Quadrature version of tabulate_tensor not yet implemented (introduced in UFC 2.0)." << std::endl;
  }

};

/// This class defines the interface for the tabulation of the cell
/// tensor corresponding to the local contribution to a form from
/// the integral over a cell.

class heat_cell_integral_1_0: public ufc::cell_integral
{
public:

  /// Constructor
  heat_cell_integral_1_0() : ufc::cell_integral()
  {
    // Do nothing
  }

  /// Destructor
  virtual ~heat_cell_integral_1_0()
  {
    // Do nothing
  }

  /// Tabulate the tensor for the contribution from a local cell
  virtual void tabulate_tensor(double* A,
                               const double * const * w,
                               const ufc::cell& c) const
  {
    // Extract vertex coordinates
    const double * const * x = c.coordinates;
    
    // 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_10 = x[1][1] - x[0][1];
    const double J_11 = x[2][1] - x[0][1];
    
    // Compute determinant of Jacobian
    const double detJ = J_00*J_11 - J_01*J_10;
    
    // Compute inverse of Jacobian
    
    // Set scale factor
    const double det = std::abs(detJ);
    
    // Cell volume.
    
    // Compute circumradius of triangle in 2D.
    
    
    // Facet Area.
    
    // Array of quadrature weights.
    static const double W3[3] = {0.16666667, 0.16666667, 0.16666667};
    // Quadrature points on the UFC reference element: (0.16666667, 0.16666667), (0.16666667, 0.66666667), (0.66666667, 0.16666667)
    
    // Value of basis functions at quadrature points.
    static const double FE0[3][3] = \
    {{0.66666667, 0.16666667, 0.16666667},
    {0.16666667, 0.16666667, 0.66666667},
    {0.16666667, 0.66666667, 0.16666667}};
    
    // Reset values in the element tensor.
    for (unsigned int r = 0; r < 3; r++)
    {
      A[r] = 0.0;
    }// end loop over 'r'
    
    // Compute element tensor using UFL quadrature representation
    // Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
    
    // Loop quadrature points for integral.
    // Number of operations to compute element tensor for following IP loop = 99
    for (unsigned int ip = 0; ip < 3; ip++)
    {
      
      // Coefficient declarations.
      double F0 = 0.0;
      double F1 = 0.0;
      
      // Total number of operations to compute function values = 12
      for (unsigned int r = 0; r < 3; r++)
      {
        F0 += FE0[ip][r]*w[0][r];
        F1 += FE0[ip][r]*w[1][r];
      }// end loop over 'r'
      
      // Number of operations for primary indices: 21
      for (unsigned int j = 0; j < 3; j++)
      {
        // Number of operations to compute entry: 7
        A[j] += (FE0[ip][j]*F0 + FE0[ip][j]*F1*w[2][0])*W3[ip]*det;
      }// end loop over 'j'
    }// end loop over 'ip'
  }

  /// Tabulate the tensor for the contribution from a local cell
  /// using the specified reference cell quadrature points/weights
  virtual void tabulate_tensor(double* A,
                               const double * const * w,
                               const ufc::cell& c,
                               std::size_t num_quadrature_points,
                               const double * const * quadrature_points,
                               const double* quadrature_weights) const
  {
    std::cerr << "*** FFC warning: " << "Quadrature version of tabulate_tensor not yet implemented (introduced in UFC 2.0)." << std::endl;
  }

};

/// This class defines the interface for the assembly of the global
/// tensor corresponding to a form with r + n arguments, that is, a
/// mapping
///
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
///
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
/// global tensor A is defined by
///
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
///
/// where each argument Vj represents the application to the
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
/// fixed functions (coefficients).

class heat_form_0: public ufc::form
{
public:

  /// Constructor
  heat_form_0() : ufc::form()
  {
    // Do nothing
  }

  /// Destructor
  virtual ~heat_form_0()
  {
    // Do nothing
  }

  /// Return a string identifying the form
  virtual const char* signature() const
  {
    return "beabc8cbfdbde2305229a06a22fc65169df459942208b7edbad5638470484817dd7b24adf340e04be94824a800269900212bbcd0e8e41c86261007b33e4bb8d6";
  }

  /// Return the rank of the global tensor (r)
  virtual std::size_t rank() const
  {
    return 2;
  }

  /// Return the number of coefficients (n)
  virtual std::size_t num_coefficients() const
  {
    return 2;
  }

  /// Return the number of cell domains
  virtual std::size_t num_cell_domains() const
  {
    return 1;
  }

  /// Return the number of exterior facet domains
  virtual std::size_t num_exterior_facet_domains() const
  {
    return 0;
  }

  /// Return the number of interior facet domains
  virtual std::size_t num_interior_facet_domains() const
  {
    return 0;
  }

  /// Return whether the form has any cell integrals
  virtual bool has_cell_integrals() const
  {
    return true;
  }

  /// Return whether the form has any exterior facet integrals
  virtual bool has_exterior_facet_integrals() const
  {
    return false;
  }

  /// Return whether the form has any interior facet integrals
  virtual bool has_interior_facet_integrals() const
  {
    return false;
  }

  /// Create a new finite element for argument function i
  virtual ufc::finite_element* create_finite_element(std::size_t i) const
  {
    switch (i)
    {
    case 0:
      {
        return new heat_finite_element_1();
        break;
      }
    case 1:
      {
        return new heat_finite_element_1();
        break;
      }
    case 2:
      {
        return new heat_finite_element_1();
        break;
      }
    case 3:
      {
        return new heat_finite_element_0();
        break;
      }
    }
    
    return 0;
  }

  /// Create a new dofmap for argument function i
  virtual ufc::dofmap* create_dofmap(std::size_t i) const
  {
    switch (i)
    {
    case 0:
      {
        return new heat_dofmap_1();
        break;
      }
    case 1:
      {
        return new heat_dofmap_1();
        break;
      }
    case 2:
      {
        return new heat_dofmap_1();
        break;
      }
    case 3:
      {
        return new heat_dofmap_0();
        break;
      }
    }
    
    return 0;
  }

  /// Create a new cell integral on sub domain i
  virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
  {
    switch (i)
    {
    case 0:
      {
        return new heat_cell_integral_0_0();
        break;
      }
    }
    
    return 0;
  }

  /// Create a new exterior facet integral on sub domain i
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
  {
    return 0;
  }

  /// Create a new interior facet integral on sub domain i
  virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
  {
    return 0;
  }

  /// Create a new cell integral on everywhere else
  virtual ufc::cell_integral* create_default_cell_integral() const
  {
    return 0;
  }

  /// Create a new exterior facet integral on everywhere else
  virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
  {
    return 0;
  }

  /// Create a new interior facet integral on everywhere else
  virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
  {
    return 0;
  }

};

/// This class defines the interface for the assembly of the global
/// tensor corresponding to a form with r + n arguments, that is, a
/// mapping
///
///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
///
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
/// global tensor A is defined by
///
///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
///
/// where each argument Vj represents the application to the
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
/// fixed functions (coefficients).

class heat_form_1: public ufc::form
{
public:

  /// Constructor
  heat_form_1() : ufc::form()
  {
    // Do nothing
  }

  /// Destructor
  virtual ~heat_form_1()
  {
    // Do nothing
  }

  /// Return a string identifying the form
  virtual const char* signature() const
  {
    return "b90e8813365e5fa06e1e0b75dbbe1a4a7729cc03a1dfd43caa3296bc2908db84d3aac4f9e1220007e11927ad07707d2e9266c2d8d4742a78eecc84fc6831ede1";
  }

  /// Return the rank of the global tensor (r)
  virtual std::size_t rank() const
  {
    return 1;
  }

  /// Return the number of coefficients (n)
  virtual std::size_t num_coefficients() const
  {
    return 3;
  }

  /// Return the number of cell domains
  virtual std::size_t num_cell_domains() const
  {
    return 1;
  }

  /// Return the number of exterior facet domains
  virtual std::size_t num_exterior_facet_domains() const
  {
    return 0;
  }

  /// Return the number of interior facet domains
  virtual std::size_t num_interior_facet_domains() const
  {
    return 0;
  }

  /// Return whether the form has any cell integrals
  virtual bool has_cell_integrals() const
  {
    return true;
  }

  /// Return whether the form has any exterior facet integrals
  virtual bool has_exterior_facet_integrals() const
  {
    return false;
  }

  /// Return whether the form has any interior facet integrals
  virtual bool has_interior_facet_integrals() const
  {
    return false;
  }

  /// Create a new finite element for argument function i
  virtual ufc::finite_element* create_finite_element(std::size_t i) const
  {
    switch (i)
    {
    case 0:
      {
        return new heat_finite_element_1();
        break;
      }
    case 1:
      {
        return new heat_finite_element_1();
        break;
      }
    case 2:
      {
        return new heat_finite_element_1();
        break;
      }
    case 3:
      {
        return new heat_finite_element_0();
        break;
      }
    }
    
    return 0;
  }

  /// Create a new dofmap for argument function i
  virtual ufc::dofmap* create_dofmap(std::size_t i) const
  {
    switch (i)
    {
    case 0:
      {
        return new heat_dofmap_1();
        break;
      }
    case 1:
      {
        return new heat_dofmap_1();
        break;
      }
    case 2:
      {
        return new heat_dofmap_1();
        break;
      }
    case 3:
      {
        return new heat_dofmap_0();
        break;
      }
    }
    
    return 0;
  }

  /// Create a new cell integral on sub domain i
  virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
  {
    switch (i)
    {
    case 0:
      {
        return new heat_cell_integral_1_0();
        break;
      }
    }
    
    return 0;
  }

  /// Create a new exterior facet integral on sub domain i
  virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
  {
    return 0;
  }

  /// Create a new interior facet integral on sub domain i
  virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
  {
    return 0;
  }

  /// Create a new cell integral on everywhere else
  virtual ufc::cell_integral* create_default_cell_integral() const
  {
    return 0;
  }

  /// Create a new exterior facet integral on everywhere else
  virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
  {
    return 0;
  }

  /// Create a new interior facet integral on everywhere else
  virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
  {
    return 0;
  }

};

#endif