// This code conforms with the UFC specification version 2018.1.0
// and was automatically generated by FFC version 2018.1.0.
//
// This code was generated with the option '-l dolfin' and
// contains DOLFIN-specific wrappers that depend on DOLFIN.
//
// This code was generated with the following parameters:
//

//  add_tabulate_tensor_timing:     False
//  convert_exceptions_to_warnings: False
//  cpp_optimize:                   True
//  cpp_optimize_flags:             '-O2'
//  epsilon:                        1e-14
//  error_control:                  False
//  external_include_dirs:          ''
//  external_includes:              ''
//  external_libraries:             ''
//  external_library_dirs:          ''
//  form_postfix:                   True
//  format:                         'dolfin'
//  generate_dummy_tabulate_tensor: False
//  max_signature_length:           0
//  optimize:                       True
//  precision:                      None
//  quadrature_degree:              None
//  quadrature_rule:                None
//  representation:                 'auto'
//  split:                          False

#ifndef __PRESSURE_H
#define __PRESSURE_H
#include <algorithm>
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <ufc.h>

class pressure_finite_element_0: public ufc::finite_element
{
public:

  pressure_finite_element_0() : ufc::finite_element()
  {
    // Do nothing
  }

  ~pressure_finite_element_0() override
  {
    // Do nothing
  }

  const char * signature() const final override
  {
    return "FiniteElement('Lagrange', triangle, 1)";
  }

  ufc::shape cell_shape() const final override
  {
    return ufc::shape::triangle;
  }

  std::size_t topological_dimension() const final override
  {
    return 2;
  }

  std::size_t geometric_dimension() const final override
  {
    return 2;
  }

  std::size_t space_dimension() const final override
  {
    return 3;
  }

  std::size_t value_rank() const final override
  {
    return 0;
  }

  std::size_t value_dimension(std::size_t i) const final override
  {
    return 1;
  }

  std::size_t value_size() const final override
  {
    return 1;
  }

  std::size_t reference_value_rank() const final override
  {
    return 0;
  }

  std::size_t reference_value_dimension(std::size_t i) const final override
  {
    return 1;
  }

  std::size_t reference_value_size() const final override
  {
    return 1;
  }

  std::size_t degree() const final override
  {
    return 1;
  }

  const char * family() const final override
  {
    return "Lagrange";
  }

  void evaluate_reference_basis(double * reference_values,
                                std::size_t num_points,
                                const double * X) const final override
  {
    static const double coefficients0[1][3] = { { 0.4714045207910317, -0.2886751345948129, -0.16666666666666666 } };
    static const double coefficients1[1][3] = { { 0.4714045207910317, 0.2886751345948129, -0.16666666666666666 } };
    static const double coefficients2[1][3] = { { 0.4714045207910316, 0.0, 0.3333333333333333 } };
    for (std::size_t k = 0; k < num_points * 3; ++k)
        reference_values[k] = 0.0;
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Map from UFC reference coordinate X to FIAT reference coordinate Y
        const double Y[2] = { 2.0 * X[ip * 2] - 1.0, 2.0 * X[ip * 2 + 1] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues1[3] = {};
        basisvalues1[0] = 1.0;
        const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
        basisvalues1[1] = tmp1_1;
        basisvalues1[2] = (0.5 + 1.5 * Y[1]) * basisvalues1[0];
        basisvalues1[0] *= std::sqrt(0.5);
        basisvalues1[2] *= std::sqrt(1.0);
        basisvalues1[1] *= std::sqrt(3.0);
        // Accumulate products of coefficients and basisvalues
        for (std::size_t r = 0; r < 3; ++r)
            reference_values[3 * ip] += coefficients0[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 3; ++r)
            reference_values[3 * ip + 1] += coefficients1[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 3; ++r)
            reference_values[3 * ip + 2] += coefficients2[0][r] * basisvalues1[r];
    }
  }

  void evaluate_reference_basis_derivatives(double * reference_values,
                                            std::size_t order,
                                            std::size_t num_points,
                                            const double * X) const final override
  {
    if (order == 0)
    {
        evaluate_reference_basis(reference_values, num_points, X);
        return;
    }
    const std::size_t num_derivatives = std::pow(2, order);
    std::fill_n(reference_values, num_points * 3 * num_derivatives, 0.0);
    if (order > 1)
        return;
    // Tables of derivatives of the polynomial base (transpose).
    alignas(32) static const double dmats0[2][3][3] =
        { { { 0.0, 0.0, 0.0 },
            { 4.8989794855663495, 0.0, 0.0 },
            { 0.0, 0.0, 0.0 } },
          { { 0.0, 0.0, 0.0 },
            { 2.449489742783182, 0.0, 0.0 },
            { 4.242640687119285, 0.0, 0.0 } } };
    static const double coefficients0[1][3] = { { 0.4714045207910317, -0.2886751345948129, -0.16666666666666666 } };
    static const double coefficients1[1][3] = { { 0.4714045207910317, 0.2886751345948129, -0.16666666666666666 } };
    static const double coefficients2[1][3] = { { 0.4714045207910316, 0.0, 0.3333333333333333 } };
    const std::size_t reference_offset[3] = {};
    const std::size_t num_components[3] = { 1, 1, 1 };
    // Precomputed combinations
    const std::size_t combinations[1][2][1] =
        { { { 0 },
            { 1 } } };
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Map from UFC reference coordinate X to FIAT reference coordinate Y
        const double Y[2] = { 2.0 * X[ip * 2] - 1.0, 2.0 * X[ip * 2 + 1] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues1[3] = {};
        basisvalues1[0] = 1.0;
        const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
        basisvalues1[1] = tmp1_1;
        basisvalues1[2] = (0.5 + 1.5 * Y[1]) * basisvalues1[0];
        basisvalues1[0] *= std::sqrt(0.5);
        basisvalues1[2] *= std::sqrt(1.0);
        basisvalues1[1] *= std::sqrt(3.0);
        // Loop over all dofs
        for (std::size_t i = 0; i < 3; ++i)
        {
            double derivatives[2] = {};
            switch (i)
            {
            case 0:
                // Compute reference derivatives for dof 0.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[3] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[3][3] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 9, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[3][3];
                        std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 9, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 3; ++s)
                        for (std::size_t t = 0; t < 3; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 3; ++s)
                        derivatives[r] += coefficients0[0][s] * aux[s];
                }
                break;
            case 1:
                // Compute reference derivatives for dof 1.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[3] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[3][3] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 9, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[3][3];
                        std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 9, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 3; ++s)
                        for (std::size_t t = 0; t < 3; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 3; ++s)
                        derivatives[r] += coefficients1[0][s] * aux[s];
                }
                break;
            case 2:
                // Compute reference derivatives for dof 2.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[3] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[3][3] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 9, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[3][3];
                        std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 9, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 3; ++s)
                        for (std::size_t t = 0; t < 3; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 3; ++s)
                        derivatives[r] += coefficients2[0][s] * aux[s];
                }
                break;
            }
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t c = 0; c < num_components[i]; ++c)
                    reference_values[3 * num_derivatives * ip + num_derivatives * i + r + (reference_offset[i] + c)] = derivatives[num_derivatives * c + r];
        }
    }
  }

  void transform_reference_basis_derivatives(double * values,
                                             std::size_t order,
                                             std::size_t num_points,
                                             const double * reference_values,
                                             const double * X,
                                             const double * J,
                                             const double * detJ,
                                             const double * K,
                                             int cell_orientation) const final override
  {
    const std::size_t num_derivatives = std::pow(2, order);
    // Precomputed combinations
    const std::size_t combinations[1][2][1] =
        { { { 0 },
            { 1 } } };
    std::fill_n(values, num_points * 3 * num_derivatives, 0.0);
    const std::size_t reference_offsets[3] = {};
    const std::size_t physical_offsets[3] = {};
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        double transform[2][2];
        for (std::size_t r = 0; r < num_derivatives; ++r)
            for (std::size_t s = 0; s < num_derivatives; ++s)
                transform[r][s] = 1.0;
        for (std::size_t r = 0; r < num_derivatives; ++r)
            for (std::size_t s = 0; s < num_derivatives; ++s)
                for (std::size_t k = 0; k < order; ++k)
                    transform[r][s] *= K[2 * 2 * ip + 2 * combinations[order - 1][s][k] + combinations[order - 1][r][k]];
        for (std::size_t d = 0; d < 3; ++d)
        {
            for (std::size_t s = 0; s < num_derivatives; ++s)
            {
                for (std::size_t i = 0; i < 1; ++i)
                {
                    // Using affine transform to map values back to the physical element.
                    const double mapped_value = reference_values[3 * num_derivatives * ip + num_derivatives * d + s + reference_offsets[d]];
                    // Mapping derivatives back to the physical element
                    for (std::size_t r = 0; r < num_derivatives; ++r)
                        values[3 * num_derivatives * ip + num_derivatives * d + r + (physical_offsets[d] + i)] += transform[r][s] * mapped_value;
                }
            }
        }
    }
  }

  void evaluate_basis(std::size_t i,
                      double * values,
                      const double * x,
                      const double * coordinate_dofs,
                      int cell_orientation,
                      const ufc::coordinate_mapping * cm=nullptr
                      ) const final override
  {
    double X[2] = {};
    double J[4];
    double detJ;
    double K[4];
    if (cm)
    {
        cm->compute_reference_geometry(X, J, &detJ, K, 1, x, coordinate_dofs, cell_orientation);
    }
    else
    {
        compute_jacobian_triangle_2d(J, coordinate_dofs);
        compute_jacobian_inverse_triangle_2d(K, detJ, J);
        // Compute constants
        const double C0 = coordinate_dofs[2] + coordinate_dofs[4];
        const double C1 = coordinate_dofs[3] + coordinate_dofs[5];
        // Get coordinates and map to the reference (FIAT) element
        double Y[2] = { (J[1] * (C1 - 2.0 * x[1]) + J[3] * (2.0 * x[0] - C0)) / detJ, (J[0] * (2.0 * x[1] - C1) + J[2] * (C0 - 2.0 * x[0])) / detJ };
        // Map to FFC reference coordinate
        for (std::size_t k = 0; k < 2; ++k)
            X[k] = (Y[k] + 1.0) / 2.0;
    }
    // Evaluate basis on reference element
    double ref_values[3];
    evaluate_reference_basis(ref_values, 1, X);
    // Push forward
    double physical_values[3];
    transform_reference_basis_derivatives(physical_values, 0, 1, ref_values, X, J, &detJ, K, cell_orientation);
    for (std::size_t k = 0; k < 1; ++k)
        values[k] = physical_values[i + k];
  }

  void evaluate_basis_all(double * values,
                          const double * x,
                          const double * coordinate_dofs,
                          int cell_orientation,
                          const ufc::coordinate_mapping * cm=nullptr
                          ) const final override
  {
    // Helper variable to hold value of a single dof.
    double dof_values = 0.0;
    // Loop dofs and call evaluate_basis
    for (std::size_t r = 0; r < 3; ++r)
    {
        evaluate_basis(r, &dof_values, x, coordinate_dofs, cell_orientation);
        values[r] = dof_values;
    }
  }

  void evaluate_basis_derivatives(std::size_t i,
                                  std::size_t n,
                                  double * values,
                                  const double * x,
                                  const double * coordinate_dofs,
                                  int cell_orientation,
                                  const ufc::coordinate_mapping * cm=nullptr
                                  ) const final override
  {
    std::size_t num_derivatives = std::pow(2, n);
    std::fill_n(values, num_derivatives, 0.0);
    // Call evaluate_basis_all if order of derivatives is equal to zero.
    if (n == 0)
    {
        evaluate_basis(i, values, x, coordinate_dofs, cell_orientation);
        return;
    }
    // If order of derivatives is greater than the maximum polynomial degree, return zeros.
    if (n > 1)
        return;
    // Compute Jacobian
    double J[4];
    compute_jacobian_triangle_2d(J, coordinate_dofs);
    // Compute Inverse Jacobian and determinant
    double K[4];
    double detJ;
    compute_jacobian_inverse_triangle_2d(K, detJ, J);
    // Compute constants
    const double C0 = coordinate_dofs[2] + coordinate_dofs[4];
    const double C1 = coordinate_dofs[3] + coordinate_dofs[5];
    // Get coordinates and map to the reference (FIAT) element
    double Y[2] = { (J[1] * (C1 - 2.0 * x[1]) + J[3] * (2.0 * x[0] - C0)) / detJ, (J[0] * (2.0 * x[1] - C1) + J[2] * (C0 - 2.0 * x[0])) / detJ };
    // Precomputed combinations
    const std::size_t combinations[1][2][1] =
        { { { 0 },
            { 1 } } };
    // Declare transformation matrix
    double transform[2][2] =
        { { 1.0, 1.0 },
          { 1.0, 1.0 } };
    // Construct transformation matrix
    for (std::size_t row = 0; row < num_derivatives; ++row)
        for (std::size_t col = 0; col < num_derivatives; ++col)
            for (std::size_t k = 0; k < n; ++k)
                transform[row][col] *= K[2 * combinations[n - 1][col][k] + combinations[n - 1][row][k]];
    switch (i)
    {
    case 0:
        {
            double basisvalues[3] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_1;
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            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.4714045207910317, -0.2886751345948129, -0.16666666666666666 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 4.8989794855663495, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0 } };
            static const double dmats1[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 2.449489742783182, 0.0, 0.0 },
                  { 4.242640687119285, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[2] = {};
            // 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 (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 9, 0.0);
                for (std::size_t t = 0; t < 3; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 9, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 3; ++s)
                    for (std::size_t t = 0; t < 3; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[r] += transform[r][s] * derivatives[s];
        }
        break;
    case 1:
        {
            double basisvalues[3] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_1;
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            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.4714045207910317, 0.2886751345948129, -0.16666666666666666 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 4.8989794855663495, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0 } };
            static const double dmats1[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 2.449489742783182, 0.0, 0.0 },
                  { 4.242640687119285, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[2] = {};
            // 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 (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 9, 0.0);
                for (std::size_t t = 0; t < 3; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 9, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 3; ++s)
                    for (std::size_t t = 0; t < 3; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[r] += transform[r][s] * derivatives[s];
        }
        break;
    case 2:
        {
            double basisvalues[3] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_1;
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            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.4714045207910316, 0.0, 0.3333333333333333 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 4.8989794855663495, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0 } };
            static const double dmats1[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 2.449489742783182, 0.0, 0.0 },
                  { 4.242640687119285, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[2] = {};
            // 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 (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 9, 0.0);
                for (std::size_t t = 0; t < 3; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 9, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 3; ++s)
                    for (std::size_t t = 0; t < 3; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[r] += transform[r][s] * derivatives[s];
        }
        break;
    }
  }

  void evaluate_basis_derivatives_all(std::size_t n,
                                      double * values,
                                      const double * x,
                                      const double * coordinate_dofs,
                                      int cell_orientation,
                                      const ufc::coordinate_mapping * cm=nullptr
                                      ) const final override
  {
    // Call evaluate_basis_all if order of derivatives is equal to zero.
    if (n == 0)
    {
        evaluate_basis_all(values, x, coordinate_dofs, cell_orientation);
        return;
    }
    unsigned int num_derivatives = std::pow(2, n);
    // Set values equal to zero.
    std::fill_n(values, num_derivatives * 3, 0.0);
    // If order of derivatives is greater than the maximum polynomial degree, return zeros.
    if (n > 1)
        return;
    // Helper variable to hold values of a single dof.
    double dof_values[2] = {};
    // Loop dofs and call evaluate_basis_derivatives.
    for (std::size_t r = 0; r < 3; ++r)
    {
        evaluate_basis_derivatives(r, n, dof_values, x, coordinate_dofs, cell_orientation);
        for (std::size_t s = 0; s < num_derivatives; ++s)
            values[num_derivatives * r + s] = dof_values[s];
    }
  }

  double evaluate_dof(std::size_t i,
                      const ufc::function& f,
                      const double * coordinate_dofs,
                      int cell_orientation,
                      const ufc::cell& c,
                      const ufc::coordinate_mapping * cm=nullptr
                      ) const final override
  {
    // Declare variables for result of evaluation
    double vals[1];
    // Declare variable for physical coordinates
    double y[2];
    switch (i)
    {
    case 0:
        {
            y[0] = coordinate_dofs[0];
            y[1] = coordinate_dofs[1];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 1:
        {
            y[0] = coordinate_dofs[2];
            y[1] = coordinate_dofs[3];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 2:
        {
            y[0] = coordinate_dofs[4];
            y[1] = coordinate_dofs[5];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    }
    return 0.0;
  }

  void evaluate_dofs(double * values,
                             const ufc::function& f,
                             const double * coordinate_dofs,
                             int cell_orientation,
                             const ufc::cell& c,
                             const ufc::coordinate_mapping * cm=nullptr
                             ) const final override
  {
    // Declare variables for result of evaluation
    double vals[1];
    // Declare variable for physical coordinates
    double y[2];
    y[0] = coordinate_dofs[0];
    y[1] = coordinate_dofs[1];
    f.evaluate(vals, y, c);
    values[0] = vals[0];
    y[0] = coordinate_dofs[2];
    y[1] = coordinate_dofs[3];
    f.evaluate(vals, y, c);
    values[1] = vals[0];
    y[0] = coordinate_dofs[4];
    y[1] = coordinate_dofs[5];
    f.evaluate(vals, y, c);
    values[2] = vals[0];
  }

  void interpolate_vertex_values(double * vertex_values,
                                 const double * dof_values,
                                 const double * coordinate_dofs,
                                 int cell_orientation,
                                 const ufc::coordinate_mapping * cm=nullptr
                                 ) const final override
  {
    // Evaluate function and change variables
    vertex_values[0] = dof_values[0];
    vertex_values[1] = dof_values[1];
    vertex_values[2] = dof_values[2];
  }

  void tabulate_dof_coordinates(double * dof_coordinates,
                                const double * coordinate_dofs,
                                const ufc::coordinate_mapping * cm=nullptr
                                ) const final override
  {
    dof_coordinates[0] = coordinate_dofs[0];
    dof_coordinates[1] = coordinate_dofs[1];
    dof_coordinates[2] = coordinate_dofs[2];
    dof_coordinates[2 + 1] = coordinate_dofs[3];
    dof_coordinates[2 * 2] = coordinate_dofs[4];
    dof_coordinates[2 * 2 + 1] = coordinate_dofs[5];
  }

  void tabulate_reference_dof_coordinates(double * reference_dof_coordinates) const final override
  {
    static const double dof_X[6] = { 0.0, 0.0, 1.0, 0.0, 0.0, 1.0 };
    std::copy_n(dof_X, 6, reference_dof_coordinates);
  }

  std::size_t num_sub_elements() const final override
  {
    return 0;
  }

  ufc::finite_element * create_sub_element(std::size_t i) const final override
  {
    return nullptr;
  }

  ufc::finite_element * create() const final override
  {
    return new pressure_finite_element_0();
  }

};


class pressure_finite_element_1: public ufc::finite_element
{
public:

  pressure_finite_element_1() : ufc::finite_element()
  {
    // Do nothing
  }

  ~pressure_finite_element_1() override
  {
    // Do nothing
  }

  const char * signature() const final override
  {
    return "VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2)";
  }

  ufc::shape cell_shape() const final override
  {
    return ufc::shape::triangle;
  }

  std::size_t topological_dimension() const final override
  {
    return 2;
  }

  std::size_t geometric_dimension() const final override
  {
    return 2;
  }

  std::size_t space_dimension() const final override
  {
    return 6;
  }

  std::size_t value_rank() const final override
  {
    return 1;
  }

  std::size_t value_dimension(std::size_t i) const final override
  {
    static const int return_values[1] = { 2 };
    if (i >= 1)
        return 1;
    return return_values[i];
  }

  std::size_t value_size() const final override
  {
    return 2;
  }

  std::size_t reference_value_rank() const final override
  {
    return 1;
  }

  std::size_t reference_value_dimension(std::size_t i) const final override
  {
    static const int return_values[1] = { 2 };
    if (i >= 1)
        return 1;
    return return_values[i];
  }

  std::size_t reference_value_size() const final override
  {
    return 2;
  }

  std::size_t degree() const final override
  {
    return 1;
  }

  const char * family() const final override
  {
    return "Lagrange";
  }

  void evaluate_reference_basis(double * reference_values,
                                std::size_t num_points,
                                const double * X) const final override
  {
    static const double coefficients0[1][3] = { { 0.4714045207910317, -0.2886751345948129, -0.16666666666666666 } };
    static const double coefficients1[1][3] = { { 0.4714045207910317, 0.2886751345948129, -0.16666666666666666 } };
    static const double coefficients2[1][3] = { { 0.4714045207910316, 0.0, 0.3333333333333333 } };
    for (std::size_t k = 0; k < num_points * 6 * 2; ++k)
        reference_values[k] = 0.0;
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Map from UFC reference coordinate X to FIAT reference coordinate Y
        const double Y[2] = { 2.0 * X[ip * 2] - 1.0, 2.0 * X[ip * 2 + 1] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues1[3] = {};
        basisvalues1[0] = 1.0;
        const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
        basisvalues1[1] = tmp1_1;
        basisvalues1[2] = (0.5 + 1.5 * Y[1]) * basisvalues1[0];
        basisvalues1[0] *= std::sqrt(0.5);
        basisvalues1[2] *= std::sqrt(1.0);
        basisvalues1[1] *= std::sqrt(3.0);
        // Accumulate products of coefficients and basisvalues
        for (std::size_t r = 0; r < 3; ++r)
            reference_values[6 * 2 * ip] += coefficients0[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 3; ++r)
            reference_values[6 * 2 * ip + 2] += coefficients1[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 3; ++r)
            reference_values[6 * 2 * ip + 2 * 2] += coefficients2[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 3; ++r)
            reference_values[6 * 2 * ip + 2 * 3 + 1] += coefficients0[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 3; ++r)
            reference_values[6 * 2 * ip + 2 * 4 + 1] += coefficients1[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 3; ++r)
            reference_values[6 * 2 * ip + 2 * 5 + 1] += coefficients2[0][r] * basisvalues1[r];
    }
  }

  void evaluate_reference_basis_derivatives(double * reference_values,
                                            std::size_t order,
                                            std::size_t num_points,
                                            const double * X) const final override
  {
    if (order == 0)
    {
        evaluate_reference_basis(reference_values, num_points, X);
        return;
    }
    const std::size_t num_derivatives = std::pow(2, order);
    std::fill_n(reference_values, num_points * 6 * num_derivatives * 2, 0.0);
    if (order > 1)
        return;
    // Tables of derivatives of the polynomial base (transpose).
    alignas(32) static const double dmats0[2][3][3] =
        { { { 0.0, 0.0, 0.0 },
            { 4.8989794855663495, 0.0, 0.0 },
            { 0.0, 0.0, 0.0 } },
          { { 0.0, 0.0, 0.0 },
            { 2.449489742783182, 0.0, 0.0 },
            { 4.242640687119285, 0.0, 0.0 } } };
    static const double coefficients0[1][3] = { { 0.4714045207910317, -0.2886751345948129, -0.16666666666666666 } };
    static const double coefficients1[1][3] = { { 0.4714045207910317, 0.2886751345948129, -0.16666666666666666 } };
    static const double coefficients2[1][3] = { { 0.4714045207910316, 0.0, 0.3333333333333333 } };
    const std::size_t reference_offset[6] = { 0, 0, 0, 1, 1, 1 };
    const std::size_t num_components[6] = { 1, 1, 1, 1, 1, 1 };
    // Precomputed combinations
    const std::size_t combinations[1][2][1] =
        { { { 0 },
            { 1 } } };
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Map from UFC reference coordinate X to FIAT reference coordinate Y
        const double Y[2] = { 2.0 * X[ip * 2] - 1.0, 2.0 * X[ip * 2 + 1] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues1[3] = {};
        basisvalues1[0] = 1.0;
        const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
        basisvalues1[1] = tmp1_1;
        basisvalues1[2] = (0.5 + 1.5 * Y[1]) * basisvalues1[0];
        basisvalues1[0] *= std::sqrt(0.5);
        basisvalues1[2] *= std::sqrt(1.0);
        basisvalues1[1] *= std::sqrt(3.0);
        // Loop over all dofs
        for (std::size_t i = 0; i < 6; ++i)
        {
            double derivatives[2] = {};
            switch (i)
            {
            case 0:
                // Compute reference derivatives for dof 0.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[3] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[3][3] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 9, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[3][3];
                        std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 9, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 3; ++s)
                        for (std::size_t t = 0; t < 3; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 3; ++s)
                        derivatives[r] += coefficients0[0][s] * aux[s];
                }
                break;
            case 1:
                // Compute reference derivatives for dof 1.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[3] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[3][3] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 9, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[3][3];
                        std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 9, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 3; ++s)
                        for (std::size_t t = 0; t < 3; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 3; ++s)
                        derivatives[r] += coefficients1[0][s] * aux[s];
                }
                break;
            case 2:
                // Compute reference derivatives for dof 2.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[3] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[3][3] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 9, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[3][3];
                        std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 9, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 3; ++s)
                        for (std::size_t t = 0; t < 3; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 3; ++s)
                        derivatives[r] += coefficients2[0][s] * aux[s];
                }
                break;
            case 3:
                // Compute reference derivatives for dof 3.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[3] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[3][3] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 9, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[3][3];
                        std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 9, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 3; ++s)
                        for (std::size_t t = 0; t < 3; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 3; ++s)
                        derivatives[r] += coefficients0[0][s] * aux[s];
                }
                break;
            case 4:
                // Compute reference derivatives for dof 4.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[3] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[3][3] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 9, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[3][3];
                        std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 9, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 3; ++s)
                        for (std::size_t t = 0; t < 3; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 3; ++s)
                        derivatives[r] += coefficients1[0][s] * aux[s];
                }
                break;
            case 5:
                // Compute reference derivatives for dof 5.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[3] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[3][3] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 9, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[3][3];
                        std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 9, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 3; ++s)
                        for (std::size_t t = 0; t < 3; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 3; ++s)
                        derivatives[r] += coefficients2[0][s] * aux[s];
                }
                break;
            }
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t c = 0; c < num_components[i]; ++c)
                    reference_values[6 * (num_derivatives * 2) * ip + num_derivatives * 2 * i + 2 * r + (reference_offset[i] + c)] = derivatives[num_derivatives * c + r];
        }
    }
  }

  void transform_reference_basis_derivatives(double * values,
                                             std::size_t order,
                                             std::size_t num_points,
                                             const double * reference_values,
                                             const double * X,
                                             const double * J,
                                             const double * detJ,
                                             const double * K,
                                             int cell_orientation) const final override
  {
    const std::size_t num_derivatives = std::pow(2, order);
    // Precomputed combinations
    const std::size_t combinations[1][2][1] =
        { { { 0 },
            { 1 } } };
    std::fill_n(values, num_points * 6 * num_derivatives * 2, 0.0);
    const std::size_t reference_offsets[6] = { 0, 0, 0, 1, 1, 1 };
    const std::size_t physical_offsets[6] = { 0, 0, 0, 1, 1, 1 };
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        double transform[2][2];
        for (std::size_t r = 0; r < num_derivatives; ++r)
            for (std::size_t s = 0; s < num_derivatives; ++s)
                transform[r][s] = 1.0;
        for (std::size_t r = 0; r < num_derivatives; ++r)
            for (std::size_t s = 0; s < num_derivatives; ++s)
                for (std::size_t k = 0; k < order; ++k)
                    transform[r][s] *= K[2 * 2 * ip + 2 * combinations[order - 1][s][k] + combinations[order - 1][r][k]];
        for (std::size_t d = 0; d < 6; ++d)
        {
            for (std::size_t s = 0; s < num_derivatives; ++s)
            {
                for (std::size_t i = 0; i < 1; ++i)
                {
                    // Using affine transform to map values back to the physical element.
                    const double mapped_value = reference_values[6 * (num_derivatives * 2) * ip + num_derivatives * 2 * d + 2 * s + reference_offsets[d]];
                    // Mapping derivatives back to the physical element
                    for (std::size_t r = 0; r < num_derivatives; ++r)
                        values[6 * (num_derivatives * 2) * ip + num_derivatives * 2 * d + 2 * r + (physical_offsets[d] + i)] += transform[r][s] * mapped_value;
                }
            }
        }
    }
  }

  void evaluate_basis(std::size_t i,
                      double * values,
                      const double * x,
                      const double * coordinate_dofs,
                      int cell_orientation,
                      const ufc::coordinate_mapping * cm=nullptr
                      ) const final override
  {
    double X[2] = {};
    double J[4];
    double detJ;
    double K[4];
    if (cm)
    {
        cm->compute_reference_geometry(X, J, &detJ, K, 1, x, coordinate_dofs, cell_orientation);
    }
    else
    {
        compute_jacobian_triangle_2d(J, coordinate_dofs);
        compute_jacobian_inverse_triangle_2d(K, detJ, J);
        // Compute constants
        const double C0 = coordinate_dofs[2] + coordinate_dofs[4];
        const double C1 = coordinate_dofs[3] + coordinate_dofs[5];
        // Get coordinates and map to the reference (FIAT) element
        double Y[2] = { (J[1] * (C1 - 2.0 * x[1]) + J[3] * (2.0 * x[0] - C0)) / detJ, (J[0] * (2.0 * x[1] - C1) + J[2] * (C0 - 2.0 * x[0])) / detJ };
        // Map to FFC reference coordinate
        for (std::size_t k = 0; k < 2; ++k)
            X[k] = (Y[k] + 1.0) / 2.0;
    }
    // Evaluate basis on reference element
    double ref_values[12];
    evaluate_reference_basis(ref_values, 1, X);
    // Push forward
    double physical_values[12];
    transform_reference_basis_derivatives(physical_values, 0, 1, ref_values, X, J, &detJ, K, cell_orientation);
    for (std::size_t k = 0; k < 2; ++k)
        values[k] = physical_values[2 * i + k];
  }

  void evaluate_basis_all(double * values,
                          const double * x,
                          const double * coordinate_dofs,
                          int cell_orientation,
                          const ufc::coordinate_mapping * cm=nullptr
                          ) const final override
  {
    // Helper variable to hold values of a single dof.
    double dof_values[2] = {};
    // Loop dofs and call evaluate_basis
    for (std::size_t r = 0; r < 6; ++r)
    {
        evaluate_basis(r, dof_values, x, coordinate_dofs, cell_orientation);
        for (std::size_t s = 0; s < 2; ++s)
            values[r * 2 + s] = dof_values[s];
    }
  }

  void evaluate_basis_derivatives(std::size_t i,
                                  std::size_t n,
                                  double * values,
                                  const double * x,
                                  const double * coordinate_dofs,
                                  int cell_orientation,
                                  const ufc::coordinate_mapping * cm=nullptr
                                  ) const final override
  {
    std::size_t num_derivatives = std::pow(2, n);
    std::fill_n(values, 2 * num_derivatives, 0.0);
    // Call evaluate_basis_all if order of derivatives is equal to zero.
    if (n == 0)
    {
        evaluate_basis(i, values, x, coordinate_dofs, cell_orientation);
        return;
    }
    // If order of derivatives is greater than the maximum polynomial degree, return zeros.
    if (n > 1)
        return;
    // Compute Jacobian
    double J[4];
    compute_jacobian_triangle_2d(J, coordinate_dofs);
    // Compute Inverse Jacobian and determinant
    double K[4];
    double detJ;
    compute_jacobian_inverse_triangle_2d(K, detJ, J);
    // Compute constants
    const double C0 = coordinate_dofs[2] + coordinate_dofs[4];
    const double C1 = coordinate_dofs[3] + coordinate_dofs[5];
    // Get coordinates and map to the reference (FIAT) element
    double Y[2] = { (J[1] * (C1 - 2.0 * x[1]) + J[3] * (2.0 * x[0] - C0)) / detJ, (J[0] * (2.0 * x[1] - C1) + J[2] * (C0 - 2.0 * x[0])) / detJ };
    // Precomputed combinations
    const std::size_t combinations[1][2][1] =
        { { { 0 },
            { 1 } } };
    // Declare transformation matrix
    double transform[2][2] =
        { { 1.0, 1.0 },
          { 1.0, 1.0 } };
    // Construct transformation matrix
    for (std::size_t row = 0; row < num_derivatives; ++row)
        for (std::size_t col = 0; col < num_derivatives; ++col)
            for (std::size_t k = 0; k < n; ++k)
                transform[row][col] *= K[2 * combinations[n - 1][col][k] + combinations[n - 1][row][k]];
    switch (i)
    {
    case 0:
        {
            double basisvalues[3] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_1;
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            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.4714045207910317, -0.2886751345948129, -0.16666666666666666 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 4.8989794855663495, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0 } };
            static const double dmats1[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 2.449489742783182, 0.0, 0.0 },
                  { 4.242640687119285, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[2] = {};
            // 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 (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 9, 0.0);
                for (std::size_t t = 0; t < 3; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 9, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 3; ++s)
                    for (std::size_t t = 0; t < 3; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[r] += transform[r][s] * derivatives[s];
        }
        break;
    case 1:
        {
            double basisvalues[3] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_1;
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            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.4714045207910317, 0.2886751345948129, -0.16666666666666666 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 4.8989794855663495, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0 } };
            static const double dmats1[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 2.449489742783182, 0.0, 0.0 },
                  { 4.242640687119285, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[2] = {};
            // 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 (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 9, 0.0);
                for (std::size_t t = 0; t < 3; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 9, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 3; ++s)
                    for (std::size_t t = 0; t < 3; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[r] += transform[r][s] * derivatives[s];
        }
        break;
    case 2:
        {
            double basisvalues[3] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_1;
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            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.4714045207910316, 0.0, 0.3333333333333333 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 4.8989794855663495, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0 } };
            static const double dmats1[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 2.449489742783182, 0.0, 0.0 },
                  { 4.242640687119285, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[2] = {};
            // 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 (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 9, 0.0);
                for (std::size_t t = 0; t < 3; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 9, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 3; ++s)
                    for (std::size_t t = 0; t < 3; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[r] += transform[r][s] * derivatives[s];
        }
        break;
    case 3:
        {
            double basisvalues[3] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_1;
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            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.4714045207910317, -0.2886751345948129, -0.16666666666666666 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 4.8989794855663495, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0 } };
            static const double dmats1[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 2.449489742783182, 0.0, 0.0 },
                  { 4.242640687119285, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[2] = {};
            // 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 (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 9, 0.0);
                for (std::size_t t = 0; t < 3; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 9, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 3; ++s)
                    for (std::size_t t = 0; t < 3; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[num_derivatives + r] += transform[r][s] * derivatives[s];
        }
        break;
    case 4:
        {
            double basisvalues[3] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_1;
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            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.4714045207910317, 0.2886751345948129, -0.16666666666666666 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 4.8989794855663495, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0 } };
            static const double dmats1[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 2.449489742783182, 0.0, 0.0 },
                  { 4.242640687119285, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[2] = {};
            // 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 (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 9, 0.0);
                for (std::size_t t = 0; t < 3; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 9, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 3; ++s)
                    for (std::size_t t = 0; t < 3; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[num_derivatives + r] += transform[r][s] * derivatives[s];
        }
        break;
    case 5:
        {
            double basisvalues[3] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_1;
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            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.4714045207910316, 0.0, 0.3333333333333333 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 4.8989794855663495, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0 } };
            static const double dmats1[3][3] =
                { { 0.0, 0.0, 0.0 },
                  { 2.449489742783182, 0.0, 0.0 },
                  { 4.242640687119285, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[2] = {};
            // 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 (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 9, 0.0);
                for (std::size_t t = 0; t < 3; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 9, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 9, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 3; ++t)
                            for (std::size_t u = 0; u < 3; ++u)
                                for (std::size_t tu = 0; tu < 3; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 3; ++s)
                    for (std::size_t t = 0; t < 3; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[num_derivatives + r] += transform[r][s] * derivatives[s];
        }
        break;
    }
  }

  void evaluate_basis_derivatives_all(std::size_t n,
                                      double * values,
                                      const double * x,
                                      const double * coordinate_dofs,
                                      int cell_orientation,
                                      const ufc::coordinate_mapping * cm=nullptr
                                      ) const final override
  {
    // Call evaluate_basis_all if order of derivatives is equal to zero.
    if (n == 0)
    {
        evaluate_basis_all(values, x, coordinate_dofs, cell_orientation);
        return;
    }
    unsigned int num_derivatives = std::pow(2, n);
    // Set values equal to zero.
    std::fill_n(values, 2 * num_derivatives * 6, 0.0);
    // If order of derivatives is greater than the maximum polynomial degree, return zeros.
    if (n > 1)
        return;
    // Helper variable to hold values of a single dof.
    double dof_values[4] = {};
    // Loop dofs and call evaluate_basis_derivatives.
    for (std::size_t r = 0; r < 6; ++r)
    {
        evaluate_basis_derivatives(r, n, dof_values, x, coordinate_dofs, cell_orientation);
        for (std::size_t s = 0; s < 2 * num_derivatives; ++s)
            values[2 * num_derivatives * r + s] = dof_values[s];
    }
  }

  double evaluate_dof(std::size_t i,
                      const ufc::function& f,
                      const double * coordinate_dofs,
                      int cell_orientation,
                      const ufc::cell& c,
                      const ufc::coordinate_mapping * cm=nullptr
                      ) const final override
  {
    // Declare variables for result of evaluation
    double vals[2];
    // Declare variable for physical coordinates
    double y[2];
    switch (i)
    {
    case 0:
        {
            y[0] = coordinate_dofs[0];
            y[1] = coordinate_dofs[1];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 1:
        {
            y[0] = coordinate_dofs[2];
            y[1] = coordinate_dofs[3];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 2:
        {
            y[0] = coordinate_dofs[4];
            y[1] = coordinate_dofs[5];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 3:
        {
            y[0] = coordinate_dofs[0];
            y[1] = coordinate_dofs[1];
            f.evaluate(vals, y, c);
            return vals[1];
        }
        break;
    case 4:
        {
            y[0] = coordinate_dofs[2];
            y[1] = coordinate_dofs[3];
            f.evaluate(vals, y, c);
            return vals[1];
        }
        break;
    case 5:
        {
            y[0] = coordinate_dofs[4];
            y[1] = coordinate_dofs[5];
            f.evaluate(vals, y, c);
            return vals[1];
        }
        break;
    }
    return 0.0;
  }

  void evaluate_dofs(double * values,
                             const ufc::function& f,
                             const double * coordinate_dofs,
                             int cell_orientation,
                             const ufc::cell& c,
                             const ufc::coordinate_mapping * cm=nullptr
                             ) const final override
  {
    // Declare variables for result of evaluation
    double vals[2];
    // Declare variable for physical coordinates
    double y[2];
    y[0] = coordinate_dofs[0];
    y[1] = coordinate_dofs[1];
    f.evaluate(vals, y, c);
    values[0] = vals[0];
    y[0] = coordinate_dofs[2];
    y[1] = coordinate_dofs[3];
    f.evaluate(vals, y, c);
    values[1] = vals[0];
    y[0] = coordinate_dofs[4];
    y[1] = coordinate_dofs[5];
    f.evaluate(vals, y, c);
    values[2] = vals[0];
    y[0] = coordinate_dofs[0];
    y[1] = coordinate_dofs[1];
    f.evaluate(vals, y, c);
    values[3] = vals[1];
    y[0] = coordinate_dofs[2];
    y[1] = coordinate_dofs[3];
    f.evaluate(vals, y, c);
    values[4] = vals[1];
    y[0] = coordinate_dofs[4];
    y[1] = coordinate_dofs[5];
    f.evaluate(vals, y, c);
    values[5] = vals[1];
  }

  void interpolate_vertex_values(double * vertex_values,
                                 const double * dof_values,
                                 const double * coordinate_dofs,
                                 int cell_orientation,
                                 const ufc::coordinate_mapping * cm=nullptr
                                 ) const final override
  {
    // Evaluate function and change variables
    vertex_values[0] = dof_values[0];
    vertex_values[2] = dof_values[1];
    vertex_values[4] = dof_values[2];
    // Evaluate function and change variables
    vertex_values[1] = dof_values[3];
    vertex_values[3] = dof_values[4];
    vertex_values[5] = dof_values[5];
  }

  void tabulate_dof_coordinates(double * dof_coordinates,
                                const double * coordinate_dofs,
                                const ufc::coordinate_mapping * cm=nullptr
                                ) const final override
  {
    dof_coordinates[0] = coordinate_dofs[0];
    dof_coordinates[1] = coordinate_dofs[1];
    dof_coordinates[2] = coordinate_dofs[2];
    dof_coordinates[2 + 1] = coordinate_dofs[3];
    dof_coordinates[2 * 2] = coordinate_dofs[4];
    dof_coordinates[2 * 2 + 1] = coordinate_dofs[5];
    dof_coordinates[2 * 3] = coordinate_dofs[0];
    dof_coordinates[2 * 3 + 1] = coordinate_dofs[1];
    dof_coordinates[2 * 4] = coordinate_dofs[2];
    dof_coordinates[2 * 4 + 1] = coordinate_dofs[3];
    dof_coordinates[2 * 5] = coordinate_dofs[4];
    dof_coordinates[2 * 5 + 1] = coordinate_dofs[5];
  }

  void tabulate_reference_dof_coordinates(double * reference_dof_coordinates) const final override
  {
    static const double dof_X[12] = { 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0 };
    std::copy_n(dof_X, 12, reference_dof_coordinates);
  }

  std::size_t num_sub_elements() const final override
  {
    return 2;
  }

  ufc::finite_element * create_sub_element(std::size_t i) const final override
  {
    switch (i)
    {
    case 0:
        return new pressure_finite_element_0();
    case 1:
        return new pressure_finite_element_0();
    default:
        return nullptr;
    }
  }

  ufc::finite_element * create() const final override
  {
    return new pressure_finite_element_1();
  }

};


class pressure_finite_element_2: public ufc::finite_element
{
public:

  pressure_finite_element_2() : ufc::finite_element()
  {
    // Do nothing
  }

  ~pressure_finite_element_2() override
  {
    // Do nothing
  }

  const char * signature() const final override
  {
    return "FiniteElement('Lagrange', triangle, 2)";
  }

  ufc::shape cell_shape() const final override
  {
    return ufc::shape::triangle;
  }

  std::size_t topological_dimension() const final override
  {
    return 2;
  }

  std::size_t geometric_dimension() const final override
  {
    return 2;
  }

  std::size_t space_dimension() const final override
  {
    return 6;
  }

  std::size_t value_rank() const final override
  {
    return 0;
  }

  std::size_t value_dimension(std::size_t i) const final override
  {
    return 1;
  }

  std::size_t value_size() const final override
  {
    return 1;
  }

  std::size_t reference_value_rank() const final override
  {
    return 0;
  }

  std::size_t reference_value_dimension(std::size_t i) const final override
  {
    return 1;
  }

  std::size_t reference_value_size() const final override
  {
    return 1;
  }

  std::size_t degree() const final override
  {
    return 2;
  }

  const char * family() const final override
  {
    return "Lagrange";
  }

  void evaluate_reference_basis(double * reference_values,
                                std::size_t num_points,
                                const double * X) const final override
  {
    static const double coefficients0[1][6] = { { 0.0, -0.17320508075688776, -0.09999999999999999, 0.12171612389003691, 0.09428090415820634, 0.05443310539518172 } };
    static const double coefficients1[1][6] = { { 0.0, 0.1732050807568877, -0.09999999999999999, 0.12171612389003691, -0.09428090415820636, 0.05443310539518175 } };
    static const double coefficients2[1][6] = { { 0.0, 0.0, 0.20000000000000004, 0.0, 0.0, 0.16329931618554522 } };
    static const double coefficients3[1][6] = { { 0.47140452079103157, 0.23094010767585027, 0.1333333333333333, 0.0, 0.18856180831641264, -0.1632993161855452 } };
    static const double coefficients4[1][6] = { { 0.47140452079103157, -0.23094010767585027, 0.13333333333333333, 0.0, -0.18856180831641264, -0.16329931618554522 } };
    static const double coefficients5[1][6] = { { 0.4714045207910316, 0.0, -0.26666666666666666, -0.24343224778007383, 0.0, 0.05443310539518172 } };
    for (std::size_t k = 0; k < num_points * 6; ++k)
        reference_values[k] = 0.0;
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Map from UFC reference coordinate X to FIAT reference coordinate Y
        const double Y[2] = { 2.0 * X[ip * 2] - 1.0, 2.0 * X[ip * 2 + 1] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues2[6] = {};
        basisvalues2[0] = 1.0;
        const double tmp1_2 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
        basisvalues2[1] = tmp1_2;
        const double tmp2_2 = (1.0 - Y[1]) / 2.0;
        const double tmp3_2 = tmp2_2 * tmp2_2;
        basisvalues2[3] = 1.5 * tmp1_2 * basisvalues2[1] - 0.5 * tmp3_2 * basisvalues2[0];
        basisvalues2[2] = (0.5 + 1.5 * Y[1]) * basisvalues2[0];
        basisvalues2[4] = (1.5 + 2.5 * Y[1]) * basisvalues2[1];
        basisvalues2[5] = (0.1111111111111111 + 1.6666666666666667 * Y[1]) * basisvalues2[2] - 0.5555555555555556 * basisvalues2[0];
        basisvalues2[0] *= std::sqrt(0.5);
        basisvalues2[2] *= std::sqrt(1.0);
        basisvalues2[5] *= std::sqrt(1.5);
        basisvalues2[1] *= std::sqrt(3.0);
        basisvalues2[4] *= std::sqrt(4.5);
        basisvalues2[3] *= std::sqrt(7.5);
        // Accumulate products of coefficients and basisvalues
        for (std::size_t r = 0; r < 6; ++r)
            reference_values[6 * ip] += coefficients0[0][r] * basisvalues2[r];
        for (std::size_t r = 0; r < 6; ++r)
            reference_values[6 * ip + 1] += coefficients1[0][r] * basisvalues2[r];
        for (std::size_t r = 0; r < 6; ++r)
            reference_values[6 * ip + 2] += coefficients2[0][r] * basisvalues2[r];
        for (std::size_t r = 0; r < 6; ++r)
            reference_values[6 * ip + 3] += coefficients3[0][r] * basisvalues2[r];
        for (std::size_t r = 0; r < 6; ++r)
            reference_values[6 * ip + 4] += coefficients4[0][r] * basisvalues2[r];
        for (std::size_t r = 0; r < 6; ++r)
            reference_values[6 * ip + 5] += coefficients5[0][r] * basisvalues2[r];
    }
  }

  void evaluate_reference_basis_derivatives(double * reference_values,
                                            std::size_t order,
                                            std::size_t num_points,
                                            const double * X) const final override
  {
    if (order == 0)
    {
        evaluate_reference_basis(reference_values, num_points, X);
        return;
    }
    const std::size_t num_derivatives = std::pow(2, order);
    std::fill_n(reference_values, num_points * 6 * num_derivatives, 0.0);
    if (order > 2)
        return;
    // Tables of derivatives of the polynomial base (transpose).
    alignas(32) static const double dmats0[2][6][6] =
        { { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
            { 4.898979485566347, 0.0, 0.0, 0.0, 0.0, 0.0 },
            { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
            { 0.0, 9.486832980505149, 0.0, 0.0, 0.0, 0.0 },
            { 3.999999999999994, 0.0, 7.0710678118654675, 0.0, 0.0, 0.0 },
            { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } },
          { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
            { 2.4494897427831805, 0.0, 0.0, 0.0, 0.0, 0.0 },
            { 4.242640687119282, 0.0, 0.0, 0.0, 0.0, 0.0 },
            { 2.581988897471618, 4.743416490252574, -0.91287092917528, 0.0, 0.0, 0.0 },
            { 1.9999999999999152, 6.123724356957939, 3.5355339059327338, 0.0, 0.0, 0.0 },
            { -2.3094010767584887, 0.0, 8.164965809277266, 0.0, 0.0, 0.0 } } };
    static const double coefficients0[1][6] = { { 0.0, -0.17320508075688776, -0.09999999999999999, 0.12171612389003691, 0.09428090415820634, 0.05443310539518172 } };
    static const double coefficients1[1][6] = { { 0.0, 0.1732050807568877, -0.09999999999999999, 0.12171612389003691, -0.09428090415820636, 0.05443310539518175 } };
    static const double coefficients2[1][6] = { { 0.0, 0.0, 0.20000000000000004, 0.0, 0.0, 0.16329931618554522 } };
    static const double coefficients3[1][6] = { { 0.47140452079103157, 0.23094010767585027, 0.1333333333333333, 0.0, 0.18856180831641264, -0.1632993161855452 } };
    static const double coefficients4[1][6] = { { 0.47140452079103157, -0.23094010767585027, 0.13333333333333333, 0.0, -0.18856180831641264, -0.16329931618554522 } };
    static const double coefficients5[1][6] = { { 0.4714045207910316, 0.0, -0.26666666666666666, -0.24343224778007383, 0.0, 0.05443310539518172 } };
    const std::size_t reference_offset[6] = {};
    const std::size_t num_components[6] = { 1, 1, 1, 1, 1, 1 };
    // Precomputed combinations
    const std::size_t combinations[2][4][2] =
        { { { 0, 0 },
            { 1, 0 },
            { 0, 0 },
            { 1, 0 } },
          { { 0, 0 },
            { 0, 1 },
            { 1, 0 },
            { 1, 1 } } };
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Map from UFC reference coordinate X to FIAT reference coordinate Y
        const double Y[2] = { 2.0 * X[ip * 2] - 1.0, 2.0 * X[ip * 2 + 1] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues2[6] = {};
        basisvalues2[0] = 1.0;
        const double tmp1_2 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
        basisvalues2[1] = tmp1_2;
        const double tmp2_2 = (1.0 - Y[1]) / 2.0;
        const double tmp3_2 = tmp2_2 * tmp2_2;
        basisvalues2[3] = 1.5 * tmp1_2 * basisvalues2[1] - 0.5 * tmp3_2 * basisvalues2[0];
        basisvalues2[2] = (0.5 + 1.5 * Y[1]) * basisvalues2[0];
        basisvalues2[4] = (1.5 + 2.5 * Y[1]) * basisvalues2[1];
        basisvalues2[5] = (0.1111111111111111 + 1.6666666666666667 * Y[1]) * basisvalues2[2] - 0.5555555555555556 * basisvalues2[0];
        basisvalues2[0] *= std::sqrt(0.5);
        basisvalues2[2] *= std::sqrt(1.0);
        basisvalues2[5] *= std::sqrt(1.5);
        basisvalues2[1] *= std::sqrt(3.0);
        basisvalues2[4] *= std::sqrt(4.5);
        basisvalues2[3] *= std::sqrt(7.5);
        // Loop over all dofs
        for (std::size_t i = 0; i < 6; ++i)
        {
            double derivatives[4] = {};
            switch (i)
            {
            case 0:
                // Compute reference derivatives for dof 0.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[6] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[6][6] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 36, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[6][6];
                        std::copy_n(&dmats[0][0], 36, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 36, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 6; ++s)
                        for (std::size_t t = 0; t < 6; ++t)
                            aux[s] += dmats[s][t] * basisvalues2[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 6; ++s)
                        derivatives[r] += coefficients0[0][s] * aux[s];
                }
                break;
            case 1:
                // Compute reference derivatives for dof 1.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[6] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[6][6] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 36, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[6][6];
                        std::copy_n(&dmats[0][0], 36, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 36, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 6; ++s)
                        for (std::size_t t = 0; t < 6; ++t)
                            aux[s] += dmats[s][t] * basisvalues2[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 6; ++s)
                        derivatives[r] += coefficients1[0][s] * aux[s];
                }
                break;
            case 2:
                // Compute reference derivatives for dof 2.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[6] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[6][6] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 36, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[6][6];
                        std::copy_n(&dmats[0][0], 36, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 36, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 6; ++s)
                        for (std::size_t t = 0; t < 6; ++t)
                            aux[s] += dmats[s][t] * basisvalues2[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 6; ++s)
                        derivatives[r] += coefficients2[0][s] * aux[s];
                }
                break;
            case 3:
                // Compute reference derivatives for dof 3.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[6] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[6][6] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 36, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[6][6];
                        std::copy_n(&dmats[0][0], 36, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 36, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 6; ++s)
                        for (std::size_t t = 0; t < 6; ++t)
                            aux[s] += dmats[s][t] * basisvalues2[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 6; ++s)
                        derivatives[r] += coefficients3[0][s] * aux[s];
                }
                break;
            case 4:
                // Compute reference derivatives for dof 4.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[6] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[6][6] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 36, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[6][6];
                        std::copy_n(&dmats[0][0], 36, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 36, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 6; ++s)
                        for (std::size_t t = 0; t < 6; ++t)
                            aux[s] += dmats[s][t] * basisvalues2[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 6; ++s)
                        derivatives[r] += coefficients4[0][s] * aux[s];
                }
                break;
            case 5:
                // Compute reference derivatives for dof 5.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[6] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[6][6] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 36, &dmats[0][0]);
                    // Looping derivative order to generate dmats.
                    for (std::size_t s = 1; s < order; ++s)
                    {
                        // Store previous dmats matrix.
                        double dmats_old[6][6];
                        std::copy_n(&dmats[0][0], 36, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 36, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 6; ++s)
                        for (std::size_t t = 0; t < 6; ++t)
                            aux[s] += dmats[s][t] * basisvalues2[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 6; ++s)
                        derivatives[r] += coefficients5[0][s] * aux[s];
                }
                break;
            }
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t c = 0; c < num_components[i]; ++c)
                    reference_values[6 * num_derivatives * ip + num_derivatives * i + r + (reference_offset[i] + c)] = derivatives[num_derivatives * c + r];
        }
    }
  }

  void transform_reference_basis_derivatives(double * values,
                                             std::size_t order,
                                             std::size_t num_points,
                                             const double * reference_values,
                                             const double * X,
                                             const double * J,
                                             const double * detJ,
                                             const double * K,
                                             int cell_orientation) const final override
  {
    const std::size_t num_derivatives = std::pow(2, order);
    // Precomputed combinations
    const std::size_t combinations[2][4][2] =
        { { { 0, 0 },
            { 1, 0 },
            { 0, 0 },
            { 1, 0 } },
          { { 0, 0 },
            { 0, 1 },
            { 1, 0 },
            { 1, 1 } } };
    std::fill_n(values, num_points * 6 * num_derivatives, 0.0);
    const std::size_t reference_offsets[6] = {};
    const std::size_t physical_offsets[6] = {};
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        double transform[4][4];
        for (std::size_t r = 0; r < num_derivatives; ++r)
            for (std::size_t s = 0; s < num_derivatives; ++s)
                transform[r][s] = 1.0;
        for (std::size_t r = 0; r < num_derivatives; ++r)
            for (std::size_t s = 0; s < num_derivatives; ++s)
                for (std::size_t k = 0; k < order; ++k)
                    transform[r][s] *= K[2 * 2 * ip + 2 * combinations[order - 1][s][k] + combinations[order - 1][r][k]];
        for (std::size_t d = 0; d < 6; ++d)
        {
            for (std::size_t s = 0; s < num_derivatives; ++s)
            {
                for (std::size_t i = 0; i < 1; ++i)
                {
                    // Using affine transform to map values back to the physical element.
                    const double mapped_value = reference_values[6 * num_derivatives * ip + num_derivatives * d + s + reference_offsets[d]];
                    // Mapping derivatives back to the physical element
                    for (std::size_t r = 0; r < num_derivatives; ++r)
                        values[6 * num_derivatives * ip + num_derivatives * d + r + (physical_offsets[d] + i)] += transform[r][s] * mapped_value;
                }
            }
        }
    }
  }

  void evaluate_basis(std::size_t i,
                      double * values,
                      const double * x,
                      const double * coordinate_dofs,
                      int cell_orientation,
                      const ufc::coordinate_mapping * cm=nullptr
                      ) const final override
  {
    double X[2] = {};
    double J[4];
    double detJ;
    double K[4];
    if (cm)
    {
        cm->compute_reference_geometry(X, J, &detJ, K, 1, x, coordinate_dofs, cell_orientation);
    }
    else
    {
        compute_jacobian_triangle_2d(J, coordinate_dofs);
        compute_jacobian_inverse_triangle_2d(K, detJ, J);
        // Compute constants
        const double C0 = coordinate_dofs[2] + coordinate_dofs[4];
        const double C1 = coordinate_dofs[3] + coordinate_dofs[5];
        // Get coordinates and map to the reference (FIAT) element
        double Y[2] = { (J[1] * (C1 - 2.0 * x[1]) + J[3] * (2.0 * x[0] - C0)) / detJ, (J[0] * (2.0 * x[1] - C1) + J[2] * (C0 - 2.0 * x[0])) / detJ };
        // Map to FFC reference coordinate
        for (std::size_t k = 0; k < 2; ++k)
            X[k] = (Y[k] + 1.0) / 2.0;
    }
    // Evaluate basis on reference element
    double ref_values[6];
    evaluate_reference_basis(ref_values, 1, X);
    // Push forward
    double physical_values[6];
    transform_reference_basis_derivatives(physical_values, 0, 1, ref_values, X, J, &detJ, K, cell_orientation);
    for (std::size_t k = 0; k < 1; ++k)
        values[k] = physical_values[i + k];
  }

  void evaluate_basis_all(double * values,
                          const double * x,
                          const double * coordinate_dofs,
                          int cell_orientation,
                          const ufc::coordinate_mapping * cm=nullptr
                          ) const final override
  {
    // Helper variable to hold value of a single dof.
    double dof_values = 0.0;
    // Loop dofs and call evaluate_basis
    for (std::size_t r = 0; r < 6; ++r)
    {
        evaluate_basis(r, &dof_values, x, coordinate_dofs, cell_orientation);
        values[r] = dof_values;
    }
  }

  void evaluate_basis_derivatives(std::size_t i,
                                  std::size_t n,
                                  double * values,
                                  const double * x,
                                  const double * coordinate_dofs,
                                  int cell_orientation,
                                  const ufc::coordinate_mapping * cm=nullptr
                                  ) const final override
  {
    std::size_t num_derivatives = std::pow(2, n);
    std::fill_n(values, num_derivatives, 0.0);
    // Call evaluate_basis_all if order of derivatives is equal to zero.
    if (n == 0)
    {
        evaluate_basis(i, values, x, coordinate_dofs, cell_orientation);
        return;
    }
    // If order of derivatives is greater than the maximum polynomial degree, return zeros.
    if (n > 2)
        return;
    // Compute Jacobian
    double J[4];
    compute_jacobian_triangle_2d(J, coordinate_dofs);
    // Compute Inverse Jacobian and determinant
    double K[4];
    double detJ;
    compute_jacobian_inverse_triangle_2d(K, detJ, J);
    // Compute constants
    const double C0 = coordinate_dofs[2] + coordinate_dofs[4];
    const double C1 = coordinate_dofs[3] + coordinate_dofs[5];
    // Get coordinates and map to the reference (FIAT) element
    double Y[2] = { (J[1] * (C1 - 2.0 * x[1]) + J[3] * (2.0 * x[0] - C0)) / detJ, (J[0] * (2.0 * x[1] - C1) + J[2] * (C0 - 2.0 * x[0])) / detJ };
    // Precomputed combinations
    const std::size_t combinations[2][4][2] =
        { { { 0, 0 },
            { 1, 0 },
            { 0, 0 },
            { 1, 0 } },
          { { 0, 0 },
            { 0, 1 },
            { 1, 0 },
            { 1, 1 } } };
    // Declare transformation matrix
    double transform[4][4] =
        { { 1.0, 1.0, 1.0, 1.0 },
          { 1.0, 1.0, 1.0, 1.0 },
          { 1.0, 1.0, 1.0, 1.0 },
          { 1.0, 1.0, 1.0, 1.0 } };
    // Construct transformation matrix
    for (std::size_t row = 0; row < num_derivatives; ++row)
        for (std::size_t col = 0; col < num_derivatives; ++col)
            for (std::size_t k = 0; k < n; ++k)
                transform[row][col] *= K[2 * combinations[n - 1][col][k] + combinations[n - 1][row][k]];
    switch (i)
    {
    case 0:
        {
            double basisvalues[6] = {};
            basisvalues[0] = 1.0;
            const double tmp1_2 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_2;
            const double tmp2_2 = (1.0 - Y[1]) / 2.0;
            const double tmp3_2 = tmp2_2 * tmp2_2;
            basisvalues[3] = 1.5 * tmp1_2 * basisvalues[1] - 0.5 * tmp3_2 * basisvalues[0];
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            basisvalues[4] = (1.5 + 2.5 * Y[1]) * basisvalues[1];
            basisvalues[5] = (0.1111111111111111 + 1.6666666666666667 * Y[1]) * basisvalues[2] - 0.5555555555555556 * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.5);
            basisvalues[2] *= std::sqrt(1.0);
            basisvalues[5] *= std::sqrt(1.5);
            basisvalues[1] *= std::sqrt(3.0);
            basisvalues[4] *= std::sqrt(4.5);
            basisvalues[3] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[6] = { 0.0, -0.17320508075688776, -0.09999999999999999, 0.12171612389003691, 0.09428090415820634, 0.05443310539518172 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[6][6] =
                { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 4.898979485566347, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 9.486832980505149, 0.0, 0.0, 0.0, 0.0 },
                  { 3.999999999999994, 0.0, 7.0710678118654675, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[6][6] =
                { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 2.4494897427831805, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 4.242640687119282, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 2.581988897471618, 4.743416490252574, -0.91287092917528, 0.0, 0.0, 0.0 },
                  { 1.9999999999999152, 6.123724356957939, 3.5355339059327338, 0.0, 0.0, 0.0 },
                  { -2.3094010767584887, 0.0, 8.164965809277266, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[4] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[6][6] =
                { { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[6][6] =
                { { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };
            // Loop possible derivatives.
            for (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 36, 0.0);
                for (std::size_t t = 0; t < 6; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 36, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 36, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 6; ++s)
                    for (std::size_t t = 0; t < 6; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[r] += transform[r][s] * derivatives[s];
        }
        break;
    case 1:
        {
            double basisvalues[6] = {};
            basisvalues[0] = 1.0;
            const double tmp1_2 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_2;
            const double tmp2_2 = (1.0 - Y[1]) / 2.0;
            const double tmp3_2 = tmp2_2 * tmp2_2;
            basisvalues[3] = 1.5 * tmp1_2 * basisvalues[1] - 0.5 * tmp3_2 * basisvalues[0];
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            basisvalues[4] = (1.5 + 2.5 * Y[1]) * basisvalues[1];
            basisvalues[5] = (0.1111111111111111 + 1.6666666666666667 * Y[1]) * basisvalues[2] - 0.5555555555555556 * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.5);
            basisvalues[2] *= std::sqrt(1.0);
            basisvalues[5] *= std::sqrt(1.5);
            basisvalues[1] *= std::sqrt(3.0);
            basisvalues[4] *= std::sqrt(4.5);
            basisvalues[3] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[6] = { 0.0, 0.1732050807568877, -0.09999999999999999, 0.12171612389003691, -0.09428090415820636, 0.05443310539518175 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[6][6] =
                { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 4.898979485566347, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 9.486832980505149, 0.0, 0.0, 0.0, 0.0 },
                  { 3.999999999999994, 0.0, 7.0710678118654675, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[6][6] =
                { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 2.4494897427831805, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 4.242640687119282, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 2.581988897471618, 4.743416490252574, -0.91287092917528, 0.0, 0.0, 0.0 },
                  { 1.9999999999999152, 6.123724356957939, 3.5355339059327338, 0.0, 0.0, 0.0 },
                  { -2.3094010767584887, 0.0, 8.164965809277266, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[4] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[6][6] =
                { { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[6][6] =
                { { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };
            // Loop possible derivatives.
            for (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 36, 0.0);
                for (std::size_t t = 0; t < 6; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 36, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 36, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 6; ++s)
                    for (std::size_t t = 0; t < 6; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[r] += transform[r][s] * derivatives[s];
        }
        break;
    case 2:
        {
            double basisvalues[6] = {};
            basisvalues[0] = 1.0;
            const double tmp1_2 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_2;
            const double tmp2_2 = (1.0 - Y[1]) / 2.0;
            const double tmp3_2 = tmp2_2 * tmp2_2;
            basisvalues[3] = 1.5 * tmp1_2 * basisvalues[1] - 0.5 * tmp3_2 * basisvalues[0];
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            basisvalues[4] = (1.5 + 2.5 * Y[1]) * basisvalues[1];
            basisvalues[5] = (0.1111111111111111 + 1.6666666666666667 * Y[1]) * basisvalues[2] - 0.5555555555555556 * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.5);
            basisvalues[2] *= std::sqrt(1.0);
            basisvalues[5] *= std::sqrt(1.5);
            basisvalues[1] *= std::sqrt(3.0);
            basisvalues[4] *= std::sqrt(4.5);
            basisvalues[3] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[6] = { 0.0, 0.0, 0.20000000000000004, 0.0, 0.0, 0.16329931618554522 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[6][6] =
                { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 4.898979485566347, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 9.486832980505149, 0.0, 0.0, 0.0, 0.0 },
                  { 3.999999999999994, 0.0, 7.0710678118654675, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[6][6] =
                { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 2.4494897427831805, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 4.242640687119282, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 2.581988897471618, 4.743416490252574, -0.91287092917528, 0.0, 0.0, 0.0 },
                  { 1.9999999999999152, 6.123724356957939, 3.5355339059327338, 0.0, 0.0, 0.0 },
                  { -2.3094010767584887, 0.0, 8.164965809277266, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[4] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[6][6] =
                { { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[6][6] =
                { { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };
            // Loop possible derivatives.
            for (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 36, 0.0);
                for (std::size_t t = 0; t < 6; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 36, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 36, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 6; ++s)
                    for (std::size_t t = 0; t < 6; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[r] += transform[r][s] * derivatives[s];
        }
        break;
    case 3:
        {
            double basisvalues[6] = {};
            basisvalues[0] = 1.0;
            const double tmp1_2 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_2;
            const double tmp2_2 = (1.0 - Y[1]) / 2.0;
            const double tmp3_2 = tmp2_2 * tmp2_2;
            basisvalues[3] = 1.5 * tmp1_2 * basisvalues[1] - 0.5 * tmp3_2 * basisvalues[0];
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            basisvalues[4] = (1.5 + 2.5 * Y[1]) * basisvalues[1];
            basisvalues[5] = (0.1111111111111111 + 1.6666666666666667 * Y[1]) * basisvalues[2] - 0.5555555555555556 * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.5);
            basisvalues[2] *= std::sqrt(1.0);
            basisvalues[5] *= std::sqrt(1.5);
            basisvalues[1] *= std::sqrt(3.0);
            basisvalues[4] *= std::sqrt(4.5);
            basisvalues[3] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[6] = { 0.47140452079103157, 0.23094010767585027, 0.1333333333333333, 0.0, 0.18856180831641264, -0.1632993161855452 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[6][6] =
                { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 4.898979485566347, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 9.486832980505149, 0.0, 0.0, 0.0, 0.0 },
                  { 3.999999999999994, 0.0, 7.0710678118654675, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[6][6] =
                { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 2.4494897427831805, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 4.242640687119282, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 2.581988897471618, 4.743416490252574, -0.91287092917528, 0.0, 0.0, 0.0 },
                  { 1.9999999999999152, 6.123724356957939, 3.5355339059327338, 0.0, 0.0, 0.0 },
                  { -2.3094010767584887, 0.0, 8.164965809277266, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[4] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[6][6] =
                { { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[6][6] =
                { { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };
            // Loop possible derivatives.
            for (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 36, 0.0);
                for (std::size_t t = 0; t < 6; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 36, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 36, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 6; ++s)
                    for (std::size_t t = 0; t < 6; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[r] += transform[r][s] * derivatives[s];
        }
        break;
    case 4:
        {
            double basisvalues[6] = {};
            basisvalues[0] = 1.0;
            const double tmp1_2 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_2;
            const double tmp2_2 = (1.0 - Y[1]) / 2.0;
            const double tmp3_2 = tmp2_2 * tmp2_2;
            basisvalues[3] = 1.5 * tmp1_2 * basisvalues[1] - 0.5 * tmp3_2 * basisvalues[0];
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            basisvalues[4] = (1.5 + 2.5 * Y[1]) * basisvalues[1];
            basisvalues[5] = (0.1111111111111111 + 1.6666666666666667 * Y[1]) * basisvalues[2] - 0.5555555555555556 * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.5);
            basisvalues[2] *= std::sqrt(1.0);
            basisvalues[5] *= std::sqrt(1.5);
            basisvalues[1] *= std::sqrt(3.0);
            basisvalues[4] *= std::sqrt(4.5);
            basisvalues[3] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[6] = { 0.47140452079103157, -0.23094010767585027, 0.13333333333333333, 0.0, -0.18856180831641264, -0.16329931618554522 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[6][6] =
                { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 4.898979485566347, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 9.486832980505149, 0.0, 0.0, 0.0, 0.0 },
                  { 3.999999999999994, 0.0, 7.0710678118654675, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[6][6] =
                { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 2.4494897427831805, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 4.242640687119282, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 2.581988897471618, 4.743416490252574, -0.91287092917528, 0.0, 0.0, 0.0 },
                  { 1.9999999999999152, 6.123724356957939, 3.5355339059327338, 0.0, 0.0, 0.0 },
                  { -2.3094010767584887, 0.0, 8.164965809277266, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[4] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[6][6] =
                { { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[6][6] =
                { { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };
            // Loop possible derivatives.
            for (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 36, 0.0);
                for (std::size_t t = 0; t < 6; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 36, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 36, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 6; ++s)
                    for (std::size_t t = 0; t < 6; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[r] += transform[r][s] * derivatives[s];
        }
        break;
    case 5:
        {
            double basisvalues[6] = {};
            basisvalues[0] = 1.0;
            const double tmp1_2 = (1.0 + 2.0 * Y[0] + Y[1]) / 2.0;
            basisvalues[1] = tmp1_2;
            const double tmp2_2 = (1.0 - Y[1]) / 2.0;
            const double tmp3_2 = tmp2_2 * tmp2_2;
            basisvalues[3] = 1.5 * tmp1_2 * basisvalues[1] - 0.5 * tmp3_2 * basisvalues[0];
            basisvalues[2] = (0.5 + 1.5 * Y[1]) * basisvalues[0];
            basisvalues[4] = (1.5 + 2.5 * Y[1]) * basisvalues[1];
            basisvalues[5] = (0.1111111111111111 + 1.6666666666666667 * Y[1]) * basisvalues[2] - 0.5555555555555556 * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.5);
            basisvalues[2] *= std::sqrt(1.0);
            basisvalues[5] *= std::sqrt(1.5);
            basisvalues[1] *= std::sqrt(3.0);
            basisvalues[4] *= std::sqrt(4.5);
            basisvalues[3] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[6] = { 0.4714045207910316, 0.0, -0.26666666666666666, -0.24343224778007383, 0.0, 0.05443310539518172 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[6][6] =
                { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 4.898979485566347, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 9.486832980505149, 0.0, 0.0, 0.0, 0.0 },
                  { 3.999999999999994, 0.0, 7.0710678118654675, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[6][6] =
                { { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 2.4494897427831805, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 4.242640687119282, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 2.581988897471618, 4.743416490252574, -0.91287092917528, 0.0, 0.0, 0.0 },
                  { 1.9999999999999152, 6.123724356957939, 3.5355339059327338, 0.0, 0.0, 0.0 },
                  { -2.3094010767584887, 0.0, 8.164965809277266, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[4] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[6][6] =
                { { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[6][6] =
                { { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 } };
            // Loop possible derivatives.
            for (std::size_t r = 0; r < num_derivatives; ++r)
            {
                // Reset dmats to identity
                std::fill_n(&dmats[0][0], 36, 0.0);
                for (std::size_t t = 0; t < 6; ++t)
                    dmats[t][t] = 1.0;
                // Looping derivative order to generate dmats.
                for (std::size_t s = 0; s < n; ++s)
                {
                    std::copy_n(&dmats[0][0], 36, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 36, 0.0);
                    // Update dmats using an inner product.
                    // _dmats_product(shape_dmats, comb[r][s], 0)
                    if (combinations[n - 1][r][s] == 0)
                    {
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 1)
                    if (combinations[n - 1][r][s] == 1)
                    {
                        for (std::size_t t = 0; t < 6; ++t)
                            for (std::size_t u = 0; u < 6; ++u)
                                for (std::size_t tu = 0; tu < 6; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 6; ++s)
                    for (std::size_t t = 0; t < 6; ++t)
                        derivatives[r] += coefficients0[s] * dmats[s][t] * basisvalues[t];
            }
            // Transform derivatives back to physical element
            for (std::size_t r = 0; r < num_derivatives; ++r)
                for (std::size_t s = 0; s < num_derivatives; ++s)
                    values[r] += transform[r][s] * derivatives[s];
        }
        break;
    }
  }

  void evaluate_basis_derivatives_all(std::size_t n,
                                      double * values,
                                      const double * x,
                                      const double * coordinate_dofs,
                                      int cell_orientation,
                                      const ufc::coordinate_mapping * cm=nullptr
                                      ) const final override
  {
    // Call evaluate_basis_all if order of derivatives is equal to zero.
    if (n == 0)
    {
        evaluate_basis_all(values, x, coordinate_dofs, cell_orientation);
        return;
    }
    unsigned int num_derivatives = std::pow(2, n);
    // Set values equal to zero.
    std::fill_n(values, num_derivatives * 6, 0.0);
    // If order of derivatives is greater than the maximum polynomial degree, return zeros.
    if (n > 2)
        return;
    // Helper variable to hold values of a single dof.
    double dof_values[4] = {};
    // Loop dofs and call evaluate_basis_derivatives.
    for (std::size_t r = 0; r < 6; ++r)
    {
        evaluate_basis_derivatives(r, n, dof_values, x, coordinate_dofs, cell_orientation);
        for (std::size_t s = 0; s < num_derivatives; ++s)
            values[num_derivatives * r + s] = dof_values[s];
    }
  }

  double evaluate_dof(std::size_t i,
                      const ufc::function& f,
                      const double * coordinate_dofs,
                      int cell_orientation,
                      const ufc::cell& c,
                      const ufc::coordinate_mapping * cm=nullptr
                      ) const final override
  {
    // Declare variables for result of evaluation
    double vals[1];
    // Declare variable for physical coordinates
    double y[2];
    switch (i)
    {
    case 0:
        {
            y[0] = coordinate_dofs[0];
            y[1] = coordinate_dofs[1];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 1:
        {
            y[0] = coordinate_dofs[2];
            y[1] = coordinate_dofs[3];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 2:
        {
            y[0] = coordinate_dofs[4];
            y[1] = coordinate_dofs[5];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 3:
        {
            y[0] = 0.5 * coordinate_dofs[2] + 0.5 * coordinate_dofs[4];
            y[1] = 0.5 * coordinate_dofs[3] + 0.5 * coordinate_dofs[5];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 4:
        {
            y[0] = 0.5 * coordinate_dofs[0] + 0.5 * coordinate_dofs[4];
            y[1] = 0.5 * coordinate_dofs[1] + 0.5 * coordinate_dofs[5];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 5:
        {
            y[0] = 0.5 * coordinate_dofs[0] + 0.5 * coordinate_dofs[2];
            y[1] = 0.5 * coordinate_dofs[1] + 0.5 * coordinate_dofs[3];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    }
    return 0.0;
  }

  void evaluate_dofs(double * values,
                             const ufc::function& f,
                             const double * coordinate_dofs,
                             int cell_orientation,
                             const ufc::cell& c,
                             const ufc::coordinate_mapping * cm=nullptr
                             ) const final override
  {
    // Declare variables for result of evaluation
    double vals[1];
    // Declare variable for physical coordinates
    double y[2];
    y[0] = coordinate_dofs[0];
    y[1] = coordinate_dofs[1];
    f.evaluate(vals, y, c);
    values[0] = vals[0];
    y[0] = coordinate_dofs[2];
    y[1] = coordinate_dofs[3];
    f.evaluate(vals, y, c);
    values[1] = vals[0];
    y[0] = coordinate_dofs[4];
    y[1] = coordinate_dofs[5];
    f.evaluate(vals, y, c);
    values[2] = vals[0];
    y[0] = 0.5 * coordinate_dofs[2] + 0.5 * coordinate_dofs[4];
    y[1] = 0.5 * coordinate_dofs[3] + 0.5 * coordinate_dofs[5];
    f.evaluate(vals, y, c);
    values[3] = vals[0];
    y[0] = 0.5 * coordinate_dofs[0] + 0.5 * coordinate_dofs[4];
    y[1] = 0.5 * coordinate_dofs[1] + 0.5 * coordinate_dofs[5];
    f.evaluate(vals, y, c);
    values[4] = vals[0];
    y[0] = 0.5 * coordinate_dofs[0] + 0.5 * coordinate_dofs[2];
    y[1] = 0.5 * coordinate_dofs[1] + 0.5 * coordinate_dofs[3];
    f.evaluate(vals, y, c);
    values[5] = vals[0];
  }

  void interpolate_vertex_values(double * vertex_values,
                                 const double * dof_values,
                                 const double * coordinate_dofs,
                                 int cell_orientation,
                                 const ufc::coordinate_mapping * cm=nullptr
                                 ) const final override
  {
    // Evaluate function and change variables
    vertex_values[0] = dof_values[0];
    vertex_values[1] = dof_values[1];
    vertex_values[2] = dof_values[2];
  }

  void tabulate_dof_coordinates(double * dof_coordinates,
                                const double * coordinate_dofs,
                                const ufc::coordinate_mapping * cm=nullptr
                                ) const final override
  {
    dof_coordinates[0] = coordinate_dofs[0];
    dof_coordinates[1] = coordinate_dofs[1];
    dof_coordinates[2] = coordinate_dofs[2];
    dof_coordinates[2 + 1] = coordinate_dofs[3];
    dof_coordinates[2 * 2] = coordinate_dofs[4];
    dof_coordinates[2 * 2 + 1] = coordinate_dofs[5];
    dof_coordinates[2 * 3] = 0.5 * coordinate_dofs[2] + 0.5 * coordinate_dofs[4];
    dof_coordinates[2 * 3 + 1] = 0.5 * coordinate_dofs[3] + 0.5 * coordinate_dofs[5];
    dof_coordinates[2 * 4] = 0.5 * coordinate_dofs[0] + 0.5 * coordinate_dofs[4];
    dof_coordinates[2 * 4 + 1] = 0.5 * coordinate_dofs[1] + 0.5 * coordinate_dofs[5];
    dof_coordinates[2 * 5] = 0.5 * coordinate_dofs[0] + 0.5 * coordinate_dofs[2];
    dof_coordinates[2 * 5 + 1] = 0.5 * coordinate_dofs[1] + 0.5 * coordinate_dofs[3];
  }

  void tabulate_reference_dof_coordinates(double * reference_dof_coordinates) const final override
  {
    static const double dof_X[12] = { 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.5, 0.0, 0.5, 0.5, 0.0 };
    std::copy_n(dof_X, 12, reference_dof_coordinates);
  }

  std::size_t num_sub_elements() const final override
  {
    return 0;
  }

  ufc::finite_element * create_sub_element(std::size_t i) const final override
  {
    return nullptr;
  }

  ufc::finite_element * create() const final override
  {
    return new pressure_finite_element_2();
  }

};


class pressure_dofmap_0: public ufc::dofmap
{
public:

  pressure_dofmap_0() : ufc::dofmap()
  {
    // Do nothing
  }

  ~pressure_dofmap_0() override
  {
    // Do nothing
  }

  const char * signature() const final override
  {
    return "FFC dofmap for FiniteElement('Lagrange', triangle, 1)";
  }

  bool needs_mesh_entities(std::size_t d) const final override
  {
    static const bool return_values[3] = { true, false, false };
    if (d >= 3)
        return false;
    return return_values[d];
  }

  std::size_t topological_dimension() const final override
  {
    return 2;
  }

  std::size_t global_dimension(const std::vector<std::size_t>&
                               num_global_entities) const final override
  {
    return num_global_entities[0];
  }

  std::size_t num_global_support_dofs() const final override
  {
    return 0;
  }

  std::size_t num_element_support_dofs() const final override
  {
    return 3;
  }

  std::size_t num_element_dofs() const final override
  {
    return 3;
  }

  std::size_t num_facet_dofs() const final override
  {
    return 2;
  }

  std::size_t num_entity_dofs(std::size_t d) const final override
  {
    static const std::size_t return_values[3] = { 1, 0, 0 };
    if (d >= 3)
        return 0;
    return return_values[d];
  }

  std::size_t num_entity_closure_dofs(std::size_t d) const final override
  {
    static const std::size_t return_values[3] = { 1, 2, 3 };
    if (d >= 3)
        return 0;
    return return_values[d];
  }

  void tabulate_dofs(std::size_t * dofs,
                     const std::vector<std::size_t>& num_global_entities,
                     const std::vector<std::vector<std::size_t>>& entity_indices) const final override
  {
    dofs[0] = entity_indices[0][0];
    dofs[1] = entity_indices[0][1];
    dofs[2] = entity_indices[0][2];
  }

  void tabulate_facet_dofs(std::size_t * dofs,
                           std::size_t facet) const final override
  {
    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;
    }
  }

  void tabulate_entity_dofs(std::size_t * dofs,
                            std::size_t d, std::size_t i) const final override
  {
    switch (d)
    {
    case 0:
        switch (i)
        {
        case 0:
            dofs[0] = 0;
            break;
        case 1:
            dofs[0] = 1;
            break;
        case 2:
            dofs[0] = 2;
            break;
        }
        break;
    }
  }

  void tabulate_entity_closure_dofs(std::size_t * dofs,
                                    std::size_t d, std::size_t i) const final override
  {
    switch (d)
    {
    case 0:
        switch (i)
        {
        case 0:
            dofs[0] = 0;
            break;
        case 1:
            dofs[0] = 1;
            break;
        case 2:
            dofs[0] = 2;
            break;
        }
        break;
    case 1:
        switch (i)
        {
        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;
        }
        break;
    case 2:
        switch (i)
        {
        case 0:
            dofs[0] = 0;
            dofs[1] = 1;
            dofs[2] = 2;
            break;
        }
        break;
    }
  }


  std::size_t num_sub_dofmaps() const final override
  {
    return 0;
  }

  ufc::dofmap * create_sub_dofmap(std::size_t i) const final override
  {
    return nullptr;
  }

  ufc::dofmap * create() const final override
  {
    return new pressure_dofmap_0();
  }

};


class pressure_dofmap_1: public ufc::dofmap
{
public:

  pressure_dofmap_1() : ufc::dofmap()
  {
    // Do nothing
  }

  ~pressure_dofmap_1() override
  {
    // Do nothing
  }

  const char * signature() const final override
  {
    return "FFC dofmap for VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2)";
  }

  bool needs_mesh_entities(std::size_t d) const final override
  {
    static const bool return_values[3] = { true, false, false };
    if (d >= 3)
        return false;
    return return_values[d];
  }

  std::size_t topological_dimension() const final override
  {
    return 2;
  }

  std::size_t global_dimension(const std::vector<std::size_t>&
                               num_global_entities) const final override
  {
    return 2 * num_global_entities[0];
  }

  std::size_t num_global_support_dofs() const final override
  {
    return 0;
  }

  std::size_t num_element_support_dofs() const final override
  {
    return 6;
  }

  std::size_t num_element_dofs() const final override
  {
    return 6;
  }

  std::size_t num_facet_dofs() const final override
  {
    return 4;
  }

  std::size_t num_entity_dofs(std::size_t d) const final override
  {
    static const std::size_t return_values[3] = { 2, 0, 0 };
    if (d >= 3)
        return 0;
    return return_values[d];
  }

  std::size_t num_entity_closure_dofs(std::size_t d) const final override
  {
    static const std::size_t return_values[3] = { 2, 4, 6 };
    if (d >= 3)
        return 0;
    return return_values[d];
  }

  void tabulate_dofs(std::size_t * dofs,
                     const std::vector<std::size_t>& num_global_entities,
                     const std::vector<std::vector<std::size_t>>& entity_indices) const final override
  {
    std::size_t offset = 0;
    dofs[0] = offset + entity_indices[0][0];
    dofs[1] = offset + entity_indices[0][1];
    dofs[2] = offset + entity_indices[0][2];
    offset += num_global_entities[0];
    dofs[3] = offset + entity_indices[0][0];
    dofs[4] = offset + entity_indices[0][1];
    dofs[5] = offset + entity_indices[0][2];
    offset += num_global_entities[0];
  }

  void tabulate_facet_dofs(std::size_t * dofs,
                           std::size_t facet) const final override
  {
    switch (facet)
    {
    case 0:
        dofs[0] = 1;
        dofs[1] = 2;
        dofs[2] = 4;
        dofs[3] = 5;
        break;
    case 1:
        dofs[0] = 0;
        dofs[1] = 2;
        dofs[2] = 3;
        dofs[3] = 5;
        break;
    case 2:
        dofs[0] = 0;
        dofs[1] = 1;
        dofs[2] = 3;
        dofs[3] = 4;
        break;
    }
  }

  void tabulate_entity_dofs(std::size_t * dofs,
                            std::size_t d, std::size_t i) const final override
  {
    switch (d)
    {
    case 0:
        switch (i)
        {
        case 0:
            dofs[0] = 0;
            dofs[1] = 3;
            break;
        case 1:
            dofs[0] = 1;
            dofs[1] = 4;
            break;
        case 2:
            dofs[0] = 2;
            dofs[1] = 5;
            break;
        }
        break;
    }
  }

  void tabulate_entity_closure_dofs(std::size_t * dofs,
                                    std::size_t d, std::size_t i) const final override
  {
    switch (d)
    {
    case 0:
        switch (i)
        {
        case 0:
            dofs[0] = 0;
            dofs[1] = 3;
            break;
        case 1:
            dofs[0] = 1;
            dofs[1] = 4;
            break;
        case 2:
            dofs[0] = 2;
            dofs[1] = 5;
            break;
        }
        break;
    case 1:
        switch (i)
        {
        case 0:
            dofs[0] = 1;
            dofs[1] = 2;
            dofs[2] = 4;
            dofs[3] = 5;
            break;
        case 1:
            dofs[0] = 0;
            dofs[1] = 2;
            dofs[2] = 3;
            dofs[3] = 5;
            break;
        case 2:
            dofs[0] = 0;
            dofs[1] = 1;
            dofs[2] = 3;
            dofs[3] = 4;
            break;
        }
        break;
    case 2:
        switch (i)
        {
        case 0:
            dofs[0] = 0;
            dofs[1] = 1;
            dofs[2] = 2;
            dofs[3] = 3;
            dofs[4] = 4;
            dofs[5] = 5;
            break;
        }
        break;
    }
  }


  std::size_t num_sub_dofmaps() const final override
  {
    return 2;
  }

  ufc::dofmap * create_sub_dofmap(std::size_t i) const final override
  {
    switch (i)
    {
    case 0:
        return new pressure_dofmap_0();
    case 1:
        return new pressure_dofmap_0();
    default:
        return nullptr;
    }
  }

  ufc::dofmap * create() const final override
  {
    return new pressure_dofmap_1();
  }

};


class pressure_dofmap_2: public ufc::dofmap
{
public:

  pressure_dofmap_2() : ufc::dofmap()
  {
    // Do nothing
  }

  ~pressure_dofmap_2() override
  {
    // Do nothing
  }

  const char * signature() const final override
  {
    return "FFC dofmap for FiniteElement('Lagrange', triangle, 2)";
  }

  bool needs_mesh_entities(std::size_t d) const final override
  {
    static const bool return_values[3] = { true, true, false };
    if (d >= 3)
        return false;
    return return_values[d];
  }

  std::size_t topological_dimension() const final override
  {
    return 2;
  }

  std::size_t global_dimension(const std::vector<std::size_t>&
                               num_global_entities) const final override
  {
    return num_global_entities[0] + num_global_entities[1];
  }

  std::size_t num_global_support_dofs() const final override
  {
    return 0;
  }

  std::size_t num_element_support_dofs() const final override
  {
    return 6;
  }

  std::size_t num_element_dofs() const final override
  {
    return 6;
  }

  std::size_t num_facet_dofs() const final override
  {
    return 3;
  }

  std::size_t num_entity_dofs(std::size_t d) const final override
  {
    static const std::size_t return_values[3] = { 1, 1, 0 };
    if (d >= 3)
        return 0;
    return return_values[d];
  }

  std::size_t num_entity_closure_dofs(std::size_t d) const final override
  {
    static const std::size_t return_values[3] = { 1, 3, 6 };
    if (d >= 3)
        return 0;
    return return_values[d];
  }

  void tabulate_dofs(std::size_t * dofs,
                     const std::vector<std::size_t>& num_global_entities,
                     const std::vector<std::vector<std::size_t>>& entity_indices) const final override
  {
    std::size_t offset = 0;
    dofs[0] = offset + entity_indices[0][0];
    dofs[1] = offset + entity_indices[0][1];
    dofs[2] = offset + entity_indices[0][2];
    offset += num_global_entities[0];
    dofs[3] = offset + entity_indices[1][0];
    dofs[4] = offset + entity_indices[1][1];
    dofs[5] = offset + entity_indices[1][2];
    offset += num_global_entities[1];
  }

  void tabulate_facet_dofs(std::size_t * dofs,
                           std::size_t facet) const final override
  {
    switch (facet)
    {
    case 0:
        dofs[0] = 1;
        dofs[1] = 2;
        dofs[2] = 3;
        break;
    case 1:
        dofs[0] = 0;
        dofs[1] = 2;
        dofs[2] = 4;
        break;
    case 2:
        dofs[0] = 0;
        dofs[1] = 1;
        dofs[2] = 5;
        break;
    }
  }

  void tabulate_entity_dofs(std::size_t * dofs,
                            std::size_t d, std::size_t i) const final override
  {
    switch (d)
    {
    case 0:
        switch (i)
        {
        case 0:
            dofs[0] = 0;
            break;
        case 1:
            dofs[0] = 1;
            break;
        case 2:
            dofs[0] = 2;
            break;
        }
        break;
    case 1:
        switch (i)
        {
        case 0:
            dofs[0] = 3;
            break;
        case 1:
            dofs[0] = 4;
            break;
        case 2:
            dofs[0] = 5;
            break;
        }
        break;
    }
  }

  void tabulate_entity_closure_dofs(std::size_t * dofs,
                                    std::size_t d, std::size_t i) const final override
  {
    switch (d)
    {
    case 0:
        switch (i)
        {
        case 0:
            dofs[0] = 0;
            break;
        case 1:
            dofs[0] = 1;
            break;
        case 2:
            dofs[0] = 2;
            break;
        }
        break;
    case 1:
        switch (i)
        {
        case 0:
            dofs[0] = 1;
            dofs[1] = 2;
            dofs[2] = 3;
            break;
        case 1:
            dofs[0] = 0;
            dofs[1] = 2;
            dofs[2] = 4;
            break;
        case 2:
            dofs[0] = 0;
            dofs[1] = 1;
            dofs[2] = 5;
            break;
        }
        break;
    case 2:
        switch (i)
        {
        case 0:
            dofs[0] = 0;
            dofs[1] = 1;
            dofs[2] = 2;
            dofs[3] = 3;
            dofs[4] = 4;
            dofs[5] = 5;
            break;
        }
        break;
    }
  }


  std::size_t num_sub_dofmaps() const final override
  {
    return 0;
  }

  ufc::dofmap * create_sub_dofmap(std::size_t i) const final override
  {
    return nullptr;
  }

  ufc::dofmap * create() const final override
  {
    return new pressure_dofmap_2();
  }

};


class pressure_coordinate_mapping_1: public ufc::coordinate_mapping
{
public:

  pressure_coordinate_mapping_1() : ufc::coordinate_mapping()
  {
    // Do nothing
  }

  ~pressure_coordinate_mapping_1() override
  {
    // Do nothing
  }

  const char * signature() const final override
  {
    return "FFC coordinate_mapping from VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2)";
  }

  ufc::coordinate_mapping * create() const final override
  {
    return new pressure_coordinate_mapping_1();
  }

  std::size_t geometric_dimension() const final override
  {
    return 2;
  }

  std::size_t topological_dimension() const final override
  {
    return 2;
  }

  ufc::shape cell_shape() const final override
  {
    return ufc::shape::triangle;
  }

  ufc::finite_element * create_coordinate_finite_element() const final override
  {
    return new pressure_finite_element_1();
  }

  ufc::dofmap * create_coordinate_dofmap() const final override
  {
    return new pressure_dofmap_1();
  }

  void compute_physical_coordinates(
      double * x, std::size_t num_points,
      const double * X,
      const double * coordinate_dofs) const final override
  {
    pressure_finite_element_0 xelement;
    double phi[3];
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Compute basis values of coordinate element
        xelement.evaluate_reference_basis(phi, 1, &X[2 * ip]);
        // Compute x
        for (std::size_t i = 0; i < 2; ++i)
            for (std::size_t d = 0; d < 3; ++d)
                x[2 * ip + i] += coordinate_dofs[2 * d + i] * phi[d];
    }
  }

  void compute_reference_coordinates(
      double * X, std::size_t num_points,
      const double * x,
      const double * coordinate_dofs, int cell_orientation) const final override
  {
    const double phi_X0[3] = { 1.0, 0.0, 0.0 };
    const double dphi_X0[2][3] =
        { { -1.0, 1.0, 0.0 },
          { -1.0, 0.0, 1.0 } };
    double J[4];
    double detJ[1];
    double K[4];
    double x0[2] = {};
    for (std::size_t i = 0; i < 2; ++i)
        for (std::size_t k = 0; k < 3; ++k)
            x0[i] += coordinate_dofs[2 * k + i] * phi_X0[k];
    for (std::size_t i = 0; i < 2; ++i)
    {
        for (std::size_t j = 0; j < 2; ++j)
        {
            J[2 * i + j] = 0.0;
            for (std::size_t k = 0; k < 3; ++k)
                J[2 * i + j] += coordinate_dofs[2 * k + i] * dphi_X0[j][k];
        }
    }
    compute_jacobian_determinants(detJ, 1, J, cell_orientation);
    compute_jacobian_inverses(K, 1, J, detJ);
    for (std::size_t ip = 0; ip < num_points; ++ip)
        for (std::size_t j = 0; j < 2; ++j)
            for (std::size_t i = 0; i < 2; ++i)
                X[2 * ip + j] += K[2 * j + i] * (x[2 * ip + i] - x0[i]);
  }

  void compute_reference_geometry(
      double * X, double * J, double * detJ, double * K, std::size_t num_points,
      const double * x,
      const double * coordinate_dofs, int cell_orientation) const final override
  {
    const double phi_X0[3] = { 1.0, 0.0, 0.0 };
    const double dphi_X0[2][3] =
        { { -1.0, 1.0, 0.0 },
          { -1.0, 0.0, 1.0 } };
    double x0[2] = {};
    for (std::size_t i = 0; i < 2; ++i)
        for (std::size_t k = 0; k < 3; ++k)
            x0[i] += coordinate_dofs[2 * k + i] * phi_X0[k];
    for (std::size_t i = 0; i < 2; ++i)
    {
        for (std::size_t j = 0; j < 2; ++j)
        {
            J[2 * i + j] = 0.0;
            for (std::size_t k = 0; k < 3; ++k)
                J[2 * i + j] += coordinate_dofs[2 * k + i] * dphi_X0[j][k];
        }
    }
    compute_jacobian_determinants(detJ, 1, J, cell_orientation);
    compute_jacobian_inverses(K, 1, J, detJ);
    for (std::size_t ip = 0; ip < num_points; ++ip)
        for (std::size_t j = 0; j < 2; ++j)
            for (std::size_t i = 0; i < 2; ++i)
                X[2 * ip + j] += K[2 * j + i] * (x[2 * ip + i] - x0[i]);
  }

  void compute_jacobians(
      double * J, std::size_t num_points,
      const double * X,
      const double * coordinate_dofs) const final override
  {
    pressure_finite_element_0 xelement;
    double dphi[6];
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Compute basis derivatives of coordinate element
        xelement.evaluate_reference_basis_derivatives(dphi, 1, 1, &X[2 * ip]);
        // Compute J
        for (std::size_t i = 0; i < 2; ++i)
            for (std::size_t j = 0; j < 2; ++j)
                for (std::size_t d = 0; d < 3; ++d)
                    J[2 * 2 * ip + 2 * i + j] += coordinate_dofs[2 * d + i] * dphi[2 * d + j];
    }
  }

  void compute_jacobian_determinants(
      double * detJ, std::size_t num_points,
      const double * J,
      int cell_orientation) const final override
  {
    for (std::size_t ip = 0; ip < num_points; ++ip)
        detJ[ip] = J[2 * 2 * ip] * J[2 * 2 * ip + 2 + 1] - J[2 * 2 * ip + 1] * J[2 * 2 * ip + 2];
  }

  void compute_jacobian_inverses(
      double * K, std::size_t num_points,
      const double * J, const double * detJ) const final override
  {
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        K[2 * 2 * ip] = J[2 * 2 * ip + 2 + 1] / detJ[ip];
        K[2 * 2 * ip + 1] = -J[2 * 2 * ip + 1] / detJ[ip];
        K[2 * 2 * ip + 2] = -J[2 * 2 * ip + 2] / detJ[ip];
        K[2 * 2 * ip + 2 + 1] = J[2 * 2 * ip] / detJ[ip];
    }
  }

  void compute_geometry(
      double * x, double * J, double * detJ, double * K, std::size_t num_points,
      const double * X,
      const double * coordinate_dofs, int cell_orientation) const final override
  {
    compute_physical_coordinates(x, num_points, X, coordinate_dofs);
    compute_jacobians(J, num_points, X, coordinate_dofs);
    compute_jacobian_determinants(detJ, num_points, J, cell_orientation);
    compute_jacobian_inverses(K, num_points, J, detJ);
  }

  void compute_midpoint_geometry(
      double * x, double * J,
      const double * coordinate_dofs) const final override
  {
    const double phi_Xm[3] = { 0.33333333333333337, 0.33333333333333337, 0.3333333333333333 };
    const double dphi_Xm[2][3] =
        { { -1.0, 1.0, 0.0 },
          { -1.0, 0.0, 1.0 } };
    // Compute x
    for (std::size_t i = 0; i < 2; ++i)
        for (std::size_t d = 0; d < 3; ++d)
            x[i] += coordinate_dofs[2 * d + i] * phi_Xm[d];
    // Compute J
    for (std::size_t i = 0; i < 2; ++i)
        for (std::size_t j = 0; j < 2; ++j)
            for (std::size_t d = 0; d < 3; ++d)
                J[2 * i + j] += coordinate_dofs[2 * d + i] * dphi_Xm[j][d];
  }

};


class pressure_cell_integral_0_otherwise: public ufc::cell_integral
{
public:

  pressure_cell_integral_0_otherwise() : ufc::cell_integral()
  {

  }

  ~pressure_cell_integral_0_otherwise() override
  {

  }

  const std::vector<bool> & enabled_coefficients() const final override
  {
static const std::vector<bool> enabled({});
return enabled;
  }

  void tabulate_tensor(double * A,
                       const double * const * w,
                       const double * coordinate_dofs,
                       int cell_orientation) const final override
  {
    // This function was generated using 'uflacs' representation
    // with the following integrals metadata:
    // 
    // num_cells:         None
    // optimize:          True
    // precision:         16
    // quadrature_degree: 2
    // quadrature_rule:   'default'
    // representation:    'uflacs'
    // 
    // and the following integral 0 metadata:
    // 
    // estimated_polynomial_degree: 2
    // optimize:                    True
    // precision:                   16
    // quadrature_degree:           2
    // quadrature_rule:             'default'
    // representation:              'uflacs'
    // Precomputed values of basis functions and precomputations
    // FE* dimensions: [entities][points][dofs]
    // PI* dimensions: [entities][dofs][dofs] or [entities][dofs]
    // PM* dimensions: [entities][dofs][dofs]
    alignas(32) static const double FE3_C0_D01_Q3[1][1][2] = { { { -1.0, 1.0 } } };
    // Unstructured piecewise computations
    const double J_c0 = coordinate_dofs[0] * FE3_C0_D01_Q3[0][0][0] + coordinate_dofs[2] * FE3_C0_D01_Q3[0][0][1];
    const double J_c3 = coordinate_dofs[1] * FE3_C0_D01_Q3[0][0][0] + coordinate_dofs[5] * FE3_C0_D01_Q3[0][0][1];
    const double J_c1 = coordinate_dofs[0] * FE3_C0_D01_Q3[0][0][0] + coordinate_dofs[4] * FE3_C0_D01_Q3[0][0][1];
    const double J_c2 = coordinate_dofs[1] * FE3_C0_D01_Q3[0][0][0] + coordinate_dofs[3] * FE3_C0_D01_Q3[0][0][1];
    alignas(32) double sp[20];
    sp[0] = J_c0 * J_c3;
    sp[1] = J_c1 * J_c2;
    sp[2] = sp[0] + -1 * sp[1];
    sp[3] = J_c0 / sp[2];
    sp[4] = -1 * J_c1 / sp[2];
    sp[5] = sp[3] * sp[3];
    sp[6] = sp[3] * sp[4];
    sp[7] = sp[4] * sp[4];
    sp[8] = J_c3 / sp[2];
    sp[9] = -1 * J_c2 / sp[2];
    sp[10] = sp[9] * sp[9];
    sp[11] = sp[8] * sp[9];
    sp[12] = sp[8] * sp[8];
    sp[13] = sp[5] + sp[10];
    sp[14] = sp[6] + sp[11];
    sp[15] = sp[12] + sp[7];
    sp[16] = std::abs(sp[2]);
    sp[17] = sp[13] * sp[16];
    sp[18] = sp[14] * sp[16];
    sp[19] = sp[15] * sp[16];
    A[0] = 0.5 * sp[19] + 0.5 * sp[18] + 0.5 * sp[18] + 0.5 * sp[17];
    A[1] = 0.1666666666666679 * sp[19] + 0.1666666666666667 * sp[18];
    A[2] = 0.1666666666666664 * sp[18] + 0.1666666666666656 * sp[17];
    A[3] = 0.0;
    A[4] = -0.6666666666666683 * sp[18] + -0.6666666666666682 * sp[17];
    A[5] = -0.6666666666666682 * sp[19] + -0.6666666666666675 * sp[18];
    A[6] = 0.1666666666666679 * sp[19] + 0.1666666666666667 * sp[18];
    A[7] = 0.5 * sp[19];
    A[8] = -0.1666666666666665 * sp[18];
    A[9] = 0.6666666666666637 * sp[18];
    A[10] = 0.0;
    A[11] = -0.6666666666666681 * sp[19] + -0.6666666666666666 * sp[18];
    A[12] = 0.1666666666666664 * sp[18] + 0.1666666666666656 * sp[17];
    A[13] = -0.1666666666666665 * sp[18];
    A[14] = 0.5 * sp[17];
    A[15] = 0.6666666666666665 * sp[18];
    A[16] = -0.6666666666666665 * sp[18] + -0.6666666666666649 * sp[17];
    std::fill(A + 17, A + 19, 0.0);
    A[19] = 0.6666666666666637 * sp[18];
    A[20] = 0.6666666666666665 * sp[18];
    A[21] = 1.333333333333328 * sp[19] + 0.6666666666666556 * sp[18] + 0.6666666666666556 * sp[18] + 1.333333333333315 * sp[17];
    A[22] = -1.333333333333328 * sp[19] + -0.6666666666666605 * sp[18] + -0.6666666666666556 * sp[18];
    A[23] = -0.6666666666666645 * sp[18] + -0.6666666666666662 * sp[18] + -1.333333333333324 * sp[17];
    A[24] = -0.6666666666666683 * sp[18] + -0.6666666666666682 * sp[17];
    A[25] = 0.0;
    A[26] = -0.6666666666666665 * sp[18] + -0.6666666666666649 * sp[17];
    A[27] = -1.333333333333328 * sp[19] + -0.6666666666666556 * sp[18] + -0.6666666666666605 * sp[18];
    A[28] = 1.333333333333328 * sp[19] + 0.6666666666666605 * sp[18] + 0.6666666666666605 * sp[18] + 1.333333333333332 * sp[17];
    A[29] = 0.6666666666666645 * sp[18] + 0.6666666666666702 * sp[18];
    A[30] = -0.6666666666666682 * sp[19] + -0.6666666666666675 * sp[18];
    A[31] = -0.6666666666666681 * sp[19] + -0.6666666666666666 * sp[18];
    A[32] = 0.0;
    A[33] = -0.6666666666666662 * sp[18] + -0.6666666666666645 * sp[18] + -1.333333333333324 * sp[17];
    A[34] = 0.6666666666666702 * sp[18] + 0.6666666666666645 * sp[18];
    A[35] = 1.333333333333336 * sp[19] + 0.6666666666666676 * sp[18] + 0.6666666666666676 * sp[18] + 1.333333333333334 * sp[17];
  }

};


class pressure_cell_integral_1_otherwise: public ufc::cell_integral
{
public:

  pressure_cell_integral_1_otherwise() : ufc::cell_integral()
  {

  }

  ~pressure_cell_integral_1_otherwise() override
  {

  }

  const std::vector<bool> & enabled_coefficients() const final override
  {
static const std::vector<bool> enabled({true, false});
return enabled;
  }

  void tabulate_tensor(double * A,
                       const double * const * w,
                       const double * coordinate_dofs,
                       int cell_orientation) const final override
  {
    // This function was generated using 'uflacs' representation
    // with the following integrals metadata:
    // 
    // num_cells:         None
    // optimize:          True
    // precision:         16
    // quadrature_degree: 4
    // quadrature_rule:   'default'
    // representation:    'uflacs'
    // 
    // and the following integral 0 metadata:
    // 
    // estimated_polynomial_degree: 4
    // optimize:                    True
    // precision:                   16
    // quadrature_degree:           4
    // quadrature_rule:             'default'
    // representation:              'uflacs'
    // Quadrature rules
    alignas(32) static const double weights6[6] = { 0.054975871827661, 0.054975871827661, 0.054975871827661, 0.1116907948390055, 0.1116907948390055, 0.1116907948390055 };
    // Precomputed values of basis functions and precomputations
    // FE* dimensions: [entities][points][dofs]
    // PI* dimensions: [entities][dofs][dofs] or [entities][dofs]
    // PM* dimensions: [entities][dofs][dofs]
    alignas(32) static const double FE3_C0_D01_Q6[1][1][2] = { { { -1.0, 1.0 } } };
    alignas(32) static const double FE6_C0_Q6[1][6][6] =
        { { { -0.07480380774819603, 0.5176323419876735, -0.07480380774819664, 0.2992152309927869, 0.03354481152314828, 0.2992152309927838 },
            { -0.07480380774819599, -0.07480380774819662, 0.5176323419876736, 0.2992152309927869, 0.2992152309927838, 0.03354481152314827 },
            { 0.5176323419876714, -0.07480380774819663, -0.07480380774819664, 0.03354481152314866, 0.2992152309927866, 0.2992152309927866 },
            { -0.04820837781551188, -0.08473049309397782, -0.04820837781551189, 0.1928335112620478, 0.795480226200906, 0.1928335112620478 },
            { -0.04820837781551192, -0.04820837781551197, -0.08473049309397784, 0.1928335112620478, 0.1928335112620478, 0.7954802262009062 },
            { -0.08473049309397777, -0.04820837781551201, -0.04820837781551189, 0.795480226200906, 0.1928335112620478, 0.1928335112620478 } } };
    // Unstructured piecewise computations
    const double J_c0 = coordinate_dofs[0] * FE3_C0_D01_Q6[0][0][0] + coordinate_dofs[2] * FE3_C0_D01_Q6[0][0][1];
    const double J_c3 = coordinate_dofs[1] * FE3_C0_D01_Q6[0][0][0] + coordinate_dofs[5] * FE3_C0_D01_Q6[0][0][1];
    const double J_c1 = coordinate_dofs[0] * FE3_C0_D01_Q6[0][0][0] + coordinate_dofs[4] * FE3_C0_D01_Q6[0][0][1];
    const double J_c2 = coordinate_dofs[1] * FE3_C0_D01_Q6[0][0][0] + coordinate_dofs[3] * FE3_C0_D01_Q6[0][0][1];
    alignas(32) double sp[4];
    sp[0] = J_c0 * J_c3;
    sp[1] = J_c1 * J_c2;
    sp[2] = sp[0] + -1 * sp[1];
    sp[3] = std::abs(sp[2]);
    alignas(32) double BF0[6] = {};
    for (int iq = 0; iq < 6; ++iq)
    {
        // Quadrature loop body setup (num_points=6)
        // Unstructured varying computations for num_points=6
        double w0 = 0.0;
        for (int ic = 0; ic < 6; ++ic)
            w0 += w[0][ic] * FE6_C0_Q6[0][iq][ic];
        alignas(32) double sv6[1];
        sv6[0] = sp[3] * w0;
        const double fw0 = sv6[0] * weights6[iq];
        for (int i = 0; i < 6; ++i)
            BF0[i] += fw0 * FE6_C0_Q6[0][iq][i];
    }
    std::fill(A, A + 6, 0.0);
    for (int i = 0; i < 6; ++i)
        A[i] += BF0[i];
  }

};


class pressure_exterior_facet_integral_1_otherwise: public ufc::exterior_facet_integral
{
public:

  pressure_exterior_facet_integral_1_otherwise() : ufc::exterior_facet_integral()
  {

  }

  ~pressure_exterior_facet_integral_1_otherwise() override
  {

  }

  const std::vector<bool> & enabled_coefficients() const final override
  {
static const std::vector<bool> enabled({false, true});
return enabled;
  }

  void tabulate_tensor(double * A,
                       const double * const * w,
                       const double * coordinate_dofs,
                       std::size_t facet,
                       int cell_orientation) const final override
  {
    // This function was generated using 'uflacs' representation
    // with the following integrals metadata:
    // 
    // num_cells:         None
    // optimize:          True
    // precision:         16
    // quadrature_degree: 4
    // quadrature_rule:   'default'
    // representation:    'uflacs'
    // 
    // and the following integral 0 metadata:
    // 
    // estimated_polynomial_degree: 4
    // optimize:                    True
    // precision:                   16
    // quadrature_degree:           4
    // quadrature_rule:             'default'
    // representation:              'uflacs'
    // Quadrature rules
    alignas(32) static const double weights3[3] = { 0.2777777777777777, 0.4444444444444444, 0.2777777777777777 };
    // Precomputed values of basis functions and precomputations
    // FE* dimensions: [entities][points][dofs]
    // PI* dimensions: [entities][dofs][dofs] or [entities][dofs]
    // PM* dimensions: [entities][dofs][dofs]
    alignas(32) static const double FE3_C0_D01_F_Q3[1][1][2] = { { { -1.0, 1.0 } } };
    alignas(32) static const double FE6_C0_F_Q3[3][3][6] =
        { { { 0.0, 0.6872983346207419, -0.08729833462074166, 0.3999999999999999, 0.0, 0.0 },
            { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
            { 0.0, -0.08729833462074163, 0.6872983346207419, 0.3999999999999998, 0.0, 0.0 } },
          { { 0.6872983346207417, 0.0, -0.08729833462074166, 0.0, 0.3999999999999997, 0.0 },
            { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
            { -0.08729833462074162, 0.0, 0.6872983346207419, 0.0, 0.3999999999999997, 0.0 } },
          { { 0.6872983346207417, -0.08729833462074163, 0.0, 0.0, 0.0, 0.3999999999999999 },
            { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 },
            { -0.08729833462074169, 0.6872983346207417, 0.0, 0.0, 0.0, 0.3999999999999999 } } };
    // Unstructured piecewise computations
    const double J_c0 = coordinate_dofs[0] * FE3_C0_D01_F_Q3[0][0][0] + coordinate_dofs[2] * FE3_C0_D01_F_Q3[0][0][1];
    const double J_c1 = coordinate_dofs[0] * FE3_C0_D01_F_Q3[0][0][0] + coordinate_dofs[4] * FE3_C0_D01_F_Q3[0][0][1];
    const double J_c2 = coordinate_dofs[1] * FE3_C0_D01_F_Q3[0][0][0] + coordinate_dofs[3] * FE3_C0_D01_F_Q3[0][0][1];
    const double J_c3 = coordinate_dofs[1] * FE3_C0_D01_F_Q3[0][0][0] + coordinate_dofs[5] * FE3_C0_D01_F_Q3[0][0][1];
    alignas(32) double sp[10];
    sp[0] = J_c0 * triangle_reference_facet_jacobian[facet][0][0];
    sp[1] = J_c1 * triangle_reference_facet_jacobian[facet][1][0];
    sp[2] = sp[0] + sp[1];
    sp[3] = sp[2] * sp[2];
    sp[4] = triangle_reference_facet_jacobian[facet][0][0] * J_c2;
    sp[5] = triangle_reference_facet_jacobian[facet][1][0] * J_c3;
    sp[6] = sp[4] + sp[5];
    sp[7] = sp[6] * sp[6];
    sp[8] = sp[3] + sp[7];
    sp[9] = std::sqrt(sp[8]);
    alignas(32) double BF0[6] = {};
    for (int iq = 0; iq < 3; ++iq)
    {
        // Quadrature loop body setup (num_points=3)
        // Unstructured varying computations for num_points=3
        double w1 = 0.0;
        for (int ic = 0; ic < 6; ++ic)
            w1 += w[1][ic] * FE6_C0_F_Q3[facet][iq][ic];
        alignas(32) double sv3[1];
        sv3[0] = sp[9] * w1;
        const double fw0 = sv3[0] * weights3[iq];
        for (int i = 0; i < 6; ++i)
            BF0[i] += fw0 * FE6_C0_F_Q3[facet][iq][i];
    }
    std::fill(A, A + 6, 0.0);
    for (int i = 0; i < 6; ++i)
        A[i] += BF0[i];
  }

};


class pressure_form_0: public ufc::form
{
public:

  pressure_form_0() : ufc::form()
  {
    // Do nothing
  }

  ~pressure_form_0() override
  {
    // Do nothing
  }

  const char * signature() const final override
  {
    return "ef65f0965f76d60bff4102f2edb06e54369b8ab46c44a14230eb040ec6ecf67f7cc8cc488e58b038ccc7ec2439e0eb39b6c3485567de3ff258d2517effdf741e";
  }

  std::size_t rank() const final override
  {
    return 2;
  }

  std::size_t num_coefficients() const final override
  {
    return 0;
  }

  std::size_t original_coefficient_position(std::size_t i) const final override
  {
    throw std::runtime_error("Invalid original coefficient index.");
    return i;
  }

  ufc::finite_element * create_coordinate_finite_element() const final override
  {
    return new pressure_finite_element_1();
  }

  ufc::dofmap * create_coordinate_dofmap() const final override
  {
    return new pressure_dofmap_1();
   }

  ufc::coordinate_mapping * create_coordinate_mapping() const final override
  {
    return new pressure_coordinate_mapping_1();
  }

  ufc::finite_element * create_finite_element(std::size_t i) const final override
  {
    switch (i)
    {
    case 0:
        return new pressure_finite_element_2();
    case 1:
        return new pressure_finite_element_2();
    default:
        return nullptr;
    }
  }

  ufc::dofmap * create_dofmap(std::size_t i) const final override
  {
    switch (i)
    {
    case 0:
        return new pressure_dofmap_2();
    case 1:
        return new pressure_dofmap_2();
    default:
        return nullptr;
    }
  }

  std::size_t max_cell_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_exterior_facet_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_interior_facet_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_vertex_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_custom_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_cutcell_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_interface_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_overlap_subdomain_id() const final override
  {
    return 0;
  }

  bool has_cell_integrals() const final override
  {
    return true;
  }

  bool has_exterior_facet_integrals() const final override
  {
    return false;
  }

  bool has_interior_facet_integrals() const final override
  {
    return false;
  }

  bool has_vertex_integrals() const final override
  {
    return false;
  }

  bool has_custom_integrals() const final override
  {
    return false;
  }

  bool has_cutcell_integrals() const final override
  {
    return false;
  }

  bool has_interface_integrals() const final override
  {
    return false;
  }

  bool has_overlap_integrals() const final override
  {
    return false;
  }

  ufc::cell_integral * create_cell_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::exterior_facet_integral * create_exterior_facet_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::interior_facet_integral * create_interior_facet_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::vertex_integral * create_vertex_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::custom_integral * create_custom_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::cutcell_integral * create_cutcell_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::interface_integral * create_interface_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::overlap_integral * create_overlap_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::cell_integral * create_default_cell_integral() const final override
  {
    return new pressure_cell_integral_0_otherwise();
  }

  ufc::exterior_facet_integral * create_default_exterior_facet_integral() const final override
  {
    return nullptr;
  }

  ufc::interior_facet_integral * create_default_interior_facet_integral() const final override
  {
    return nullptr;
  }

  ufc::vertex_integral * create_default_vertex_integral() const final override
  {
    return nullptr;
  }

  ufc::custom_integral * create_default_custom_integral() const final override
  {
    return nullptr;
  }

  ufc::cutcell_integral * create_default_cutcell_integral() const final override
  {
    return nullptr;
  }

  ufc::interface_integral * create_default_interface_integral() const final override
  {
    return nullptr;
  }

  ufc::overlap_integral * create_default_overlap_integral() const final override
  {
    return nullptr;
  }

};


class pressure_form_1: public ufc::form
{
public:

  pressure_form_1() : ufc::form()
  {
    // Do nothing
  }

  ~pressure_form_1() override
  {
    // Do nothing
  }

  const char * signature() const final override
  {
    return "d4ffd37393e7b7769ae23f96dc4c6d8c54f15db1a15ebe547a1fcd4d496abde3fd92ee14a8997a15ccae8e8e00606d6869fe13cf2c6911218824bd99b0f04f8c";
  }

  std::size_t rank() const final override
  {
    return 1;
  }

  std::size_t num_coefficients() const final override
  {
    return 2;
  }

  std::size_t original_coefficient_position(std::size_t i) const final override
  {
    if (i >= 2)
    {
        throw std::runtime_error("Invalid original coefficient index.");
    }
    static const std::vector<std::size_t> position = {0, 1};
    return position[i];
  }

  ufc::finite_element * create_coordinate_finite_element() const final override
  {
    return new pressure_finite_element_1();
  }

  ufc::dofmap * create_coordinate_dofmap() const final override
  {
    return new pressure_dofmap_1();
   }

  ufc::coordinate_mapping * create_coordinate_mapping() const final override
  {
    return new pressure_coordinate_mapping_1();
  }

  ufc::finite_element * create_finite_element(std::size_t i) const final override
  {
    switch (i)
    {
    case 0:
        return new pressure_finite_element_2();
    case 1:
        return new pressure_finite_element_2();
    case 2:
        return new pressure_finite_element_2();
    default:
        return nullptr;
    }
  }

  ufc::dofmap * create_dofmap(std::size_t i) const final override
  {
    switch (i)
    {
    case 0:
        return new pressure_dofmap_2();
    case 1:
        return new pressure_dofmap_2();
    case 2:
        return new pressure_dofmap_2();
    default:
        return nullptr;
    }
  }

  std::size_t max_cell_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_exterior_facet_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_interior_facet_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_vertex_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_custom_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_cutcell_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_interface_subdomain_id() const final override
  {
    return 0;
  }

  std::size_t max_overlap_subdomain_id() const final override
  {
    return 0;
  }

  bool has_cell_integrals() const final override
  {
    return true;
  }

  bool has_exterior_facet_integrals() const final override
  {
    return true;
  }

  bool has_interior_facet_integrals() const final override
  {
    return false;
  }

  bool has_vertex_integrals() const final override
  {
    return false;
  }

  bool has_custom_integrals() const final override
  {
    return false;
  }

  bool has_cutcell_integrals() const final override
  {
    return false;
  }

  bool has_interface_integrals() const final override
  {
    return false;
  }

  bool has_overlap_integrals() const final override
  {
    return false;
  }

  ufc::cell_integral * create_cell_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::exterior_facet_integral * create_exterior_facet_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::interior_facet_integral * create_interior_facet_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::vertex_integral * create_vertex_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::custom_integral * create_custom_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::cutcell_integral * create_cutcell_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::interface_integral * create_interface_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::overlap_integral * create_overlap_integral(std::size_t subdomain_id) const final override
  {
    return nullptr;
  }

  ufc::cell_integral * create_default_cell_integral() const final override
  {
    return new pressure_cell_integral_1_otherwise();
  }

  ufc::exterior_facet_integral * create_default_exterior_facet_integral() const final override
  {
    return new pressure_exterior_facet_integral_1_otherwise();
  }

  ufc::interior_facet_integral * create_default_interior_facet_integral() const final override
  {
    return nullptr;
  }

  ufc::vertex_integral * create_default_vertex_integral() const final override
  {
    return nullptr;
  }

  ufc::custom_integral * create_default_custom_integral() const final override
  {
    return nullptr;
  }

  ufc::cutcell_integral * create_default_cutcell_integral() const final override
  {
    return nullptr;
  }

  ufc::interface_integral * create_default_interface_integral() const final override
  {
    return nullptr;
  }

  ufc::overlap_integral * create_default_overlap_integral() const final override
  {
    return nullptr;
  }

};

// DOLFIN wrappers

// Standard library includes
#include <string>

// DOLFIN includes
#include <dolfin/common/NoDeleter.h>
#include <dolfin/mesh/Mesh.h>
#include <dolfin/mesh/MultiMesh.h>
#include <dolfin/fem/FiniteElement.h>
#include <dolfin/fem/DofMap.h>
#include <dolfin/fem/Form.h>
#include <dolfin/fem/MultiMeshForm.h>
#include <dolfin/function/FunctionSpace.h>
#include <dolfin/function/MultiMeshFunctionSpace.h>
#include <dolfin/function/GenericFunction.h>
#include <dolfin/function/CoefficientAssigner.h>
#include <dolfin/function/MultiMeshCoefficientAssigner.h>
#include <dolfin/adaptivity/ErrorControl.h>
#include <dolfin/adaptivity/GoalFunctional.h>
#include <dolfin/la/GenericVector.h>

namespace pressure
{

class CoefficientSpace_f: public dolfin::FunctionSpace
{
public:

  // Constructor for standard function space
  CoefficientSpace_f(std::shared_ptr<const dolfin::Mesh> mesh):
    dolfin::FunctionSpace(mesh,
                          std::make_shared<const dolfin::FiniteElement>(std::make_shared<pressure_finite_element_2>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<pressure_dofmap_2>(), *mesh))
  {
    // Do nothing
  }

  // Constructor for constrained function space
  CoefficientSpace_f(std::shared_ptr<const dolfin::Mesh> mesh, std::shared_ptr<const dolfin::SubDomain> constrained_domain):
    dolfin::FunctionSpace(mesh,
                          std::make_shared<const dolfin::FiniteElement>(std::make_shared<pressure_finite_element_2>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<pressure_dofmap_2>(), *mesh, constrained_domain))
  {
    // Do nothing
  }

};

class CoefficientSpace_g: public dolfin::FunctionSpace
{
public:

  // Constructor for standard function space
  CoefficientSpace_g(std::shared_ptr<const dolfin::Mesh> mesh):
    dolfin::FunctionSpace(mesh,
                          std::make_shared<const dolfin::FiniteElement>(std::make_shared<pressure_finite_element_2>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<pressure_dofmap_2>(), *mesh))
  {
    // Do nothing
  }

  // Constructor for constrained function space
  CoefficientSpace_g(std::shared_ptr<const dolfin::Mesh> mesh, std::shared_ptr<const dolfin::SubDomain> constrained_domain):
    dolfin::FunctionSpace(mesh,
                          std::make_shared<const dolfin::FiniteElement>(std::make_shared<pressure_finite_element_2>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<pressure_dofmap_2>(), *mesh, constrained_domain))
  {
    // Do nothing
  }

};

class Form_a_FunctionSpace_0: public dolfin::FunctionSpace
{
public:

  // Constructor for standard function space
  Form_a_FunctionSpace_0(std::shared_ptr<const dolfin::Mesh> mesh):
    dolfin::FunctionSpace(mesh,
                          std::make_shared<const dolfin::FiniteElement>(std::make_shared<pressure_finite_element_2>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<pressure_dofmap_2>(), *mesh))
  {
    // Do nothing
  }

  // Constructor for constrained function space
  Form_a_FunctionSpace_0(std::shared_ptr<const dolfin::Mesh> mesh, std::shared_ptr<const dolfin::SubDomain> constrained_domain):
    dolfin::FunctionSpace(mesh,
                          std::make_shared<const dolfin::FiniteElement>(std::make_shared<pressure_finite_element_2>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<pressure_dofmap_2>(), *mesh, constrained_domain))
  {
    // Do nothing
  }

};

class Form_a_FunctionSpace_1: public dolfin::FunctionSpace
{
public:

  // Constructor for standard function space
  Form_a_FunctionSpace_1(std::shared_ptr<const dolfin::Mesh> mesh):
    dolfin::FunctionSpace(mesh,
                          std::make_shared<const dolfin::FiniteElement>(std::make_shared<pressure_finite_element_2>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<pressure_dofmap_2>(), *mesh))
  {
    // Do nothing
  }

  // Constructor for constrained function space
  Form_a_FunctionSpace_1(std::shared_ptr<const dolfin::Mesh> mesh, std::shared_ptr<const dolfin::SubDomain> constrained_domain):
    dolfin::FunctionSpace(mesh,
                          std::make_shared<const dolfin::FiniteElement>(std::make_shared<pressure_finite_element_2>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<pressure_dofmap_2>(), *mesh, constrained_domain))
  {
    // Do nothing
  }

};

class Form_a_MultiMeshFunctionSpace_0: public dolfin::MultiMeshFunctionSpace
{
public:

  // Constructor for multimesh function space
  Form_a_MultiMeshFunctionSpace_0(std::shared_ptr<const dolfin::MultiMesh> multimesh): dolfin::MultiMeshFunctionSpace(multimesh)
  {
    // Create and add standard function spaces
    for (std::size_t part = 0; part < multimesh->num_parts(); part++)
    {
      std::shared_ptr<const dolfin::FunctionSpace> V(new Form_a_FunctionSpace_0(multimesh->part(part)));
      add(V);
    }

    // Build multimesh function space
    build();
  }

};

class Form_a_MultiMeshFunctionSpace_1: public dolfin::MultiMeshFunctionSpace
{
public:

  // Constructor for multimesh function space
  Form_a_MultiMeshFunctionSpace_1(std::shared_ptr<const dolfin::MultiMesh> multimesh): dolfin::MultiMeshFunctionSpace(multimesh)
  {
    // Create and add standard function spaces
    for (std::size_t part = 0; part < multimesh->num_parts(); part++)
    {
      std::shared_ptr<const dolfin::FunctionSpace> V(new Form_a_FunctionSpace_1(multimesh->part(part)));
      add(V);
    }

    // Build multimesh function space
    build();
  }

};

class Form_a: public dolfin::Form
{
public:

  // Constructor
  Form_a(std::shared_ptr<const dolfin::FunctionSpace> V1, std::shared_ptr<const dolfin::FunctionSpace> V0):
    dolfin::Form(2, 0)
  {
    _function_spaces[0] = V0;
    _function_spaces[1] = V1;

    _ufc_form = std::make_shared<const pressure_form_0>();
  }

  // Destructor
  ~Form_a()
  {}

  /// Return the number of the coefficient with this name
  virtual std::size_t coefficient_number(const std::string& name) const
  {

    dolfin::dolfin_error("generated code for class Form",
                         "access coefficient data",
                         "There are no coefficients");
    return 0;
  }

  /// Return the name of the coefficient with this number
  virtual std::string coefficient_name(std::size_t i) const
  {

    dolfin::dolfin_error("generated code for class Form",
                         "access coefficient data",
                         "There are no coefficients");
    return "unnamed";
  }

  // Typedefs
  typedef Form_a_FunctionSpace_0 TestSpace;
  typedef Form_a_FunctionSpace_1 TrialSpace;
  typedef Form_a_MultiMeshFunctionSpace_0 MultiMeshTestSpace;
  typedef Form_a_MultiMeshFunctionSpace_1 MultiMeshTrialSpace;

  // Coefficients
};

class MultiMeshForm_a: public dolfin::MultiMeshForm
{
public:

  // Constructor
  MultiMeshForm_a(std::shared_ptr<const dolfin::MultiMeshFunctionSpace> V1, std::shared_ptr<const dolfin::MultiMeshFunctionSpace> V0):
    dolfin::MultiMeshForm(V1, V0)
  {
    // Create and add standard forms
    std::size_t num_parts = V0->num_parts(); // assume all equal and pick first
    for (std::size_t part = 0; part < num_parts; part++)
    {
      std::shared_ptr<dolfin::Form> a(new Form_a(V1->part(part), V0->part(part)));
    add(a);

    }
    // Build multimesh form
    build();

    /// Assign coefficients

  }

  // Destructor
  ~MultiMeshForm_a()
  {}

  /// Return the number of the coefficient with this name
  virtual std::size_t coefficient_number(const std::string& name) const
  {

    dolfin::dolfin_error("generated code for class Form",
                         "access coefficient data",
                         "There are no coefficients");
    return 0;
  }

  /// Return the name of the coefficient with this number
  virtual std::string coefficient_name(std::size_t i) const
  {

    dolfin::dolfin_error("generated code for class Form",
                         "access coefficient data",
                         "There are no coefficients");
    return "unnamed";
  }

  // Typedefs
  typedef Form_a_FunctionSpace_0 TestSpace;
  typedef Form_a_FunctionSpace_1 TrialSpace;
  typedef Form_a_MultiMeshFunctionSpace_0 MultiMeshTestSpace;
  typedef Form_a_MultiMeshFunctionSpace_1 MultiMeshTrialSpace;

  // Coefficients
};

class Form_L_FunctionSpace_0: public dolfin::FunctionSpace
{
public:

  // Constructor for standard function space
  Form_L_FunctionSpace_0(std::shared_ptr<const dolfin::Mesh> mesh):
    dolfin::FunctionSpace(mesh,
                          std::make_shared<const dolfin::FiniteElement>(std::make_shared<pressure_finite_element_2>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<pressure_dofmap_2>(), *mesh))
  {
    // Do nothing
  }

  // Constructor for constrained function space
  Form_L_FunctionSpace_0(std::shared_ptr<const dolfin::Mesh> mesh, std::shared_ptr<const dolfin::SubDomain> constrained_domain):
    dolfin::FunctionSpace(mesh,
                          std::make_shared<const dolfin::FiniteElement>(std::make_shared<pressure_finite_element_2>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<pressure_dofmap_2>(), *mesh, constrained_domain))
  {
    // Do nothing
  }

};

class Form_L_MultiMeshFunctionSpace_0: public dolfin::MultiMeshFunctionSpace
{
public:

  // Constructor for multimesh function space
  Form_L_MultiMeshFunctionSpace_0(std::shared_ptr<const dolfin::MultiMesh> multimesh): dolfin::MultiMeshFunctionSpace(multimesh)
  {
    // Create and add standard function spaces
    for (std::size_t part = 0; part < multimesh->num_parts(); part++)
    {
      std::shared_ptr<const dolfin::FunctionSpace> V(new Form_L_FunctionSpace_0(multimesh->part(part)));
      add(V);
    }

    // Build multimesh function space
    build();
  }

};

typedef CoefficientSpace_f Form_L_FunctionSpace_1;

typedef CoefficientSpace_g Form_L_FunctionSpace_2;

class Form_L: public dolfin::Form
{
public:

  // Constructor
  Form_L(std::shared_ptr<const dolfin::FunctionSpace> V0):
    dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
  {
    _function_spaces[0] = V0;

    _ufc_form = std::make_shared<const pressure_form_1>();
  }

  // Constructor
  Form_L(std::shared_ptr<const dolfin::FunctionSpace> V0, std::shared_ptr<const dolfin::GenericFunction> f, std::shared_ptr<const dolfin::GenericFunction> g):
    dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
  {
    _function_spaces[0] = V0;

    this->f = f;
    this->g = g;

    _ufc_form = std::make_shared<const pressure_form_1>();
  }

  // Destructor
  ~Form_L()
  {}

  /// Return the number of the coefficient with this name
  virtual std::size_t coefficient_number(const std::string& name) const
  {
    if (name == "f")
      return 0;
    else if (name == "g")
      return 1;

    dolfin::dolfin_error("generated code for class Form",
                         "access coefficient data",
                         "Invalid coefficient");
    return 0;
  }

  /// Return the name of the coefficient with this number
  virtual std::string coefficient_name(std::size_t i) const
  {
    switch (i)
    {
    case 0:
      return "f";
    case 1:
      return "g";
    }

    dolfin::dolfin_error("generated code for class Form",
                         "access coefficient data",
                         "Invalid coefficient");
    return "unnamed";
  }

  // Typedefs
  typedef Form_L_FunctionSpace_0 TestSpace;
  typedef Form_L_MultiMeshFunctionSpace_0 MultiMeshTestSpace;
  typedef Form_L_FunctionSpace_1 CoefficientSpace_f;
  typedef Form_L_FunctionSpace_2 CoefficientSpace_g;

  // Coefficients
  dolfin::CoefficientAssigner f;
  dolfin::CoefficientAssigner g;
};

class MultiMeshForm_L: public dolfin::MultiMeshForm
{
public:

  // Constructor
  MultiMeshForm_L(std::shared_ptr<const dolfin::MultiMeshFunctionSpace> V0):
    dolfin::MultiMeshForm(V0), f(*this, 0), g(*this, 1)
  {
    // Create and add standard forms
    std::size_t num_parts = V0->num_parts(); // assume all equal and pick first
    for (std::size_t part = 0; part < num_parts; part++)
    {
      std::shared_ptr<dolfin::Form> a(new Form_L(V0->part(part)));
    add(a);

    }
    // Build multimesh form
    build();

    /// Assign coefficients

  }

  // Constructor
  MultiMeshForm_L(std::shared_ptr<const dolfin::MultiMeshFunctionSpace> V0, std::shared_ptr<const dolfin::GenericFunction> f, std::shared_ptr<const dolfin::GenericFunction> g):
    dolfin::MultiMeshForm(V0), f(*this, 0), g(*this, 1)
  {
    // Create and add standard forms
    std::size_t num_parts = V0->num_parts(); // assume all equal and pick first
    for (std::size_t part = 0; part < num_parts; part++)
    {
      std::shared_ptr<dolfin::Form> a(new Form_L(V0->part(part)));
    add(a);

    }
    // Build multimesh form
    build();

    /// Assign coefficients
    this->f = f;
    this->g = g;

  }

  // Destructor
  ~MultiMeshForm_L()
  {}

  /// Return the number of the coefficient with this name
  virtual std::size_t coefficient_number(const std::string& name) const
  {
    if (name == "f")
      return 0;
    else if (name == "g")
      return 1;

    dolfin::dolfin_error("generated code for class Form",
                         "access coefficient data",
                         "Invalid coefficient");
    return 0;
  }

  /// Return the name of the coefficient with this number
  virtual std::string coefficient_name(std::size_t i) const
  {
    switch (i)
    {
    case 0:
      return "f";
    case 1:
      return "g";
    }

    dolfin::dolfin_error("generated code for class Form",
                         "access coefficient data",
                         "Invalid coefficient");
    return "unnamed";
  }

  // Typedefs
  typedef Form_L_FunctionSpace_0 TestSpace;
  typedef Form_L_MultiMeshFunctionSpace_0 MultiMeshTestSpace;
  typedef Form_L_FunctionSpace_1 CoefficientSpace_f;
  typedef Form_L_FunctionSpace_2 CoefficientSpace_g;

  // Coefficients
  dolfin::MultiMeshCoefficientAssigner f;
  dolfin::MultiMeshCoefficientAssigner g;
};

// Class typedefs
typedef Form_a BilinearForm;
typedef MultiMeshForm_a MultiMeshBilinearForm;
typedef Form_a JacobianForm;
typedef MultiMeshForm_a MultiMeshJacobianForm;
typedef Form_L LinearForm;
typedef MultiMeshForm_L MultiMeshLinearForm;
typedef Form_L ResidualForm;
typedef MultiMeshForm_L MultiMeshResidualForm;
typedef Form_a::TestSpace FunctionSpace;
typedef Form_a::MultiMeshTestSpace MultiMeshFunctionSpace;

}

#endif
