// This code conforms with the UFC specification version 2.2.0+
// and was automatically generated by FFC version 1.2.0+.
// 
// This code was generated with the following parameters:
// 
//   cache_dir:                      ''
//   convert_exceptions_to_warnings: False
//   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:                      15
//   quadrature_degree:              'auto'
//   quadrature_rule:                'auto'
//   representation:                 'auto'
//   restrict_keyword:               ''
//   split:                          False
//   swig_binary:                    'swig'
//   swig_path:                      ''

#ifndef __LAPLACIAN_H
#define __LAPLACIAN_H

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

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

class laplacian_finite_element_0: public ufc::finite_element
{
public:

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

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

  /// Return a string identifying the finite element
  virtual const char* signature() const
  {
    return "FiniteElement('Lagrange', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 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 3;
  }

  /// 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 x in cell
  virtual void evaluate_basis(std::size_t i,
                              double* values,
                              const double* x,
                              const double* vertex_coordinates,
                              int cell_orientation) const
  {
    // Compute Jacobian
    double J[6];
    compute_jacobian_triangle_3d(J, vertex_coordinates);
    
    // Compute Jacobian inverse and determinant
    double K[6];
    double detJ;
    compute_jacobian_inverse_triangle_3d(K, detJ, J);
    
    
    const double b0 = vertex_coordinates[0];
    const double b1 = vertex_coordinates[1];
    const double b2 = vertex_coordinates[2];
    
    // P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
    double X = 2*(K[0]*(x[0] - b0) + K[1]*(x[1] - b1) + K[2]*(x[2] - b2)) - 1.0;
    double Y = 2*(K[3]*(x[0] - b0) + K[4]*(x[1] - b1) + K[5]*(x[2] - b2)) - 1.0;
    
    
    // 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.471404520791032, -0.288675134594813, -0.166666666666667};
      
      // 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.471404520791032, 0.288675134594813, -0.166666666666667};
      
      // 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.471404520791032, 0.0, 0.333333333333333};
      
      // 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 x in cell
  virtual void evaluate_basis_all(double* values,
                                  const double* x,
                                  const double* vertex_coordinates,
                                  int cell_orientation) 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, x, vertex_coordinates, cell_orientation);
      values[r] = dof_values;
    }// end loop over 'r'
  }

  /// Evaluate order n derivatives of basis function i at given point x in cell
  virtual void evaluate_basis_derivatives(std::size_t i,
                                          std::size_t n,
                                          double* values,
                                          const double* x,
                                          const double* vertex_coordinates,
                                          int cell_orientation) const
  {
    // Compute Jacobian
    double J[6];
    compute_jacobian_triangle_3d(J, vertex_coordinates);
    
    // Compute Jacobian inverse and determinant
    double K[6];
    double detJ;
    compute_jacobian_inverse_triangle_3d(K, detJ, J);
    
    
    const double b0 = vertex_coordinates[0];
    const double b1 = vertex_coordinates[1];
    const double b2 = vertex_coordinates[2];
    
    // P_FFC = J^dag (p - b), P_FIAT = 2*P_FFC - (1, 1)
    double X = 2*(K[0]*(x[0] - b0) + K[1]*(x[1] - b1) + K[2]*(x[2] - b2)) - 1.0;
    double Y = 2*(K[3]*(x[0] - b0) + K[4]*(x[1] - b1) + K[5]*(x[2] - b2)) - 1.0;
    
    
    // Compute number of derivatives.
    unsigned int num_derivatives_t = 1;
    for (unsigned int r = 0; r < n; r++)
    {
      num_derivatives_t *= 2;
    }// end loop over 'r'
    
    // Compute number of derivatives.
    unsigned int num_derivatives_g = 1;
    for (unsigned int r = 0; r < n; r++)
    {
      num_derivatives_g *= 3;
    }// end loop over 'r'
    
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
    unsigned int **combinations_t = new unsigned int *[num_derivatives_t];
    for (unsigned int row = 0; row < num_derivatives_t; row++)
    {
      combinations_t[row] = new unsigned int [n];
      for (unsigned int col = 0; col < n; col++)
        combinations_t[row][col] = 0;
    }
    
    // Generate combinations of derivatives
    for (unsigned int row = 1; row < num_derivatives_t; row++)
    {
      for (unsigned int num = 0; num < row; num++)
      {
        for (unsigned int col = n-1; col+1 > 0; col--)
        {
          if (combinations_t[row][col] + 1 > 1)
            combinations_t[row][col] = 0;
          else
          {
            combinations_t[row][col] += 1;
            break;
          }
        }
      }
    }
    
    // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
    unsigned int **combinations_g = new unsigned int *[num_derivatives_g];
    for (unsigned int row = 0; row < num_derivatives_g; row++)
    {
      combinations_g[row] = new unsigned int [n];
      for (unsigned int col = 0; col < n; col++)
        combinations_g[row][col] = 0;
    }
    
    // Generate combinations of derivatives
    for (unsigned int row = 1; row < num_derivatives_g; row++)
    {
      for (unsigned int num = 0; num < row; num++)
      {
        for (unsigned int col = n-1; col+1 > 0; col--)
        {
          if (combinations_g[row][col] + 1 > 2)
            combinations_g[row][col] = 0;
          else
          {
            combinations_g[row][col] += 1;
            break;
          }
        }
      }
    }
    
    // Compute inverse of Jacobian
    const double Jinv[2][3] = {{K[0], K[1], K[2]}, {K[3], K[4], K[5]}};
    
    // Declare transformation matrix
    // Declare pointer to two dimensional array and initialise
    double **transform = new double *[num_derivatives_g];
    
    for (unsigned int j = 0; j < num_derivatives_g; j++)
    {
      transform[j] = new double [num_derivatives_t];
      for (unsigned int k = 0; k < num_derivatives_t; k++)
        transform[j][k] = 1;
    }
    
    // Construct transformation matrix
    for (unsigned int row = 0; row < num_derivatives_g; row++)
    {
      for (unsigned int col = 0; col < num_derivatives_t; col++)
      {
        for (unsigned int k = 0; k < n; k++)
          transform[row][col] *= Jinv[combinations_t[col][k]][combinations_g[row][k]];
      }
    }
    
    // Reset values. Assuming that values is always an array.
    for (unsigned int r = 0; r < num_derivatives_g; 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.471404520791032, -0.288675134594813, -0.166666666666667};
      
      // Tables of derivatives of the polynomial base (transpose).
      static const double dmats0[3][3] = \
      {{0.0, 0.0, 0.0},
      {4.89897948556636, 0.0, 0.0},
      {0.0, 0.0, 0.0}};
      
      static const double dmats1[3][3] = \
      {{0.0, 0.0, 0.0},
      {2.44948974278318, 0.0, 0.0},
      {4.24264068711928, 0.0, 0.0}};
      
      // Compute reference derivatives.
      // Declare pointer to array of derivatives on FIAT element.
      double *derivatives = new double[num_derivatives_t];
      for (unsigned int r = 0; r < num_derivatives_t; 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_t; 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_t[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_t[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_g; r++)
      {
        for (unsigned int s = 0; s < num_derivatives_t; 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_t; r++)
      {
        delete [] combinations_t[r];
      }// end loop over 'r'
      delete [] combinations_t;
      for (unsigned int r = 0; r < num_derivatives_g; r++)
      {
        delete [] combinations_g[r];
      }// end loop over 'r'
      delete [] combinations_g;
      for (unsigned int r = 0; r < num_derivatives_g; 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.471404520791032, 0.288675134594813, -0.166666666666667};
      
      // Tables of derivatives of the polynomial base (transpose).
      static const double dmats0[3][3] = \
      {{0.0, 0.0, 0.0},
      {4.89897948556636, 0.0, 0.0},
      {0.0, 0.0, 0.0}};
      
      static const double dmats1[3][3] = \
      {{0.0, 0.0, 0.0},
      {2.44948974278318, 0.0, 0.0},
      {4.24264068711928, 0.0, 0.0}};
      
      // Compute reference derivatives.
      // Declare pointer to array of derivatives on FIAT element.
      double *derivatives = new double[num_derivatives_t];
      for (unsigned int r = 0; r < num_derivatives_t; 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_t; 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_t[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_t[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_g; r++)
      {
        for (unsigned int s = 0; s < num_derivatives_t; 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_t; r++)
      {
        delete [] combinations_t[r];
      }// end loop over 'r'
      delete [] combinations_t;
      for (unsigned int r = 0; r < num_derivatives_g; r++)
      {
        delete [] combinations_g[r];
      }// end loop over 'r'
      delete [] combinations_g;
      for (unsigned int r = 0; r < num_derivatives_g; 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.471404520791032, 0.0, 0.333333333333333};
      
      // Tables of derivatives of the polynomial base (transpose).
      static const double dmats0[3][3] = \
      {{0.0, 0.0, 0.0},
      {4.89897948556636, 0.0, 0.0},
      {0.0, 0.0, 0.0}};
      
      static const double dmats1[3][3] = \
      {{0.0, 0.0, 0.0},
      {2.44948974278318, 0.0, 0.0},
      {4.24264068711928, 0.0, 0.0}};
      
      // Compute reference derivatives.
      // Declare pointer to array of derivatives on FIAT element.
      double *derivatives = new double[num_derivatives_t];
      for (unsigned int r = 0; r < num_derivatives_t; 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_t; 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_t[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_t[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_g; r++)
      {
        for (unsigned int s = 0; s < num_derivatives_t; 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_t; r++)
      {
        delete [] combinations_t[r];
      }// end loop over 'r'
      delete [] combinations_t;
      for (unsigned int r = 0; r < num_derivatives_g; r++)
      {
        delete [] combinations_g[r];
      }// end loop over 'r'
      delete [] combinations_g;
      for (unsigned int r = 0; r < num_derivatives_g; r++)
      {
        delete [] transform[r];
      }// end loop over 'r'
      delete [] transform;
        break;
      }
    }
    
  }

  /// Evaluate order n derivatives of all basis functions at given point x in cell
  virtual void evaluate_basis_derivatives_all(std::size_t n,
                                              double* values,
                                              const double* x,
                                              const double* vertex_coordinates,
                                              int cell_orientation) const
  {
    // Compute number of derivatives.
    unsigned int num_derivatives_g = 1;
    for (unsigned int r = 0; r < n; r++)
    {
      num_derivatives_g *= 3;
    }// end loop over 'r'
    
    // Helper variable to hold values of a single dof.
    double *dof_values = new double[num_derivatives_g];
    for (unsigned int r = 0; r < num_derivatives_g; 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, x, vertex_coordinates, cell_orientation);
      for (unsigned int s = 0; s < num_derivatives_g; s++)
      {
        values[r*num_derivatives_g + 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 double* vertex_coordinates,
                              int cell_orientation,
                              const ufc::cell& c) const
  {
    // Declare variables for result of evaluation
    double vals[1];
    
    // Declare variable for physical coordinates
    double y[3];
    switch (i)
    {
    case 0:
      {
        y[0] = vertex_coordinates[0];
      y[1] = vertex_coordinates[1];
      y[2] = vertex_coordinates[2];
      f.evaluate(vals, y, c);
      return vals[0];
        break;
      }
    case 1:
      {
        y[0] = vertex_coordinates[3];
      y[1] = vertex_coordinates[4];
      y[2] = vertex_coordinates[5];
      f.evaluate(vals, y, c);
      return vals[0];
        break;
      }
    case 2:
      {
        y[0] = vertex_coordinates[6];
      y[1] = vertex_coordinates[7];
      y[2] = vertex_coordinates[8];
      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 double* vertex_coordinates,
                             int cell_orientation,
                             const ufc::cell& c) const
  {
    // Declare variables for result of evaluation
    double vals[1];
    
    // Declare variable for physical coordinates
    double y[3];
    y[0] = vertex_coordinates[0];
    y[1] = vertex_coordinates[1];
    y[2] = vertex_coordinates[2];
    f.evaluate(vals, y, c);
    values[0] = vals[0];
    y[0] = vertex_coordinates[3];
    y[1] = vertex_coordinates[4];
    y[2] = vertex_coordinates[5];
    f.evaluate(vals, y, c);
    values[1] = vals[0];
    y[0] = vertex_coordinates[6];
    y[1] = vertex_coordinates[7];
    y[2] = vertex_coordinates[8];
    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 double* vertex_coordinates,
                                         int cell_orientation,
                                         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
  {
    throw std::runtime_error("map_from_reference_cell not yet implemented.");
  }

  /// 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
  {
    throw std::runtime_error("map_to_reference_cell not yet implemented.");
  }

  /// 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 laplacian_finite_element_0();
  }

};

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

class laplacian_dofmap_0: public ufc::dofmap
{
public:

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

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

  /// Return a string identifying the dofmap
  virtual const char* signature() const
  {
    return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 3), 'triangle_multiverse', 3, 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 3;
  }

  /// 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
  {
    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)
    {
    throw std::runtime_error("d is larger than dimension (2)");
    }
    
    switch (d)
    {
    case 0:
      {
        if (i > 2)
      {
      throw std::runtime_error("i is larger than number of entities (2)");
      }
      
      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** dof_coordinates,
                                    const double* vertex_coordinates) const
  {
    dof_coordinates[0][0] = vertex_coordinates[0];
    dof_coordinates[0][1] = vertex_coordinates[1];
    dof_coordinates[0][2] = vertex_coordinates[2];
    dof_coordinates[1][0] = vertex_coordinates[3];
    dof_coordinates[1][1] = vertex_coordinates[4];
    dof_coordinates[1][2] = vertex_coordinates[5];
    dof_coordinates[2][0] = vertex_coordinates[6];
    dof_coordinates[2][1] = vertex_coordinates[7];
    dof_coordinates[2][2] = vertex_coordinates[8];
  }

  /// 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 laplacian_dofmap_0();
  }

};

/// 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 laplacian_cell_integral_0_otherwise: public ufc::cell_integral
{
public:

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

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

  /// Tabulate the tensor for the contribution from a local cell
  virtual void tabulate_tensor(double*  A,
                               const double * const *  w,
                               const double*  vertex_coordinates,
                               int cell_orientation) const
  {
    // Number of operations (multiply-add pairs) for Jacobian data:      3
    // Number of operations (multiply-add pairs) for geometry tensor:    12
    // Number of operations (multiply-add pairs) for tensor contraction: 11
    // Total number of operations (multiply-add pairs):                  26
    
    // Compute Jacobian
    double J[6];
    compute_jacobian_triangle_3d(J, vertex_coordinates);
    
    // Compute Jacobian inverse and determinant
    double K[6];
    double detJ;
    compute_jacobian_inverse_triangle_3d(K, detJ, J);
    
    // Set scale factor
    const double det = std::abs(detJ);
    
    // Compute geometry tensor
    const double G0_0_0 = det*(K[0]*K[0] + K[1]*K[1] + K[2]*K[2]);
    const double G0_0_1 = det*(K[0]*K[3] + K[1]*K[4] + K[2]*K[5]);
    const double G0_1_0 = det*(K[3]*K[0] + K[4]*K[1] + K[5]*K[2]);
    const double G0_1_1 = det*(K[3]*K[3] + K[4]*K[4] + K[5]*K[5]);
    
    // Compute element tensor
    A[0] = 0.5*G0_0_0 + 0.5*G0_0_1 + 0.5*G0_1_0 + 0.5*G0_1_1;
    A[1] = -0.5*G0_0_0 - 0.5*G0_1_0;
    A[2] = -0.5*G0_0_1 - 0.5*G0_1_1;
    A[3] = -0.5*G0_0_0 - 0.5*G0_0_1;
    A[4] = 0.5*G0_0_0;
    A[5] = 0.5*G0_0_1;
    A[6] = -0.5*G0_1_0 - 0.5*G0_1_1;
    A[7] = 0.5*G0_1_0;
    A[8] = 0.5*G0_1_1;
  }

};

/// 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 laplacian_form_0: public ufc::form
{
public:

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

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

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

  /// 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 0;
  }

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

  /// 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 the number of point domains
  virtual std::size_t num_point_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;
  }

  /// Return whether the form has any point integrals
  virtual bool has_point_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 laplacian_finite_element_0();
        break;
      }
    case 1:
      {
        return new laplacian_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 laplacian_dofmap_0();
        break;
      }
    case 1:
      {
        return new laplacian_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
  {
    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 point integral on sub domain i
  virtual ufc::point_integral* create_point_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 new laplacian_cell_integral_0_otherwise();
  }

  /// 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;
  }

  /// Create a new point integral on everywhere else
  virtual ufc::point_integral* create_default_point_integral() const
  {
    return 0;
  }

};

#endif
