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

class ad1dss_finite_element_0: public ufc::finite_element
{
public:

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

  ~ad1dss_finite_element_0() override
  {
    // Do nothing
  }

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

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

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

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

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

  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][2] = { { 0.7071067811865475, -0.4082482904638631 } };
    static const double coefficients1[1][2] = { { 0.7071067811865475, 0.4082482904638631 } };
    for (std::size_t k = 0; k < num_points * 2; ++k)
        reference_values[k] = 0.0;
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Map from UFC reference coordinate X to FIAT reference coordinate Y
        const double Y[1] = { 2.0 * X[ip] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues1[2] = {};
        basisvalues1[0] = 1.0;
        basisvalues1[1] = Y[0];
        for (std::size_t p = 0; p < 2; ++p)
            basisvalues1[p] *= std::sqrt(0.5 + p);
        // Accumulate products of coefficients and basisvalues
        for (std::size_t r = 0; r < 2; ++r)
            reference_values[2 * ip] += coefficients0[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 2; ++r)
            reference_values[2 * ip + 1] += coefficients1[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(1, order);
    std::fill_n(reference_values, num_points * 2 * num_derivatives, 0.0);
    if (order > 1)
        return;
    // Tables of derivatives of the polynomial base (transpose).
    alignas(32) static const double dmats0[1][2][2] =
        { { { 0.0, 0.0 },
            { 3.464101615137754, 0.0 } } };
    static const double coefficients0[1][2] = { { 0.7071067811865475, -0.4082482904638631 } };
    static const double coefficients1[1][2] = { { 0.7071067811865475, 0.4082482904638631 } };
    const std::size_t reference_offset[2] = {};
    const std::size_t num_components[2] = { 1, 1 };
    // Precomputed combinations
    const std::size_t combinations[1][1][1] = {};
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Map from UFC reference coordinate X to FIAT reference coordinate Y
        const double Y[1] = { 2.0 * X[ip] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues1[2] = {};
        basisvalues1[0] = 1.0;
        basisvalues1[1] = Y[0];
        for (std::size_t p = 0; p < 2; ++p)
            basisvalues1[p] *= std::sqrt(0.5 + p);
        // Loop over all dofs
        for (std::size_t i = 0; i < 2; ++i)
        {
            double derivatives[1] = {};
            switch (i)
            {
            case 0:
                // Compute reference derivatives for dof 0.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[2] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[2][2] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 4, &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[2][2];
                        std::copy_n(&dmats[0][0], 4, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 4, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 2; ++t)
                            for (std::size_t u = 0; u < 2; ++u)
                                for (std::size_t tu = 0; tu < 2; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 2; ++s)
                        for (std::size_t t = 0; t < 2; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 2; ++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[2] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[2][2] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 4, &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[2][2];
                        std::copy_n(&dmats[0][0], 4, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 4, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 2; ++t)
                            for (std::size_t u = 0; u < 2; ++u)
                                for (std::size_t tu = 0; tu < 2; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 2; ++s)
                        for (std::size_t t = 0; t < 2; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 2; ++s)
                        derivatives[r] += coefficients1[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[2 * 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(1, order);
    // Precomputed combinations
    const std::size_t combinations[1][1][1] = {};
    std::fill_n(values, num_points * 2 * num_derivatives, 0.0);
    const std::size_t reference_offsets[2] = {};
    const std::size_t physical_offsets[2] = {};
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        double transform[1][1];
        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[ip + combinations[order - 1][s][k] + combinations[order - 1][r][k]];
        for (std::size_t d = 0; d < 2; ++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[2 * 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[2 * 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[1] = {};
    double J[1];
    double detJ;
    double K[1];
    if (cm)
    {
        cm->compute_reference_geometry(X, J, &detJ, K, 1, x, coordinate_dofs, cell_orientation);
    }
    else
    {
        compute_jacobian_interval_1d(J, coordinate_dofs);
        compute_jacobian_inverse_interval_1d(K, detJ, J);
        // Get coordinates and map to the reference (FIAT) element
        double Y[1] = { (2 * x[0] - coordinate_dofs[0] - coordinate_dofs[1]) / J[0] };
        // Map to FFC reference coordinate
        for (std::size_t k = 0; k < 1; ++k)
            X[k] = (Y[k] + 1.0) / 2.0;
    }
    // Evaluate basis on reference element
    double ref_values[2];
    evaluate_reference_basis(ref_values, 1, X);
    // Push forward
    double physical_values[2];
    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 < 2; ++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(1, 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[1];
    compute_jacobian_interval_1d(J, coordinate_dofs);
    // Compute Inverse Jacobian and determinant
    double K[1];
    double detJ;
    compute_jacobian_inverse_interval_1d(K, detJ, J);
    // Get coordinates and map to the reference (FIAT) element
    double Y[1] = { (2 * x[0] - coordinate_dofs[0] - coordinate_dofs[1]) / J[0] };
    // Precomputed combinations
    const std::size_t combinations[1][1][1] = {};
    // Declare transformation matrix
    double transform[1][1] = { { 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[combinations[n - 1][col][k] + combinations[n - 1][row][k]];
    switch (i)
    {
    case 0:
        {
            double basisvalues[2] = {};
            basisvalues[0] = 1.0;
            basisvalues[1] = Y[0];
            for (std::size_t p = 0; p < 2; ++p)
                basisvalues[p] *= std::sqrt(0.5 + p);
            // Table(s) of coefficients
            static const double coefficients0[2] = { 0.7071067811865475, -0.4082482904638631 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[2][2] =
                { { 0.0, 0.0 },
                  { 3.464101615137754, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[1] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[2][2] =
                { { 1.0, 0.0 },
                  { 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[2][2] =
                { { 1.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], 4, 0.0);
                for (std::size_t t = 0; t < 2; ++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], 4, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 4, 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 < 2; ++t)
                            for (std::size_t u = 0; u < 2; ++u)
                                for (std::size_t tu = 0; tu < 2; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 2; ++s)
                    for (std::size_t t = 0; t < 2; ++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[2] = {};
            basisvalues[0] = 1.0;
            basisvalues[1] = Y[0];
            for (std::size_t p = 0; p < 2; ++p)
                basisvalues[p] *= std::sqrt(0.5 + p);
            // Table(s) of coefficients
            static const double coefficients0[2] = { 0.7071067811865475, 0.4082482904638631 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[2][2] =
                { { 0.0, 0.0 },
                  { 3.464101615137754, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[1] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[2][2] =
                { { 1.0, 0.0 },
                  { 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[2][2] =
                { { 1.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], 4, 0.0);
                for (std::size_t t = 0; t < 2; ++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], 4, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 4, 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 < 2; ++t)
                            for (std::size_t u = 0; u < 2; ++u)
                                for (std::size_t tu = 0; tu < 2; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 2; ++s)
                    for (std::size_t t = 0; t < 2; ++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(1, n);
    // Set values equal to zero.
    std::fill_n(values, num_derivatives * 2, 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[1] = {};
    // Loop dofs and call evaluate_basis_derivatives.
    for (std::size_t r = 0; r < 2; ++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[1];
    switch (i)
    {
    case 0:
        {
            y[0] = coordinate_dofs[0];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 1:
        {
            y[0] = coordinate_dofs[1];
            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[1];
    y[0] = coordinate_dofs[0];
    f.evaluate(vals, y, c);
    values[0] = vals[0];
    y[0] = coordinate_dofs[1];
    f.evaluate(vals, y, c);
    values[1] = 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];
  }

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

  void tabulate_reference_dof_coordinates(double * reference_dof_coordinates) const final override
  {
    static const double dof_X[2] = { 0.0, 1.0 };
    std::copy_n(dof_X, 2, 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 ad1dss_finite_element_0();
  }

};


class ad1dss_finite_element_1: public ufc::finite_element
{
public:

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

  ~ad1dss_finite_element_1() override
  {
    // Do nothing
  }

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

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

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

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

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

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

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

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

  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][2] = { { 0.7071067811865475, -0.4082482904638631 } };
    static const double coefficients1[1][2] = { { 0.7071067811865475, 0.4082482904638631 } };
    for (std::size_t k = 0; k < num_points * 2; ++k)
        reference_values[k] = 0.0;
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Map from UFC reference coordinate X to FIAT reference coordinate Y
        const double Y[1] = { 2.0 * X[ip] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues1[2] = {};
        basisvalues1[0] = 1.0;
        basisvalues1[1] = Y[0];
        for (std::size_t p = 0; p < 2; ++p)
            basisvalues1[p] *= std::sqrt(0.5 + p);
        // Accumulate products of coefficients and basisvalues
        for (std::size_t r = 0; r < 2; ++r)
            reference_values[2 * ip] += coefficients0[0][r] * basisvalues1[r];
        for (std::size_t r = 0; r < 2; ++r)
            reference_values[2 * ip + 1] += coefficients1[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(1, order);
    std::fill_n(reference_values, num_points * 2 * num_derivatives, 0.0);
    if (order > 1)
        return;
    // Tables of derivatives of the polynomial base (transpose).
    alignas(32) static const double dmats0[1][2][2] =
        { { { 0.0, 0.0 },
            { 3.464101615137754, 0.0 } } };
    static const double coefficients0[1][2] = { { 0.7071067811865475, -0.4082482904638631 } };
    static const double coefficients1[1][2] = { { 0.7071067811865475, 0.4082482904638631 } };
    const std::size_t reference_offset[2] = {};
    const std::size_t num_components[2] = { 1, 1 };
    // Precomputed combinations
    const std::size_t combinations[1][1][1] = {};
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Map from UFC reference coordinate X to FIAT reference coordinate Y
        const double Y[1] = { 2.0 * X[ip] - 1.0 };
        // Compute basisvalues for each relevant embedded degree
        double basisvalues1[2] = {};
        basisvalues1[0] = 1.0;
        basisvalues1[1] = Y[0];
        for (std::size_t p = 0; p < 2; ++p)
            basisvalues1[p] *= std::sqrt(0.5 + p);
        // Loop over all dofs
        for (std::size_t i = 0; i < 2; ++i)
        {
            double derivatives[1] = {};
            switch (i)
            {
            case 0:
                // Compute reference derivatives for dof 0.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {
                    double aux[2] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[2][2] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 4, &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[2][2];
                        std::copy_n(&dmats[0][0], 4, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 4, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 2; ++t)
                            for (std::size_t u = 0; u < 2; ++u)
                                for (std::size_t tu = 0; tu < 2; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 2; ++s)
                        for (std::size_t t = 0; t < 2; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 2; ++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[2] = {};
                    // Declare derivative matrix (of polynomial basis).
                    double dmats[2][2] = {};
                    // Initialize dmats.
                    std::size_t comb = combinations[order - 1][r][0];
                    std::copy_n(&dmats0[comb][0][0], 4, &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[2][2];
                        std::copy_n(&dmats[0][0], 4, &dmats_old[0][0]);
                        // Resetting dmats.
                        std::fill_n(&dmats[0][0], 4, 0.0);
                        // Update dmats using an inner product.
                        comb = combinations[order - 1][r][s];
                        for (std::size_t t = 0; t < 2; ++t)
                            for (std::size_t u = 0; u < 2; ++u)
                                for (std::size_t tu = 0; tu < 2; ++tu)
                                    dmats[t][u] += dmats0[comb][t][tu] * dmats_old[tu][u];
                    }
                    for (std::size_t s = 0; s < 2; ++s)
                        for (std::size_t t = 0; t < 2; ++t)
                            aux[s] += dmats[s][t] * basisvalues1[t];
                    derivatives[r] = 0.0;
                    for (std::size_t s = 0; s < 2; ++s)
                        derivatives[r] += coefficients1[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[2 * 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(1, order);
    // Precomputed combinations
    const std::size_t combinations[1][1][1] = {};
    std::fill_n(values, num_points * 2 * num_derivatives, 0.0);
    const std::size_t reference_offsets[2] = {};
    const std::size_t physical_offsets[2] = {};
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        double transform[1][1];
        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[ip + combinations[order - 1][s][k] + combinations[order - 1][r][k]];
        for (std::size_t d = 0; d < 2; ++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[2 * 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[2 * 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[1] = {};
    double J[1];
    double detJ;
    double K[1];
    if (cm)
    {
        cm->compute_reference_geometry(X, J, &detJ, K, 1, x, coordinate_dofs, cell_orientation);
    }
    else
    {
        compute_jacobian_interval_1d(J, coordinate_dofs);
        compute_jacobian_inverse_interval_1d(K, detJ, J);
        // Get coordinates and map to the reference (FIAT) element
        double Y[1] = { (2 * x[0] - coordinate_dofs[0] - coordinate_dofs[1]) / J[0] };
        // Map to FFC reference coordinate
        for (std::size_t k = 0; k < 1; ++k)
            X[k] = (Y[k] + 1.0) / 2.0;
    }
    // Evaluate basis on reference element
    double ref_values[2];
    evaluate_reference_basis(ref_values, 1, X);
    // Push forward
    double physical_values[2];
    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 < 2; ++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(1, 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[1];
    compute_jacobian_interval_1d(J, coordinate_dofs);
    // Compute Inverse Jacobian and determinant
    double K[1];
    double detJ;
    compute_jacobian_inverse_interval_1d(K, detJ, J);
    // Get coordinates and map to the reference (FIAT) element
    double Y[1] = { (2 * x[0] - coordinate_dofs[0] - coordinate_dofs[1]) / J[0] };
    // Precomputed combinations
    const std::size_t combinations[1][1][1] = {};
    // Declare transformation matrix
    double transform[1][1] = { { 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[combinations[n - 1][col][k] + combinations[n - 1][row][k]];
    switch (i)
    {
    case 0:
        {
            double basisvalues[2] = {};
            basisvalues[0] = 1.0;
            basisvalues[1] = Y[0];
            for (std::size_t p = 0; p < 2; ++p)
                basisvalues[p] *= std::sqrt(0.5 + p);
            // Table(s) of coefficients
            static const double coefficients0[2] = { 0.7071067811865475, -0.4082482904638631 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[2][2] =
                { { 0.0, 0.0 },
                  { 3.464101615137754, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[1] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[2][2] =
                { { 1.0, 0.0 },
                  { 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[2][2] =
                { { 1.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], 4, 0.0);
                for (std::size_t t = 0; t < 2; ++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], 4, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 4, 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 < 2; ++t)
                            for (std::size_t u = 0; u < 2; ++u)
                                for (std::size_t tu = 0; tu < 2; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 2; ++s)
                    for (std::size_t t = 0; t < 2; ++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[2] = {};
            basisvalues[0] = 1.0;
            basisvalues[1] = Y[0];
            for (std::size_t p = 0; p < 2; ++p)
                basisvalues[p] *= std::sqrt(0.5 + p);
            // Table(s) of coefficients
            static const double coefficients0[2] = { 0.7071067811865475, 0.4082482904638631 };
            // Tables of derivatives of the polynomial base (transpose).
            static const double dmats0[2][2] =
                { { 0.0, 0.0 },
                  { 3.464101615137754, 0.0 } };
            // Compute reference derivatives.
            // Declare array of derivatives on FIAT element.
            double derivatives[1] = {};
            // Declare derivative matrix (of polynomial basis).
            double dmats[2][2] =
                { { 1.0, 0.0 },
                  { 0.0, 1.0 } };
            // Declare (auxiliary) derivative matrix (of polynomial basis).
            double dmats_old[2][2] =
                { { 1.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], 4, 0.0);
                for (std::size_t t = 0; t < 2; ++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], 4, &dmats_old[0][0]);
                    std::fill_n(&dmats[0][0], 4, 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 < 2; ++t)
                            for (std::size_t u = 0; u < 2; ++u)
                                for (std::size_t tu = 0; tu < 2; ++tu)
                                    dmats[t][u] += dmats_old[tu][u] * dmats0[t][tu];
                    }
                }
                for (std::size_t s = 0; s < 2; ++s)
                    for (std::size_t t = 0; t < 2; ++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(1, n);
    // Set values equal to zero.
    std::fill_n(values, num_derivatives * 2, 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[1] = {};
    // Loop dofs and call evaluate_basis_derivatives.
    for (std::size_t r = 0; r < 2; ++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[1];
    switch (i)
    {
    case 0:
        {
            y[0] = coordinate_dofs[0];
            f.evaluate(vals, y, c);
            return vals[0];
        }
        break;
    case 1:
        {
            y[0] = coordinate_dofs[1];
            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[1];
    y[0] = coordinate_dofs[0];
    f.evaluate(vals, y, c);
    values[0] = vals[0];
    y[0] = coordinate_dofs[1];
    f.evaluate(vals, y, c);
    values[1] = 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];
  }

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

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

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

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

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

};


class ad1dss_finite_element_2: public ufc::finite_element
{
public:

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

  ~ad1dss_finite_element_2() override
  {
    // Do nothing
  }

  const char * signature() const final override
  {
    return "FiniteElement('Real', interval, 0)";
  }

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

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

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

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

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

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

  void evaluate_reference_basis(double * reference_values,
                                std::size_t num_points,
                                const double * X) const final override
  {
    static const double coefficients0[1][1] = { { 1.414213562373095 } };
    for (std::size_t k = 0; k < num_points; ++k)
        reference_values[k] = 0.0;
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        // Compute basisvalues for each relevant embedded degree
        double basisvalues0[1] = {};
        basisvalues0[0] = 1.0;
        for (std::size_t p = 0; p < 1; ++p)
            basisvalues0[p] *= std::sqrt(0.5 + p);
        // Accumulate products of coefficients and basisvalues
        reference_values[ip] += coefficients0[0][0] * basisvalues0[0];
    }
  }

  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(1, order);
    std::fill_n(reference_values, num_points * num_derivatives, 0.0);
    if (order > 0)
        return;
  }

  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
  {
    std::fill_n(values, num_points, 0.0);
    const std::size_t reference_offsets[1] = {};
    const std::size_t physical_offsets[1] = {};
    for (std::size_t ip = 0; ip < num_points; ++ip)
    {
        double transform[1][1];
        for (std::size_t r = 0; r < 1; ++r)
            for (std::size_t s = 0; s < 1; ++s)
                transform[r][s] = 1.0;
        for (std::size_t d = 0; d < 1; ++d)
        {
            for (std::size_t s = 0; s < 1; ++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[ip + d + s + reference_offsets[d]];
                    // Mapping derivatives back to the physical element
                    for (std::size_t r = 0; r < 1; ++r)
                        values[ip + 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[1] = {};
    double J[1];
    double detJ;
    double K[1];
    if (cm)
    {
        cm->compute_reference_geometry(X, J, &detJ, K, 1, x, coordinate_dofs, cell_orientation);
    }
    else
    {
        compute_jacobian_interval_1d(J, coordinate_dofs);
        compute_jacobian_inverse_interval_1d(K, detJ, J);
    }
    // Evaluate basis on reference element
    double ref_values[1];
    evaluate_reference_basis(ref_values, 1, X);
    // Push forward
    double physical_values[1];
    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
  {
    // Element is constant, calling evaluate_basis.
    evaluate_basis(0, values, x, coordinate_dofs, cell_orientation);
  }

  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(1, 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;
    }
  }

  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
  {
    // Element is constant, calling evaluate_basis_derivatives.
    evaluate_basis_derivatives(0, n, values, x, coordinate_dofs, cell_orientation);
  }

  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[1];
    switch (i)
    {
    case 0:
        {
            y[0] = 0.5 * coordinate_dofs[0] + 0.5 * coordinate_dofs[1];
            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[1];
    y[0] = 0.5 * coordinate_dofs[0] + 0.5 * coordinate_dofs[1];
    f.evaluate(vals, y, c);
    values[0] = 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[0];
  }

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

  void tabulate_reference_dof_coordinates(double * reference_dof_coordinates) const final override
  {
    static const double dof_X[1] = { 0.5 };
    std::copy_n(dof_X, 1, 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 ad1dss_finite_element_2();
  }

};


class ad1dss_dofmap_0: public ufc::dofmap
{
public:

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

  ~ad1dss_dofmap_0() override
  {
    // Do nothing
  }

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

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

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

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

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

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

  std::size_t num_entity_dofs(std::size_t d) const final override
  {
    static const std::size_t return_values[2] = { 1, 0 };
    if (d >= 2)
        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[2] = { 1, 2 };
    if (d >= 2)
        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];
  }

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

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

};


class ad1dss_dofmap_1: public ufc::dofmap
{
public:

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

  ~ad1dss_dofmap_1() override
  {
    // Do nothing
  }

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

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

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

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

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

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

  std::size_t num_entity_dofs(std::size_t d) const final override
  {
    static const std::size_t return_values[2] = { 1, 0 };
    if (d >= 2)
        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[2] = { 1, 2 };
    if (d >= 2)
        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];
  }

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

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


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

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

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

};


class ad1dss_dofmap_2: public ufc::dofmap
{
public:

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

  ~ad1dss_dofmap_2() override
  {
    // Do nothing
  }

  const char * signature() const final override
  {
    return "FFC dofmap for FiniteElement('Real', interval, 0)";
  }

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

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

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

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

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

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

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

  std::size_t num_entity_dofs(std::size_t d) const final override
  {
    static const std::size_t return_values[2] = { 0, 1 };
    if (d >= 2)
        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[2] = { 0, 1 };
    if (d >= 2)
        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] = 0;
  }

  void tabulate_facet_dofs(std::size_t * dofs,
                           std::size_t facet) const final override
  {
    // Do nothing
  }

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

};


class ad1dss_coordinate_mapping_1: public ufc::coordinate_mapping
{
public:

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

  ~ad1dss_coordinate_mapping_1() override
  {
    // Do nothing
  }

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

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

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

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

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

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

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

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

  void compute_jacobians(
      double * J, std::size_t num_points,
      const double * X,
      const double * coordinate_dofs) const final override
  {
    ad1dss_finite_element_0 xelement;
    double dphi[2];
    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[ip]);
        // Compute J
        for (std::size_t i = 0; i < 1; ++i)
            for (std::size_t j = 0; j < 1; ++j)
                for (std::size_t d = 0; d < 2; ++d)
                    J[ip + i + j] += coordinate_dofs[d + i] * dphi[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[ip];
  }

  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[ip] = 1.0 / J[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[2] = { 0.5, 0.5 };
    const double dphi_Xm[1][2] = { { -1.0, 1.0 } };
    // Compute x
    for (std::size_t i = 0; i < 1; ++i)
        for (std::size_t d = 0; d < 2; ++d)
            x[i] += coordinate_dofs[d + i] * phi_Xm[d];
    // Compute J
    for (std::size_t i = 0; i < 1; ++i)
        for (std::size_t j = 0; j < 1; ++j)
            for (std::size_t d = 0; d < 2; ++d)
                J[i + j] += coordinate_dofs[d + i] * dphi_Xm[j][d];
  }

};


class ad1dss_cell_integral_0_otherwise: public ufc::cell_integral
{
public:

  ad1dss_cell_integral_0_otherwise() : ufc::cell_integral()
  {

  }

  ~ad1dss_cell_integral_0_otherwise() override
  {

  }

  const std::vector<bool> & enabled_coefficients() const final override
  {
static const std::vector<bool> enabled({true, 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: 1
    // quadrature_rule:   'default'
    // representation:    'uflacs'
    // 
    // and the following integral 0 metadata:
    // 
    // estimated_polynomial_degree: 1
    // optimize:                    True
    // precision:                   16
    // quadrature_degree:           1
    // 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 FE0_C0_D1_Q1[1][1][2] = { { { -1.0, 1.0 } } };
    // Unstructured piecewise computations
    const double J_c0 = coordinate_dofs[0] * FE0_C0_D1_Q1[0][0][0] + coordinate_dofs[1] * FE0_C0_D1_Q1[0][0][1];
    alignas(32) double sp[6];
    sp[0] = 1.0 / J_c0 * w[0][0];
    sp[1] = sp[0] * (1.0 / J_c0);
    sp[2] = 1.0 / J_c0 * w[1][0];
    sp[3] = std::abs(J_c0);
    sp[4] = sp[1] * sp[3];
    sp[5] = sp[2] * sp[3];
    A[0] = -0.5 * sp[5] + sp[4];
    A[1] = 0.5 * sp[5] - sp[4];
    A[2] = -0.5 * sp[5] - sp[4];
    A[3] = 0.5 * sp[5] + sp[4];
  }

};


class ad1dss_cell_integral_1_otherwise: public ufc::cell_integral
{
public:

  ad1dss_cell_integral_1_otherwise() : ufc::cell_integral()
  {

  }

  ~ad1dss_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 weights2[2] = { 0.5, 0.5 };
    // 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 FE0_C0_D1_Q2[1][1][2] = { { { -1.0, 1.0 } } };
    alignas(32) static const double FE0_C0_Q2[1][2][2] =
        { { { 0.7886751345948129, 0.2113248654051871 },
            { 0.2113248654051871, 0.7886751345948129 } } };
    // Unstructured piecewise computations
    const double J_c0 = coordinate_dofs[0] * FE0_C0_D1_Q2[0][0][0] + coordinate_dofs[1] * FE0_C0_D1_Q2[0][0][1];
    alignas(32) double sp[1];
    sp[0] = std::abs(J_c0);
    alignas(32) double BF0[2] = {};
    for (int iq = 0; iq < 2; ++iq)
    {
        // Quadrature loop body setup (num_points=2)
        // Unstructured varying computations for num_points=2
        double w0 = 0.0;
        for (int ic = 0; ic < 2; ++ic)
            w0 += w[0][ic] * FE0_C0_Q2[0][iq][ic];
        alignas(32) double sv2[1];
        sv2[0] = -1 * w0 * sp[0];
        const double fw0 = sv2[0] * weights2[iq];
        for (int i = 0; i < 2; ++i)
            BF0[i] += fw0 * FE0_C0_Q2[0][iq][i];
    }
    std::fill(A, A + 2, 0.0);
    for (int i = 0; i < 2; ++i)
        A[i] += BF0[i];
  }

};


class ad1dss_form_0: public ufc::form
{
public:

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

  ~ad1dss_form_0() override
  {
    // Do nothing
  }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  ufc::cell_integral * create_default_cell_integral() const final override
  {
    return new ad1dss_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 ad1dss_form_1: public ufc::form
{
public:

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

  ~ad1dss_form_1() override
  {
    // Do nothing
  }

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

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

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

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

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

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

class CoefficientSpace_D: public dolfin::FunctionSpace
{
public:

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

  // Constructor for constrained function space
  CoefficientSpace_D(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<ad1dss_finite_element_2>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<ad1dss_dofmap_2>(), *mesh, constrained_domain))
  {
    // Do nothing
  }

};

class CoefficientSpace_f: public dolfin::FunctionSpace
{
public:

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

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

};

class CoefficientSpace_v: public dolfin::FunctionSpace
{
public:

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

  // Constructor for constrained function space
  CoefficientSpace_v(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<ad1dss_finite_element_2>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<ad1dss_dofmap_2>(), *mesh, constrained_domain))
  {
    // Do nothing
  }

};

class Form_a_FunctionSpace_0: public dolfin::FunctionSpace
{
public:

  // Constructor for standard function space
  Form_a_FunctionSpace_0(std::shared_ptr<const dolfin::Mesh> mesh):
    dolfin::FunctionSpace(mesh,
                          std::make_shared<const dolfin::FiniteElement>(std::make_shared<ad1dss_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<ad1dss_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<ad1dss_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<ad1dss_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<ad1dss_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<ad1dss_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<ad1dss_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<ad1dss_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();
  }

};

typedef CoefficientSpace_D Form_a_FunctionSpace_2;

typedef CoefficientSpace_v Form_a_FunctionSpace_3;

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, 2), D(*this, 0), v(*this, 1)
  {
    _function_spaces[0] = V0;
    _function_spaces[1] = V1;

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

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

    this->D = D;
    this->v = v;

    _ufc_form = std::make_shared<const ad1dss_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
  {
    if (name == "D")
      return 0;
    else if (name == "v")
      return 1;

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

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

    dolfin::dolfin_error("generated code for class Form",
                         "access coefficient data",
                         "Invalid coefficient");
    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;
  typedef Form_a_FunctionSpace_2 CoefficientSpace_D;
  typedef Form_a_FunctionSpace_3 CoefficientSpace_v;

  // Coefficients
  dolfin::CoefficientAssigner D;
  dolfin::CoefficientAssigner v;
};

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), D(*this, 0), v(*this, 1)
  {
    // Create and add standard forms
    std::size_t num_parts = V0->num_parts(); // assume all equal and pick first
    for (std::size_t part = 0; part < num_parts; part++)
    {
      std::shared_ptr<dolfin::Form> a(new Form_a(V1->part(part), V0->part(part)));
    add(a);

    }
    // Build multimesh form
    build();

    /// Assign coefficients

  }

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

    }
    // Build multimesh form
    build();

    /// Assign coefficients
    this->D = D;
    this->v = v;

  }

  // Destructor
  ~MultiMeshForm_a()
  {}

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

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

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

    dolfin::dolfin_error("generated code for class Form",
                         "access coefficient data",
                         "Invalid coefficient");
    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;
  typedef Form_a_FunctionSpace_2 CoefficientSpace_D;
  typedef Form_a_FunctionSpace_3 CoefficientSpace_v;

  // Coefficients
  dolfin::MultiMeshCoefficientAssigner D;
  dolfin::MultiMeshCoefficientAssigner v;
};

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<ad1dss_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<ad1dss_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<ad1dss_finite_element_0>()),
                          std::make_shared<const dolfin::DofMap>(std::make_shared<ad1dss_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_f 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), f(*this, 0)
  {
    _function_spaces[0] = V0;

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

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

    this->f = f;

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

  // Destructor
  ~Form_L()
  {}

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

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

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

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

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

  // Coefficients
  dolfin::CoefficientAssigner f;
};

class MultiMeshForm_L: public dolfin::MultiMeshForm
{
public:

  // Constructor
  MultiMeshForm_L(std::shared_ptr<const dolfin::MultiMeshFunctionSpace> V0):
    dolfin::MultiMeshForm(V0), f(*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> f):
    dolfin::MultiMeshForm(V0), f(*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->f = f;

  }

  // Destructor
  ~MultiMeshForm_L()
  {}

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

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

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

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

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

  // Coefficients
  dolfin::MultiMeshCoefficientAssigner f;
};

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