MyLagrange.h

// This code conforms with the UFC specification version 2018.1.0
// and was automatically generated by FFC version 2019.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 __MYLAGRANGE_H
#define __MYLAGRANGE_H
#include <algorithm>
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <ufc.h>

class mylagrange_finite_element_0: public ufc::finite_element
{
public:

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

  ~mylagrange_finite_element_0() override
  {
    // Do nothing
  }

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

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

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

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

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

  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][4] = { { 0.288675134594813, -0.18257418583505536, -0.10540925533894598, -0.07453559924999298 } };
    static const double coefficients1[1][4] = { { 0.2886751345948129, 0.18257418583505536, -0.10540925533894596, -0.07453559924999296 } };
    static const double coefficients2[1][4] = { { 0.2886751345948129, 0.0, 0.21081851067789195, -0.07453559924999299 } };
    static const double coefficients3[1][4] = { { 0.2886751345948129, 0.0, 0.0, 0.22360679774997896 } };
    for (std::size_t k = 0; k < num_points * 4; ++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[3] = { 2.0 * X[ip * 3] - 1.0, 2.0 * X[ip * 3 + 1] - 1.0, 2.0 * X[ip * 3 + 2] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues1[4] = {};
        basisvalues1[0] = 1.0;
        const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
        basisvalues1[1] = tmp1_1;
        basisvalues1[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues1[0];
        basisvalues1[3] = (2.0 * Y[2] + 1.0) * basisvalues1[0];
        basisvalues1[0] *= std::sqrt(0.75);
        basisvalues1[3] *= std::sqrt(1.25);
        basisvalues1[2] *= std::sqrt(2.5);
        basisvalues1[1] *= std::sqrt(7.5);
        // Accumulate products of coefficients and basisvalues
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[4 * ip] += coefficients0[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[4 * ip + 1] += coefficients1[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[4 * ip + 2] += coefficients2[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[4 * ip + 3] += coefficients3[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(3, order);
    std::fill_n(reference_values, num_points * 4 * num_derivatives, 0.0);
    if (order > 1)
        return;
    // Tables of derivatives of the polynomial base (transpose).
    alignas(32) static const double dmats0[3][4][4] =
        { { { 0.0, 0.0, 0.0, 0.0 },
            { 6.324555320336757, 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 },
            { 3.1622776601683786, 0.0, 0.0, 0.0 },
            { 5.477225575051662, 0.0, 0.0, 0.0 },
            { 0.0, 0.0, 0.0, 0.0 } },
          { { 0.0, 0.0, 0.0, 0.0 },
            { 3.1622776601683786, 0.0, 0.0, 0.0 },
            { 1.825741858350554, 0.0, 0.0, 0.0 },
            { 5.163977794943222, 0.0, 0.0, 0.0 } } };
    static const double coefficients0[1][4] = { { 0.288675134594813, -0.18257418583505536, -0.10540925533894598, -0.07453559924999298 } };
    static const double coefficients1[1][4] = { { 0.2886751345948129, 0.18257418583505536, -0.10540925533894596, -0.07453559924999296 } };
    static const double coefficients2[1][4] = { { 0.2886751345948129, 0.0, 0.21081851067789195, -0.07453559924999299 } };
    static const double coefficients3[1][4] = { { 0.2886751345948129, 0.0, 0.0, 0.22360679774997896 } };
    const std::size_t reference_offset[4] = {};
    const std::size_t num_components[4] = { 1, 1, 1, 1 };
    // Precomputed combinations
    const std::size_t combinations[1][3][1] =
        { { { 0 },
            { 1 },
            { 2 } } };
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Map from UFC reference coordinate X to FIAT reference coordinate Y
        const double Y[3] = { 2.0 * X[ip * 3] - 1.0, 2.0 * X[ip * 3 + 1] - 1.0, 2.0 * X[ip * 3 + 2] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues1[4] = {};
        basisvalues1[0] = 1.0;
        const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
        basisvalues1[1] = tmp1_1;
        basisvalues1[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues1[0];
        basisvalues1[3] = (2.0 * Y[2] + 1.0) * basisvalues1[0];
        basisvalues1[0] *= std::sqrt(0.75);
        basisvalues1[3] *= std::sqrt(1.25);
        basisvalues1[2] *= std::sqrt(2.5);
        basisvalues1[1] *= std::sqrt(7.5);
        // Loop over all dofs
        for (std::size_t i = 0; i < 4; ++i)
        {
            double derivatives[3] = {};
            switch (i)
            {
            case 0:
                // Compute reference derivatives for dof 0.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++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[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++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[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++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[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++s)
                        derivatives[r] += coefficients3[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[4 * 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(3, order);
    // Precomputed combinations
    const std::size_t combinations[1][3][1] =
        { { { 0 },
            { 1 },
            { 2 } } };
    std::fill_n(values, num_points * 4 * num_derivatives, 0.0);
    const std::size_t reference_offsets[4] = {};
    const std::size_t physical_offsets[4] = {};
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        double transform[3][3];
        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[3 * 3 * ip + 3 * combinations[order - 1][s][k] + combinations[order - 1][r][k]];
        for (std::size_t d = 0; d < 4; ++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[4 * 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[4 * 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[3] = {};
    double J[9];
    double detJ;
    double K[9];
    if (cm)
    {
        cm->compute_reference_geometry(X, J, &detJ, K, 1, x, coordinate_dofs, cell_orientation);
    }
    else
    {
        compute_jacobian_tetrahedron_3d(J, coordinate_dofs);
        compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
        // Compute constants
        const double C0 = coordinate_dofs[9] + coordinate_dofs[6] + coordinate_dofs[3] - coordinate_dofs[0];
        const double C1 = coordinate_dofs[10] + coordinate_dofs[7] + coordinate_dofs[4] - coordinate_dofs[1];
        const double C2 = coordinate_dofs[11] + coordinate_dofs[8] + coordinate_dofs[5] - coordinate_dofs[2];
        // Compute subdeterminants
        const double d[9] = { J[4] * J[8] - J[5] * J[7], J[5] * J[6] - J[3] * J[8], J[3] * J[7] - J[4] * J[6], J[2] * J[7] - J[1] * J[8], J[0] * J[8] - J[2] * J[6], J[1] * J[6] - J[0] * J[7], J[1] * J[5] - J[2] * J[4], J[2] * J[3] - J[0] * J[5], J[0] * J[4] - J[1] * J[3] };
        // Get coordinates and map to the reference (FIAT) element
        double Y[3] = { (d[0] * (2.0 * x[0] - C0) + d[3] * (2.0 * x[1] - C1) + d[6] * (2.0 * x[2] - C2)) / detJ, (d[1] * (2.0 * x[0] - C0) + d[4] * (2.0 * x[1] - C1) + d[7] * (2.0 * x[2] - C2)) / detJ, (d[2] * (2.0 * x[0] - C0) + d[5] * (2.0 * x[1] - C1) + d[8] * (2.0 * x[2] - C2)) / detJ };
        // Map to FFC reference coordinate
        for (std::size_t k = 0; k < 3; ++k)
            X[k] = (Y[k] + 1.0) / 2.0;
    }
    // Evaluate basis on reference element
    double ref_values[4];
    evaluate_reference_basis(ref_values, 1, X);
    // Push forward
    double physical_values[4];
    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 < 4; ++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(3, 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[9];
    compute_jacobian_tetrahedron_3d(J, coordinate_dofs);
    // Compute Inverse Jacobian and determinant
    double K[9];
    double detJ;
    compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
    // Compute constants
    const double C0 = coordinate_dofs[9] + coordinate_dofs[6] + coordinate_dofs[3] - coordinate_dofs[0];
    const double C1 = coordinate_dofs[10] + coordinate_dofs[7] + coordinate_dofs[4] - coordinate_dofs[1];
    const double C2 = coordinate_dofs[11] + coordinate_dofs[8] + coordinate_dofs[5] - coordinate_dofs[2];
    // Compute subdeterminants
    const double d[9] = { J[4] * J[8] - J[5] * J[7], J[5] * J[6] - J[3] * J[8], J[3] * J[7] - J[4] * J[6], J[2] * J[7] - J[1] * J[8], J[0] * J[8] - J[2] * J[6], J[1] * J[6] - J[0] * J[7], J[1] * J[5] - J[2] * J[4], J[2] * J[3] - J[0] * J[5], J[0] * J[4] - J[1] * J[3] };
    // Get coordinates and map to the reference (FIAT) element
    double Y[3] = { (d[0] * (2.0 * x[0] - C0) + d[3] * (2.0 * x[1] - C1) + d[6] * (2.0 * x[2] - C2)) / detJ, (d[1] * (2.0 * x[0] - C0) + d[4] * (2.0 * x[1] - C1) + d[7] * (2.0 * x[2] - C2)) / detJ, (d[2] * (2.0 * x[0] - C0) + d[5] * (2.0 * x[1] - C1) + d[8] * (2.0 * x[2] - C2)) / detJ };
    // Precomputed combinations
    const std::size_t combinations[1][3][1] =
        { { { 0 },
            { 1 },
            { 2 } } };
    // Declare transformation matrix
    double transform[3][3] =
        { { 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[3 * combinations[n - 1][col][k] + combinations[n - 1][row][k]];
    switch (i)
    {
    case 0:
        {
            double basisvalues[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.288675134594813, -0.18257418583505536, -0.10540925533894598, -0.07453559924999298 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.2886751345948129, 0.18257418583505536, -0.10540925533894596, -0.07453559924999296 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.2886751345948129, 0.0, 0.21081851067789195, -0.07453559924999299 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.2886751345948129, 0.0, 0.0, 0.22360679774997896 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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(3, n);
    // Set values equal to zero.
    std::fill_n(values, num_derivatives * 4, 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[3] = {};
    // Loop dofs and call evaluate_basis_derivatives.
    for (std::size_t r = 0; r < 4; ++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[3];
    switch (i)
    {
    case 0:
        {
            y[0] = coordinate_dofs[0];
            y[1] = coordinate_dofs[1];
            y[2] = coordinate_dofs[2];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 1:
        {
            y[0] = coordinate_dofs[3];
            y[1] = coordinate_dofs[4];
            y[2] = coordinate_dofs[5];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 2:
        {
            y[0] = coordinate_dofs[6];
            y[1] = coordinate_dofs[7];
            y[2] = coordinate_dofs[8];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 3:
        {
            y[0] = coordinate_dofs[9];
            y[1] = coordinate_dofs[10];
            y[2] = coordinate_dofs[11];
            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[3];
    y[0] = coordinate_dofs[0];
    y[1] = coordinate_dofs[1];
    y[2] = coordinate_dofs[2];
    f.evaluate(vals, y, c);
    values[0] = vals[0];
    y[0] = coordinate_dofs[3];
    y[1] = coordinate_dofs[4];
    y[2] = coordinate_dofs[5];
    f.evaluate(vals, y, c);
    values[1] = vals[0];
    y[0] = coordinate_dofs[6];
    y[1] = coordinate_dofs[7];
    y[2] = coordinate_dofs[8];
    f.evaluate(vals, y, c);
    values[2] = vals[0];
    y[0] = coordinate_dofs[9];
    y[1] = coordinate_dofs[10];
    y[2] = coordinate_dofs[11];
    f.evaluate(vals, y, c);
    values[3] = 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];
    vertex_values[3] = dof_values[3];
  }

  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[3] = coordinate_dofs[3];
    dof_coordinates[3 + 1] = coordinate_dofs[4];
    dof_coordinates[3 + 2] = coordinate_dofs[5];
    dof_coordinates[3 * 2] = coordinate_dofs[6];
    dof_coordinates[3 * 2 + 1] = coordinate_dofs[7];
    dof_coordinates[3 * 2 + 2] = coordinate_dofs[8];
    dof_coordinates[3 * 3] = coordinate_dofs[9];
    dof_coordinates[3 * 3 + 1] = coordinate_dofs[10];
    dof_coordinates[3 * 3 + 2] = coordinate_dofs[11];
  }

  void tabulate_reference_dof_coordinates(double * reference_dof_coordinates) const final override
  {
    static const double dof_X[12] = { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.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 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 mylagrange_finite_element_0();
  }

};


class mylagrange_finite_element_1: public ufc::finite_element
{
public:

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

  ~mylagrange_finite_element_1() override
  {
    // Do nothing
  }

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

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

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

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

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

  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] = { 3 };
    if (i >= 1)
        return 1;
    return return_values[i];
  }

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

  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] = { 3 };
    if (i >= 1)
        return 1;
    return return_values[i];
  }

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

  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][4] = { { 0.288675134594813, -0.18257418583505536, -0.10540925533894598, -0.07453559924999298 } };
    static const double coefficients1[1][4] = { { 0.2886751345948129, 0.18257418583505536, -0.10540925533894596, -0.07453559924999296 } };
    static const double coefficients2[1][4] = { { 0.2886751345948129, 0.0, 0.21081851067789195, -0.07453559924999299 } };
    static const double coefficients3[1][4] = { { 0.2886751345948129, 0.0, 0.0, 0.22360679774997896 } };
    for (std::size_t k = 0; k < num_points * 12 * 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[3] = { 2.0 * X[ip * 3] - 1.0, 2.0 * X[ip * 3 + 1] - 1.0, 2.0 * X[ip * 3 + 2] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues1[4] = {};
        basisvalues1[0] = 1.0;
        const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
        basisvalues1[1] = tmp1_1;
        basisvalues1[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues1[0];
        basisvalues1[3] = (2.0 * Y[2] + 1.0) * basisvalues1[0];
        basisvalues1[0] *= std::sqrt(0.75);
        basisvalues1[3] *= std::sqrt(1.25);
        basisvalues1[2] *= std::sqrt(2.5);
        basisvalues1[1] *= std::sqrt(7.5);
        // Accumulate products of coefficients and basisvalues
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[12 * 3 * ip] += coefficients0[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[12 * 3 * ip + 3] += coefficients1[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[12 * 3 * ip + 3 * 2] += coefficients2[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[12 * 3 * ip + 3 * 3] += coefficients3[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[12 * 3 * ip + 3 * 4 + 1] += coefficients0[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[12 * 3 * ip + 3 * 5 + 1] += coefficients1[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[12 * 3 * ip + 3 * 6 + 1] += coefficients2[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[12 * 3 * ip + 3 * 7 + 1] += coefficients3[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[12 * 3 * ip + 3 * 8 + 2] += coefficients0[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[12 * 3 * ip + 3 * 9 + 2] += coefficients1[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[12 * 3 * ip + 3 * 10 + 2] += coefficients2[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 4; ++r)
            reference_values[12 * 3 * ip + 3 * 11 + 2] += coefficients3[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(3, order);
    std::fill_n(reference_values, num_points * 12 * num_derivatives * 3, 0.0);
    if (order > 1)
        return;
    // Tables of derivatives of the polynomial base (transpose).
    alignas(32) static const double dmats0[3][4][4] =
        { { { 0.0, 0.0, 0.0, 0.0 },
            { 6.324555320336757, 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 },
            { 3.1622776601683786, 0.0, 0.0, 0.0 },
            { 5.477225575051662, 0.0, 0.0, 0.0 },
            { 0.0, 0.0, 0.0, 0.0 } },
          { { 0.0, 0.0, 0.0, 0.0 },
            { 3.1622776601683786, 0.0, 0.0, 0.0 },
            { 1.825741858350554, 0.0, 0.0, 0.0 },
            { 5.163977794943222, 0.0, 0.0, 0.0 } } };
    static const double coefficients0[1][4] = { { 0.288675134594813, -0.18257418583505536, -0.10540925533894598, -0.07453559924999298 } };
    static const double coefficients1[1][4] = { { 0.2886751345948129, 0.18257418583505536, -0.10540925533894596, -0.07453559924999296 } };
    static const double coefficients2[1][4] = { { 0.2886751345948129, 0.0, 0.21081851067789195, -0.07453559924999299 } };
    static const double coefficients3[1][4] = { { 0.2886751345948129, 0.0, 0.0, 0.22360679774997896 } };
    const std::size_t reference_offset[12] = { 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2 };
    const std::size_t num_components[12] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
    // Precomputed combinations
    const std::size_t combinations[1][3][1] =
        { { { 0 },
            { 1 },
            { 2 } } };
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Map from UFC reference coordinate X to FIAT reference coordinate Y
        const double Y[3] = { 2.0 * X[ip * 3] - 1.0, 2.0 * X[ip * 3 + 1] - 1.0, 2.0 * X[ip * 3 + 2] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues1[4] = {};
        basisvalues1[0] = 1.0;
        const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
        basisvalues1[1] = tmp1_1;
        basisvalues1[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues1[0];
        basisvalues1[3] = (2.0 * Y[2] + 1.0) * basisvalues1[0];
        basisvalues1[0] *= std::sqrt(0.75);
        basisvalues1[3] *= std::sqrt(1.25);
        basisvalues1[2] *= std::sqrt(2.5);
        basisvalues1[1] *= std::sqrt(7.5);
        // Loop over all dofs
        for (std::size_t i = 0; i < 12; ++i)
        {
            double derivatives[3] = {};
            switch (i)
            {
            case 0:
                // Compute reference derivatives for dof 0.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++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[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++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[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++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[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++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[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++s)
                        derivatives[r] += coefficients0[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[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++s)
                        derivatives[r] += coefficients1[0][s] * aux[s];
                }
                break;
            case 6:
                // Compute reference derivatives for dof 6.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++s)
                        derivatives[r] += coefficients2[0][s] * aux[s];
                }
                break;
            case 7:
                // Compute reference derivatives for dof 7.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++s)
                        derivatives[r] += coefficients3[0][s] * aux[s];
                }
                break;
            case 8:
                // Compute reference derivatives for dof 8.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++s)
                        derivatives[r] += coefficients0[0][s] * aux[s];
                }
                break;
            case 9:
                // Compute reference derivatives for dof 9.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++s)
                        derivatives[r] += coefficients1[0][s] * aux[s];
                }
                break;
            case 10:
                // Compute reference derivatives for dof 10.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++s)
                        derivatives[r] += coefficients2[0][s] * aux[s];
                }
                break;
            case 11:
                // Compute reference derivatives for dof 11.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[4] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[4][4] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 16, &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[4][4];
                        std::copy_n(&dmats[0][0], 16, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 16, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 4; ++s)
                        for (std::size_t t = 0; t < 4; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 4; ++s)
                        derivatives[r] += coefficients3[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[12 * (num_derivatives * 3) * ip + num_derivatives * 3 * i + 3 * 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(3, order);
    // Precomputed combinations
    const std::size_t combinations[1][3][1] =
        { { { 0 },
            { 1 },
            { 2 } } };
    std::fill_n(values, num_points * 12 * num_derivatives * 3, 0.0);
    const std::size_t reference_offsets[12] = { 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2 };
    const std::size_t physical_offsets[12] = { 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2 };
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        double transform[3][3];
        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[3 * 3 * ip + 3 * combinations[order - 1][s][k] + combinations[order - 1][r][k]];
        for (std::size_t d = 0; d < 12; ++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[12 * (num_derivatives * 3) * ip + num_derivatives * 3 * d + 3 * s + reference_offsets[d]];
                    // Mapping derivatives back to the physical element
                    for (std::size_t r = 0; r < num_derivatives; ++r)
                        values[12 * (num_derivatives * 3) * ip + num_derivatives * 3 * d + 3 * 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[3] = {};
    double J[9];
    double detJ;
    double K[9];
    if (cm)
    {
        cm->compute_reference_geometry(X, J, &detJ, K, 1, x, coordinate_dofs, cell_orientation);
    }
    else
    {
        compute_jacobian_tetrahedron_3d(J, coordinate_dofs);
        compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
        // Compute constants
        const double C0 = coordinate_dofs[9] + coordinate_dofs[6] + coordinate_dofs[3] - coordinate_dofs[0];
        const double C1 = coordinate_dofs[10] + coordinate_dofs[7] + coordinate_dofs[4] - coordinate_dofs[1];
        const double C2 = coordinate_dofs[11] + coordinate_dofs[8] + coordinate_dofs[5] - coordinate_dofs[2];
        // Compute subdeterminants
        const double d[9] = { J[4] * J[8] - J[5] * J[7], J[5] * J[6] - J[3] * J[8], J[3] * J[7] - J[4] * J[6], J[2] * J[7] - J[1] * J[8], J[0] * J[8] - J[2] * J[6], J[1] * J[6] - J[0] * J[7], J[1] * J[5] - J[2] * J[4], J[2] * J[3] - J[0] * J[5], J[0] * J[4] - J[1] * J[3] };
        // Get coordinates and map to the reference (FIAT) element
        double Y[3] = { (d[0] * (2.0 * x[0] - C0) + d[3] * (2.0 * x[1] - C1) + d[6] * (2.0 * x[2] - C2)) / detJ, (d[1] * (2.0 * x[0] - C0) + d[4] * (2.0 * x[1] - C1) + d[7] * (2.0 * x[2] - C2)) / detJ, (d[2] * (2.0 * x[0] - C0) + d[5] * (2.0 * x[1] - C1) + d[8] * (2.0 * x[2] - C2)) / detJ };
        // Map to FFC reference coordinate
        for (std::size_t k = 0; k < 3; ++k)
            X[k] = (Y[k] + 1.0) / 2.0;
    }
    // Evaluate basis on reference element
    double ref_values[36];
    evaluate_reference_basis(ref_values, 1, X);
    // Push forward
    double physical_values[36];
    transform_reference_basis_derivatives(physical_values, 0, 1, ref_values, X, J, &detJ, K, cell_orientation);
    for (std::size_t k = 0; k < 3; ++k)
        values[k] = physical_values[3 * 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[3] = {};
    // Loop dofs and call evaluate_basis
    for (std::size_t r = 0; r < 12; ++r)
    {
        evaluate_basis(r, dof_values, x, coordinate_dofs, cell_orientation);
        for (std::size_t s = 0; s < 3; ++s)
            values[r * 3 + 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(3, n);
    std::fill_n(values, 3 * 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[9];
    compute_jacobian_tetrahedron_3d(J, coordinate_dofs);
    // Compute Inverse Jacobian and determinant
    double K[9];
    double detJ;
    compute_jacobian_inverse_tetrahedron_3d(K, detJ, J);
    // Compute constants
    const double C0 = coordinate_dofs[9] + coordinate_dofs[6] + coordinate_dofs[3] - coordinate_dofs[0];
    const double C1 = coordinate_dofs[10] + coordinate_dofs[7] + coordinate_dofs[4] - coordinate_dofs[1];
    const double C2 = coordinate_dofs[11] + coordinate_dofs[8] + coordinate_dofs[5] - coordinate_dofs[2];
    // Compute subdeterminants
    const double d[9] = { J[4] * J[8] - J[5] * J[7], J[5] * J[6] - J[3] * J[8], J[3] * J[7] - J[4] * J[6], J[2] * J[7] - J[1] * J[8], J[0] * J[8] - J[2] * J[6], J[1] * J[6] - J[0] * J[7], J[1] * J[5] - J[2] * J[4], J[2] * J[3] - J[0] * J[5], J[0] * J[4] - J[1] * J[3] };
    // Get coordinates and map to the reference (FIAT) element
    double Y[3] = { (d[0] * (2.0 * x[0] - C0) + d[3] * (2.0 * x[1] - C1) + d[6] * (2.0 * x[2] - C2)) / detJ, (d[1] * (2.0 * x[0] - C0) + d[4] * (2.0 * x[1] - C1) + d[7] * (2.0 * x[2] - C2)) / detJ, (d[2] * (2.0 * x[0] - C0) + d[5] * (2.0 * x[1] - C1) + d[8] * (2.0 * x[2] - C2)) / detJ };
    // Precomputed combinations
    const std::size_t combinations[1][3][1] =
        { { { 0 },
            { 1 },
            { 2 } } };
    // Declare transformation matrix
    double transform[3][3] =
        { { 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[3 * combinations[n - 1][col][k] + combinations[n - 1][row][k]];
    switch (i)
    {
    case 0:
        {
            double basisvalues[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.288675134594813, -0.18257418583505536, -0.10540925533894598, -0.07453559924999298 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.2886751345948129, 0.18257418583505536, -0.10540925533894596, -0.07453559924999296 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.2886751345948129, 0.0, 0.21081851067789195, -0.07453559924999299 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.2886751345948129, 0.0, 0.0, 0.22360679774997896 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.288675134594813, -0.18257418583505536, -0.10540925533894598, -0.07453559924999298 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.2886751345948129, 0.18257418583505536, -0.10540925533894596, -0.07453559924999296 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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 6:
        {
            double basisvalues[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.2886751345948129, 0.0, 0.21081851067789195, -0.07453559924999299 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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 7:
        {
            double basisvalues[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.2886751345948129, 0.0, 0.0, 0.22360679774997896 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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 8:
        {
            double basisvalues[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.288675134594813, -0.18257418583505536, -0.10540925533894598, -0.07453559924999298 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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[2 * num_derivatives + r] += transform[r][s] * derivatives[s];
        }
        break;
    case 9:
        {
            double basisvalues[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.2886751345948129, 0.18257418583505536, -0.10540925533894596, -0.07453559924999296 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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[2 * num_derivatives + r] += transform[r][s] * derivatives[s];
        }
        break;
    case 10:
        {
            double basisvalues[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.2886751345948129, 0.0, 0.21081851067789195, -0.07453559924999299 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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[2 * num_derivatives + r] += transform[r][s] * derivatives[s];
        }
        break;
    case 11:
        {
            double basisvalues[4] = {};
            basisvalues[0] = 1.0;
            const double tmp1_1 = 0.5 * (2.0 + 2.0 * Y[0] + Y[1] + Y[2]);
            basisvalues[1] = tmp1_1;
            basisvalues[2] = 0.5 * (2.0 + 3.0 * Y[1] + Y[2]) * basisvalues[0];
            basisvalues[3] = (2.0 * Y[2] + 1.0) * basisvalues[0];
            basisvalues[0] *= std::sqrt(0.75);
            basisvalues[3] *= std::sqrt(1.25);
            basisvalues[2] *= std::sqrt(2.5);
            basisvalues[1] *= std::sqrt(7.5);
            // Table(s) of coefficients
            static const double coefficients0[4] = { 0.2886751345948129, 0.0, 0.0, 0.22360679774997896 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 6.324555320336757, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats1[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 5.477225575051662, 0.0, 0.0, 0.0 },
                  { 0.0, 0.0, 0.0, 0.0 } };
            static const double dmats2[4][4] =
                { { 0.0, 0.0, 0.0, 0.0 },
                  { 3.1622776601683786, 0.0, 0.0, 0.0 },
                  { 1.825741858350554, 0.0, 0.0, 0.0 },
                  { 5.163977794943222, 0.0, 0.0, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[3] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.0, 0.0 },
                  { 0.0, 0.0, 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[4][4] =
                { { 1.0, 0.0, 0.0, 0.0 },
                  { 0.0, 1.0, 0.0, 0.0 },
                  { 0.0, 0.0, 1.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], 16, 0.0);
                for (std::size_t t = 0; t < 4; ++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], 16, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 16, 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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++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 < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats1[t][tu];
                    }
                    // _dmats_product(shape_dmats, comb[r][s], 2)
                    if (combinations[n - 1][r][s] == 2)
                    {
                        for (std::size_t t = 0; t < 4; ++t)
                            for (std::size_t u = 0; u < 4; ++u)
                                for (std::size_t tu = 0; tu < 4; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats2[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 4; ++s)
                    for (std::size_t t = 0; t < 4; ++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[2 * 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(3, n);
    // Set values equal to zero.
    std::fill_n(values, 3 * num_derivatives * 12, 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[9] = {};
    // Loop dofs and call evaluate_basis_derivatives.
    for (std::size_t r = 0; r < 12; ++r)
    {
        evaluate_basis_derivatives(r, n, dof_values, x, coordinate_dofs, cell_orientation);
        for (std::size_t s = 0; s < 3 * num_derivatives; ++s)
            values[3 * 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[3];
    // Declare variable for physical coordinates
    double y[3];
    switch (i)
    {
    case 0:
        {
            y[0] = coordinate_dofs[0];
            y[1] = coordinate_dofs[1];
            y[2] = coordinate_dofs[2];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 1:
        {
            y[0] = coordinate_dofs[3];
            y[1] = coordinate_dofs[4];
            y[2] = coordinate_dofs[5];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 2:
        {
            y[0] = coordinate_dofs[6];
            y[1] = coordinate_dofs[7];
            y[2] = coordinate_dofs[8];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 3:
        {
            y[0] = coordinate_dofs[9];
            y[1] = coordinate_dofs[10];
            y[2] = coordinate_dofs[11];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 4:
        {
            y[0] = coordinate_dofs[0];
            y[1] = coordinate_dofs[1];
            y[2] = coordinate_dofs[2];
            f.evaluate(vals, y, c);
            return vals[1];
        }
        break;
    case 5:
        {
            y[0] = coordinate_dofs[3];
            y[1] = coordinate_dofs[4];
            y[2] = coordinate_dofs[5];
            f.evaluate(vals, y, c);
            return vals[1];
        }
        break;
    case 6:
        {
            y[0] = coordinate_dofs[6];
            y[1] = coordinate_dofs[7];
            y[2] = coordinate_dofs[8];
            f.evaluate(vals, y, c);
            return vals[1];
        }
        break;
    case 7:
        {
            y[0] = coordinate_dofs[9];
            y[1] = coordinate_dofs[10];
            y[2] = coordinate_dofs[11];
            f.evaluate(vals, y, c);
            return vals[1];
        }
        break;
    case 8:
        {
            y[0] = coordinate_dofs[0];
            y[1] = coordinate_dofs[1];
            y[2] = coordinate_dofs[2];
            f.evaluate(vals, y, c);
            return vals[2];
        }
        break;
    case 9:
        {
            y[0] = coordinate_dofs[3];
            y[1] = coordinate_dofs[4];
            y[2] = coordinate_dofs[5];
            f.evaluate(vals, y, c);
            return vals[2];
        }
        break;
    case 10:
        {
            y[0] = coordinate_dofs[6];
            y[1] = coordinate_dofs[7];
            y[2] = coordinate_dofs[8];
            f.evaluate(vals, y, c);
            return vals[2];
        }
        break;
    case 11:
        {
            y[0] = coordinate_dofs[9];
            y[1] = coordinate_dofs[10];
            y[2] = coordinate_dofs[11];
            f.evaluate(vals, y, c);
            return vals[2];
        }
        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[3];
    // Declare variable for physical coordinates
    double y[3];
    y[0] = coordinate_dofs[0];
    y[1] = coordinate_dofs[1];
    y[2] = coordinate_dofs[2];
    f.evaluate(vals, y, c);
    values[0] = vals[0];
    y[0] = coordinate_dofs[3];
    y[1] = coordinate_dofs[4];
    y[2] = coordinate_dofs[5];
    f.evaluate(vals, y, c);
    values[1] = vals[0];
    y[0] = coordinate_dofs[6];
    y[1] = coordinate_dofs[7];
    y[2] = coordinate_dofs[8];
    f.evaluate(vals, y, c);
    values[2] = vals[0];
    y[0] = coordinate_dofs[9];
    y[1] = coordinate_dofs[10];
    y[2] = coordinate_dofs[11];
    f.evaluate(vals, y, c);
    values[3] = vals[0];
    y[0] = coordinate_dofs[0];
    y[1] = coordinate_dofs[1];
    y[2] = coordinate_dofs[2];
    f.evaluate(vals, y, c);
    values[4] = vals[1];
    y[0] = coordinate_dofs[3];
    y[1] = coordinate_dofs[4];
    y[2] = coordinate_dofs[5];
    f.evaluate(vals, y, c);
    values[5] = vals[1];
    y[0] = coordinate_dofs[6];
    y[1] = coordinate_dofs[7];
    y[2] = coordinate_dofs[8];
    f.evaluate(vals, y, c);
    values[6] = vals[1];
    y[0] = coordinate_dofs[9];
    y[1] = coordinate_dofs[10];
    y[2] = coordinate_dofs[11];
    f.evaluate(vals, y, c);
    values[7] = vals[1];
    y[0] = coordinate_dofs[0];
    y[1] = coordinate_dofs[1];
    y[2] = coordinate_dofs[2];
    f.evaluate(vals, y, c);
    values[8] = vals[2];
    y[0] = coordinate_dofs[3];
    y[1] = coordinate_dofs[4];
    y[2] = coordinate_dofs[5];
    f.evaluate(vals, y, c);
    values[9] = vals[2];
    y[0] = coordinate_dofs[6];
    y[1] = coordinate_dofs[7];
    y[2] = coordinate_dofs[8];
    f.evaluate(vals, y, c);
    values[10] = vals[2];
    y[0] = coordinate_dofs[9];
    y[1] = coordinate_dofs[10];
    y[2] = coordinate_dofs[11];
    f.evaluate(vals, y, c);
    values[11] = vals[2];
  }

  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[3] = dof_values[1];
    vertex_values[6] = dof_values[2];
    vertex_values[9] = dof_values[3];
    // Evaluate function and change variables
    vertex_values[1] = dof_values[4];
    vertex_values[4] = dof_values[5];
    vertex_values[7] = dof_values[6];
    vertex_values[10] = dof_values[7];
    // Evaluate function and change variables
    vertex_values[2] = dof_values[8];
    vertex_values[5] = dof_values[9];
    vertex_values[8] = dof_values[10];
    vertex_values[11] = dof_values[11];
  }

  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[3] = coordinate_dofs[3];
    dof_coordinates[3 + 1] = coordinate_dofs[4];
    dof_coordinates[3 + 2] = coordinate_dofs[5];
    dof_coordinates[3 * 2] = coordinate_dofs[6];
    dof_coordinates[3 * 2 + 1] = coordinate_dofs[7];
    dof_coordinates[3 * 2 + 2] = coordinate_dofs[8];
    dof_coordinates[3 * 3] = coordinate_dofs[9];
    dof_coordinates[3 * 3 + 1] = coordinate_dofs[10];
    dof_coordinates[3 * 3 + 2] = coordinate_dofs[11];
    dof_coordinates[3 * 4] = coordinate_dofs[0];
    dof_coordinates[3 * 4 + 1] = coordinate_dofs[1];
    dof_coordinates[3 * 4 + 2] = coordinate_dofs[2];
    dof_coordinates[3 * 5] = coordinate_dofs[3];
    dof_coordinates[3 * 5 + 1] = coordinate_dofs[4];
    dof_coordinates[3 * 5 + 2] = coordinate_dofs[5];
    dof_coordinates[3 * 6] = coordinate_dofs[6];
    dof_coordinates[3 * 6 + 1] = coordinate_dofs[7];
    dof_coordinates[3 * 6 + 2] = coordinate_dofs[8];
    dof_coordinates[3 * 7] = coordinate_dofs[9];
    dof_coordinates[3 * 7 + 1] = coordinate_dofs[10];
    dof_coordinates[3 * 7 + 2] = coordinate_dofs[11];
    dof_coordinates[3 * 8] = coordinate_dofs[0];
    dof_coordinates[3 * 8 + 1] = coordinate_dofs[1];
    dof_coordinates[3 * 8 + 2] = coordinate_dofs[2];
    dof_coordinates[3 * 9] = coordinate_dofs[3];
    dof_coordinates[3 * 9 + 1] = coordinate_dofs[4];
    dof_coordinates[3 * 9 + 2] = coordinate_dofs[5];
    dof_coordinates[3 * 10] = coordinate_dofs[6];
    dof_coordinates[3 * 10 + 1] = coordinate_dofs[7];
    dof_coordinates[3 * 10 + 2] = coordinate_dofs[8];
    dof_coordinates[3 * 11] = coordinate_dofs[9];
    dof_coordinates[3 * 11 + 1] = coordinate_dofs[10];
    dof_coordinates[3 * 11 + 2] = coordinate_dofs[11];
  }

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

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

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

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

};


class mylagrange_dofmap_0: public ufc::dofmap
{
public:

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

  ~mylagrange_dofmap_0() override
  {
    // Do nothing
  }

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

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

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

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

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

  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[4] = { 1, 0, 0, 0 };
    if (d >= 4)
        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[4] = { 1, 2, 3, 4 };
    if (d >= 4)
        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];
    dofs[3] = entity_indices[0][3];
  }

  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] = 3;
        break;
    case 2:
        dofs[0] = 0;
        dofs[1] = 1;
        dofs[2] = 3;
        break;
    case 3:
        dofs[0] = 0;
        dofs[1] = 1;
        dofs[2] = 2;
        break;
    }
  }

  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;
        case 3:
            dofs[0] = 3;
            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;
        case 3:
            dofs[0] = 3;
            break;
        }
        break;
    case 1:
        switch (i)
        {
        case 0:
            dofs[0] = 2;
            dofs[1] = 3;
            break;
        case 1:
            dofs[0] = 1;
            dofs[1] = 3;
            break;
        case 2:
            dofs[0] = 1;
            dofs[1] = 2;
            break;
        case 3:
            dofs[0] = 0;
            dofs[1] = 3;
            break;
        case 4:
            dofs[0] = 0;
            dofs[1] = 2;
            break;
        case 5:
            dofs[0] = 0;
            dofs[1] = 1;
            break;
        }
        break;
    case 2:
        switch (i)
        {
        case 0:
            dofs[0] = 1;
            dofs[1] = 2;
            dofs[2] = 3;
            break;
        case 1:
            dofs[0] = 0;
            dofs[1] = 2;
            dofs[2] = 3;
            break;
        case 2:
            dofs[0] = 0;
            dofs[1] = 1;
            dofs[2] = 3;
            break;
        case 3:
            dofs[0] = 0;
            dofs[1] = 1;
            dofs[2] = 2;
            break;
        }
        break;
    case 3:
        switch (i)
        {
        case 0:
            dofs[0] = 0;
            dofs[1] = 1;
            dofs[2] = 2;
            dofs[3] = 3;
            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 mylagrange_dofmap_0();
  }

};


class mylagrange_dofmap_1: public ufc::dofmap
{
public:

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

  ~mylagrange_dofmap_1() override
  {
    // Do nothing
  }

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

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

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

  std::size_t global_dimension(const std::vector<std::size_t>&
                               num_global_entities) const final override
  {
    return 3 * 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 12;
  }

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

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

  std::size_t num_entity_dofs(std::size_t d) const final override
  {
    static const std::size_t return_values[4] = { 3, 0, 0, 0 };
    if (d >= 4)
        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[4] = { 3, 6, 9, 12 };
    if (d >= 4)
        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];
    dofs[3] = offset + entity_indices[0][3];
    offset += num_global_entities[0];
    dofs[4] = offset + entity_indices[0][0];
    dofs[5] = offset + entity_indices[0][1];
    dofs[6] = offset + entity_indices[0][2];
    dofs[7] = offset + entity_indices[0][3];
    offset += num_global_entities[0];
    dofs[8] = offset + entity_indices[0][0];
    dofs[9] = offset + entity_indices[0][1];
    dofs[10] = offset + entity_indices[0][2];
    dofs[11] = offset + entity_indices[0][3];
    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] = 3;
        dofs[3] = 5;
        dofs[4] = 6;
        dofs[5] = 7;
        dofs[6] = 9;
        dofs[7] = 10;
        dofs[8] = 11;
        break;
    case 1:
        dofs[0] = 0;
        dofs[1] = 2;
        dofs[2] = 3;
        dofs[3] = 4;
        dofs[4] = 6;
        dofs[5] = 7;
        dofs[6] = 8;
        dofs[7] = 10;
        dofs[8] = 11;
        break;
    case 2:
        dofs[0] = 0;
        dofs[1] = 1;
        dofs[2] = 3;
        dofs[3] = 4;
        dofs[4] = 5;
        dofs[5] = 7;
        dofs[6] = 8;
        dofs[7] = 9;
        dofs[8] = 11;
        break;
    case 3:
        dofs[0] = 0;
        dofs[1] = 1;
        dofs[2] = 2;
        dofs[3] = 4;
        dofs[4] = 5;
        dofs[5] = 6;
        dofs[6] = 8;
        dofs[7] = 9;
        dofs[8] = 10;
        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] = 4;
            dofs[2] = 8;
            break;
        case 1:
            dofs[0] = 1;
            dofs[1] = 5;
            dofs[2] = 9;
            break;
        case 2:
            dofs[0] = 2;
            dofs[1] = 6;
            dofs[2] = 10;
            break;
        case 3:
            dofs[0] = 3;
            dofs[1] = 7;
            dofs[2] = 11;
            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] = 4;
            dofs[2] = 8;
            break;
        case 1:
            dofs[0] = 1;
            dofs[1] = 5;
            dofs[2] = 9;
            break;
        case 2:
            dofs[0] = 2;
            dofs[1] = 6;
            dofs[2] = 10;
            break;
        case 3:
            dofs[0] = 3;
            dofs[1] = 7;
            dofs[2] = 11;
            break;
        }
        break;
    case 1:
        switch (i)
        {
        case 0:
            dofs[0] = 2;
            dofs[1] = 3;
            dofs[2] = 6;
            dofs[3] = 7;
            dofs[4] = 10;
            dofs[5] = 11;
            break;
        case 1:
            dofs[0] = 1;
            dofs[1] = 3;
            dofs[2] = 5;
            dofs[3] = 7;
            dofs[4] = 9;
            dofs[5] = 11;
            break;
        case 2:
            dofs[0] = 1;
            dofs[1] = 2;
            dofs[2] = 5;
            dofs[3] = 6;
            dofs[4] = 9;
            dofs[5] = 10;
            break;
        case 3:
            dofs[0] = 0;
            dofs[1] = 3;
            dofs[2] = 4;
            dofs[3] = 7;
            dofs[4] = 8;
            dofs[5] = 11;
            break;
        case 4:
            dofs[0] = 0;
            dofs[1] = 2;
            dofs[2] = 4;
            dofs[3] = 6;
            dofs[4] = 8;
            dofs[5] = 10;
            break;
        case 5:
            dofs[0] = 0;
            dofs[1] = 1;
            dofs[2] = 4;
            dofs[3] = 5;
            dofs[4] = 8;
            dofs[5] = 9;
            break;
        }
        break;
    case 2:
        switch (i)
        {
        case 0:
            dofs[0] = 1;
            dofs[1] = 2;
            dofs[2] = 3;
            dofs[3] = 5;
            dofs[4] = 6;
            dofs[5] = 7;
            dofs[6] = 9;
            dofs[7] = 10;
            dofs[8] = 11;
            break;
        case 1:
            dofs[0] = 0;
            dofs[1] = 2;
            dofs[2] = 3;
            dofs[3] = 4;
            dofs[4] = 6;
            dofs[5] = 7;
            dofs[6] = 8;
            dofs[7] = 10;
            dofs[8] = 11;
            break;
        case 2:
            dofs[0] = 0;
            dofs[1] = 1;
            dofs[2] = 3;
            dofs[3] = 4;
            dofs[4] = 5;
            dofs[5] = 7;
            dofs[6] = 8;
            dofs[7] = 9;
            dofs[8] = 11;
            break;
        case 3:
            dofs[0] = 0;
            dofs[1] = 1;
            dofs[2] = 2;
            dofs[3] = 4;
            dofs[4] = 5;
            dofs[5] = 6;
            dofs[6] = 8;
            dofs[7] = 9;
            dofs[8] = 10;
            break;
        }
        break;
    case 3:
        switch (i)
        {
        case 0:
            dofs[0] = 0;
            dofs[1] = 1;
            dofs[2] = 2;
            dofs[3] = 3;
            dofs[4] = 4;
            dofs[5] = 5;
            dofs[6] = 6;
            dofs[7] = 7;
            dofs[8] = 8;
            dofs[9] = 9;
            dofs[10] = 10;
            dofs[11] = 11;
            break;
        }
        break;
    }
  }


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

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

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

};


class mylagrange_coordinate_mapping_1: public ufc::coordinate_mapping
{
public:

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

  ~mylagrange_coordinate_mapping_1() override
  {
    // Do nothing
  }

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

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

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

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

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

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

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

  void compute_physical_coordinates(
      double * x, std::size_t num_points,
      const double * X,
      const double * coordinate_dofs) const final override
  {
    mylagrange_finite_element_0 xelement;
    double phi[4];
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Compute basis values of coordinate element
        xelement.evaluate_reference_basis(phi, 1, &X[3 * ip]);
        // Compute x
        for (std::size_t i = 0; i < 3; ++i)
            for (std::size_t d = 0; d < 4; ++d)
                x[3 * ip + i] += coordinate_dofs[3 * 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[4] = { 1.0, 0.0, 0.0, 0.0 };
    const double dphi_X0[3][4] =
        { { -1.0, 1.0, 0.0, 0.0 },
          { -1.0, 0.0, 1.0, 0.0 },
          { -1.0, 0.0, 0.0, 1.0 } };
    double J[9];
    double detJ[1];
    double K[9];
    double x0[3] = {};
    for (std::size_t i = 0; i < 3; ++i)
        for (std::size_t k = 0; k < 4; ++k)
            x0[i] += coordinate_dofs[3 * k + i] * phi_X0[k];
    for (std::size_t i = 0; i < 3; ++i)
    {
        for (std::size_t j = 0; j < 3; ++j)
        {
            J[3 * i + j] = 0.0;
            for (std::size_t k = 0; k < 4; ++k)
                J[3 * i + j] += coordinate_dofs[3 * 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 < 3; ++j)
            for (std::size_t i = 0; i < 3; ++i)
                X[3 * ip + j] += K[3 * j + i] * (x[3 * 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[4] = { 1.0, 0.0, 0.0, 0.0 };
    const double dphi_X0[3][4] =
        { { -1.0, 1.0, 0.0, 0.0 },
          { -1.0, 0.0, 1.0, 0.0 },
          { -1.0, 0.0, 0.0, 1.0 } };
    double x0[3] = {};
    for (std::size_t i = 0; i < 3; ++i)
        for (std::size_t k = 0; k < 4; ++k)
            x0[i] += coordinate_dofs[3 * k + i] * phi_X0[k];
    for (std::size_t i = 0; i < 3; ++i)
    {
        for (std::size_t j = 0; j < 3; ++j)
        {
            J[3 * i + j] = 0.0;
            for (std::size_t k = 0; k < 4; ++k)
                J[3 * i + j] += coordinate_dofs[3 * 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 < 3; ++j)
            for (std::size_t i = 0; i < 3; ++i)
                X[3 * ip + j] += K[3 * j + i] * (x[3 * ip + i] - x0[i]);
  }

  void compute_jacobians(
      double * J, std::size_t num_points,
      const double * X,
      const double * coordinate_dofs) const final override
  {
    mylagrange_finite_element_0 xelement;
    double dphi[12];
    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[3 * ip]);
        // Compute J
        for (std::size_t i = 0; i < 3; ++i)
            for (std::size_t j = 0; j < 3; ++j)
                for (std::size_t d = 0; d < 4; ++d)
                    J[3 * 3 * ip + 3 * i + j] += coordinate_dofs[3 * d + i] * dphi[3 * 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[3 * 3 * ip] * (J[3 * 3 * ip + 3 + 1] * J[3 * 3 * ip + 3 * 2 + 2] - J[3 * 3 * ip + 3 + 2] * J[3 * 3 * ip + 3 * 2 + 1]) + J[3 * 3 * ip + 1] * (J[3 * 3 * ip + 3 + 2] * J[3 * 3 * ip + 3 * 2] - J[3 * 3 * ip + 3] * J[3 * 3 * ip + 3 * 2 + 2]) + J[3 * 3 * ip + 2] * (J[3 * 3 * ip + 3] * J[3 * 3 * ip + 3 * 2 + 1] - J[3 * 3 * ip + 3 + 1] * J[3 * 3 * ip + 3 * 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[3 * 3 * ip] = (J[3 * 3 * ip + 3 * 2 + 2] * J[3 * 3 * ip + 3 + 1] - J[3 * 3 * ip + 3 + 2] * J[3 * 3 * ip + 3 * 2 + 1]) / detJ[ip];
        K[3 * 3 * ip + 1] = (J[3 * 3 * ip + 2] * J[3 * 3 * ip + 3 * 2 + 1] - J[3 * 3 * ip + 1] * J[3 * 3 * ip + 3 * 2 + 2]) / detJ[ip];
        K[3 * 3 * ip + 2] = (J[3 * 3 * ip + 1] * J[3 * 3 * ip + 3 + 2] - J[3 * 3 * ip + 2] * J[3 * 3 * ip + 3 + 1]) / detJ[ip];
        K[3 * 3 * ip + 3] = (J[3 * 3 * ip + 3 + 2] * J[3 * 3 * ip + 3 * 2] - J[3 * 3 * ip + 3 * 2 + 2] * J[3 * 3 * ip + 3]) / detJ[ip];
        K[3 * 3 * ip + 3 + 1] = (J[3 * 3 * ip + 3 * 2 + 2] * J[3 * 3 * ip] - J[3 * 3 * ip + 2] * J[3 * 3 * ip + 3 * 2]) / detJ[ip];
        K[3 * 3 * ip + 3 + 2] = (J[3 * 3 * ip + 2] * J[3 * 3 * ip + 3] - J[3 * 3 * ip + 3 + 2] * J[3 * 3 * ip]) / detJ[ip];
        K[3 * 3 * ip + 3 * 2] = (J[3 * 3 * ip + 3] * J[3 * 3 * ip + 3 * 2 + 1] - J[3 * 3 * ip + 3 * 2] * J[3 * 3 * ip + 3 + 1]) / detJ[ip];
        K[3 * 3 * ip + 3 * 2 + 1] = (J[3 * 3 * ip + 1] * J[3 * 3 * ip + 3 * 2] - J[3 * 3 * ip] * J[3 * 3 * ip + 3 * 2 + 1]) / detJ[ip];
        K[3 * 3 * ip + 3 * 2 + 2] = (J[3 * 3 * ip] * J[3 * 3 * ip + 3 + 1] - J[3 * 3 * ip + 1] * J[3 * 3 * ip + 3]) / 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[4] = { 0.25000000000000006, 0.25, 0.25, 0.25 };
    const double dphi_Xm[3][4] =
        { { -1.0, 1.0, 0.0, 0.0 },
          { -1.0, 0.0, 1.0, 0.0 },
          { -1.0, 0.0, 0.0, 1.0 } };
    // Compute x
    for (std::size_t i = 0; i < 3; ++i)
        for (std::size_t d = 0; d < 4; ++d)
            x[i] += coordinate_dofs[3 * d + i] * phi_Xm[d];
    // Compute J
    for (std::size_t i = 0; i < 3; ++i)
        for (std::size_t j = 0; j < 3; ++j)
            for (std::size_t d = 0; d < 4; ++d)
                J[3 * i + j] += coordinate_dofs[3 * d + i] * dphi_Xm[j][d];
  }

};


class mylagrange_cell_integral_0_otherwise: public ufc::cell_integral
{
public:

  mylagrange_cell_integral_0_otherwise() : ufc::cell_integral()
  {

  }

  ~mylagrange_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 FE8_C0_D001_Q4[1][1][2] = { { { -1.0, 1.0 } } };
    // Unstructured piecewise computations
    const double J_c0 = coordinate_dofs[0] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[3] * FE8_C0_D001_Q4[0][0][1];
    const double J_c4 = coordinate_dofs[1] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[7] * FE8_C0_D001_Q4[0][0][1];
    const double J_c8 = coordinate_dofs[2] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[11] * FE8_C0_D001_Q4[0][0][1];
    const double J_c5 = coordinate_dofs[1] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[10] * FE8_C0_D001_Q4[0][0][1];
    const double J_c7 = coordinate_dofs[2] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[8] * FE8_C0_D001_Q4[0][0][1];
    const double J_c1 = coordinate_dofs[0] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[6] * FE8_C0_D001_Q4[0][0][1];
    const double J_c6 = coordinate_dofs[2] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[5] * FE8_C0_D001_Q4[0][0][1];
    const double J_c3 = coordinate_dofs[1] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[4] * FE8_C0_D001_Q4[0][0][1];
    const double J_c2 = coordinate_dofs[0] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[9] * FE8_C0_D001_Q4[0][0][1];
    alignas(32) double sp[15];
    sp[0] = J_c4 * J_c8;
    sp[1] = J_c5 * J_c7;
    sp[2] = sp[0] + -1 * sp[1];
    sp[3] = J_c0 * sp[2];
    sp[4] = J_c5 * J_c6;
    sp[5] = J_c3 * J_c8;
    sp[6] = sp[4] + -1 * sp[5];
    sp[7] = J_c1 * sp[6];
    sp[8] = sp[3] + sp[7];
    sp[9] = J_c3 * J_c7;
    sp[10] = J_c4 * J_c6;
    sp[11] = sp[9] + -1 * sp[10];
    sp[12] = J_c2 * sp[11];
    sp[13] = sp[8] + sp[12];
    sp[14] = std::abs(sp[13]);
    A[0] = 0.01666666666666656 * sp[14];
    A[1] = 0.008333333333333295 * sp[14];
    A[2] = 0.008333333333333297 * sp[14];
    A[3] = 0.008333333333333297 * sp[14];
    A[4] = 0.008333333333333295 * sp[14];
    A[5] = 0.01666666666666671 * sp[14];
    A[6] = 0.008333333333333373 * sp[14];
    A[7] = 0.008333333333333373 * sp[14];
    A[8] = 0.008333333333333297 * sp[14];
    A[9] = 0.008333333333333373 * sp[14];
    A[10] = 0.01666666666666671 * sp[14];
    A[11] = 0.008333333333333373 * sp[14];
    A[12] = 0.008333333333333297 * sp[14];
    A[13] = 0.008333333333333373 * sp[14];
    A[14] = 0.008333333333333373 * sp[14];
    A[15] = 0.01666666666666671 * sp[14];
  }

};


class mylagrange_cell_integral_1_otherwise: public ufc::cell_integral
{
public:

  mylagrange_cell_integral_1_otherwise() : ufc::cell_integral()
  {

  }

  ~mylagrange_cell_integral_1_otherwise() override
  {

  }

  const std::vector<bool> & enabled_coefficients() const final override
  {
static const std::vector<bool> enabled({true});
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'
    // Quadrature rules
    alignas(32) static const double weights4[4] = { 0.04166666666666666, 0.04166666666666666, 0.04166666666666666, 0.04166666666666666 };
    // 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 FE8_C0_D001_Q4[1][1][2] = { { { -1.0, 1.0 } } };
    alignas(32) static const double FE8_C0_Q4[1][4][4] =
        { { { 0.1381966011250091, 0.585410196624969, 0.138196601125011, 0.138196601125011 },
            { 0.1381966011250091, 0.138196601125011, 0.585410196624969, 0.138196601125011 },
            { 0.1381966011250091, 0.138196601125011, 0.138196601125011, 0.585410196624969 },
            { 0.5854101966249672, 0.1381966011250109, 0.138196601125011, 0.138196601125011 } } };
    // Unstructured piecewise computations
    const double J_c0 = coordinate_dofs[0] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[3] * FE8_C0_D001_Q4[0][0][1];
    const double J_c4 = coordinate_dofs[1] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[7] * FE8_C0_D001_Q4[0][0][1];
    const double J_c8 = coordinate_dofs[2] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[11] * FE8_C0_D001_Q4[0][0][1];
    const double J_c5 = coordinate_dofs[1] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[10] * FE8_C0_D001_Q4[0][0][1];
    const double J_c7 = coordinate_dofs[2] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[8] * FE8_C0_D001_Q4[0][0][1];
    const double J_c1 = coordinate_dofs[0] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[6] * FE8_C0_D001_Q4[0][0][1];
    const double J_c6 = coordinate_dofs[2] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[5] * FE8_C0_D001_Q4[0][0][1];
    const double J_c3 = coordinate_dofs[1] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[4] * FE8_C0_D001_Q4[0][0][1];
    const double J_c2 = coordinate_dofs[0] * FE8_C0_D001_Q4[0][0][0] + coordinate_dofs[9] * FE8_C0_D001_Q4[0][0][1];
    alignas(32) double sp[15];
    sp[0] = J_c4 * J_c8;
    sp[1] = J_c5 * J_c7;
    sp[2] = sp[0] + -1 * sp[1];
    sp[3] = J_c0 * sp[2];
    sp[4] = J_c5 * J_c6;
    sp[5] = J_c3 * J_c8;
    sp[6] = sp[4] + -1 * sp[5];
    sp[7] = J_c1 * sp[6];
    sp[8] = sp[3] + sp[7];
    sp[9] = J_c3 * J_c7;
    sp[10] = J_c4 * J_c6;
    sp[11] = sp[9] + -1 * sp[10];
    sp[12] = J_c2 * sp[11];
    sp[13] = sp[8] + sp[12];
    sp[14] = std::abs(sp[13]);
    alignas(32) double BF0[4] = {};
    for (int iq = 0; iq < 4; ++iq)
    {
        // Quadrature loop body setup (num_points=4)
        // Unstructured varying computations for num_points=4
        double w0 = 0.0;
        for (int ic = 0; ic < 4; ++ic)
            w0 += w[0][ic] * FE8_C0_Q4[0][iq][ic];
        alignas(32) double sv4[1];
        sv4[0] = sp[14] * w0;
        const double fw0 = sv4[0] * weights4[iq];
        for (int i = 0; i < 4; ++i)
            BF0[i] += fw0 * FE8_C0_Q4[0][iq][i];
    }
    std::fill(A, A + 4, 0.0);
    for (int i = 0; i < 4; ++i)
        A[i] += BF0[i];
  }

};


class mylagrange_form_0: public ufc::form
{
public:

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

  ~mylagrange_form_0() override
  {
    // Do nothing
  }

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

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

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

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

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

  ufc::dofmap * create_dofmap(std::size_t i) const final override
  {
    switch (i)
    {
    case 0:
        return new mylagrange_dofmap_0();
    case 1:
        return new mylagrange_dofmap_0();
    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 mylagrange_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 mylagrange_form_1: public ufc::form
{
public:

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

  ~mylagrange_form_1() override
  {
    // Do nothing
  }

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

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

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

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

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

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

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

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

  ufc::dofmap * create_dofmap(std::size_t i) const final override
  {
    switch (i)
    {
    case 0:
        return new mylagrange_dofmap_0();
    case 1:
        return new mylagrange_dofmap_0();
    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 mylagrange_cell_integral_1_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;
  }

};

// 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 MyLagrange
{

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<mylagrange_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<mylagrange_dofmap_0>(), *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<mylagrange_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<mylagrange_dofmap_0>(), *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<mylagrange_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<mylagrange_dofmap_0>(), *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<mylagrange_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<mylagrange_dofmap_0>(), *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<mylagrange_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<mylagrange_dofmap_0>(), *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<mylagrange_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<mylagrange_dofmap_0>(), *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 mylagrange_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<mylagrange_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<mylagrange_dofmap_0>(), *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<mylagrange_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<mylagrange_dofmap_0>(), *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_g Form_L_FunctionSpace_1;

class Form_L: public dolfin::Form
{
public:

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

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

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

    this->g = g;

    _ufc_form = std::make_shared<const mylagrange_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 == "g")
      return 0;

    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 "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_g;

  // Coefficients
  dolfin::CoefficientAssigner g;
};

class MultiMeshForm_L: public dolfin::MultiMeshForm
{
public:

  // Constructor
  MultiMeshForm_L(std::shared_ptr<const dolfin::MultiMeshFunctionSpace> V0):
    dolfin::MultiMeshForm(V0), g(*this, 0)
  {
    // 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> g):
    dolfin::MultiMeshForm(V0), g(*this, 0)
  {
    // 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->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 == "g")
      return 0;

    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 "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_g;

  // Coefficients
  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