Subject: [Ffc] evaluate_integrand

From: Anders Logg <[email-masked]>
Date: Mon, 12 Apr 2010 14:47:28 +0200

We are doing some work where we need to do run-time quadrature over
arbitrary polyhedra. For this reason, we would like to generate code
for evaluating the integrand of a form at an arbitrary point within
the cell.

I would propose adding something like this to each of the *_integral
classes in the UFC interface:

  /// Evaluate integrand at given point
  virtual void evaluate_integrand(double* A,
                                  const double * const * w,
                                  const cell& c,
                                  const double* coordinates) const = 0;

I suspect/hope this should be easy to add to FFC since the current
quadrature code must do something like this, but in addition iterate
over the points and compute the weighted sum. This would let a user
handle the quadrature manually in cases where that is necessary.

Would it be possible to add?

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 12 Apr 2010 21:19:50 +0800

On 12/04/10 20:47, Anders Logg wrote:
> arbitrary polyhedra.

We (Mehdi and I) do this already (using UFC), so I don't see why a new 
function is required. Can you explain why evaluate_tensor is not enough?

Garth

> For this reason, we would like to generate code
> _______________________________________________
> Mailing list: https://launchpad.net/~ffc
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~ffc
> More help   : https://help.launchpad.net/ListHelp

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 12 Apr 2010 21:21:32 +0800

> function is required. Can you explain why evaluate_tensor is not enough?

I meant 'tabulate_tensor'.

Garth

---

From: Anders Logg <[email-masked]>
Date: Mon, 12 Apr 2010 15:29:47 +0200

> I meant 'tabulate_tensor'.

Which function do you call for evaluating the integrand?

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 12 Apr 2010 21:34:38 +0800

> Which function do you call for evaluating the integrand?

We evaluate it inside ufc::tabulate_tensor. We construct our forms with 
an extra argument, say an object "CutCellIntegrator", which can provide 
quadrature schemes which depend on the considered cell.

Garth

---

From: Anders Logg <[email-masked]>
Date: Mon, 12 Apr 2010 15:49:46 +0200

On Mon, Apr 12, 2010 at 09:34:38PM +0800, Garth N. Wells wrote:
> We evaluate it inside ufc::tabulate_tensor. We construct our forms
> with an extra argument, say an object "CutCellIntegrator", which can
> provide quadrature schemes which depend on the considered cell.

That would require a special purpose code generator. Having
evaluate_integrand would allow more flexibility for users to implement
their own special quadrature scheme.

I also don't understand how it can work at all since the quadrature
points are not known at compile-time. Or do you a fixed set of reference
polyhedra that you map to?

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 12 Apr 2010 22:20:13 +0800

On 12/04/10 21:49, Anders Logg wrote:
> That would require a special purpose code generator.

What's wrong with that? FFC won't (and shouldn't) be able to do 
everything. Just adding a function to UFC won't make FFC do what we do 
now. We reuse FFC (import modules) and add special purpose extensions.

> Having

We make "CutCellIntegrator" an abstract base class, so the user has 
*complete* freedom to define the quadrature scheme and the generated 
code does not depend on the scheme, since the scheme may depend on 
things like how the cell 'cut' is represented.

> polyhedra that you map to?

The "CutCellIntegrator" object is asked for the quadrature points, and 
it returns them. One way to do it is to sub-triangulate polyhedra (I 
think that CGAL can do this). "CutCellIntegrator" would usually be aware 
of intersecting surfaces, etc. It can also return schemes on surfaces.

I would like a quadrature generator on polyhedra to eventually be part 
of DOLFIN.

Garth

---

From: Anders Logg <[email-masked]>
Date: Mon, 12 Apr 2010 17:35:05 +0200

On Mon, Apr 12, 2010 at 10:20:13PM +0800, Garth N. Wells wrote:
> everything. Just adding a function to UFC won't make FFC do what we
> do now. We reuse FFC (import modules) and add special purpose
> extensions.

Exactly, it won't make FFC do what we need, but we could *use* FFC to
do what we need (without adding a special-purpose code generator).

> things like how the cell 'cut' is represented.

Then it sounds to me like that generated code is not at all special,
but instead general purpose and should be added to UFC/FFC.

And the most general interface would (I think) be an interface for
evaluating the integrand at a given point. We already have the same
for basis functions (evaluate_basis_function) so it is a natural
extension.

> The "CutCellIntegrator" object is asked for the quadrature points,
> and it returns them. One way to do it is to sub-triangulate
> polyhedra (I think that CGAL can do this). "CutCellIntegrator" would
> usually be aware of intersecting surfaces, etc. It can also return
> schemes on surfaces.
>
> I would like a quadrature generator on polyhedra to eventually be
> part of DOLFIN.

It is already. See the BarycenterQuadrature class in the Nitsche
repository. I'm currently merging it. It computes the simplest
possible quadrature (exact for linear polynomials) for arbitrary
polyhedra.

--
Anders

---

From: Andre Massing <[email-masked]>
Date: Mon, 12 Apr 2010 18:16:31 +0200

Hi!

Anders and I are working on an implementation for Nitsche's method, where the 
question of integration on cut cell and facets also (beside XFEM for you) came 
up. 

On Monday 12. April 2010 16.20.13 Garth N. Wells wrote:
> >>>>> We (Mehdi and I) do this already (using UFC), so I don't see why a
> >>>>> new function is required. Can you explain why evaluate_tensor is not
> >>>>> enough?
> >>>> 

I think code generator is cool :) But in this case and maybe also other use-
cases for DOLFIN it presents a high barrier to try out new FEM schemes.  
And for example adding an evaluate integrand function would *DOLFIN* make to 
do what we want to to with Nitsche (see below).

> things like how the cell 'cut' is represented.

That's sounds like a pretty cool and sophisticated framework but to me it is 
not clear how our implementation would easily fit into this framework. Although 
there are quite a lot similiarities as integration on cut cells and facets my 
understanding of the pum code is not deep enough to see 
how we could *easily* reuse your code in our case. Maybe you could give a hint 
of how it would possible?  For example how tight is your ffc code generation 
implementation for XFEM connected with the use of the CutCellIntegrator class?
 
Since Nitsche methods can be understand as a kind of domain decompositon we 
have *several* (overlapping) meshes. 
Then we want is to integrate forms (possibly different equation) on each (not 
covered part of the) mesh and some additonal interface integrals along the 
interfaces. Basically we have almost everything in place. We have an 
intersection map, mapping each cell/facet to the covering cells/facets and we 
have at least a barycenter quadrature rule for arbitrary cut cells and facet 
at hand.  

So for trying out another FEM scheme in DOLFIN all we need is to assemble it 
and merging the domain specific matrices. Everything can already be done in 
DOLFIN with resonable effort except the actual calculation of the integrals. 

To summarize, this only one evaluate_integrand  function would help to make 
DOLFIN/FENICS flexible enough for an average programmer like me to implement  
new, semi-automated :) FEM schemes  with a reasonable amount of time.

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 13 Apr 2010 07:17:28 +0800

> for basis functions (evaluate_basis_function) so it is a natural

I still don't see the need for 'evaluate_integrand' unless you plan to 
call it directly from the assembler side (i.e. DOLFIN). Is that the 
case? Perhaps you can give me a concrete example of how you plan to use it.

> possible quadrature (exact for linear polynomials) for arbitrary

Nice. We'll take a look.

Garth

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 13 Apr 2010 07:41:01 +0800

> implementation for XFEM connected with the use of the CutCellIntegrator class?

You wouldn't be able to re-use all the code, but we could share some 
common base classes, for example for locating quadrature points, which 
in turn would depend on how the surface is represented, and for 
integrating on surfaces.

> DOLFIN with resonable effort except the actual calculation of the integrals.

We do this now more-or-less in tabulate_tensor with our PUM object 
providing the necessary data for taking care of the merging.

Garth

---

From: Anders Logg <[email-masked]>
Date: Tue, 13 Apr 2010 09:45:32 +0200

On Tue, Apr 13, 2010 at 07:17:28AM +0800, Garth N. Wells wrote:
> I still don't see the need for 'evaluate_integrand' unless you plan
> to call it directly from the assembler side (i.e. DOLFIN). Is that
> the case? Perhaps you can give me a concrete example of how you plan
> to use it.

Yes, that's the plan. In pseudo-code, this is what we want to do:

  for polyhedron in intersection.cut_cells:

    quadrature_rule = QuadratureRule(polyhedron)
    AK = 0

    for (x, w) in quadrature_rule:

      AK += w * evaluate_integrand(x)

    A += AK

All data representing the geometry, the polyhedra, mappings from those
polyhedra to the original cells etc is in the Intersection class (in
the sandbox):

  intersection = Intersection(mesh0, mesh1)

Eventually, we might want to move part of the functionality into
either DOLFIN or FFC, but having access to and evaluate_integrand
function makes it possible to experiment (from the C++ side) without
the need for building complex abstractions at this point.

--
Anders

---

From: Mehdi Nikbakht <[email-masked]>
Date: Tue, 13 Apr 2010 10:31:51 +0200

> the need for building complex abstractions at this point.

As far as I understood you want to compute the integration points for
polyhedrons inside Dolfin and evaluate_integrands will just compute that
integrand in that specific integration points. If it is the case how
would you determine what is the order of quadrature rule that you want
to use?

Since to evaluate integrands you need to tabulate basis functions and
their derivatives on arbitrary points. How do you want to tabulate basis
functions and their derivatives inside evaluate_integrands?

I think it is better not modifying UFC at this point. You can derive a
class from ufc::foo_integrals which you can pass information on the
geometry of polyhedra to compute the integration points. 

Mehdi

---

From: Anders Logg <[email-masked]>
Date: Tue, 13 Apr 2010 10:43:40 +0200

> would you determine what is the order of quadrature rule that you want

The quadrature rule would be a simple option for the user to
set. Currently we only have one general rule implemented which is
barycentric quadrature.

> functions and their derivatives inside evaluate_integrands?

That's a good point. We would need to evaluate the basis functions and
their derivatives at a given arbitrary point which is not known at
compile-time.

> geometry of polyhedra to compute the integration points.

I don't understand how that helps. Don't you run into the same
problem, that the location of the quadrature points is not known at
compile-time?

--
Anders

---

From: Anders Logg <[email-masked]>
Date: Tue, 13 Apr 2010 11:59:06 +0200

On Tue, Apr 13, 2010 at 11:31:38AM +0200, Mehdi Nikbakht wrote:
> On Tue, 2010-04-13 at 10:43 +0200, Anders Logg wrote:
> Then if you use higher order elements, you need to update your
> quadrature rule manually. I think it would be nice to compute your own
> quadrature rule inside tabualte_tensor by using the standard quadrature
> rule.

Yes, but the problem is that the computation of the quadrature rule is
nontrivial and it might be better to do it from C++ than as part of
the generated code, especially when the code depends on external
libraries like CGAL. See BarycenterQuadrature.cpp.

> Yes, you need to use the tabulate_basis* functions implemented by
> Kristian. Then I am not the only one who is using them. Good for
> Kristian. ;)

ok, good. Then there is no principal problem of generating code for
evaluate_integrand (if it is allowed to call tabaulate_basis).

I don't see what the problem is of adding evaluate_integrand. It is a
natural extension (we have evaluate_basis already), it would be
"simple" to implement (Kristian can correct me if I'm wrong) and it
makes generated code useful for assembly over cut meshes (without the
need for writing a special-purpose code generator).

--
Anders

---

From: Mehdi Nikbakht <[email-masked]>
Date: Tue, 13 Apr 2010 12:51:29 +0200

> libraries like CGAL. See BarycenterQuadrature.cpp.

Where can I find this file?

> need for writing a special-purpose code generator).

I don't think it would be straightforward in the current ffc modules.
You need to tabulate basis functions and their derivatives in the
compile time which differs from what is already in the place. You should
have your own quadrature transformer to generates entries required for
this case. It is similar to what we have done for XFEM. 

Please have a look on the tabulate_Tensor function inside generated code
for Poisson demo in the attached files.

Mehdi

---

From: Anders Logg <[email-masked]>
Date: Tue, 13 Apr 2010 13:10:52 +0200

> Where can I find this file?

In the directory dolfin/quadrature/ in DOLFIN.

> this case. It is similar to what we have done for XFEM.

What is so complicated? The current quadrature code does this:

  1. Tabulate basis functions at quadrature points (precomputed values)

  2. Iterate over quadrature points

  3. Evaluate integrand at each quadrature points

The new function would need to

  1. Evaluate basis functions at given quadrature point

  2. Evaluate integrand at each quadrature points

> for Poisson demo in the attached files.

It looks too complex for me. What I propose is a very simple extension
of UFC (a function that is simpler than what we have now) which would
allow users (such as myself) to handle the complexity elsewhere.

--
Anders

> // This code conforms with the UFC specification version 1.4
> // and was automatically generated by FFC version 0.9.2+.
> //
> // 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:
> //
> //   cache_dir:                      ''
> //   convert_exceptions_to_warnings: False
> //   cpp_optimize:                   False
> //   epsilon:                        1e-14
> //   form_postfix:                   True
> //   format:                         'dolfin'
> //   log_level:                      20
> //   log_prefix:                     ''
> //   optimize:                       False
> //   output_dir:                     '.'
> //   precision:                      15
> //   quadrature_degree:              'auto'
> //   quadrature_rule:                'auto'
> //   representation:                 'quadrature'
> //   split:                          False
>
> #ifndef __POISSON_H
> #define __POISSON_H
>
> #include <cmath>
> #include <algorithm>
> #include <stdexcept>
> #include <fstream>
> #include <boost/assign/list_of.hpp>
> #include <ufc.h>
> #include <pum/GenericPUM.h>
>
> /// This class defines the interface for a finite element.
>
> class poisson_finite_element_0: public ufc::finite_element
> {
> public:
>
>   /// Constructor
>   poisson_finite_element_0() : ufc::finite_element()
>   {
>     // Do nothing
>   }
>
>   /// Destructor
>   virtual ~poisson_finite_element_0()
>   {
>     // Do nothing
>   }
>
>   /// Return a string identifying the finite element
>   virtual const char* signature() const
>   {
>     return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
>   }
>
>   /// Return the cell shape
>   virtual ufc::shape cell_shape() const
>   {
>     return ufc::triangle;
>   }
>
>   /// Return the dimension of the finite element function space
>   virtual unsigned int space_dimension() const
>   {
>     return 3;
>   }
>
>   /// Return the rank of the value space
>   virtual unsigned int value_rank() const
>   {
>     return 0;
>   }
>
>   /// Return the dimension of the value space for axis i
>   virtual unsigned int value_dimension(unsigned int i) const
>   {
>     return 1;
>   }
>
>   /// Evaluate basis function i at given point in cell
>   virtual void evaluate_basis(unsigned int i,
>                               double* values,
>                               const double* coordinates,
>                               const ufc::cell& c) const
>   {
>     // Extract vertex coordinates
>     const double * const * x = c.coordinates;
>
>     // Compute Jacobian of affine map from reference cell
>     const double J_00 = x[1][0] - x[0][0];
>     const double J_01 = x[2][0] - x[0][0];
>     const double J_10 = x[1][1] - x[0][1];
>     const double J_11 = x[2][1] - x[0][1];
>
>     // Compute determinant of Jacobian
>     double detJ = J_00*J_11 - J_01*J_10;
>
>     // Compute inverse of Jacobian
>
>     // Compute constants
>     const double C0 = x[1][0] + x[2][0];
>     const double C1 = x[1][1] + x[2][1];
>
>     // Get coordinates and map to the reference (FIAT) element
>     double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
>     double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
>
>     // Reset values.
>     *values = 0.000000000000000;
>
>     // Map degree of freedom to element degree of freedom
>     const unsigned int dof = i;
>
>     // Array of basisvalues.
>     double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
>
>     // Declare helper variables.
>     unsigned int rr = 0;
>     unsigned int ss = 0;
>     double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
>
>     // Compute basisvalues.
>     basisvalues[0] = 1.000000000000000;
>     basisvalues[1] = tmp0;
>     for (unsigned int r = 0; r < 1; r++)
>     {
>       rr = (r + 1)*(r + 1 + 1)/2 + 1;
>       ss = r*(r + 1)/2;
>       basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
>     }// end loop over 'r'
>     for (unsigned int r = 0; r < 2; r++)
>     {
>       for (unsigned int s = 0; s < 2 - r; s++)
>       {
>         rr = (r + s)*(r + s + 1)/2 + s;
>         basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
>       }// end loop over 's'
>     }// end loop over 'r'
>
>     // Table(s) of coefficients.
>     static const double coefficients0[3][3] = \
>     {{0.471404520791032, -0.288675134594813, -0.166666666666667},
>     {0.471404520791032, 0.288675134594813, -0.166666666666667},
>     {0.471404520791032, 0.000000000000000, 0.333333333333333}};
>
>     // Compute value(s).
>     for (unsigned int r = 0; r < 3; r++)
>     {
>       *values += coefficients0[dof][r]*basisvalues[r];
>     }// end loop over 'r'
>   }
>
>   /// Evaluate all basis functions at given point in cell
>   virtual void evaluate_basis_all(double* values,
>                                   const double* coordinates,
>                                   const ufc::cell& c) const
>   {
>     // Helper variable to hold values of a single dof.
>     double dof_values = 0.000000000000000;
>
>     // Loop dofs and call evaluate_basis.
>     for (unsigned int r = 0; r < 3; r++)
>     {
>       evaluate_basis(r, &dof_values, coordinates, c);
>       values[r] = dof_values;
>     }// end loop over 'r'
>   }
>
>   /// Evaluate order n derivatives of basis function i at given point in cell
>   virtual void evaluate_basis_derivatives(unsigned int i,
>                                           unsigned int n,
>                                           double* values,
>                                           const double* coordinates,
>                                           const ufc::cell& c) const
>   {
>     // Extract vertex coordinates
>     const double * const * x = c.coordinates;
>
>     // Compute Jacobian of affine map from reference cell
>     const double J_00 = x[1][0] - x[0][0];
>     const double J_01 = x[2][0] - x[0][0];
>     const double J_10 = x[1][1] - x[0][1];
>     const double J_11 = x[2][1] - x[0][1];
>
>     // Compute determinant of Jacobian
>     double detJ = J_00*J_11 - J_01*J_10;
>
>     // Compute inverse of Jacobian
>     const double K_00 =  J_11 / detJ;
>     const double K_01 = -J_01 / detJ;
>     const double K_10 = -J_10 / detJ;
>     const double K_11 =  J_00 / detJ;
>
>     // Compute constants
>     const double C0 = x[1][0] + x[2][0];
>     const double C1 = x[1][1] + x[2][1];
>
>     // Get coordinates and map to the reference (FIAT) element
>     double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
>     double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
>
>     // Compute number of derivatives.
>     unsigned int num_derivatives = 1;
>     for (unsigned int r = 0; r < n; r++)
>     {
>       num_derivatives *= 2;
>     }// end loop over 'r'
>
>     // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
>     unsigned int **combinations = new unsigned int *[num_derivatives];
>     for (unsigned int row = 0; row < num_derivatives; row++)
>     {
>       combinations[row] = new unsigned int [n];
>       for (unsigned int col = 0; col < n; col++)
>         combinations[row][col] = 0;
>     }
>
>     // Generate combinations of derivatives
>     for (unsigned int row = 1; row < num_derivatives; row++)
>     {
>       for (unsigned int num = 0; num < row; num++)
>       {
>         for (unsigned int col = n-1; col+1 > 0; col--)
>         {
>           if (combinations[row][col] + 1 > 1)
>             combinations[row][col] = 0;
>           else
>           {
>             combinations[row][col] += 1;
>             break;
>           }
>         }
>       }
>     }
>
>     // Compute inverse of Jacobian
>     const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
>
>     // Declare transformation matrix
>     // Declare pointer to two dimensional array and initialise
>     double **transform = new double *[num_derivatives];
>
>     for (unsigned int j = 0; j < num_derivatives; j++)
>     {
>       transform[j] = new double [num_derivatives];
>       for (unsigned int k = 0; k < num_derivatives; k++)
>         transform[j][k] = 1;
>     }
>
>     // Construct transformation matrix
>     for (unsigned int row = 0; row < num_derivatives; row++)
>     {
>       for (unsigned int col = 0; col < num_derivatives; col++)
>       {
>         for (unsigned int k = 0; k < n; k++)
>           transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
>       }
>     }
>
>     // Reset values. Assuming that values is always an array.
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       values[r] = 0.000000000000000;
>     }// end loop over 'r'
>
>     // Map degree of freedom to element degree of freedom
>     const unsigned int dof = i;
>
>     // Array of basisvalues.
>     double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
>
>     // Declare helper variables.
>     unsigned int rr = 0;
>     unsigned int ss = 0;
>     double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
>
>     // Compute basisvalues.
>     basisvalues[0] = 1.000000000000000;
>     basisvalues[1] = tmp0;
>     for (unsigned int r = 0; r < 1; r++)
>     {
>       rr = (r + 1)*(r + 1 + 1)/2 + 1;
>       ss = r*(r + 1)/2;
>       basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
>     }// end loop over 'r'
>     for (unsigned int r = 0; r < 2; r++)
>     {
>       for (unsigned int s = 0; s < 2 - r; s++)
>       {
>         rr = (r + s)*(r + s + 1)/2 + s;
>         basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
>       }// end loop over 's'
>     }// end loop over 'r'
>
>     // Table(s) of coefficients.
>     static const double coefficients0[3][3] = \
>     {{0.471404520791032, -0.288675134594813, -0.166666666666667},
>     {0.471404520791032, 0.288675134594813, -0.166666666666667},
>     {0.471404520791032, 0.000000000000000, 0.333333333333333}};
>
>     // Tables of derivatives of the polynomial base (transpose).
>     static const double dmats0[3][3] = \
>     {{0.000000000000000, 0.000000000000000, 0.000000000000000},
>     {4.898979485566356, 0.000000000000000, 0.000000000000000},
>     {0.000000000000000, 0.000000000000000, 0.000000000000000}};
>
>     static const double dmats1[3][3] = \
>     {{0.000000000000000, 0.000000000000000, 0.000000000000000},
>     {2.449489742783178, 0.000000000000000, 0.000000000000000},
>     {4.242640687119285, 0.000000000000000, 0.000000000000000}};
>
>     // Compute reference derivatives.
>     // Declare pointer to array of derivatives on FIAT element.
>     double *derivatives = new double[num_derivatives];
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       derivatives[r] = 0.000000000000000;
>     }// end loop over 'r'
>
>     // Declare derivative matrix (of polynomial basis).
>     double dmats[3][3] = \
>     {{1.000000000000000, 0.000000000000000, 0.000000000000000},
>     {0.000000000000000, 1.000000000000000, 0.000000000000000},
>     {0.000000000000000, 0.000000000000000, 1.000000000000000}};
>
>     // Declare (auxiliary) derivative matrix (of polynomial basis).
>     double dmats_old[3][3] = \
>     {{1.000000000000000, 0.000000000000000, 0.000000000000000},
>     {0.000000000000000, 1.000000000000000, 0.000000000000000},
>     {0.000000000000000, 0.000000000000000, 1.000000000000000}};
>
>     // Loop possible derivatives.
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       // Resetting dmats values to compute next derivative.
>       for (unsigned int t = 0; t < 3; t++)
>       {
>         for (unsigned int u = 0; u < 3; u++)
>         {
>           dmats[t][u] = 0.000000000000000;
>           if (t == u)
>           {
>           dmats[t][u] = 1.000000000000000;
>           }
>
>         }// end loop over 'u'
>       }// end loop over 't'
>
>       // Looping derivative order to generate dmats.
>       for (unsigned int s = 0; s < n; s++)
>       {
>         // Updating dmats_old with new values and resetting dmats.
>         for (unsigned int t = 0; t < 3; t++)
>         {
>           for (unsigned int u = 0; u < 3; u++)
>           {
>             dmats_old[t][u] = dmats[t][u];
>             dmats[t][u] = 0.000000000000000;
>           }// end loop over 'u'
>         }// end loop over 't'
>
>         // Update dmats using an inner product.
>         if (combinations[r][s] == 0)
>         {
>         for (unsigned int t = 0; t < 3; t++)
>         {
>           for (unsigned int u = 0; u < 3; u++)
>           {
>             for (unsigned int tu = 0; tu < 3; tu++)
>             {
>               dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
>             }// end loop over 'tu'
>           }// end loop over 'u'
>         }// end loop over 't'
>         }
>
>         if (combinations[r][s] == 1)
>         {
>         for (unsigned int t = 0; t < 3; t++)
>         {
>           for (unsigned int u = 0; u < 3; u++)
>           {
>             for (unsigned int tu = 0; tu < 3; tu++)
>             {
>               dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
>             }// end loop over 'tu'
>           }// end loop over 'u'
>         }// end loop over 't'
>         }
>
>       }// end loop over 's'
>       for (unsigned int s = 0; s < 3; s++)
>       {
>         for (unsigned int t = 0; t < 3; t++)
>         {
>           derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
>         }// end loop over 't'
>       }// end loop over 's'
>     }// end loop over 'r'
>
>     // Transform derivatives back to physical element
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       for (unsigned int s = 0; s < num_derivatives; s++)
>       {
>         values[r] += transform[r][s]*derivatives[s];
>       }// end loop over 's'
>     }// end loop over 'r'
>
>     // Delete pointer to array of derivatives on FIAT element
>     delete [] derivatives;
>
>     // Delete pointer to array of combinations of derivatives and transform
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       delete [] combinations[r];
>     }// end loop over 'r'
>     delete [] combinations;
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       delete [] transform[r];
>     }// end loop over 'r'
>     delete [] transform;
>   }
>
>   /// Evaluate order n derivatives of all basis functions at given point in cell
>   virtual void evaluate_basis_derivatives_all(unsigned int n,
>                                               double* values,
>                                               const double* coordinates,
>                                               const ufc::cell& c) const
>   {
>     // Compute number of derivatives.
>     unsigned int num_derivatives = 1;
>     for (unsigned int r = 0; r < n; r++)
>     {
>       num_derivatives *= 2;
>     }// end loop over 'r'
>
>     // Helper variable to hold values of a single dof.
>     double *dof_values = new double[num_derivatives];
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       dof_values[r] = 0.000000000000000;
>     }// end loop over 'r'
>
>     // Loop dofs and call evaluate_basis_derivatives.
>     for (unsigned int r = 0; r < 3; r++)
>     {
>       evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
>       for (unsigned int s = 0; s < num_derivatives; s++)
>       {
>         values[r*num_derivatives + s] = dof_values[s];
>       }// end loop over 's'
>     }// end loop over 'r'
>
>     // Delete pointer.
>     delete [] dof_values;
>   }
>
>   /// Evaluate linear functional for dof i on the function f
>   virtual double evaluate_dof(unsigned int i,
>                               const ufc::function& f,
>                               const ufc::cell& c) const
>   {
>     // Declare variables for result of evaluation.
>     double vals[1];
>
>     // Declare variable for physical coordinates.
>     double y[2];
>     const double * const * x = c.coordinates;
>     switch (i)
>     {
>     case 0:
>       {
>         y[0] = x[0][0];
>       y[1] = x[0][1];
>       f.evaluate(vals, y, c);
>       return vals[0];
>         break;
>       }
>     case 1:
>       {
>         y[0] = x[1][0];
>       y[1] = x[1][1];
>       f.evaluate(vals, y, c);
>       return vals[0];
>         break;
>       }
>     case 2:
>       {
>         y[0] = x[2][0];
>       y[1] = x[2][1];
>       f.evaluate(vals, y, c);
>       return vals[0];
>         break;
>       }
>     }
>
>     return 0.000000000000000;
>   }
>
>   /// Evaluate linear functionals for all dofs on the function f
>   virtual void evaluate_dofs(double* values,
>                              const ufc::function& f,
>                              const ufc::cell& c) const
>   {
>     // Declare variables for result of evaluation.
>     double vals[1];
>
>     // Declare variable for physical coordinates.
>     double y[2];
>     const double * const * x = c.coordinates;
>     y[0] = x[0][0];
>     y[1] = x[0][1];
>     f.evaluate(vals, y, c);
>     values[0] = vals[0];
>     y[0] = x[1][0];
>     y[1] = x[1][1];
>     f.evaluate(vals, y, c);
>     values[1] = vals[0];
>     y[0] = x[2][0];
>     y[1] = x[2][1];
>     f.evaluate(vals, y, c);
>     values[2] = vals[0];
>   }
>
>   /// Interpolate vertex values from dof values
>   virtual void interpolate_vertex_values(double* vertex_values,
>                                          const double* dof_values,
>                                          const ufc::cell& c) const
>   {
>     // Evaluate function and change variables
>     vertex_values[0] = dof_values[0];
>     vertex_values[1] = dof_values[1];
>     vertex_values[2] = dof_values[2];
>   }
>
>   /// Return the number of sub elements (for a mixed element)
>   virtual unsigned int num_sub_elements() const
>   {
>     return 0;
>   }
>
>   /// Create a new finite element for sub element i (for a mixed element)
>   virtual ufc::finite_element* create_sub_element(unsigned int i) const
>   {
>     return 0;
>   }
>
> };
>
> /// This class defines the interface for a finite element.
>
> class poisson_finite_element_1: public ufc::finite_element
> {
> public:
>
>   /// Constructor
>   poisson_finite_element_1() : ufc::finite_element()
>   {
>     // Do nothing
>   }
>
>   /// Destructor
>   virtual ~poisson_finite_element_1()
>   {
>     // Do nothing
>   }
>
>   /// Return a string identifying the finite element
>   virtual const char* signature() const
>   {
>     return "EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None)))";
>   }
>
>   /// Return the cell shape
>   virtual ufc::shape cell_shape() const
>   {
>     return ufc::triangle;
>   }
>
>   /// Return the dimension of the finite element function space
>   virtual unsigned int space_dimension() const
>   {
>     return 6;
>   }
>
>   /// Return the rank of the value space
>   virtual unsigned int value_rank() const
>   {
>     return 0;
>   }
>
>   /// Return the dimension of the value space for axis i
>   virtual unsigned int value_dimension(unsigned int i) const
>   {
>     return 1;
>   }
>
>   /// Evaluate basis function i at given point in cell
>   virtual void evaluate_basis(unsigned int i,
>                               double* values,
>                               const double* coordinates,
>                               const ufc::cell& c) const
>   {
>     // Extract vertex coordinates
>     const double * const * x = c.coordinates;
>
>     // Compute Jacobian of affine map from reference cell
>     const double J_00 = x[1][0] - x[0][0];
>     const double J_01 = x[2][0] - x[0][0];
>     const double J_10 = x[1][1] - x[0][1];
>     const double J_11 = x[2][1] - x[0][1];
>
>     // Compute determinant of Jacobian
>     double detJ = J_00*J_11 - J_01*J_10;
>
>     // Compute inverse of Jacobian
>
>     // Compute constants
>     const double C0 = x[1][0] + x[2][0];
>     const double C1 = x[1][1] + x[2][1];
>
>     // Get coordinates and map to the reference (FIAT) element
>     double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
>     double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
>
>     // Reset values.
>     values[0] = 0.000000000000000;
>     values[1] = 0.000000000000000;
>     if (0 <= i && i <= 2)
>     {
>       // Map degree of freedom to element degree of freedom
>       const unsigned int dof = i;
>
>       // Array of basisvalues.
>       double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
>
>       // Declare helper variables.
>       unsigned int rr = 0;
>       unsigned int ss = 0;
>       double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
>
>       // Compute basisvalues.
>       basisvalues[0] = 1.000000000000000;
>       basisvalues[1] = tmp0;
>       for (unsigned int r = 0; r < 1; r++)
>       {
>         rr = (r + 1)*(r + 1 + 1)/2 + 1;
>         ss = r*(r + 1)/2;
>         basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
>       }// end loop over 'r'
>       for (unsigned int r = 0; r < 2; r++)
>       {
>         for (unsigned int s = 0; s < 2 - r; s++)
>         {
>           rr = (r + s)*(r + s + 1)/2 + s;
>           basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
>         }// end loop over 's'
>       }// end loop over 'r'
>
>       // Table(s) of coefficients.
>       static const double coefficients0[3][3] = \
>       {{0.471404520791032, -0.288675134594813, -0.166666666666667},
>       {0.471404520791032, 0.288675134594813, -0.166666666666667},
>       {0.471404520791032, 0.000000000000000, 0.333333333333333}};
>
>       // Compute value(s).
>       for (unsigned int r = 0; r < 3; r++)
>       {
>         values[0] += coefficients0[dof][r]*basisvalues[r];
>       }// end loop over 'r'
>     }
>
>     if (3 <= i && i <= 5)
>     {
>       // Map degree of freedom to element degree of freedom
>       const unsigned int dof = i - 3;
>
>       // Array of basisvalues.
>       double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
>
>       // Declare helper variables.
>       unsigned int rr = 0;
>       unsigned int ss = 0;
>       double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
>
>       // Compute basisvalues.
>       basisvalues[0] = 1.000000000000000;
>       basisvalues[1] = tmp0;
>       for (unsigned int r = 0; r < 1; r++)
>       {
>         rr = (r + 1)*(r + 1 + 1)/2 + 1;
>         ss = r*(r + 1)/2;
>         basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
>       }// end loop over 'r'
>       for (unsigned int r = 0; r < 2; r++)
>       {
>         for (unsigned int s = 0; s < 2 - r; s++)
>         {
>           rr = (r + s)*(r + s + 1)/2 + s;
>           basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
>         }// end loop over 's'
>       }// end loop over 'r'
>
>       // Table(s) of coefficients.
>       static const double coefficients0[3][3] = \
>       {{0.471404520791032, -0.288675134594813, -0.166666666666667},
>       {0.471404520791032, 0.288675134594813, -0.166666666666667},
>       {0.471404520791032, 0.000000000000000, 0.333333333333333}};
>
>       // Compute value(s).
>       for (unsigned int r = 0; r < 3; r++)
>       {
>         values[1] += coefficients0[dof][r]*basisvalues[r];
>       }// end loop over 'r'
>     }
>
>   }
>
>   /// Evaluate all basis functions at given point in cell
>   virtual void evaluate_basis_all(double* values,
>                                   const double* coordinates,
>                                   const ufc::cell& c) const
>   {
>     // Helper variable to hold values of a single dof.
>     double dof_values[2] = {0.000000000000000, 0.000000000000000};
>
>     // Loop dofs and call evaluate_basis.
>     for (unsigned int r = 0; r < 6; r++)
>     {
>       evaluate_basis(r, dof_values, coordinates, c);
>       for (unsigned int s = 0; s < 2; s++)
>       {
>         values[r*2 + s] = dof_values[s];
>       }// end loop over 's'
>     }// end loop over 'r'
>   }
>
>   /// Evaluate order n derivatives of basis function i at given point in cell
>   virtual void evaluate_basis_derivatives(unsigned int i,
>                                           unsigned int n,
>                                           double* values,
>                                           const double* coordinates,
>                                           const ufc::cell& c) const
>   {
>     // Extract vertex coordinates
>     const double * const * x = c.coordinates;
>
>     // Compute Jacobian of affine map from reference cell
>     const double J_00 = x[1][0] - x[0][0];
>     const double J_01 = x[2][0] - x[0][0];
>     const double J_10 = x[1][1] - x[0][1];
>     const double J_11 = x[2][1] - x[0][1];
>
>     // Compute determinant of Jacobian
>     double detJ = J_00*J_11 - J_01*J_10;
>
>     // Compute inverse of Jacobian
>     const double K_00 =  J_11 / detJ;
>     const double K_01 = -J_01 / detJ;
>     const double K_10 = -J_10 / detJ;
>     const double K_11 =  J_00 / detJ;
>
>     // Compute constants
>     const double C0 = x[1][0] + x[2][0];
>     const double C1 = x[1][1] + x[2][1];
>
>     // Get coordinates and map to the reference (FIAT) element
>     double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
>     double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
>
>     // Compute number of derivatives.
>     unsigned int num_derivatives = 1;
>     for (unsigned int r = 0; r < n; r++)
>     {
>       num_derivatives *= 2;
>     }// end loop over 'r'
>
>     // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
>     unsigned int **combinations = new unsigned int *[num_derivatives];
>     for (unsigned int row = 0; row < num_derivatives; row++)
>     {
>       combinations[row] = new unsigned int [n];
>       for (unsigned int col = 0; col < n; col++)
>         combinations[row][col] = 0;
>     }
>
>     // Generate combinations of derivatives
>     for (unsigned int row = 1; row < num_derivatives; row++)
>     {
>       for (unsigned int num = 0; num < row; num++)
>       {
>         for (unsigned int col = n-1; col+1 > 0; col--)
>         {
>           if (combinations[row][col] + 1 > 1)
>             combinations[row][col] = 0;
>           else
>           {
>             combinations[row][col] += 1;
>             break;
>           }
>         }
>       }
>     }
>
>     // Compute inverse of Jacobian
>     const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
>
>     // Declare transformation matrix
>     // Declare pointer to two dimensional array and initialise
>     double **transform = new double *[num_derivatives];
>
>     for (unsigned int j = 0; j < num_derivatives; j++)
>     {
>       transform[j] = new double [num_derivatives];
>       for (unsigned int k = 0; k < num_derivatives; k++)
>         transform[j][k] = 1;
>     }
>
>     // Construct transformation matrix
>     for (unsigned int row = 0; row < num_derivatives; row++)
>     {
>       for (unsigned int col = 0; col < num_derivatives; col++)
>       {
>         for (unsigned int k = 0; k < n; k++)
>           transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
>       }
>     }
>
>     // Reset values. Assuming that values is always an array.
>     for (unsigned int r = 0; r < 2*num_derivatives; r++)
>     {
>       values[r] = 0.000000000000000;
>     }// end loop over 'r'
>
>     if (0 <= i && i <= 2)
>     {
>       // Map degree of freedom to element degree of freedom
>       const unsigned int dof = i;
>
>       // Array of basisvalues.
>       double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
>
>       // Declare helper variables.
>       unsigned int rr = 0;
>       unsigned int ss = 0;
>       double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
>
>       // Compute basisvalues.
>       basisvalues[0] = 1.000000000000000;
>       basisvalues[1] = tmp0;
>       for (unsigned int r = 0; r < 1; r++)
>       {
>         rr = (r + 1)*(r + 1 + 1)/2 + 1;
>         ss = r*(r + 1)/2;
>         basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
>       }// end loop over 'r'
>       for (unsigned int r = 0; r < 2; r++)
>       {
>         for (unsigned int s = 0; s < 2 - r; s++)
>         {
>           rr = (r + s)*(r + s + 1)/2 + s;
>           basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
>         }// end loop over 's'
>       }// end loop over 'r'
>
>       // Table(s) of coefficients.
>       static const double coefficients0[3][3] = \
>       {{0.471404520791032, -0.288675134594813, -0.166666666666667},
>       {0.471404520791032, 0.288675134594813, -0.166666666666667},
>       {0.471404520791032, 0.000000000000000, 0.333333333333333}};
>
>       // Tables of derivatives of the polynomial base (transpose).
>       static const double dmats0[3][3] = \
>       {{0.000000000000000, 0.000000000000000, 0.000000000000000},
>       {4.898979485566356, 0.000000000000000, 0.000000000000000},
>       {0.000000000000000, 0.000000000000000, 0.000000000000000}};
>
>       static const double dmats1[3][3] = \
>       {{0.000000000000000, 0.000000000000000, 0.000000000000000},
>       {2.449489742783178, 0.000000000000000, 0.000000000000000},
>       {4.242640687119285, 0.000000000000000, 0.000000000000000}};
>
>       // Compute reference derivatives.
>       // Declare pointer to array of derivatives on FIAT element.
>       double *derivatives = new double[num_derivatives];
>       for (unsigned int r = 0; r < num_derivatives; r++)
>       {
>         derivatives[r] = 0.000000000000000;
>       }// end loop over 'r'
>
>       // Declare derivative matrix (of polynomial basis).
>       double dmats[3][3] = \
>       {{1.000000000000000, 0.000000000000000, 0.000000000000000},
>       {0.000000000000000, 1.000000000000000, 0.000000000000000},
>       {0.000000000000000, 0.000000000000000, 1.000000000000000}};
>
>       // Declare (auxiliary) derivative matrix (of polynomial basis).
>       double dmats_old[3][3] = \
>       {{1.000000000000000, 0.000000000000000, 0.000000000000000},
>       {0.000000000000000, 1.000000000000000, 0.000000000000000},
>       {0.000000000000000, 0.000000000000000, 1.000000000000000}};
>
>       // Loop possible derivatives.
>       for (unsigned int r = 0; r < num_derivatives; r++)
>       {
>         // Resetting dmats values to compute next derivative.
>         for (unsigned int t = 0; t < 3; t++)
>         {
>           for (unsigned int u = 0; u < 3; u++)
>           {
>             dmats[t][u] = 0.000000000000000;
>             if (t == u)
>             {
>             dmats[t][u] = 1.000000000000000;
>             }
>
>           }// end loop over 'u'
>         }// end loop over 't'
>
>         // Looping derivative order to generate dmats.
>         for (unsigned int s = 0; s < n; s++)
>         {
>           // Updating dmats_old with new values and resetting dmats.
>           for (unsigned int t = 0; t < 3; t++)
>           {
>             for (unsigned int u = 0; u < 3; u++)
>             {
>               dmats_old[t][u] = dmats[t][u];
>               dmats[t][u] = 0.000000000000000;
>             }// end loop over 'u'
>           }// end loop over 't'
>
>           // Update dmats using an inner product.
>           if (combinations[r][s] == 0)
>           {
>           for (unsigned int t = 0; t < 3; t++)
>           {
>             for (unsigned int u = 0; u < 3; u++)
>             {
>               for (unsigned int tu = 0; tu < 3; tu++)
>               {
>                 dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
>               }// end loop over 'tu'
>             }// end loop over 'u'
>           }// end loop over 't'
>           }
>
>           if (combinations[r][s] == 1)
>           {
>           for (unsigned int t = 0; t < 3; t++)
>           {
>             for (unsigned int u = 0; u < 3; u++)
>             {
>               for (unsigned int tu = 0; tu < 3; tu++)
>               {
>                 dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
>               }// end loop over 'tu'
>             }// end loop over 'u'
>           }// end loop over 't'
>           }
>
>         }// end loop over 's'
>         for (unsigned int s = 0; s < 3; s++)
>         {
>           for (unsigned int t = 0; t < 3; t++)
>           {
>             derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
>           }// end loop over 't'
>         }// end loop over 's'
>       }// end loop over 'r'
>
>       // Transform derivatives back to physical element
>       for (unsigned int r = 0; r < num_derivatives; r++)
>       {
>         for (unsigned int s = 0; s < num_derivatives; s++)
>         {
>           values[r] += transform[r][s]*derivatives[s];
>         }// end loop over 's'
>       }// end loop over 'r'
>
>       // Delete pointer to array of derivatives on FIAT element
>       delete [] derivatives;
>
>       // Delete pointer to array of combinations of derivatives and transform
>       for (unsigned int r = 0; r < num_derivatives; r++)
>       {
>         delete [] combinations[r];
>       }// end loop over 'r'
>       delete [] combinations;
>       for (unsigned int r = 0; r < num_derivatives; r++)
>       {
>         delete [] transform[r];
>       }// end loop over 'r'
>       delete [] transform;
>     }
>
>     if (3 <= i && i <= 5)
>     {
>       // Map degree of freedom to element degree of freedom
>       const unsigned int dof = i - 3;
>
>       // Array of basisvalues.
>       double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
>
>       // Declare helper variables.
>       unsigned int rr = 0;
>       unsigned int ss = 0;
>       double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
>
>       // Compute basisvalues.
>       basisvalues[0] = 1.000000000000000;
>       basisvalues[1] = tmp0;
>       for (unsigned int r = 0; r < 1; r++)
>       {
>         rr = (r + 1)*(r + 1 + 1)/2 + 1;
>         ss = r*(r + 1)/2;
>         basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
>       }// end loop over 'r'
>       for (unsigned int r = 0; r < 2; r++)
>       {
>         for (unsigned int s = 0; s < 2 - r; s++)
>         {
>           rr = (r + s)*(r + s + 1)/2 + s;
>           basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
>         }// end loop over 's'
>       }// end loop over 'r'
>
>       // Table(s) of coefficients.
>       static const double coefficients0[3][3] = \
>       {{0.471404520791032, -0.288675134594813, -0.166666666666667},
>       {0.471404520791032, 0.288675134594813, -0.166666666666667},
>       {0.471404520791032, 0.000000000000000, 0.333333333333333}};
>
>       // Tables of derivatives of the polynomial base (transpose).
>       static const double dmats0[3][3] = \
>       {{0.000000000000000, 0.000000000000000, 0.000000000000000},
>       {4.898979485566356, 0.000000000000000, 0.000000000000000},
>       {0.000000000000000, 0.000000000000000, 0.000000000000000}};
>
>       static const double dmats1[3][3] = \
>       {{0.000000000000000, 0.000000000000000, 0.000000000000000},
>       {2.449489742783178, 0.000000000000000, 0.000000000000000},
>       {4.242640687119285, 0.000000000000000, 0.000000000000000}};
>
>       // Compute reference derivatives.
>       // Declare pointer to array of derivatives on FIAT element.
>       double *derivatives = new double[num_derivatives];
>       for (unsigned int r = 0; r < num_derivatives; r++)
>       {
>         derivatives[r] = 0.000000000000000;
>       }// end loop over 'r'
>
>       // Declare derivative matrix (of polynomial basis).
>       double dmats[3][3] = \
>       {{1.000000000000000, 0.000000000000000, 0.000000000000000},
>       {0.000000000000000, 1.000000000000000, 0.000000000000000},
>       {0.000000000000000, 0.000000000000000, 1.000000000000000}};
>
>       // Declare (auxiliary) derivative matrix (of polynomial basis).
>       double dmats_old[3][3] = \
>       {{1.000000000000000, 0.000000000000000, 0.000000000000000},
>       {0.000000000000000, 1.000000000000000, 0.000000000000000},
>       {0.000000000000000, 0.000000000000000, 1.000000000000000}};
>
>       // Loop possible derivatives.
>       for (unsigned int r = 0; r < num_derivatives; r++)
>       {
>         // Resetting dmats values to compute next derivative.
>         for (unsigned int t = 0; t < 3; t++)
>         {
>           for (unsigned int u = 0; u < 3; u++)
>           {
>             dmats[t][u] = 0.000000000000000;
>             if (t == u)
>             {
>             dmats[t][u] = 1.000000000000000;
>             }
>
>           }// end loop over 'u'
>         }// end loop over 't'
>
>         // Looping derivative order to generate dmats.
>         for (unsigned int s = 0; s < n; s++)
>         {
>           // Updating dmats_old with new values and resetting dmats.
>           for (unsigned int t = 0; t < 3; t++)
>           {
>             for (unsigned int u = 0; u < 3; u++)
>             {
>               dmats_old[t][u] = dmats[t][u];
>               dmats[t][u] = 0.000000000000000;
>             }// end loop over 'u'
>           }// end loop over 't'
>
>           // Update dmats using an inner product.
>           if (combinations[r][s] == 0)
>           {
>           for (unsigned int t = 0; t < 3; t++)
>           {
>             for (unsigned int u = 0; u < 3; u++)
>             {
>               for (unsigned int tu = 0; tu < 3; tu++)
>               {
>                 dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
>               }// end loop over 'tu'
>             }// end loop over 'u'
>           }// end loop over 't'
>           }
>
>           if (combinations[r][s] == 1)
>           {
>           for (unsigned int t = 0; t < 3; t++)
>           {
>             for (unsigned int u = 0; u < 3; u++)
>             {
>               for (unsigned int tu = 0; tu < 3; tu++)
>               {
>                 dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
>               }// end loop over 'tu'
>             }// end loop over 'u'
>           }// end loop over 't'
>           }
>
>         }// end loop over 's'
>         for (unsigned int s = 0; s < 3; s++)
>         {
>           for (unsigned int t = 0; t < 3; t++)
>           {
>             derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
>           }// end loop over 't'
>         }// end loop over 's'
>       }// end loop over 'r'
>
>       // Transform derivatives back to physical element
>       for (unsigned int r = 0; r < num_derivatives; r++)
>       {
>         for (unsigned int s = 0; s < num_derivatives; s++)
>         {
>           values[num_derivatives + r] += transform[r][s]*derivatives[s];
>         }// end loop over 's'
>       }// end loop over 'r'
>
>       // Delete pointer to array of derivatives on FIAT element
>       delete [] derivatives;
>
>       // Delete pointer to array of combinations of derivatives and transform
>       for (unsigned int r = 0; r < num_derivatives; r++)
>       {
>         delete [] combinations[r];
>       }// end loop over 'r'
>       delete [] combinations;
>       for (unsigned int r = 0; r < num_derivatives; r++)
>       {
>         delete [] transform[r];
>       }// end loop over 'r'
>       delete [] transform;
>     }
>
>   }
>
>   /// Evaluate order n derivatives of all basis functions at given point in cell
>   virtual void evaluate_basis_derivatives_all(unsigned int n,
>                                               double* values,
>                                               const double* coordinates,
>                                               const ufc::cell& c) const
>   {
>     // Compute number of derivatives.
>     unsigned int num_derivatives = 1;
>     for (unsigned int r = 0; r < n; r++)
>     {
>       num_derivatives *= 2;
>     }// end loop over 'r'
>
>     // Helper variable to hold values of a single dof.
>     double *dof_values = new double[2*num_derivatives];
>     for (unsigned int r = 0; r < 2*num_derivatives; r++)
>     {
>       dof_values[r] = 0.000000000000000;
>     }// end loop over 'r'
>
>     // Loop dofs and call evaluate_basis_derivatives.
>     for (unsigned int r = 0; r < 6; r++)
>     {
>       evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
>       for (unsigned int s = 0; s < 2*num_derivatives; s++)
>       {
>         values[r*2*num_derivatives + s] = dof_values[s];
>       }// end loop over 's'
>     }// end loop over 'r'
>
>     // Delete pointer.
>     delete [] dof_values;
>   }
>
>   /// Evaluate linear functional for dof i on the function f
>   virtual double evaluate_dof(unsigned int i,
>                               const ufc::function& f,
>                               const ufc::cell& c) const
>   {
>     // Declare variables for result of evaluation.
>     double vals[1];
>
>     // Declare variable for physical coordinates.
>     double y[2];
>     const double * const * x = c.coordinates;
>     switch (i)
>     {
>     case 0:
>       {
>         y[0] = x[0][0];
>       y[1] = x[0][1];
>       f.evaluate(vals, y, c);
>       return vals[0];
>         break;
>       }
>     case 1:
>       {
>         y[0] = x[1][0];
>       y[1] = x[1][1];
>       f.evaluate(vals, y, c);
>       return vals[0];
>         break;
>       }
>     case 2:
>       {
>         y[0] = x[2][0];
>       y[1] = x[2][1];
>       f.evaluate(vals, y, c);
>       return vals[0];
>         break;
>       }
>     case 3:
>       {
>         y[0] = x[0][0];
>       y[1] = x[0][1];
>       f.evaluate(vals, y, c);
>       return vals[0];
>         break;
>       }
>     case 4:
>       {
>         y[0] = x[1][0];
>       y[1] = x[1][1];
>       f.evaluate(vals, y, c);
>       return vals[0];
>         break;
>       }
>     case 5:
>       {
>         y[0] = x[2][0];
>       y[1] = x[2][1];
>       f.evaluate(vals, y, c);
>       return vals[0];
>         break;
>       }
>     }
>
>     return 0.000000000000000;
>   }
>
>   /// Evaluate linear functionals for all dofs on the function f
>   virtual void evaluate_dofs(double* values,
>                              const ufc::function& f,
>                              const ufc::cell& c) const
>   {
>     // Declare variables for result of evaluation.
>     double vals[1];
>
>     // Declare variable for physical coordinates.
>     double y[2];
>     const double * const * x = c.coordinates;
>     y[0] = x[0][0];
>     y[1] = x[0][1];
>     f.evaluate(vals, y, c);
>     values[0] = vals[0];
>     y[0] = x[1][0];
>     y[1] = x[1][1];
>     f.evaluate(vals, y, c);
>     values[1] = vals[0];
>     y[0] = x[2][0];
>     y[1] = x[2][1];
>     f.evaluate(vals, y, c);
>     values[2] = vals[0];
>     y[0] = x[0][0];
>     y[1] = x[0][1];
>     f.evaluate(vals, y, c);
>     values[3] = vals[0];
>     y[0] = x[1][0];
>     y[1] = x[1][1];
>     f.evaluate(vals, y, c);
>     values[4] = vals[0];
>     y[0] = x[2][0];
>     y[1] = x[2][1];
>     f.evaluate(vals, y, c);
>     values[5] = vals[0];
>   }
>
>   /// Interpolate vertex values from dof values
>   virtual void interpolate_vertex_values(double* vertex_values,
>                                          const double* dof_values,
>                                          const ufc::cell& c) const
>   {
>     // Evaluate function and change variables
>     vertex_values[0] = dof_values[0] + dof_values[3];
>     vertex_values[1] = dof_values[1] + dof_values[4];
>     vertex_values[2] = dof_values[2] + dof_values[5];
>   }
>
>   /// Return the number of sub elements (for a mixed element)
>   virtual unsigned int num_sub_elements() const
>   {
>     return 0;
>   }
>
>   /// Create a new finite element for sub element i (for a mixed element)
>   virtual ufc::finite_element* create_sub_element(unsigned int i) const
>   {
>     return 0;
>   }
>
> };
>
> /// This class defines the interface for a local-to-global mapping of
> /// degrees of freedom (dofs).
>
>
> class poisson_dof_map_0: public ufc::dof_map
> {
> private:
>
>   unsigned int _global_dimension;
>
> public:
>
>   /// Constructor
>   poisson_dof_map_0(    ) :ufc::dof_map()
>   {
>     _global_dimension = 0;
>   }
>   /// Destructor
>   virtual ~poisson_dof_map_0()
>   {
>     // Do nothing
>   }
>
>   /// Return a string identifying the dof map
>   virtual const char* signature() const
>   {
>     return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
>   }
>
>   /// Return true iff mesh entities of topological dimension d are needed
>   virtual bool needs_mesh_entities(unsigned int d) const
>   {
>     switch (d)
>     {
>     case 0:
>       {
>         return true;
>         break;
>       }
>     case 1:
>       {
>         return false;
>         break;
>       }
>     case 2:
>       {
>         return false;
>         break;
>       }
>     }
>
>     return false;
>   }
>
>   /// Initialize dof map for mesh (return true iff init_cell() is needed)
>   virtual bool init_mesh(const ufc::mesh& m)
>   {
>     _global_dimension = m.num_entities[0];
>     return false;
>   }
>
>   /// Initialize dof map for given cell
>   virtual void init_cell(const ufc::mesh& m,
>                          const ufc::cell& c)
>   {
>     // Do nothing
>   }
>
>   /// Finish initialization of dof map for cells
>   virtual void init_cell_finalize()
>   {
>     // Do nothing
>   }
>
>   /// Return the dimension of the global finite element function space
>   virtual unsigned int global_dimension() const
>   {
>     return _global_dimension;
>   }
>
>   /// Return the dimension of the local finite element function space for a cell
>   virtual unsigned int local_dimension(const ufc::cell& c) const
>   {
>     return 3;
>   }
>
>   /// Return the maximum dimension of the local finite element function space
>   virtual unsigned int max_local_dimension() const
>   {
>     return 3;
>   }
>
>
>   // Return the geometric dimension of the coordinates this dof map provides
>   virtual unsigned int geometric_dimension() const
>   {
>     return 2;
>   }
>
>   /// Return the number of dofs on each cell facet
>   virtual unsigned int num_facet_dofs() const
>   {
>     return 2;
>   }
>
>   /// Return the number of dofs associated with each cell entity of dimension d
>   virtual unsigned int num_entity_dofs(unsigned int d) const
>   {
>     switch (d)
>     {
>     case 0:
>       {
>         return 1;
>         break;
>       }
>     case 1:
>       {
>         return 0;
>         break;
>       }
>     case 2:
>       {
>         return 0;
>         break;
>       }
>     }
>
>     return 0;
>   }
>
>   /// Tabulate the local-to-global mapping of dofs on a cell
>   virtual void tabulate_dofs(unsigned int* dofs,
>                              const ufc::mesh& m,
>                              const ufc::cell& c) const
>   {
>     dofs[0] = c.entity_indices[0][0];
>     dofs[1] = c.entity_indices[0][1];
>     dofs[2] = c.entity_indices[0][2];
>   }
>
>   /// Tabulate the local-to-local mapping from facet dofs to cell dofs
>   virtual void tabulate_facet_dofs(unsigned int* dofs,
>                                    unsigned int facet) const
>   {
>     switch (facet)
>     {
>     case 0:
>       {
>         dofs[0] = 1;
>       dofs[1] = 2;
>         break;
>       }
>     case 1:
>       {
>         dofs[0] = 0;
>       dofs[1] = 2;
>         break;
>       }
>     case 2:
>       {
>         dofs[0] = 0;
>       dofs[1] = 1;
>         break;
>       }
>     }
>
>   }
>
>   /// Tabulate the local-to-local mapping of dofs on entity (d, i)
>   virtual void tabulate_entity_dofs(unsigned int* dofs,
>                                     unsigned int d, unsigned int i) const
>   {
>     if (d > 2)
>     {
>     throw std::runtime_error("d is larger than dimension (2)");
>     }
>
>     switch (d)
>     {
>     case 0:
>       {
>         if (i > 2)
>       {
>       throw std::runtime_error("i is larger than number of entities (2)");
>       }
>
>       switch (i)
>       {
>       case 0:
>         {
>           dofs[0] = 0;
>           break;
>         }
>       case 1:
>         {
>           dofs[0] = 1;
>           break;
>         }
>       case 2:
>         {
>           dofs[0] = 2;
>           break;
>         }
>       }
>
>         break;
>       }
>     case 1:
>       {
>
>         break;
>       }
>     case 2:
>       {
>
>         break;
>       }
>     }
>
>   }
>
>   /// Tabulate the coordinates of all dofs on a cell
>   virtual void tabulate_coordinates(double** coordinates,
>                                     const ufc::cell& c) const
>   {
>     const double * const * x = c.coordinates;
>
>     coordinates[0][0] = x[0][0];
>     coordinates[0][1] = x[0][1];
>     coordinates[1][0] = x[1][0];
>     coordinates[1][1] = x[1][1];
>     coordinates[2][0] = x[2][0];
>     coordinates[2][1] = x[2][1];
>   }
>
>   /// Return the number of sub dof maps (for a mixed element)
>   virtual unsigned int num_sub_dof_maps() const
>   {
>     return 0;
>   }
>
>   /// Create a new dof_map for sub dof map i (for a mixed element)
>   virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
>   {
>     return 0;
>   }
>
> };
>
> /// This class defines the interface for a local-to-global mapping of
> /// degrees of freedom (dofs).
>
>
> class poisson_dof_map_1: public ufc::dof_map
> {
> private:
>
>   unsigned int _global_dimension;
>   const std::vector<const pum::GenericPUM*>& pums;
> public:
>
>   /// Constructor
>   poisson_dof_map_1(    const std::vector<const pum::GenericPUM*>& pums) :ufc::dof_map()    , pums(pums)
>   {
>     _global_dimension = 0;
>   }
>   /// Destructor
>   virtual ~poisson_dof_map_1()
>   {
>     // Do nothing
>   }
>
>   /// Return a string identifying the dof map
>   virtual const char* signature() const
>   {
>     return "FFC dofmap for EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None)))";
>   }
>
>   /// Return true iff mesh entities of topological dimension d are needed
>   virtual bool needs_mesh_entities(unsigned int d) const
>   {
>     switch (d)
>     {
>     case 0:
>       {
>         return true;
>         break;
>       }
>     case 1:
>       {
>         return false;
>         break;
>       }
>     case 2:
>       {
>         return false;
>         break;
>       }
>     }
>
>     return false;
>   }
>
>   /// Initialize dof map for mesh (return true iff init_cell() is needed)
>   virtual bool init_mesh(const ufc::mesh& m)
>   {
>     _global_dimension = m.num_entities[0];
>     return false;
>   }
>
>   /// Initialize dof map for given cell
>   virtual void init_cell(const ufc::mesh& m,
>                          const ufc::cell& c)
>   {
>     // Do nothing
>   }
>
>   /// Finish initialization of dof map for cells
>   virtual void init_cell_finalize()
>   {
>     // Do nothing
>   }
>
>   /// Return the dimension of the global finite element function space
>   virtual unsigned int global_dimension() const
>   {
>     return _global_dimension + pums[0]->enriched_global_dimension();
>   }
>
>   /// Return the dimension of the local finite element function space for a cell
>   virtual unsigned int local_dimension(const ufc::cell& c) const
>   {
>     return 3 + pums[0]->enriched_local_dimension(c);
>   }
>
>   /// Return the maximum dimension of the local finite element function space
>   virtual unsigned int max_local_dimension() const
>   {
>     return 3 + pums[0]->enriched_max_local_dimension();
>   }
>
>
>   // Return the geometric dimension of the coordinates this dof map provides
>   virtual unsigned int geometric_dimension() const
>   {
>     return 2;
>   }
>
>   /// Return the number of dofs on each cell facet
>   virtual unsigned int num_facet_dofs() const
>   {
>     return 2;
>   }
>
>   /// Return the number of dofs associated with each cell entity of dimension d
>   virtual unsigned int num_entity_dofs(unsigned int d) const
>   {
>     switch (d)
>     {
>     case 0:
>       {
>         return 1;
>         break;
>       }
>     case 1:
>       {
>         return 0;
>         break;
>       }
>     case 2:
>       {
>         return 0;
>         break;
>       }
>     }
>
>     return 0;
>   }
>
>   /// Tabulate the local-to-global mapping of dofs on a cell
>   virtual void tabulate_dofs(unsigned int* dofs,
>                              const ufc::mesh& m,
>                              const ufc::cell& c) const
>   {
>     dofs[0] = c.entity_indices[0][0];
>     dofs[1] = c.entity_indices[0][1];
>     dofs[2] = c.entity_indices[0][2];
>     \
>
>     // Generate code for tabulating extra degrees of freedom.
>     unsigned int local_offset = 3;
>     unsigned int global_offset = m.num_entities[0];
>
>     // Calculate local-to-global mapping for the enriched dofs related to the discontinuous field 0
>     pums[0]->tabulate_enriched_dofs(dofs, c, local_offset, global_offset);
>   }
>
>   /// Tabulate the local-to-local mapping from facet dofs to cell dofs
>   virtual void tabulate_facet_dofs(unsigned int* dofs,
>                                    unsigned int facet) const
>   {
>     switch (facet)
>     {
>     case 0:
>       {
>         dofs[0] = 1;
>       dofs[1] = 2;
>         break;
>       }
>     case 1:
>       {
>         dofs[0] = 0;
>       dofs[1] = 2;
>         break;
>       }
>     case 2:
>       {
>         dofs[0] = 0;
>       dofs[1] = 1;
>         break;
>       }
>     }
>
>   }
>
>   /// Tabulate the local-to-local mapping of dofs on entity (d, i)
>   virtual void tabulate_entity_dofs(unsigned int* dofs,
>                                     unsigned int d, unsigned int i) const
>   {
>     if (d > 2)
>     {
>     throw std::runtime_error("d is larger than dimension (2)");
>     }
>
>     switch (d)
>     {
>     case 0:
>       {
>         if (i > 2)
>       {
>       throw std::runtime_error("i is larger than number of entities (2)");
>       }
>
>       switch (i)
>       {
>       case 0:
>         {
>           dofs[0] = 0;
>           break;
>         }
>       case 1:
>         {
>           dofs[0] = 1;
>           break;
>         }
>       case 2:
>         {
>           dofs[0] = 2;
>           break;
>         }
>       }
>
>         break;
>       }
>     case 1:
>       {
>
>         break;
>       }
>     case 2:
>       {
>
>         break;
>       }
>     }
>
>   }
>
>   /// Tabulate the coordinates of all dofs on a cell
>   virtual void tabulate_coordinates(double** coordinates,
>                                     const ufc::cell& c) const
>   {
>     const double * const * x = c.coordinates;
>
>     coordinates[0][0] = x[0][0];
>     coordinates[0][1] = x[0][1];
>     coordinates[1][0] = x[1][0];
>     coordinates[1][1] = x[1][1];
>     coordinates[2][0] = x[2][0];
>     coordinates[2][1] = x[2][1];
>   }
>
>   /// Return the number of sub dof maps (for a mixed element)
>   virtual unsigned int num_sub_dof_maps() const
>   {
>     return 0;
>   }
>
>   /// Create a new dof_map for sub dof map i (for a mixed element)
>   virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
>   {
>     return 0;
>   }
>
> };
>
> /// This class defines the interface for the tabulation of the cell
> /// tensor corresponding to the local contribution to a form from
> /// the integral over a cell.
>
> class poisson_cell_integral_0_0: public ufc::cell_integral
> {
>   const std::vector<const pum::GenericPUM*>& pums;
>   mutable std::vector <double> Aa;
>
>   /// Tabulate the regular entities of tensor for the contribution from a local cell
>   virtual void tabulate_regular_tensor(double* A,
>                                        const ufc::cell& c) const
>   {
>     // Extract vertex coordinates
>     const double * const * x = c.coordinates;
>
>     // Compute Jacobian of affine map from reference cell
>     const double J_00 = x[1][0] - x[0][0];
>     const double J_01 = x[2][0] - x[0][0];
>     const double J_10 = x[1][1] - x[0][1];
>     const double J_11 = x[2][1] - x[0][1];
>
>     // Compute determinant of Jacobian
>     double detJ = J_00*J_11 - J_01*J_10;
>
>     // Compute inverse of Jacobian
>     const double K_00 =  J_11 / detJ;
>     const double K_01 = -J_01 / detJ;
>     const double K_10 = -J_10 / detJ;
>     const double K_11 =  J_00 / detJ;
>
>     // Set scale factor
>     const double det = std::abs(detJ);
>
>     // Array of quadrature weights.
>     static const double W1 = 0.500000000000000;
>     // Quadrature points on the UFC reference element: (0.333333333333333, 0.333333333333333)
>
>     // Value of basis functions at quadrature points.
>     static const double FE0_D01[1][6] = \
>     {{-1.000000000000000, 0.000000000000000, 1.000000000000000, -1.000000000000000, 0.000000000000000, 1.000000000000000}};
>
>     static const double FE0_D10[1][6] = \
>     {{-1.000000000000000, 1.000000000000000, 0.000000000000000, -1.000000000000000, 1.000000000000000, 0.000000000000000}};
>
>     static const double FE1[1][3] = \
>     {{0.333333333333333, 0.333333333333333, 0.333333333333333}};
>
>
>     // enriched local dimension of the current cell
>     unsigned int offset = pums[0]->enriched_local_dimension(c);
>
>     // Reset values in the element tensor.
>     const unsigned int num_entries = (3 + offset)*(3 + offset);
>
>     for (unsigned int r = 0; r < num_entries; r++)
>     {
>       A[r] = 0.000000000000000;
>     }// end loop over 'r'
>
>     // Compute element tensor using UFL quadrature representation
>     // Optimisations: ('optimisation', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False), ('ignore zero tables', False)
>
>     // Loop quadrature points for integral.
>     // Number of operations to compute element tensor for following IP loop = 25
>     // Only 1 integration point, omitting IP loop.
>     // Coefficient declarations.
>     double F0 = 0.000000000000000;
>
>     // Total number of operations to compute function values = 6
>     for (unsigned int r = 0; r < 3; r++)
>     {
>       F0 += FE1[0][r]*w[0][r];
>     }// end loop over 'r'
>     unsigned int m = 0;
>
>     // Number of operations for primary indices = 19
>     for (unsigned int j = 0; j < 6; j++)
>     {
>       for (unsigned int k = 0; k < 6; k++)
>       {
>         // Number of operations to compute entry: 19
>         if (((((0 <= j && j < 3)) && ((0 <= k && k < 3)))))
>         {
>           A[m] += ((((K_00*FE0_D10[0][j] + K_10*FE0_D01[0][j]))*((K_00*FE0_D10[0][k] + K_10*FE0_D01[0][k])) + ((K_01*FE0_D10[0][j] + K_11*FE0_D01[0][j]))*((K_01*FE0_D10[0][k] + K_11*FE0_D01[0][k]))))*F0*W1*det;
>         ++m;
>         }
>
>       }// end loop over 'k'
>
>       // Offset the entries corresponding to enriched terms
>       if  ((((0 <= j && j < 3))))
>         m += offset;
>       }// end loop over 'j'
>   }
>
> public:
>
>   /// Constructor
>   poisson_cell_integral_0_0(    const std::vector<const pum::GenericPUM*>& pums) : ufc::cell_integral()    , pums(pums)
>   {
>     // Do nothing
>   }
>
>   /// Destructor
>   virtual ~poisson_cell_integral_0_0()
>   {
>     // Do nothing
>   }
>
>
>   /// Tabulate the tensor for the contribution from a local cell
>   virtual void tabulate_tensor(double* A,
>                                const ufc::cell& c) const
>   {
>     // Tabulate regular entires of element tensor
>     tabulate_regular_tensor(A, w, c);
>
>
>     // enriched local dimension of the current cell(s)
>     unsigned int offset = pums[0]->enriched_local_dimension(c);
>
>     if (offset == 0)
>     {
>       return;
>     }
>
>
>     // Extract vertex coordinates
>     const double * const * x = c.coordinates;
>
>     // Compute Jacobian of affine map from reference cell
>     const double J_00 = x[1][0] - x[0][0];
>     const double J_01 = x[2][0] - x[0][0];
>     const double J_10 = x[1][1] - x[0][1];
>     const double J_11 = x[2][1] - x[0][1];
>
>     // Compute determinant of Jacobian
>     double detJ = J_00*J_11 - J_01*J_10;
>
>
>     // Set scale factor
>     const double det = std::abs(detJ);
>
>     // FIXME: It will crash for multiple discontinuities, if we don't have at least one cell which all dofs are enriched
>     const unsigned int min_entries = 36;
>     const unsigned int num_entries = std::max((offset + 3)*(offset + 3), min_entries);
>
>     // Resizing and reseting auxiliary tensors
>     Aa.resize(num_entries);
>     std::fill(Aa.begin(), Aa.end(), 0.0);
>
>     // Define an array to save current quadrature point
>     double coordinate[2];
>
>     // Define ufc::finite_element object(s) to evalaute shape functions or their derivatives on runtime
>     poisson_finite_element_0  element_0;
>     poisson_finite_element_1  element_1;
>
>     // Array of quadrature weights.
>     static const double W1[1] = {0.500000000000000};
>     // Quadrature points on the UFC reference element: (0.333333333333333, 0.333333333333333)
>
>
>     // Array of quadrature points.
>     static const double P1[2] = \
>     {0.333333333333333, 0.333333333333333};
>
>     // Define vectors for quadrature points and weights(note that the sizes are determined in compile time)
>     std::vector<double> Wn1;
>     std::vector<double> Pn1;
>
>
>     // Check whether there is any need to use modified integration scheme
>     if (pums[0]->modified_quadrature(c))
>     {
>
>       const std::vector<double> weight1(W1, W1 + 1);
>       const std::vector<double> point1(P1, P1 + 2);
>
>       ConstQuadratureRule standard_gauss = std::make_pair(point1, weight1);
>       QuadratureRule modified_gauss;
>
>       pums[0]->cell_quadrature_rule(modified_gauss, standard_gauss, c);
>
>       Pn1 = modified_gauss.first;
>       Wn1 = modified_gauss.second;
>
>     }
>     else
>     {
>       // Map quadrature points from the reference cell to the physical cell
>       Wn1.resize(1);
>       Pn1.resize(2);
>
>
>       for (unsigned int i = 0; i < 1; i++)
>       {
>         Wn1[i] = W1[i];
>         for (unsigned int j = 0; j < 2; j++)
>           Pn1[2*i + j] = x[0][j]*(1.0 - P1[2*i] - P1[2*i + 1]) + x[1][j]*P1[2*i + 1] + x[2][j]*P1[2*i];
>       }
>     }
>
>
>     // Return the values of enriched function at the quadrature points
>     std::vector<double>  enriched_values_1;
>     pums[0]->tabulate_enriched_basis(enriched_values_1, Pn1, c);
>
>     // Define auxilary indices: m, n
>     unsigned int m = 0;
>     unsigned int n = 0;
>
>
>     // Loop quadrature points for integral.
>     for (unsigned int ip = 0; ip < Wn1.size(); ip++)
>     {
>       // Evalaute tables and entries in the element tensor, if the enhanced value at this quadrature point is non-zero
>       if (enriched_values_1[ip] != 0)
>       {
>         // Pick up the coordinates of the current quadrature point
>         coordinate[0] = Pn1[2*ip];
>         coordinate[1] = Pn1[2*ip + 1];
>
>
>         // Creating a table to save the values of derivatives order 1 at the current guass point for <<CG1 on a <triangle of degree 1>> + <<CG1 on a <triangle of degree 1>>>|_{dc0}>
>         double value_0[2];
>         double table_1_D1[6][2];
>         for (unsigned int j = 0; j < 6; j++)
>         {
>           element_1.evaluate_basis_derivatives(j, 1, value_0, coordinate, c);
>           for (unsigned int k = 0; k < 2; k++)
>             table_1_D1[j][k] = value_0[k];
>         }
>
>
>         // Creating a table to save the values of shape functions at the current guass point for <CG1 on a <triangle of degree 1>>
>         double value_1[1];
>         double table_0_D0[3][1];
>         for (unsigned int j = 0; j < 3; j++)
>         {
>           element_0.evaluate_basis(j, value_1, coordinate, c);
>           for (unsigned int k = 0; k < 1; k++)
>             table_0_D0[j][k] = value_1[k];
>         }
>
>
>         // Creating a table to save the values of shape functions at the current guass point for <<CG1 on a <triangle of degree 1>> + <<CG1 on a <triangle of degree 1>>>|_{dc0}>
>         double value_2[1];
>         double table_1_D0[6][1];
>         for (unsigned int j = 0; j < 6; j++)
>         {
>           element_1.evaluate_basis(j, value_2, coordinate, c);
>           for (unsigned int k = 0; k < 1; k++)
>             table_1_D0[j][k] = value_2[k];
>         }
>         // Coefficient declarations.
>         double F0 = 0.000000000000000;
>
>         // Total number of operations to compute function values = 6
>       for (unsigned int r = 0; r < 3; r++)
>       {
>         F0 += table_0_D0[r][0]*w[0][r];
>       }// end loop over 'r'
>
>         // Number of operations for primary indices: 7
>         for (unsigned int j = 0; j < 6; j++)
>         {
>           for (unsigned int k = 0; k < 6; k++)
>           {
>             if (!(((0 <= j && j < 3)) && ((0 <= k && k < 3))))
>             {
>               // Move the indices of discontinuous spaces to the end of mixed space
>               if ((0 <= j && j < 3) && (3 <= k && k < 6))
>               {
>                 m = j;
>                 n = k;
>               }
>               else if ((3 <= j && j < 6) && (0 <= k && k < 3))
>               {
>                 m = j;
>                 n = k;
>               }
>               else if ((3 <= j && j < 6) && (3 <= k && k < 6))
>               {
>                 m = j;
>                 n = k;
>               }
>
>               // Number of operations to compute entry: 7
>               Aa[m*6 + n] += ((table_1_D1[j][1]*table_1_D1[k][1] + table_1_D1[j][0]*table_1_D1[k][0]))*F0*Wn1[ip]*det;
>             }// end check for enriched entiries
>           }// end loop over 'k'
>         }// end loop over 'j'
>       }
>     }// end loop over 'ip'
>
>
>     // Pick up entries from the total element tensor for the nodes active in the enrichment
>
>     // Determine a vector that contains the local numbering of enriched degrees of freedom in ufc::cell c for the field 0
>     std::vector<unsigned int> active_dofs_0;
>     pums[0]->tabulate_enriched_local_dofs(active_dofs_0, c);
>     std::vector<unsigned int>::iterator   it_0_0, it_0_1;
>
>
>     m = 0;
>     for (unsigned int j = 0; j < 6; j++)
>       for (unsigned int k = 0; k < 6; k++)
>         if ((0 <= j && j < 3) && (0 <= k && k < 3))
>           ++m;
>         else
>         {
>           it_0_0 = find(active_dofs_0.begin(), active_dofs_0.end(), j - 3);
>           it_0_1 = find(active_dofs_0.begin(), active_dofs_0.end(), k - 3);
>
>
>           // Check whether the entry is coressponding to the active enriched node
>           if  (it_0_0 != active_dofs_0.end() || it_0_1 != active_dofs_0.end())
>             if (((0 <= j && j < 3)) || ((0 <= k && k < 3)) || (it_0_0 != active_dofs_0.end() && it_0_1 != active_dofs_0.end()))
>             {
>               A[m] = Aa[j*6 + k];
>               ++m;
>             }
>         }
>   }
>
> };
>
> /// This class defines the interface for the tabulation of the cell
> /// tensor corresponding to the local contribution to a form from
> /// the integral over a cell.
>
> class poisson_cell_integral_1_0: public ufc::cell_integral
> {
>   const std::vector<const pum::GenericPUM*>& pums;
>   mutable std::vector <double> Aa;
>
>   /// Tabulate the regular entities of tensor for the contribution from a local cell
>   virtual void tabulate_regular_tensor(double* A,
>                                        const ufc::cell& c) const
>   {
>     // Extract vertex coordinates
>     const double * const * x = c.coordinates;
>
>     // Compute Jacobian of affine map from reference cell
>     const double J_00 = x[1][0] - x[0][0];
>     const double J_01 = x[2][0] - x[0][0];
>     const double J_10 = x[1][1] - x[0][1];
>     const double J_11 = x[2][1] - x[0][1];
>
>     // Compute determinant of Jacobian
>     double detJ = J_00*J_11 - J_01*J_10;
>
>     // Compute inverse of Jacobian
>
>     // Set scale factor
>     const double det = std::abs(detJ);
>
>     // Array of quadrature weights.
>     static const double W4[4] = {0.159020690871988, 0.090979309128011, 0.159020690871988, 0.090979309128011};
>     // Quadrature points on the UFC reference element: (0.178558728263616, 0.155051025721682), (0.075031110222608, 0.644948974278318), (0.666390246014701, 0.155051025721682), (0.280019915499074, 0.644948974278318)
>
>     // Value of basis functions at quadrature points.
>     static const double FE0[4][6] = \
>     {{0.666390246014701, 0.178558728263616, 0.155051025721682, 0.666390246014701, 0.178558728263616, 0.155051025721682},
>     {0.280019915499074, 0.075031110222608, 0.644948974278318, 0.280019915499074, 0.075031110222608, 0.644948974278318},
>     {0.178558728263616, 0.666390246014701, 0.155051025721682, 0.178558728263616, 0.666390246014701, 0.155051025721682},
>     {0.075031110222608, 0.280019915499074, 0.644948974278318, 0.075031110222608, 0.280019915499074, 0.644948974278318}};
>
>     static const double FE1[4][3] = \
>     {{0.666390246014701, 0.178558728263616, 0.155051025721682},
>     {0.280019915499074, 0.075031110222608, 0.644948974278318},
>     {0.178558728263616, 0.666390246014701, 0.155051025721682},
>     {0.075031110222608, 0.280019915499074, 0.644948974278318}};
>
>
>     // enriched local dimension of the current cell
>     unsigned int offset = pums[0]->enriched_local_dimension(c);
>
>     // Reset values in the element tensor.
>     const unsigned int num_entries = (3 + offset);
>
>     for (unsigned int r = 0; r < num_entries; r++)
>     {
>       A[r] = 0.000000000000000;
>     }// end loop over 'r'
>
>     // Compute element tensor using UFL quadrature representation
>     // Optimisations: ('optimisation', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False), ('ignore zero tables', False)
>
>     // Loop quadrature points for integral.
>     // Number of operations to compute element tensor for following IP loop = 40
>     for (unsigned int ip = 0; ip < 4; ip++)
>     {
>       // Coefficient declarations.
>       double F0 = 0.000000000000000;
>
>       // Total number of operations to compute function values = 6
>       for (unsigned int r = 0; r < 3; r++)
>       {
>         F0 += FE1[ip][r]*w[0][r];
>       }// end loop over 'r'
>       unsigned int m = 0;
>
>       // Number of operations for primary indices = 4
>       for (unsigned int j = 0; j < 6; j++)
>       {
>         // Number of operations to compute entry: 4
>         if (((((0 <= j && j < 3)))))
>         {
>           A[m] += FE0[ip][j]*F0*W4[ip]*det;
>         ++m;
>         }
>
>       }// end loop over 'j'
>     }// end loop over 'ip'
>   }
>
> public:
>
>   /// Constructor
>   poisson_cell_integral_1_0(    const std::vector<const pum::GenericPUM*>& pums) : ufc::cell_integral()    , pums(pums)
>   {
>     // Do nothing
>   }
>
>   /// Destructor
>   virtual ~poisson_cell_integral_1_0()
>   {
>     // Do nothing
>   }
>
>
>   /// Tabulate the tensor for the contribution from a local cell
>   virtual void tabulate_tensor(double* A,
>                                const ufc::cell& c) const
>   {
>     // Tabulate regular entires of element tensor
>     tabulate_regular_tensor(A, w, c);
>
>
>     // enriched local dimension of the current cell(s)
>     unsigned int offset = pums[0]->enriched_local_dimension(c);
>
>     if (offset == 0)
>     {
>       return;
>     }
>
>
>     // Extract vertex coordinates
>     const double * const * x = c.coordinates;
>
>     // Compute Jacobian of affine map from reference cell
>     const double J_00 = x[1][0] - x[0][0];
>     const double J_01 = x[2][0] - x[0][0];
>     const double J_10 = x[1][1] - x[0][1];
>     const double J_11 = x[2][1] - x[0][1];
>
>     // Compute determinant of Jacobian
>     double detJ = J_00*J_11 - J_01*J_10;
>
>
>     // Set scale factor
>     const double det = std::abs(detJ);
>
>     // FIXME: It will crash for multiple discontinuities, if we don't have at least one cell which all dofs are enriched
>     const unsigned int min_entries = 6;
>     const unsigned int num_entries = std::max((offset + 3), min_entries);
>
>     // Resizing and reseting auxiliary tensors
>     Aa.resize(num_entries);
>     std::fill(Aa.begin(), Aa.end(), 0.0);
>
>     // Define an array to save current quadrature point
>     double coordinate[2];
>
>     // Define ufc::finite_element object(s) to evalaute shape functions or their derivatives on runtime
>     poisson_finite_element_0  element_0;
>     poisson_finite_element_1  element_1;
>
>     // Array of quadrature weights.
>     static const double W4[4] = {0.159020690871988, 0.090979309128011, 0.159020690871988, 0.090979309128011};
>     // Quadrature points on the UFC reference element: (0.178558728263616, 0.155051025721682), (0.075031110222608, 0.644948974278318), (0.666390246014701, 0.155051025721682), (0.280019915499074, 0.644948974278318)
>
>
>     // Array of quadrature points.
>     static const double P4[8] = \
>     {0.178558728263616, 0.155051025721682,
>     0.075031110222608, 0.644948974278318,
>     0.666390246014701, 0.155051025721682,
>     0.280019915499074, 0.644948974278318};
>
>     // Define vectors for quadrature points and weights(note that the sizes are determined in compile time)
>     std::vector<double> Wn4;
>     std::vector<double> Pn4;
>
>
>     // Check whether there is any need to use modified integration scheme
>     if (pums[0]->modified_quadrature(c))
>     {
>
>       const std::vector<double> weight4(W4, W4 + 4);
>       const std::vector<double> point4(P4, P4 + 8);
>
>       ConstQuadratureRule standard_gauss = std::make_pair(point4, weight4);
>       QuadratureRule modified_gauss;
>
>       pums[0]->cell_quadrature_rule(modified_gauss, standard_gauss, c);
>
>       Pn4 = modified_gauss.first;
>       Wn4 = modified_gauss.second;
>
>     }
>     else
>     {
>       // Map quadrature points from the reference cell to the physical cell
>       Wn4.resize(4);
>       Pn4.resize(8);
>
>
>       for (unsigned int i = 0; i < 4; i++)
>       {
>         Wn4[i] = W4[i];
>         for (unsigned int j = 0; j < 2; j++)
>           Pn4[2*i + j] = x[0][j]*(1.0 - P4[2*i] - P4[2*i + 1]) + x[1][j]*P4[2*i + 1] + x[2][j]*P4[2*i];
>       }
>     }
>
>
>     // Return the values of enriched function at the quadrature points
>     std::vector<double>  enriched_values_4;
>     pums[0]->tabulate_enriched_basis(enriched_values_4, Pn4, c);
>
>     // Define an auxilary index: m
>     unsigned int m = 0;
>
>
>     // Loop quadrature points for integral.
>     for (unsigned int ip = 0; ip < Wn4.size(); ip++)
>     {
>       // Evalaute tables and entries in the element tensor, if the enhanced value at this quadrature point is non-zero
>       if (enriched_values_4[ip] != 0)
>       {
>         // Pick up the coordinates of the current quadrature point
>         coordinate[0] = Pn4[2*ip];
>         coordinate[1] = Pn4[2*ip + 1];
>
>
>         // Creating a table to save the values of shape functions at the current guass point for <CG1 on a <triangle of degree 1>>
>         double value_0[1];
>         double table_0_D0[3][1];
>         for (unsigned int j = 0; j < 3; j++)
>         {
>           element_0.evaluate_basis(j, value_0, coordinate, c);
>           for (unsigned int k = 0; k < 1; k++)
>             table_0_D0[j][k] = value_0[k];
>         }
>
>
>         // Creating a table to save the values of shape functions at the current guass point for <<CG1 on a <triangle of degree 1>> + <<CG1 on a <triangle of degree 1>>>|_{dc0}>
>         double value_1[1];
>         double table_1_D0[6][1];
>         for (unsigned int j = 0; j < 6; j++)
>         {
>           element_1.evaluate_basis(j, value_1, coordinate, c);
>           for (unsigned int k = 0; k < 1; k++)
>             table_1_D0[j][k] = value_1[k];
>         }
>         // Coefficient declarations.
>         double F0 = 0.000000000000000;
>
>         // Total number of operations to compute function values = 6
>       for (unsigned int r = 0; r < 3; r++)
>       {
>         F0 += table_0_D0[r][0]*w[0][r];
>       }// end loop over 'r'
>
>         // Number of operations for primary indices: 4
>         for (unsigned int j = 0; j < 6; j++)
>         {
>           if (!(((0 <= j && j < 3))))
>           {
>             // Move the indices of discontinuous spaces to the end of mixed space
>             if ((3 <= j && j < 6))
>             {
>               m = j;
>             }
>
>             // Number of operations to compute entry: 4
>             Aa[m] += table_1_D0[j][0]*F0*Wn4[ip]*det;
>           }// end check for enriched entiries
>         }// end loop over 'j'
>       }
>     }// end loop over 'ip'
>
>
>     // Pick up entries from the total element tensor for the nodes active in the enrichment
>
>     // Determine a vector that contains the local numbering of enriched degrees of freedom in ufc::cell c for the field 0
>     std::vector<unsigned int> active_dofs_0;
>     pums[0]->tabulate_enriched_local_dofs(active_dofs_0, c);
>     std::vector<unsigned int>::iterator   it_0_0;
>
>
>     m = 0;
>     for (unsigned int j = 0; j < 6; j++)
>       if ((0 <= j && j < 3))
>         ++m;
>       else
>       {
>         it_0_0 = find(active_dofs_0.begin(), active_dofs_0.end(), j - 3);
>
>
>         // Check whether the entry is coressponding to the active enriched node
>         if  (it_0_0 != active_dofs_0.end())
>         {
>           A[m] = Aa[j];
>           ++m;
>         }
>       }
>   }
>
> };
>
> /// This class defines the interface for the tabulation of the
> /// exterior facet tensor corresponding to the local contribution to
> /// a form from the integral over an exterior facet.
>
> class poisson_exterior_facet_integral_1_0: public ufc::exterior_facet_integral
> {
>   const std::vector<const pum::GenericPUM*>& pums;
>   mutable std::vector <double> Aa;
>
>   /// Tabulate the regular entities of the tensor for the contribution from a local exterior facet
>   virtual void tabulate_regular_tensor(double* A,
>                                        const ufc::cell& c,
>                                        unsigned int facet) const
>   {
>     // Extract vertex coordinates
>     const double * const * x = c.coordinates;
>
>     // Compute Jacobian of affine map from reference cell
>
>     // Compute determinant of Jacobian
>
>     // Compute inverse of Jacobian
>
>     // Get vertices on edge
>     static unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};
>     const unsigned int v0 = edge_vertices[facet][0];
>     const unsigned int v1 = edge_vertices[facet][1];
>
>     // Compute scale factor (length of edge scaled by length of reference interval)
>     const double dx0 = x[v1][0] - x[v0][0];
>     const double dx1 = x[v1][1] - x[v0][1];
>     const double det = std::sqrt(dx0*dx0 + dx1*dx1);
>
>
>     // Array of quadrature weights.
>     static const double W2[2] = {0.500000000000000, 0.500000000000000};
>     // Quadrature points on the UFC reference element: (0.211324865405187), (0.788675134594813)
>
>     // Value of basis functions at quadrature points.
>     static const double FE0_f0[2][6] = \
>     {{0.000000000000000, 0.788675134594813, 0.211324865405187, 0.000000000000000, 0.788675134594813, 0.211324865405187},
>     {0.000000000000000, 0.211324865405187, 0.788675134594813, 0.000000000000000, 0.211324865405187, 0.788675134594813}};
>
>     static const double FE0_f1[2][6] = \
>     {{0.788675134594813, 0.000000000000000, 0.211324865405187, 0.788675134594813, 0.000000000000000, 0.211324865405187},
>     {0.211324865405187, 0.000000000000000, 0.788675134594813, 0.211324865405187, 0.000000000000000, 0.788675134594813}};
>
>     static const double FE0_f2[2][6] = \
>     {{0.788675134594813, 0.211324865405187, 0.000000000000000, 0.788675134594813, 0.211324865405187, 0.000000000000000},
>     {0.211324865405187, 0.788675134594813, 0.000000000000000, 0.211324865405187, 0.788675134594813, 0.000000000000000}};
>
>     static const double FE1_f0[2][3] = \
>     {{0.000000000000000, 0.788675134594813, 0.211324865405187},
>     {0.000000000000000, 0.211324865405187, 0.788675134594813}};
>
>     static const double FE1_f1[2][3] = \
>     {{0.788675134594813, 0.000000000000000, 0.211324865405187},
>     {0.211324865405187, 0.000000000000000, 0.788675134594813}};
>
>     static const double FE1_f2[2][3] = \
>     {{0.788675134594813, 0.211324865405187, 0.000000000000000},
>     {0.211324865405187, 0.788675134594813, 0.000000000000000}};
>
>
>     // enriched local dimension of the current cell
>     unsigned int offset = pums[0]->enriched_local_dimension(c);
>
>     // Reset values in the element tensor.
>     const unsigned int num_entries = (3 + offset);
>
>     for (unsigned int r = 0; r < num_entries; r++)
>     {
>       A[r] = 0.000000000000000;
>     }// end loop over 'r'
>
>     // Compute element tensor using UFL quadrature representation
>     // Optimisations: ('optimisation', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False), ('ignore zero tables', False)
>     switch (facet)
>     {
>     case 0:
>       {
>         // Total number of operations to compute element tensor (from this point): 22
>
>       // Loop quadrature points for integral.
>       // Number of operations to compute element tensor for following IP loop = 22
>       for (unsigned int ip = 0; ip < 2; ip++)
>       {
>         // Coefficient declarations.
>         double F0 = 0.000000000000000;
>
>         // Total number of operations to compute function values = 6
>         for (unsigned int r = 0; r < 3; r++)
>         {
>           F0 += FE1_f0[ip][r]*w[1][r];
>         }// end loop over 'r'
>         unsigned int m = 0;
>
>         // Number of operations for primary indices = 5
>         for (unsigned int j = 0; j < 6; j++)
>         {
>           // Number of operations to compute entry: 5
>           if (((((0 <= j && j < 3)))))
>           {
>             A[m] += FE0_f0[ip][j]*F0*(-1.000000000000000)*W2[ip]*det;
>           ++m;
>           }
>
>         }// end loop over 'j'
>       }// end loop over 'ip'
>         break;
>       }
>     case 1:
>       {
>         // Total number of operations to compute element tensor (from this point): 22
>
>       // Loop quadrature points for integral.
>       // Number of operations to compute element tensor for following IP loop = 22
>       for (unsigned int ip = 0; ip < 2; ip++)
>       {
>         // Coefficient declarations.
>         double F0 = 0.000000000000000;
>
>         // Total number of operations to compute function values = 6
>         for (unsigned int r = 0; r < 3; r++)
>         {
>           F0 += FE1_f1[ip][r]*w[1][r];
>         }// end loop over 'r'
>         unsigned int m = 0;
>
>         // Number of operations for primary indices = 5
>         for (unsigned int j = 0; j < 6; j++)
>         {
>           // Number of operations to compute entry: 5
>           if (((((0 <= j && j < 3)))))
>           {
>             A[m] += FE0_f1[ip][j]*F0*(-1.000000000000000)*W2[ip]*det;
>           ++m;
>           }
>
>         }// end loop over 'j'
>       }// end loop over 'ip'
>         break;
>       }
>     case 2:
>       {
>         // Total number of operations to compute element tensor (from this point): 22
>
>       // Loop quadrature points for integral.
>       // Number of operations to compute element tensor for following IP loop = 22
>       for (unsigned int ip = 0; ip < 2; ip++)
>       {
>         // Coefficient declarations.
>         double F0 = 0.000000000000000;
>
>         // Total number of operations to compute function values = 6
>         for (unsigned int r = 0; r < 3; r++)
>         {
>           F0 += FE1_f2[ip][r]*w[1][r];
>         }// end loop over 'r'
>         unsigned int m = 0;
>
>         // Number of operations for primary indices = 5
>         for (unsigned int j = 0; j < 6; j++)
>         {
>           // Number of operations to compute entry: 5
>           if (((((0 <= j && j < 3)))))
>           {
>             A[m] += FE0_f2[ip][j]*F0*(-1.000000000000000)*W2[ip]*det;
>           ++m;
>           }
>
>         }// end loop over 'j'
>       }// end loop over 'ip'
>         break;
>       }
>     }
>
>   }
>
> public:
>
>   /// Constructor
>   poisson_exterior_facet_integral_1_0(    const std::vector<const pum::GenericPUM*>& pums) : ufc::exterior_facet_integral()    , pums(pums)
>   {
>     // Do nothing
>   }
>
>   /// Destructor
>   virtual ~poisson_exterior_facet_integral_1_0()
>   {
>     // Do nothing
>   }
>
>   /// Tabulate the tensor for the contribution from a local exterior facet
>   virtual void tabulate_tensor(double* A,
>                                const ufc::cell& c,
>                                unsigned int facet) const
>   {
>     // Tabulate regular entires of element tensor
>     tabulate_regular_tensor(A, w, c, facet);
>
>
>     // enriched local dimension of the current cell(s)
>     unsigned int offset = pums[0]->enriched_local_dimension(c);
>
>     if (offset == 0)
>     {
>       return;
>     }
>
>
>     // Extract vertex coordinates
>     const double * const * x = c.coordinates;
>
>     // Get vertices on edge
>     static unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};
>     const unsigned int v0 = edge_vertices[facet][0];
>     const unsigned int v1 = edge_vertices[facet][1];
>
>     // Compute scale factor (length of edge scaled by length of reference interval)
>     const double dx0 = x[v1][0] - x[v0][0];
>     const double dx1 = x[v1][1] - x[v0][1];
>     const double det = std::sqrt(dx0*dx0 + dx1*dx1);
>
>     // FIXME: It will crash for multiple discontinuities, if we don't have at least one cell which all dofs are enriched
>     const unsigned int min_entries = 6;
>     const unsigned int num_entries = std::max((offset + 3), min_entries);
>
>     // Resizing and reseting auxiliary tensors
>     Aa.resize(num_entries);
>     std::fill(Aa.begin(), Aa.end(), 0.0);
>
>     // Define an array to save current quadrature point
>     double coordinate[2];
>
>     // Define ufc::finite_element object(s) to evalaute shape functions or their derivatives on runtime
>     poisson_finite_element_0  element_0;
>     poisson_finite_element_1  element_1;
>
>     // Array of quadrature weights.
>     static const double W2[2] = {0.500000000000000, 0.500000000000000};
>     // Quadrature points on the UFC reference element: (0.211324865405187), (0.788675134594813)
>
>
>     // Array of quadrature points.
>     static const double P2[2] = \
>     {0.211324865405187,
>     0.788675134594813};
>
>     // Define vectors for quadrature points and weights(note that the sizes are determined in compile time)
>     std::vector<double> Wn2;
>     std::vector<double> Pn2;
>
>
>     // Check whether there is any need to use modified integration scheme
>     if (pums[0]->modified_quadrature(c, facet))
>     {
>
>       const std::vector<double> weight2(W2, W2 + 2);
>       const std::vector<double> point2(P2, P2 + 2);
>
>       ConstQuadratureRule standard_gauss = std::make_pair(point2, weight2);
>       QuadratureRule modified_gauss;
>
>       pums[0]->facet_quadrature_rule(modified_gauss, standard_gauss, c, facet);
>
>       Pn2 = modified_gauss.first;
>       Wn2 = modified_gauss.second;
>
>     }
>     else
>     {
>       // Map quadrature points from the reference cell to the physical cell
>       Wn2.resize(2);
>       Pn2.resize(4);
>
>
>       for (unsigned int i = 0; i < 2; i++)
>       {
>         Wn2[i] = W2[i];
>         for (unsigned int j = 0; j < 2; j++)
>           Pn2[2*i + j] = x[v0][j]*(1.0 - P2[i]) + x[v1][j]*P2[i];
>       }
>     }
>
>
>     // Return the values of enriched function at the quadrature points
>     std::vector<double>  enriched_values_2;
>     pums[0]->tabulate_enriched_basis(enriched_values_2, Pn2, c);
>
>     // Define an auxilary index: m
>     unsigned int m = 0;
>
>
>     // Loop quadrature points for integral.
>     for (unsigned int ip = 0; ip < Wn2.size(); ip++)
>     {
>       // Evalaute tables and entries in the element tensor, if the enhanced value at this quadrature point is non-zero
>       if (enriched_values_2[ip] != 0)
>       {
>         // Pick up the coordinates of the current quadrature point
>         coordinate[0] = Pn2[2*ip];
>         coordinate[1] = Pn2[2*ip + 1];
>
>
>         // Creating a table to save the values of shape functions at the current guass point for <CG1 on a <triangle of degree 1>>
>         double value_0[1];
>         double table_0_D0[3][1];
>         for (unsigned int j = 0; j < 3; j++)
>         {
>           element_0.evaluate_basis(j, value_0, coordinate, c);
>           for (unsigned int k = 0; k < 1; k++)
>             table_0_D0[j][k] = value_0[k];
>         }
>
>
>         // Creating a table to save the values of shape functions at the current guass point for <<CG1 on a <triangle of degree 1>> + <<CG1 on a <triangle of degree 1>>>|_{dc0}>
>         double value_1[1];
>         double table_1_D0[6][1];
>         for (unsigned int j = 0; j < 6; j++)
>         {
>           element_1.evaluate_basis(j, value_1, coordinate, c);
>           for (unsigned int k = 0; k < 1; k++)
>             table_1_D0[j][k] = value_1[k];
>         }
>         // Coefficient declarations.
>         double F0 = 0.000000000000000;
>
>         // Total number of operations to compute function values = 6
>       for (unsigned int r = 0; r < 3; r++)
>       {
>         F0 += table_0_D0[r][0]*w[1][r];
>       }// end loop over 'r'
>
>         // Number of operations for primary indices: 5
>         for (unsigned int j = 0; j < 6; j++)
>         {
>           if (!(((0 <= j && j < 3))))
>           {
>             // Move the indices of discontinuous spaces to the end of mixed space
>             if ((3 <= j && j < 6))
>             {
>               m = j;
>             }
>
>             // Number of operations to compute entry: 5
>             Aa[m] += table_1_D0[j][0]*F0*(-1.000000000000000)*Wn2[ip]*det;
>           }// end check for enriched entiries
>         }// end loop over 'j'
>       }
>     }// end loop over 'ip'
>
>
>     // Pick up entries from the total element tensor for the nodes active in the enrichment
>
>     // Determine a vector that contains the local numbering of enriched degrees of freedom in ufc::cell c for the field 0
>     std::vector<unsigned int> active_dofs_0;
>     pums[0]->tabulate_enriched_local_dofs(active_dofs_0, c);
>     std::vector<unsigned int>::iterator   it_0_0;
>
>
>     m = 0;
>     for (unsigned int j = 0; j < 6; j++)
>       if ((0 <= j && j < 3))
>         ++m;
>       else
>       {
>         it_0_0 = find(active_dofs_0.begin(), active_dofs_0.end(), j - 3);
>
>
>         // Check whether the entry is coressponding to the active enriched node
>         if  (it_0_0 != active_dofs_0.end())
>         {
>           A[m] = Aa[j];
>           ++m;
>         }
>       }
>   }
>
> };
>
> /// This class defines the interface for the assembly of the global
> /// tensor corresponding to a form with r + n arguments, that is, a
> /// mapping
> ///
> ///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
> ///
> /// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
> /// global tensor A is defined by
> ///
> ///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
> ///
> /// where each argument Vj represents the application to the
> /// sequence of basis functions of Vj and w1, w2, ..., wn are given
> /// fixed functions (coefficients).
>
> class poisson_form_0: public ufc::form
> {
>   const std::vector<const pum::GenericPUM*>& pums;
> public:
>
>   /// Constructor
>   poisson_form_0(    const std::vector<const pum::GenericPUM*>& pums) : ufc::form()    , pums(pums)
>   {
>     // Do nothing
>   }
>
>   /// Destructor
>   virtual ~poisson_form_0()
>   {
>     // Do nothing
>   }
>
>   /// Return a string identifying the form
>   virtual const char* signature() const
>   {
>     return "Form([Integral(Product(Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), IndexSum(Product(Indexed(ComponentTensor(SpatialDerivative(Argument(EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None))), 0), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ComponentTensor(SpatialDerivative(Argument(EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None))), 1), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2}))), Measure('cell', 0, None))])";
>   }
>
>   /// Return the rank of the global tensor (r)
>   virtual unsigned int rank() const
>   {
>     return 2;
>   }
>
>   /// Return the number of coefficients (n)
>   virtual unsigned int num_coefficients() const
>   {
>     return 1;
>   }
>
>   /// Return the number of cell integrals
>   virtual unsigned int num_cell_integrals() const
>   {
>     return 1;
>   }
>
>   /// Return the number of exterior facet integrals
>   virtual unsigned int num_exterior_facet_integrals() const
>   {
>     return 0;
>   }
>
>   /// Return the number of interior facet integrals
>   virtual unsigned int num_interior_facet_integrals() const
>   {
>     return 0;
>   }
>
>   /// Create a new finite element for argument function i
>   virtual ufc::finite_element* create_finite_element(unsigned int i) const
>   {
>     switch (i)
>     {
>     case 0:
>       {
>         return new poisson_finite_element_1();
>         break;
>       }
>     case 1:
>       {
>         return new poisson_finite_element_1();
>         break;
>       }
>     case 2:
>       {
>         return new poisson_finite_element_0();
>         break;
>       }
>     }
>
>     return 0;
>   }
>
>   /// Create a new dof map for argument function i
>   virtual ufc::dof_map* create_dof_map(unsigned int i) const
>   {
>     switch (i)
>     {
>     case 0:
>       {
>         return new poisson_dof_map_1(pums);
>         break;
>       }
>     case 1:
>       {
>         return new poisson_dof_map_1(pums);
>         break;
>       }
>     case 2:
>       {
>         return new poisson_dof_map_0();
>         break;
>       }
>     }
>
>     return 0;
>   }
>
>   /// Create a new cell integral on sub domain i
>   virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
>   {
>     switch (i)
>     {
>     case 0:
>       {
>         return new poisson_cell_integral_0_0(pums);
>         break;
>       }
>     }
>
>     return 0;
>   }
>
>   /// Create a new exterior facet integral on sub domain i
>   virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
>   {
>     return 0;
>   }
>
>   /// Create a new interior facet integral on sub domain i
>   virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
>   {
>     return 0;
>   }
>
> };
>
> /// This class defines the interface for the assembly of the global
> /// tensor corresponding to a form with r + n arguments, that is, a
> /// mapping
> ///
> ///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
> ///
> /// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
> /// global tensor A is defined by
> ///
> ///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
> ///
> /// where each argument Vj represents the application to the
> /// sequence of basis functions of Vj and w1, w2, ..., wn are given
> /// fixed functions (coefficients).
>
> class poisson_form_1: public ufc::form
> {
>   const std::vector<const pum::GenericPUM*>& pums;
> public:
>
>   /// Constructor
>   poisson_form_1(    const std::vector<const pum::GenericPUM*>& pums) : ufc::form()    , pums(pums)
>   {
>     // Do nothing
>   }
>
>   /// Destructor
>   virtual ~poisson_form_1()
>   {
>     // Do nothing
>   }
>
>   /// Return a string identifying the form
>   virtual const char* signature() const
>   {
>     return "Form([Integral(Product(Argument(EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None))), 0), Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0)), Measure('cell', 0, None)), Integral(Product(IntValue(-1, (), (), {}), Product(Argument(EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None))), 0), Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 1))), Measure('exterior_facet', 0, None))])";
>   }
>
>   /// Return the rank of the global tensor (r)
>   virtual unsigned int rank() const
>   {
>     return 1;
>   }
>
>   /// Return the number of coefficients (n)
>   virtual unsigned int num_coefficients() const
>   {
>     return 2;
>   }
>
>   /// Return the number of cell integrals
>   virtual unsigned int num_cell_integrals() const
>   {
>     return 1;
>   }
>
>   /// Return the number of exterior facet integrals
>   virtual unsigned int num_exterior_facet_integrals() const
>   {
>     return 1;
>   }
>
>   /// Return the number of interior facet integrals
>   virtual unsigned int num_interior_facet_integrals() const
>   {
>     return 0;
>   }
>
>   /// Create a new finite element for argument function i
>   virtual ufc::finite_element* create_finite_element(unsigned int i) const
>   {
>     switch (i)
>     {
>     case 0:
>       {
>         return new poisson_finite_element_1();
>         break;
>       }
>     case 1:
>       {
>         return new poisson_finite_element_0();
>         break;
>       }
>     case 2:
>       {
>         return new poisson_finite_element_0();
>         break;
>       }
>     }
>
>     return 0;
>   }
>
>   /// Create a new dof map for argument function i
>   virtual ufc::dof_map* create_dof_map(unsigned int i) const
>   {
>     switch (i)
>     {
>     case 0:
>       {
>         return new poisson_dof_map_1(pums);
>         break;
>       }
>     case 1:
>       {
>         return new poisson_dof_map_0();
>         break;
>       }
>     case 2:
>       {
>         return new poisson_dof_map_0();
>         break;
>       }
>     }
>
>     return 0;
>   }
>
>   /// Create a new cell integral on sub domain i
>   virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
>   {
>     switch (i)
>     {
>     case 0:
>       {
>         return new poisson_cell_integral_1_0(pums);
>         break;
>       }
>     }
>
>     return 0;
>   }
>
>   /// Create a new exterior facet integral on sub domain i
>   virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
>   {
>     switch (i)
>     {
>     case 0:
>       {
>         return new poisson_exterior_facet_integral_1_0(pums);
>         break;
>       }
>     }
>
>     return 0;
>   }
>
>   /// Create a new interior facet integral on sub domain i
>   virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
>   {
>     return 0;
>   }
>
> };
>
> /// This class defines the interface for a finite element.
>
> class poisson_auxiliary_finite_element_0: public ufc::finite_element
> {
> public:
>
>   /// Constructor
>   poisson_auxiliary_finite_element_0() : ufc::finite_element()
>   {
>     // Do nothing
>   }
>
>   /// Destructor
>   virtual ~poisson_auxiliary_finite_element_0()
>   {
>     // Do nothing
>   }
>
>   /// Return a string identifying the finite element
>   virtual const char* signature() const
>   {
>     return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
>   }
>
>   /// Return the cell shape
>   virtual ufc::shape cell_shape() const
>   {
>     return ufc::triangle;
>   }
>
>   /// Return the dimension of the finite element function space
>   virtual unsigned int space_dimension() const
>   {
>     return 3;
>   }
>
>   /// Return the rank of the value space
>   virtual unsigned int value_rank() const
>   {
>     return 0;
>   }
>
>   /// Return the dimension of the value space for axis i
>   virtual unsigned int value_dimension(unsigned int i) const
>   {
>     return 1;
>   }
>
>   /// Evaluate basis function i at given point in cell
>   virtual void evaluate_basis(unsigned int i,
>                               double* values,
>                               const double* coordinates,
>                               const ufc::cell& c) const
>   {
>     // Extract vertex coordinates
>     const double * const * x = c.coordinates;
>
>     // Compute Jacobian of affine map from reference cell
>     const double J_00 = x[1][0] - x[0][0];
>     const double J_01 = x[2][0] - x[0][0];
>     const double J_10 = x[1][1] - x[0][1];
>     const double J_11 = x[2][1] - x[0][1];
>
>     // Compute determinant of Jacobian
>     double detJ = J_00*J_11 - J_01*J_10;
>
>     // Compute inverse of Jacobian
>
>     // Compute constants
>     const double C0 = x[1][0] + x[2][0];
>     const double C1 = x[1][1] + x[2][1];
>
>     // Get coordinates and map to the reference (FIAT) element
>     double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
>     double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
>
>     // Reset values.
>     *values = 0.000000000000000;
>
>     // Map degree of freedom to element degree of freedom
>     const unsigned int dof = i;
>
>     // Array of basisvalues.
>     double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
>
>     // Declare helper variables.
>     unsigned int rr = 0;
>     unsigned int ss = 0;
>     double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
>
>     // Compute basisvalues.
>     basisvalues[0] = 1.000000000000000;
>     basisvalues[1] = tmp0;
>     for (unsigned int r = 0; r < 1; r++)
>     {
>       rr = (r + 1)*(r + 1 + 1)/2 + 1;
>       ss = r*(r + 1)/2;
>       basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
>     }// end loop over 'r'
>     for (unsigned int r = 0; r < 2; r++)
>     {
>       for (unsigned int s = 0; s < 2 - r; s++)
>       {
>         rr = (r + s)*(r + s + 1)/2 + s;
>         basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
>       }// end loop over 's'
>     }// end loop over 'r'
>
>     // Table(s) of coefficients.
>     static const double coefficients0[3][3] = \
>     {{0.471404520791032, -0.288675134594813, -0.166666666666667},
>     {0.471404520791032, 0.288675134594813, -0.166666666666667},
>     {0.471404520791032, 0.000000000000000, 0.333333333333333}};
>
>     // Compute value(s).
>     for (unsigned int r = 0; r < 3; r++)
>     {
>       *values += coefficients0[dof][r]*basisvalues[r];
>     }// end loop over 'r'
>   }
>
>   /// Evaluate all basis functions at given point in cell
>   virtual void evaluate_basis_all(double* values,
>                                   const double* coordinates,
>                                   const ufc::cell& c) const
>   {
>     // Helper variable to hold values of a single dof.
>     double dof_values = 0.000000000000000;
>
>     // Loop dofs and call evaluate_basis.
>     for (unsigned int r = 0; r < 3; r++)
>     {
>       evaluate_basis(r, &dof_values, coordinates, c);
>       values[r] = dof_values;
>     }// end loop over 'r'
>   }
>
>   /// Evaluate order n derivatives of basis function i at given point in cell
>   virtual void evaluate_basis_derivatives(unsigned int i,
>                                           unsigned int n,
>                                           double* values,
>                                           const double* coordinates,
>                                           const ufc::cell& c) const
>   {
>     // Extract vertex coordinates
>     const double * const * x = c.coordinates;
>
>     // Compute Jacobian of affine map from reference cell
>     const double J_00 = x[1][0] - x[0][0];
>     const double J_01 = x[2][0] - x[0][0];
>     const double J_10 = x[1][1] - x[0][1];
>     const double J_11 = x[2][1] - x[0][1];
>
>     // Compute determinant of Jacobian
>     double detJ = J_00*J_11 - J_01*J_10;
>
>     // Compute inverse of Jacobian
>     const double K_00 =  J_11 / detJ;
>     const double K_01 = -J_01 / detJ;
>     const double K_10 = -J_10 / detJ;
>     const double K_11 =  J_00 / detJ;
>
>     // Compute constants
>     const double C0 = x[1][0] + x[2][0];
>     const double C1 = x[1][1] + x[2][1];
>
>     // Get coordinates and map to the reference (FIAT) element
>     double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
>     double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
>
>     // Compute number of derivatives.
>     unsigned int num_derivatives = 1;
>     for (unsigned int r = 0; r < n; r++)
>     {
>       num_derivatives *= 2;
>     }// end loop over 'r'
>
>     // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
>     unsigned int **combinations = new unsigned int *[num_derivatives];
>     for (unsigned int row = 0; row < num_derivatives; row++)
>     {
>       combinations[row] = new unsigned int [n];
>       for (unsigned int col = 0; col < n; col++)
>         combinations[row][col] = 0;
>     }
>
>     // Generate combinations of derivatives
>     for (unsigned int row = 1; row < num_derivatives; row++)
>     {
>       for (unsigned int num = 0; num < row; num++)
>       {
>         for (unsigned int col = n-1; col+1 > 0; col--)
>         {
>           if (combinations[row][col] + 1 > 1)
>             combinations[row][col] = 0;
>           else
>           {
>             combinations[row][col] += 1;
>             break;
>           }
>         }
>       }
>     }
>
>     // Compute inverse of Jacobian
>     const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
>
>     // Declare transformation matrix
>     // Declare pointer to two dimensional array and initialise
>     double **transform = new double *[num_derivatives];
>
>     for (unsigned int j = 0; j < num_derivatives; j++)
>     {
>       transform[j] = new double [num_derivatives];
>       for (unsigned int k = 0; k < num_derivatives; k++)
>         transform[j][k] = 1;
>     }
>
>     // Construct transformation matrix
>     for (unsigned int row = 0; row < num_derivatives; row++)
>     {
>       for (unsigned int col = 0; col < num_derivatives; col++)
>       {
>         for (unsigned int k = 0; k < n; k++)
>           transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
>       }
>     }
>
>     // Reset values. Assuming that values is always an array.
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       values[r] = 0.000000000000000;
>     }// end loop over 'r'
>
>     // Map degree of freedom to element degree of freedom
>     const unsigned int dof = i;
>
>     // Array of basisvalues.
>     double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
>
>     // Declare helper variables.
>     unsigned int rr = 0;
>     unsigned int ss = 0;
>     double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
>
>     // Compute basisvalues.
>     basisvalues[0] = 1.000000000000000;
>     basisvalues[1] = tmp0;
>     for (unsigned int r = 0; r < 1; r++)
>     {
>       rr = (r + 1)*(r + 1 + 1)/2 + 1;
>       ss = r*(r + 1)/2;
>       basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
>     }// end loop over 'r'
>     for (unsigned int r = 0; r < 2; r++)
>     {
>       for (unsigned int s = 0; s < 2 - r; s++)
>       {
>         rr = (r + s)*(r + s + 1)/2 + s;
>         basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
>       }// end loop over 's'
>     }// end loop over 'r'
>
>     // Table(s) of coefficients.
>     static const double coefficients0[3][3] = \
>     {{0.471404520791032, -0.288675134594813, -0.166666666666667},
>     {0.471404520791032, 0.288675134594813, -0.166666666666667},
>     {0.471404520791032, 0.000000000000000, 0.333333333333333}};
>
>     // Tables of derivatives of the polynomial base (transpose).
>     static const double dmats0[3][3] = \
>     {{0.000000000000000, 0.000000000000000, 0.000000000000000},
>     {4.898979485566356, 0.000000000000000, 0.000000000000000},
>     {0.000000000000000, 0.000000000000000, 0.000000000000000}};
>
>     static const double dmats1[3][3] = \
>     {{0.000000000000000, 0.000000000000000, 0.000000000000000},
>     {2.449489742783178, 0.000000000000000, 0.000000000000000},
>     {4.242640687119285, 0.000000000000000, 0.000000000000000}};
>
>     // Compute reference derivatives.
>     // Declare pointer to array of derivatives on FIAT element.
>     double *derivatives = new double[num_derivatives];
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       derivatives[r] = 0.000000000000000;
>     }// end loop over 'r'
>
>     // Declare derivative matrix (of polynomial basis).
>     double dmats[3][3] = \
>     {{1.000000000000000, 0.000000000000000, 0.000000000000000},
>     {0.000000000000000, 1.000000000000000, 0.000000000000000},
>     {0.000000000000000, 0.000000000000000, 1.000000000000000}};
>
>     // Declare (auxiliary) derivative matrix (of polynomial basis).
>     double dmats_old[3][3] = \
>     {{1.000000000000000, 0.000000000000000, 0.000000000000000},
>     {0.000000000000000, 1.000000000000000, 0.000000000000000},
>     {0.000000000000000, 0.000000000000000, 1.000000000000000}};
>
>     // Loop possible derivatives.
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       // Resetting dmats values to compute next derivative.
>       for (unsigned int t = 0; t < 3; t++)
>       {
>         for (unsigned int u = 0; u < 3; u++)
>         {
>           dmats[t][u] = 0.000000000000000;
>           if (t == u)
>           {
>           dmats[t][u] = 1.000000000000000;
>           }
>
>         }// end loop over 'u'
>       }// end loop over 't'
>
>       // Looping derivative order to generate dmats.
>       for (unsigned int s = 0; s < n; s++)
>       {
>         // Updating dmats_old with new values and resetting dmats.
>         for (unsigned int t = 0; t < 3; t++)
>         {
>           for (unsigned int u = 0; u < 3; u++)
>           {
>             dmats_old[t][u] = dmats[t][u];
>             dmats[t][u] = 0.000000000000000;
>           }// end loop over 'u'
>         }// end loop over 't'
>
>         // Update dmats using an inner product.
>         if (combinations[r][s] == 0)
>         {
>         for (unsigned int t = 0; t < 3; t++)
>         {
>           for (unsigned int u = 0; u < 3; u++)
>           {
>             for (unsigned int tu = 0; tu < 3; tu++)
>             {
>               dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
>             }// end loop over 'tu'
>           }// end loop over 'u'
>         }// end loop over 't'
>         }
>
>         if (combinations[r][s] == 1)
>         {
>         for (unsigned int t = 0; t < 3; t++)
>         {
>           for (unsigned int u = 0; u < 3; u++)
>           {
>             for (unsigned int tu = 0; tu < 3; tu++)
>             {
>               dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
>             }// end loop over 'tu'
>           }// end loop over 'u'
>         }// end loop over 't'
>         }
>
>       }// end loop over 's'
>       for (unsigned int s = 0; s < 3; s++)
>       {
>         for (unsigned int t = 0; t < 3; t++)
>         {
>           derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
>         }// end loop over 't'
>       }// end loop over 's'
>     }// end loop over 'r'
>
>     // Transform derivatives back to physical element
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       for (unsigned int s = 0; s < num_derivatives; s++)
>       {
>         values[r] += transform[r][s]*derivatives[s];
>       }// end loop over 's'
>     }// end loop over 'r'
>
>     // Delete pointer to array of derivatives on FIAT element
>     delete [] derivatives;
>
>     // Delete pointer to array of combinations of derivatives and transform
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       delete [] combinations[r];
>     }// end loop over 'r'
>     delete [] combinations;
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       delete [] transform[r];
>     }// end loop over 'r'
>     delete [] transform;
>   }
>
>   /// Evaluate order n derivatives of all basis functions at given point in cell
>   virtual void evaluate_basis_derivatives_all(unsigned int n,
>                                               double* values,
>                                               const double* coordinates,
>                                               const ufc::cell& c) const
>   {
>     // Compute number of derivatives.
>     unsigned int num_derivatives = 1;
>     for (unsigned int r = 0; r < n; r++)
>     {
>       num_derivatives *= 2;
>     }// end loop over 'r'
>
>     // Helper variable to hold values of a single dof.
>     double *dof_values = new double[num_derivatives];
>     for (unsigned int r = 0; r < num_derivatives; r++)
>     {
>       dof_values[r] = 0.000000000000000;
>     }// end loop over 'r'
>
>     // Loop dofs and call evaluate_basis_derivatives.
>     for (unsigned int r = 0; r < 3; r++)
>     {
>       evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
>       for (unsigned int s = 0; s < num_derivatives; s++)
>       {
>         values[r*num_derivatives + s] = dof_values[s];
>       }// end loop over 's'
>     }// end loop over 'r'
>
>     // Delete pointer.
>     delete [] dof_values;
>   }
>
>   /// Evaluate linear functional for dof i on the function f
>   virtual double evaluate_dof(unsigned int i,
>                               const ufc::function& f,
>                               const ufc::cell& c) const
>   {
>     // Declare variables for result of evaluation.
>     double vals[1];
>
>     // Declare variable for physical coordinates.
>     double y[2];
>     const double * const * x = c.coordinates;
>     switch (i)
>     {
>     case 0:
>       {
>         y[0] = x[0][0];
>       y[1] = x[0][1];
>       f.evaluate(vals, y, c);
>       return vals[0];
>         break;
>       }
>     case 1:
>       {
>         y[0] = x[1][0];
>       y[1] = x[1][1];
>       f.evaluate(vals, y, c);
>       return vals[0];
>         break;
>       }
>     case 2:
>       {
>         y[0] = x[2][0];
>       y[1] = x[2][1];
>       f.evaluate(vals, y, c);
>       return vals[0];
>         break;
>       }
>     }
>
>     return 0.000000000000000;
>   }
>
>   /// Evaluate linear functionals for all dofs on the function f
>   virtual void evaluate_dofs(double* values,
>                              const ufc::function& f,
>                              const ufc::cell& c) const
>   {
>     // Declare variables for result of evaluation.
>     double vals[1];
>
>     // Declare variable for physical coordinates.
>     double y[2];
>     const double * const * x = c.coordinates;
>     y[0] = x[0][0];
>     y[1] = x[0][1];
>     f.evaluate(vals, y, c);
>     values[0] = vals[0];
>     y[0] = x[1][0];
>     y[1] = x[1][1];
>     f.evaluate(vals, y, c);
>     values[1] = vals[0];
>     y[0] = x[2][0];
>     y[1] = x[2][1];
>     f.evaluate(vals, y, c);
>     values[2] = vals[0];
>   }
>
>   /// Interpolate vertex values from dof values
>   virtual void interpolate_vertex_values(double* vertex_values,
>                                          const double* dof_values,
>                                          const ufc::cell& c) const
>   {
>     // Evaluate function and change variables
>     vertex_values[0] = dof_values[0];
>     vertex_values[1] = dof_values[1];
>     vertex_values[2] = dof_values[2];
>   }
>
>   /// Return the number of sub elements (for a mixed element)
>   virtual unsigned int num_sub_elements() const
>   {
>     return 0;
>   }
>
>   /// Create a new finite element for sub element i (for a mixed element)
>   virtual ufc::finite_element* create_sub_element(unsigned int i) const
>   {
> return 0;
>   }
>
> };
>
> /// This class defines the interface for a local-to-global mapping of
> /// degrees of freedom (dofs).
>
>
> class poisson_auxiliary_dof_map_0: public ufc::dof_map
> {
> private:
>
>   unsigned int _global_dimension;
> public:
>
>   /// Constructor
>   poisson_auxiliary_dof_map_0() :ufc::dof_map()
>   {
>     _global_dimension = 0;
>   }
>   /// Destructor
>   virtual ~poisson_auxiliary_dof_map_0()
>   {
>     // Do nothing
>   }
>
>   /// Return a string identifying the dof map
>   virtual const char* signature() const
>   {
>     return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
>   }
>
>   /// Return true iff mesh entities of topological dimension d are needed
>   virtual bool needs_mesh_entities(unsigned int d) const
>   {
>     switch (d)
>     {
>     case 0:
>       {
>         return true;
>         break;
>       }
>     case 1:
>       {
>         return false;
>         break;
>       }
>     case 2:
>       {
>         return false;
>         break;
>       }
>     }
>
>     return false;
>   }
>
>   /// Initialize dof map for mesh (return true iff init_cell() is needed)
>   virtual bool init_mesh(const ufc::mesh& m)
>   {
>     _global_dimension = m.num_entities[0];
>     return false;
>   }
>
>   /// Initialize dof map for given cell
>   virtual void init_cell(const ufc::mesh& m,
>                          const ufc::cell& c)
>   {
>     // Do nothing
>   }
>
>   /// Finish initialization of dof map for cells
>   virtual void init_cell_finalize()
>   {
>     // Do nothing
>   }
>
>   /// Return the dimension of the global finite element function space
>   virtual unsigned int global_dimension() const
>   {
>     return _global_dimension;
>   }
>
>   /// Return the dimension of the local finite element function space for a cell
>   virtual unsigned int local_dimension(const ufc::cell& c) const
>   {
>     return 3;
>   }
>
>   /// Return the maximum dimension of the local finite element function space
>   virtual unsigned int max_local_dimension() const
>   {
>     return 3;
>   }
>
>
>   // Return the geometric dimension of the coordinates this dof map provides
>   virtual unsigned int geometric_dimension() const
>   {
>     return 2;
>   }
>
>   /// Return the number of dofs on each cell facet
>   virtual unsigned int num_facet_dofs() const
>   {
>     return 2;
>   }
>
>   /// Return the number of dofs associated with each cell entity of dimension d
>   virtual unsigned int num_entity_dofs(unsigned int d) const
>   {
>     switch (d)
>     {
>     case 0:
>       {
>         return 1;
>         break;
>       }
>     case 1:
>       {
>         return 0;
>         break;
>       }
>     case 2:
>       {
>         return 0;
>         break;
>       }
>     }
>
>     return 0;
>   }
>
>   /// Tabulate the local-to-global mapping of dofs on a cell
>   virtual void tabulate_dofs(unsigned int* dofs,
>                              const ufc::mesh& m,
>                              const ufc::cell& c) const
>   {
>     dofs[0] = c.entity_indices[0][0];
>     dofs[1] = c.entity_indices[0][1];
>     dofs[2] = c.entity_indices[0][2];
>   }
>
>   /// Tabulate the local-to-local mapping from facet dofs to cell dofs
>   virtual void tabulate_facet_dofs(unsigned int* dofs,
>                                    unsigned int facet) const
>   {
>     switch (facet)
>     {
>     case 0:
>       {
>         dofs[0] = 1;
>       dofs[1] = 2;
>         break;
>       }
>     case 1:
>       {
>         dofs[0] = 0;
>       dofs[1] = 2;
>         break;
>       }
>     case 2:
>       {
>         dofs[0] = 0;
>       dofs[1] = 1;
>         break;
>       }
>     }
>
>   }
>
>   /// Tabulate the local-to-local mapping of dofs on entity (d, i)
>   virtual void tabulate_entity_dofs(unsigned int* dofs,
>                                     unsigned int d, unsigned int i) const
>   {
>     if (d > 2)
>     {
>     throw std::runtime_error("d is larger than dimension (2)");
>     }
>
>     switch (d)
>     {
>     case 0:
>       {
>         if (i > 2)
>       {
>       throw std::runtime_error("i is larger than number of entities (2)");
>       }
>
>       switch (i)
>       {
>       case 0:
>         {
>           dofs[0] = 0;
>           break;
>         }
>       case 1:
>         {
>           dofs[0] = 1;
>           break;
>         }
>       case 2:
>         {
>           dofs[0] = 2;
>           break;
>         }
>       }
>
>         break;
>       }
>     case 1:
>       {
>
>         break;
>       }
>     case 2:
>       {
>
>         break;
>       }
>     }
>
>   }
>
>   /// Tabulate the coordinates of all dofs on a cell
>   virtual void tabulate_coordinates(double** coordinates,
>                                     const ufc::cell& c) const
>   {
>     const double * const * x = c.coordinates;
>
>     coordinates[0][0] = x[0][0];
>     coordinates[0][1] = x[0][1];
>     coordinates[1][0] = x[1][0];
>     coordinates[1][1] = x[1][1];
>     coordinates[2][0] = x[2][0];
>     coordinates[2][1] = x[2][1];
>   }
>
>   /// Return the number of sub dof maps (for a mixed element)
>   virtual unsigned int num_sub_dof_maps() const
>   {
>     return 0;
>   }
>
>   /// Create a new dof_map for sub dof map i (for a mixed element)
>   virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
>   {
> return 0;
>   }
>
> };
>
> /// This class defines the interface for the assembly of the global
> /// tensor corresponding to a form with r + n arguments, that is, a
> /// mapping
> ///
> ///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
> ///
> /// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
> /// global tensor A is defined by
> ///
> ///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
> ///
> /// where each argument Vj represents the application to the
> /// sequence of basis functions of Vj and w1, w2, ..., wn are given
> /// fixed functions (coefficients).
>
> class poisson_auxiliary_form_0: public ufc::form
> {
>
> public:
>
>   /// Constructor
>   poisson_auxiliary_form_0() : ufc::form()
>   {
>     // Do nothing
>   }
>
>   /// Destructor
>   virtual ~poisson_auxiliary_form_0()
>   {
>     // Do nothing
>   }
>
>   /// Return a string identifying the form
>   virtual const char* signature() const
>   {
>     return "auxiliary form";
>   }
>
>   /// Return the rank of the global tensor (r)
>   virtual unsigned int rank() const
>   {
>     return 0;
>   }
>
>   /// Return the number of coefficients (n)
>   virtual unsigned int num_coefficients() const
>   {
>     return 0;
>   }
>
>   /// Return the number of cell integrals
>   virtual unsigned int num_cell_integrals() const
>   {
>     return 0;
>   }
>
>   /// Return the number of exterior facet integrals
>   virtual unsigned int num_exterior_facet_integrals() const
>   {
>     return 0;
>   }
>
>   /// Return the number of interior facet integrals
>   virtual unsigned int num_interior_facet_integrals() const
>   {
>     return 0;
>   }
>
>   /// Create a new finite element for argument function i
>   virtual ufc::finite_element* create_finite_element(unsigned int i) const
>   {
>     switch (i)
>     {
>     case 0:
>       {
>         return new poisson_finite_element_0();
>         break;
>       }
>     }
>
>     return 0;
>   }
>
>   /// Create a new dof map for argument function i
>   virtual ufc::dof_map* create_dof_map(unsigned int i) const
>   {
> switch (i)
> {
> case 0:
>   {
>     return new poisson_auxiliary_dof_map_0();
>     break;
>   }
> }
>
> return 0;
>   }
>
>   /// Create a new cell integral on sub domain i
>   virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
>   {
>     return 0;
>   }
>
>   /// Create a new exterior facet integral on sub domain i
>   virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
>   {
>     return 0;
>   }
>
>   /// Create a new interior facet integral on sub domain i
>   virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
>   {
>     return 0;
>   }
>
> };
>
> /// This class defines the interface for the assembly of the global
> /// tensor corresponding to a form with r + n arguments, that is, a
> /// mapping
> ///
> ///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
> ///
> /// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
> /// global tensor A is defined by
> ///
> ///     A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
> ///
> /// where each argument Vj represents the application to the
> /// sequence of basis functions of Vj and w1, w2, ..., wn are given
> /// fixed functions (coefficients).
>
> class poisson_auxiliary_form_1: public ufc::form
> {
>
> public:
>
>   /// Constructor
>   poisson_auxiliary_form_1() : ufc::form()
>   {
>     // Do nothing
>   }
>
>   /// Destructor
>   virtual ~poisson_auxiliary_form_1()
>   {
>     // Do nothing
>   }
>
>   /// Return a string identifying the form
>   virtual const char* signature() const
>   {
>     return "auxiliary form";
>   }
>
>   /// Return the rank of the global tensor (r)
>   virtual unsigned int rank() const
>   {
>     return 0;
>   }
>
>   /// Return the number of coefficients (n)
>   virtual unsigned int num_coefficients() const
>   {
>     return 0;
>   }
>
>   /// Return the number of cell integrals
>   virtual unsigned int num_cell_integrals() const
>   {
>     return 0;
>   }
>
>   /// Return the number of exterior facet integrals
>   virtual unsigned int num_exterior_facet_integrals() const
>   {
>     return 0;
>   }
>
>   /// Return the number of interior facet integrals
>   virtual unsigned int num_interior_facet_integrals() const
>   {
>     return 0;
>   }
>
>   /// Create a new finite element for argument function i
>   virtual ufc::finite_element* create_finite_element(unsigned int i) const
>   {
>     switch (i)
>     {
>     case 0:
>       {
>         return new poisson_finite_element_0();
>         break;
>       }
>     }
>
>     return 0;
>   }
>
>   /// Create a new dof map for argument function i
>   virtual ufc::dof_map* create_dof_map(unsigned int i) const
>   {
> switch (i)
> {
> case 0:
>   {
>     return new poisson_auxiliary_dof_map_0();
>     break;
>   }
> }
>
> return 0;
>   }
>
>   /// Create a new cell integral on sub domain i
>   virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
>   {
>     return 0;
>   }
>
>   /// Create a new exterior facet integral on sub domain i
>   virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
>   {
>     return 0;
>   }
>
>   /// Create a new interior facet integral on sub domain i
>   virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
>   {
>     return 0;
>   }
>
> };
>
> /// This class defines the interface for post-processing on vector x
> /// to obtain x0, u and j where,
> ///
> /// - x is the solution vector containing standard and enriched degrees of freedom
> /// defined on continuous/discontinuous space
> /// - u is the standard part of solution vector defined on continuous space
> /// - j is the enriched part pf solution vector defined on continuous space
> /// - x0 is the result vector, equall to u + j, defined on continuous space
> /// by considering enrichement function
> //
>
> // Dolfin includes
> #include <dolfin/common/NoDeleter.h>
> #include <dolfin/mesh/Mesh.h>
> #include <dolfin/fem/DofMap.h>
> #include <dolfin/la/GenericVector.h>
>
> // PartitionOfUnity includes
> #include <pum/PostProcess.h>
> #include <pum/FunctionSpace.h>
> #include <pum/PUM.h>
> #include <pum/GenericSurface.h>
>
> namespace Poisson
> {
>   class PostProcess:  public pum::PostProcess
>   {
>     const dolfin::Mesh& mesh;
>     pum::FunctionSpace& function_space;
>
>   public:
>
>     /// Constructor
>     PostProcess(pum::FunctionSpace& function_space):
>                 pum::PostProcess(function_space.mesh()), mesh(function_space.mesh()),
>                 function_space(function_space)
>     {
>       // Do nothing
>     }
>
>     /// Destructor
>     ~PostProcess()
>     {
>        // Do nothing
>     }
>
>     /// Return a string identifying the underling element
>     const char* signature() const
>     {
>       return "Interpolating results to the continuous space of EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None)))";
>     }
>
>     /// Obtain result vector 'x0' from solution vector 'x'
>     void interpolate(const dolfin::GenericVector& x, dolfin::GenericVector& x0) const
>     {
>
>       // Compute number of standard dofs for field 0
>       dolfin::DofMap dof_map_0(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh);
>       unsigned int num_standard_dofs_0 = dof_map_0.global_dimension();
>
>       double value;
>       double h;
>
>       /// Selecting standard degrees of freedom related to field 0 from the solution vector
>       double* values_0 = new double[num_standard_dofs_0];
>       unsigned int* positions_0 = new unsigned int [num_standard_dofs_0];
>
>       for (unsigned int i = 0; i < num_standard_dofs_0; ++i)
>         positions_0[i] = i;
>
>       x.get(values_0, num_standard_dofs_0, positions_0);
>       x0.set(values_0, num_standard_dofs_0, positions_0);
>
>       x0.apply("insert");
>
>       const std::vector<const pum::GenericPUM*>& pums = function_space.pums;
>
>       /// Selecting enriched degrees of freedom related to field 0 from the solution vector
>       std::vector <std::vector<unsigned int> > enhanced_dof_maps_0;
>       enhanced_dof_maps_0.resize(num_standard_dofs_0);
>
>       std::vector<unsigned int> enhanced_dof_values_0;
>       enhanced_dof_values_0.resize(num_standard_dofs_0);
>
>       compute_enhanced_dof_maps(*pums[0], dof_map_0, enhanced_dof_maps_0);
>       compute_enhanced_dof_values(*pums[0], dof_map_0, enhanced_dof_values_0);
>
>       //compute_enhanced_dof_maps(pum_0, dof_map_0, enhanced_dof_maps_0);
>       //compute_enhanced_dof_values(pum_0, dof_map_0, enhanced_dof_values_0);
>
>       for (unsigned int i = 0; i != num_standard_dofs_0; ++i)
>       {
>         unsigned int pos = i;
>
>         for (std::vector<unsigned int>::const_iterator it = enhanced_dof_maps_0[i].begin();
>                              it != enhanced_dof_maps_0[i].end(); ++it)
>         {
>           h = enhanced_dof_values_0[i];
>           unsigned int pos_n = *it + num_standard_dofs_0;
>
>           x.get(&value, 1, &pos_n);
>           value *= h;
>           x0.add(&value, 1, &pos);
>         }
>       }
>
>
>       // memory clean up
>
>       delete[] values_0;
>       delete[] positions_0;
>
>       x0.apply("add");
>     }
>
>     /// Obtain continuous u and discontinuous j parts of solution vector 'x'
>     void interpolate(const dolfin::GenericVector& x, dolfin::GenericVector& u, dolfin::GenericVector& j) const
>     {
>
>       // Compute number of standard dofs for field 0
>       dolfin::DofMap dof_map_0(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh);
>       unsigned int num_standard_dofs_0 = dof_map_0.global_dimension();
>
>       double value;
>
>       /// Selecting standard degrees of freedom related to field 0 from the solution vector
>       double* values_0 = new double[num_standard_dofs_0];
>       unsigned int* positions_0 = new unsigned int [num_standard_dofs_0];
>
>       for (unsigned int i = 0; i < num_standard_dofs_0; ++i)
>         positions_0[i] = i;
>
>       x.get(values_0, num_standard_dofs_0, positions_0);
>       u.set(values_0, num_standard_dofs_0, positions_0);
>
>
>       const std::vector<const pum::GenericPUM*>& pums = function_space.pums;
>
>       /// Selecting enriched degrees of freedom related to field 0 from the solution vector
>       std::vector <std::vector<unsigned int> > enhanced_dof_maps_0;
>       enhanced_dof_maps_0.resize(num_standard_dofs_0);
>
>       // Compute enhanced dof maps
>       compute_enhanced_dof_maps(*pums[0], dof_map_0, enhanced_dof_maps_0);
>
>
>       for (unsigned int i = 0; i != num_standard_dofs_0; ++i)
>       {
>         unsigned int pos = i ;
>
>         for (std::vector<unsigned int>::const_iterator it = enhanced_dof_maps_0[i].begin();
>                            it != enhanced_dof_maps_0[i].end(); ++it)
>         {
>           unsigned int pos_n = *it + num_standard_dofs_0;
>
>           x.get(&value, 1, &pos_n);
>           j.set(&value, 1, &pos);
>         }
>       }
>
>
>       // memory clean up
>
>       delete[] values_0;
>       delete[] positions_0;
>
>       u.apply("insert");
>       j.apply("insert");
>     }
>
>   };
> }
>
> // DOLFIN wrappers
>
> // Standard library includes
> #include <string>
>
> // DOLFIN includes
> #include <dolfin/common/NoDeleter.h>
> #include <dolfin/fem/FiniteElement.h>
> #include <dolfin/fem/DofMap.h>
> #include <dolfin/fem/Form.h>
> #include <dolfin/function/FunctionSpace.h>
> #include <dolfin/function/Function.h>
> #include <dolfin/function/GenericFunction.h>
> #include <dolfin/function/CoefficientAssigner.h>
>
> namespace Poisson
> {
>
> class CoefficientSpace_f: public dolfin::FunctionSpace
> {
> public:
>
>
>   CoefficientSpace_f(const dolfin::Mesh & mesh):
>       dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
>                             (boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
>                             (new poisson_dof_map_0()), mesh)))
>   {
>     // Do nothing
>   }
>
>   CoefficientSpace_f(dolfin::Mesh & mesh):
>     dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                           boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
>                           boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh)))
>   {
>     // Do nothing
>   }
>
>   CoefficientSpace_f(boost::shared_ptr<dolfin::Mesh> mesh):
>       dolfin::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), *mesh)))
>   {
>       // Do nothing
>   }
>
>   CoefficientSpace_f(boost::shared_ptr<const dolfin::Mesh> mesh):
>       dolfin::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), *mesh)))
>   {
>       // Do nothing
>   }
>
>
>   ~CoefficientSpace_f()
>   {
>   }
>
> };
>
> class CoefficientSpace_g: public dolfin::FunctionSpace
> {
> public:
>
>
>   CoefficientSpace_g(const dolfin::Mesh & mesh):
>       dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
>                             (boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
>                             (new poisson_dof_map_0()), mesh)))
>   {
>     // Do nothing
>   }
>
>   CoefficientSpace_g(dolfin::Mesh & mesh):
>     dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                           boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
>                           boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh)))
>   {
>     // Do nothing
>   }
>
>   CoefficientSpace_g(boost::shared_ptr<dolfin::Mesh> mesh):
>       dolfin::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), *mesh)))
>   {
>       // Do nothing
>   }
>
>   CoefficientSpace_g(boost::shared_ptr<const dolfin::Mesh> mesh):
>       dolfin::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), *mesh)))
>   {
>       // Do nothing
>   }
>
>
>   ~CoefficientSpace_g()
>   {
>   }
>
> };
>
> class CoefficientSpace_w: public dolfin::FunctionSpace
> {
> public:
>
>
>   CoefficientSpace_w(const dolfin::Mesh & mesh):
>       dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
>                             (boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
>                             (new poisson_dof_map_0()), mesh)))
>   {
>     // Do nothing
>   }
>
>   CoefficientSpace_w(dolfin::Mesh & mesh):
>     dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                           boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
>                           boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh)))
>   {
>     // Do nothing
>   }
>
>   CoefficientSpace_w(boost::shared_ptr<dolfin::Mesh> mesh):
>       dolfin::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), *mesh)))
>   {
>       // Do nothing
>   }
>
>   CoefficientSpace_w(boost::shared_ptr<const dolfin::Mesh> mesh):
>       dolfin::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), *mesh)))
>   {
>       // Do nothing
>   }
>
>
>   ~CoefficientSpace_w()
>   {
>   }
>
> };
>
> class Form_0_FunctionSpace_0: public pum::FunctionSpace
> {
> public:
>
>
>   Form_0_FunctionSpace_0(const dolfin::Mesh & mesh ,const std::vector<const pum::GenericPUM*>& pums):
>       pum::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
>                             (boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
>                             (new poisson_dof_map_1(pums)), mesh)) ,pums)
>   {
>     // Do nothing
>   }
>
>   Form_0_FunctionSpace_0(dolfin::Mesh & mesh ,const std::vector<const pum::GenericPUM*>& pums):
>     pum::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                           boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
>                           boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), mesh)) ,pums)
>   {
>     // Do nothing
>   }
>
>   Form_0_FunctionSpace_0(boost::shared_ptr<dolfin::Mesh> mesh ,const std::vector<const pum::GenericPUM*>& pums):
>       pum::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), *mesh)) ,pums)
>   {
>       // Do nothing
>   }
>
>   Form_0_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh ,const std::vector<const pum::GenericPUM*>& pums):
>       pum::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), *mesh)) ,pums)
>   {
>       // Do nothing
>   }
>
>
>   ~Form_0_FunctionSpace_0()
>   {
>   }
>
> };
>
> class Form_0_FunctionSpace_1: public pum::FunctionSpace
> {
> public:
>
>
>   Form_0_FunctionSpace_1(const dolfin::Mesh & mesh ,const std::vector<const pum::GenericPUM*>& pums):
>       pum::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
>                             (boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
>                             (new poisson_dof_map_1(pums)), mesh)) ,pums)
>   {
>     // Do nothing
>   }
>
>   Form_0_FunctionSpace_1(dolfin::Mesh & mesh ,const std::vector<const pum::GenericPUM*>& pums):
>     pum::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                           boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
>                           boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), mesh)) ,pums)
>   {
>     // Do nothing
>   }
>
>   Form_0_FunctionSpace_1(boost::shared_ptr<dolfin::Mesh> mesh ,const std::vector<const pum::GenericPUM*>& pums):
>       pum::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), *mesh)) ,pums)
>   {
>       // Do nothing
>   }
>
>   Form_0_FunctionSpace_1(boost::shared_ptr<const dolfin::Mesh> mesh ,const std::vector<const pum::GenericPUM*>& pums):
>       pum::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), *mesh)) ,pums)
>   {
>       // Do nothing
>   }
>
>
>   ~Form_0_FunctionSpace_1()
>   {
>   }
>
> };
>
> typedef CoefficientSpace_w Form_0_FunctionSpace_2;
>
> class Form_0: public dolfin::Form
> {
> public:
>
>   // Constructor
>   Form_0(const pum::FunctionSpace& V0, const pum::FunctionSpace& V1):
>     dolfin::Form(2, 1), w(*this, 0)
>   {
>     _function_spaces[0] = reference_to_no_delete_pointer(V0);
>     _function_spaces[1] = reference_to_no_delete_pointer(V1);
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_0(V0.pums));
>   }
>
>   // Constructor
>   Form_0(const pum::FunctionSpace& V0, const pum::FunctionSpace& V1, dolfin::GenericFunction & w):
>     dolfin::Form(2, 1), w(*this, 0)
>   {
>     _function_spaces[0] = reference_to_no_delete_pointer(V0);
>     _function_spaces[1] = reference_to_no_delete_pointer(V1);
>
>     this->w = w;
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_0(V0.pums));
>   }
>
>   // Constructor
>   Form_0(const pum::FunctionSpace& V0, const pum::FunctionSpace& V1, boost::shared_ptr<dolfin::GenericFunction> w):
>     dolfin::Form(2, 1), w(*this, 0)
>   {
>     _function_spaces[0] = reference_to_no_delete_pointer(V0);
>     _function_spaces[1] = reference_to_no_delete_pointer(V1);
>
>     this->w = *w;
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_0(V0.pums));
>   }
>
>   // Constructor
>   Form_0(boost::shared_ptr<const pum::FunctionSpace> V0, boost::shared_ptr<const pum::FunctionSpace> V1):
>     dolfin::Form(2, 1), w(*this, 0)
>   {
>     _function_spaces[0] = V0;
>     _function_spaces[1] = V1;
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_0(V0->pums));
>   }
>
>   // Constructor
>   Form_0(boost::shared_ptr<const pum::FunctionSpace> V0, boost::shared_ptr<const pum::FunctionSpace> V1, dolfin::GenericFunction & w):
>     dolfin::Form(2, 1), w(*this, 0)
>   {
>     _function_spaces[0] = V0;
>     _function_spaces[1] = V1;
>
>     this->w = w;
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_0(V0->pums));
>   }
>
>   // Constructor
>   Form_0(boost::shared_ptr<const pum::FunctionSpace> V0, boost::shared_ptr<const pum::FunctionSpace> V1, boost::shared_ptr<dolfin::GenericFunction> w):
>     dolfin::Form(2, 1), w(*this, 0)
>   {
>     _function_spaces[0] = V0;
>     _function_spaces[1] = V1;
>
>     this->w = *w;
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_0(V0->pums));
>   }
>
>   // Destructor
>   ~Form_0()
>   {}
>
>   /// Return the number of the coefficient with this name
>   virtual dolfin::uint coefficient_number(const std::string& name) const
>   {
>     if(name == "w") return 0;
>     dolfin::error("Invalid coefficient.");
>     return 0;
>   }
>
>   /// Return the name of the coefficient with this number
>   virtual std::string coefficient_name(dolfin::uint i) const
>   {
>     switch(i)
>     {
>       case 0: return "w";
>     }
>     dolfin::error("Invalid coefficient.");
>     return "unnamed";
>   }
>
>   // Typedefs
>   typedef Form_0_FunctionSpace_0 TestSpace;
>   typedef Form_0_FunctionSpace_1 TrialSpace;
>   typedef Form_0_FunctionSpace_2 CoefficientSpace_w;
>
>   // Coefficients
>   dolfin::CoefficientAssigner w;
> };
>
> class Form_1_FunctionSpace_0: public pum::FunctionSpace
> {
> public:
>
>
>   Form_1_FunctionSpace_0(const dolfin::Mesh & mesh ,const std::vector<const pum::GenericPUM*>& pums):
>       pum::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
>                             (boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
>                             (new poisson_dof_map_1(pums)), mesh)) ,pums)
>   {
>     // Do nothing
>   }
>
>   Form_1_FunctionSpace_0(dolfin::Mesh & mesh ,const std::vector<const pum::GenericPUM*>& pums):
>     pum::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                           boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
>                           boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), mesh)) ,pums)
>   {
>     // Do nothing
>   }
>
>   Form_1_FunctionSpace_0(boost::shared_ptr<dolfin::Mesh> mesh ,const std::vector<const pum::GenericPUM*>& pums):
>       pum::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), *mesh)) ,pums)
>   {
>       // Do nothing
>   }
>
>   Form_1_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh ,const std::vector<const pum::GenericPUM*>& pums):
>       pum::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), *mesh)) ,pums)
>   {
>       // Do nothing
>   }
>
>
>   ~Form_1_FunctionSpace_0()
>   {
>   }
>
> };
>
> typedef CoefficientSpace_f Form_1_FunctionSpace_1;
>
> typedef CoefficientSpace_g Form_1_FunctionSpace_2;
>
> class Form_1: public dolfin::Form
> {
> public:
>
>   // Constructor
>   Form_1(const pum::FunctionSpace& V0):
>     dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
>   {
>     _function_spaces[0] = reference_to_no_delete_pointer(V0);
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_1(V0.pums));
>   }
>
>   // Constructor
>   Form_1(const pum::FunctionSpace& V0, dolfin::GenericFunction & f, dolfin::GenericFunction & g):
>     dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
>   {
>     _function_spaces[0] = reference_to_no_delete_pointer(V0);
>
>     this->f = f;
>     this->g = g;
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_1(V0.pums));
>   }
>
>   // Constructor
>   Form_1(const pum::FunctionSpace& V0, boost::shared_ptr<dolfin::GenericFunction> f, boost::shared_ptr<dolfin::GenericFunction> g):
>     dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
>   {
>     _function_spaces[0] = reference_to_no_delete_pointer(V0);
>
>     this->f = *f;
>     this->g = *g;
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_1(V0.pums));
>   }
>
>   // Constructor
>   Form_1(boost::shared_ptr<const pum::FunctionSpace> V0):
>     dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
>   {
>     _function_spaces[0] = V0;
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_1(V0->pums));
>   }
>
>   // Constructor
>   Form_1(boost::shared_ptr<const pum::FunctionSpace> V0, dolfin::GenericFunction & f, dolfin::GenericFunction & g):
>     dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
>   {
>     _function_spaces[0] = V0;
>
>     this->f = f;
>     this->g = g;
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_1(V0->pums));
>   }
>
>   // Constructor
>   Form_1(boost::shared_ptr<const pum::FunctionSpace> V0, boost::shared_ptr<dolfin::GenericFunction> f, boost::shared_ptr<dolfin::GenericFunction> g):
>     dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
>   {
>     _function_spaces[0] = V0;
>
>     this->f = *f;
>     this->g = *g;
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_1(V0->pums));
>   }
>
>   // Destructor
>   ~Form_1()
>   {}
>
>   /// Return the number of the coefficient with this name
>   virtual dolfin::uint coefficient_number(const std::string& name) const
>   {
>     if(name == "f") return 0;
>     else if(name == "g") return 1;
>     dolfin::error("Invalid coefficient.");
>     return 0;
>   }
>
>   /// Return the name of the coefficient with this number
>   virtual std::string coefficient_name(dolfin::uint i) const
>   {
>     switch(i)
>     {
>       case 0: return "f";
>       case 1: return "g";
>     }
>     dolfin::error("Invalid coefficient.");
>     return "unnamed";
>   }
>
>   // Typedefs
>   typedef Form_1_FunctionSpace_0 TestSpace;
>   typedef Form_1_FunctionSpace_1 CoefficientSpace_f;
>   typedef Form_1_FunctionSpace_2 CoefficientSpace_g;
>
>   // Coefficients
>   dolfin::CoefficientAssigner f;
>   dolfin::CoefficientAssigner g;
> };
>
> // Class typedefs
> typedef Form_0 BilinearForm;
> typedef Form_1 LinearForm;
> typedef Form_0::TestSpace FunctionSpace;
>
> } // namespace Poisson
>
> // DOLFIN wrappers
>
> // Standard library includes
> #include <string>
>
> // DOLFIN includes
> #include <dolfin/common/NoDeleter.h>
> #include <dolfin/fem/FiniteElement.h>
> #include <dolfin/fem/DofMap.h>
> #include <dolfin/fem/Form.h>
> #include <dolfin/function/FunctionSpace.h>
> #include <dolfin/function/Function.h>
> #include <dolfin/function/GenericFunction.h>
> #include <dolfin/function/CoefficientAssigner.h>
>
> namespace Poisson
> {
>
> class Form_auxiliary_0_FunctionSpace_auxiliary_0: public dolfin::FunctionSpace
> {
> public:
>
>
>   Form_auxiliary_0_FunctionSpace_auxiliary_0(const dolfin::Mesh & mesh):
>       dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
>                             (boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
>                             (new poisson_auxiliary_dof_map_0()), mesh)))
>   {
>     // Do nothing
>   }
>
>   Form_auxiliary_0_FunctionSpace_auxiliary_0(dolfin::Mesh & mesh):
>     dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                           boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
>                           boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_auxiliary_dof_map_0()), mesh)))
>   {
>     // Do nothing
>   }
>
>   Form_auxiliary_0_FunctionSpace_auxiliary_0(boost::shared_ptr<dolfin::Mesh> mesh):
>       dolfin::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_auxiliary_dof_map_0()), *mesh)))
>   {
>       // Do nothing
>   }
>
>   Form_auxiliary_0_FunctionSpace_auxiliary_0(boost::shared_ptr<const dolfin::Mesh> mesh):
>       dolfin::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_auxiliary_dof_map_0()), *mesh)))
>   {
>       // Do nothing
>   }
>
>
>   ~Form_auxiliary_0_FunctionSpace_auxiliary_0()
>   {
>   }
>
> };
>
> class Form_auxiliary_0: public dolfin::Form
> {
> public:
>
>   // Constructor
>   Form_auxiliary_0(const dolfin::FunctionSpace& V0):
>     dolfin::Form(1, 0)
>   {
>     _function_spaces[0] = reference_to_no_delete_pointer(V0);
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_auxiliary_form_0());
>   }
>
>   // Constructor
>   Form_auxiliary_0(boost::shared_ptr<const dolfin::FunctionSpace> V0):
>     dolfin::Form(1, 0)
>   {
>     _function_spaces[0] = V0;
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_auxiliary_form_0());
>   }
>
>   // Destructor
>   ~Form_auxiliary_0()
>   {}
>
>   /// Return the number of the coefficient with this name
>   virtual dolfin::uint coefficient_number(const std::string& name) const
>   {
>     dolfin::error("No coefficients.");
>     return 0;
>   }
>
>   /// Return the name of the coefficient with this number
>   virtual std::string coefficient_name(dolfin::uint i) const
>   {
>     dolfin::error("No coefficients.");
>     return "unnamed";
>   }
>
>   // Typedefs
>   typedef Form_auxiliary_0_FunctionSpace_auxiliary_0 TestSpace;
>
>   // Coefficients
> };
>
> class Form_auxiliary_1_FunctionSpace_auxiliary_0: public dolfin::FunctionSpace
> {
> public:
>
>
>   Form_auxiliary_1_FunctionSpace_auxiliary_0(const dolfin::Mesh & mesh):
>       dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
>                             (boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
>                             (new poisson_auxiliary_dof_map_0()), mesh)))
>   {
>     // Do nothing
>   }
>
>   Form_auxiliary_1_FunctionSpace_auxiliary_0(dolfin::Mesh & mesh):
>     dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
>                           boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
>                           boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_auxiliary_dof_map_0()), mesh)))
>   {
>     // Do nothing
>   }
>
>   Form_auxiliary_1_FunctionSpace_auxiliary_0(boost::shared_ptr<dolfin::Mesh> mesh):
>       dolfin::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_auxiliary_dof_map_0()), *mesh)))
>   {
>       // Do nothing
>   }
>
>   Form_auxiliary_1_FunctionSpace_auxiliary_0(boost::shared_ptr<const dolfin::Mesh> mesh):
>       dolfin::FunctionSpace(mesh,
>                             boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
>                             boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_auxiliary_dof_map_0()), *mesh)))
>   {
>       // Do nothing
>   }
>
>
>   ~Form_auxiliary_1_FunctionSpace_auxiliary_0()
>   {
>   }
>
> };
>
> class Form_auxiliary_1: public dolfin::Form
> {
> public:
>
>   // Constructor
>   Form_auxiliary_1(const dolfin::FunctionSpace& V0):
>     dolfin::Form(1, 0)
>   {
>     _function_spaces[0] = reference_to_no_delete_pointer(V0);
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_auxiliary_form_1());
>   }
>
>   // Constructor
>   Form_auxiliary_1(boost::shared_ptr<const dolfin::FunctionSpace> V0):
>     dolfin::Form(1, 0)
>   {
>     _function_spaces[0] = V0;
>
>     _ufc_form = boost::shared_ptr<const ufc::form>(new poisson_auxiliary_form_1());
>   }
>
>   // Destructor
>   ~Form_auxiliary_1()
>   {}
>
>   /// Return the number of the coefficient with this name
>   virtual dolfin::uint coefficient_number(const std::string& name) const
>   {
>     dolfin::error("No coefficients.");
>     return 0;
>   }
>
>   /// Return the name of the coefficient with this number
>   virtual std::string coefficient_name(dolfin::uint i) const
>   {
>     dolfin::error("No coefficients.");
>     return "unnamed";
>   }
>
>   // Typedefs
>   typedef Form_auxiliary_1_FunctionSpace_auxiliary_0 TestSpace;
>
>   // Coefficients
> };
>
> // Class typedefs
> typedef Form_auxiliary_0::TestSpace FunctionSpace_auxiliary;
>
> } // namespace Poisson
>
> /// This class defines the interface for computing enriched function space
> /// and auxiliary function spaces
> ///
> //
>
> // Dolfin includes
> #include <dolfin/common/NoDeleter.h>
> #include <dolfin/mesh/Mesh.h>
> #include <dolfin/fem/DofMap.h>
>
> // PartitionOfUnity includes
> #include <pum/FunctionSpace.h>
> #include <pum/PUM.h>
> #include <pum/GenericSurface.h>
>
> namespace Poisson
> {
>   class FunctionSpaces
>   {
>     dolfin::Mesh& mesh;
>     std::vector<const pum::GenericSurface*>& surfaces;
>     std::vector<const pum::GenericPUM*> pums;
>
>
>     pum::PUM* pum_0;
>     dolfin::DofMap* dof_map_0;
>
>     /// Build PUM objects to intialize enriched FunctionSpace
>     void init()
>     {
>
>       // Create DofMap instance for the field 0
>       //dolfin::DofMap dof_map_0(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh);
>       dof_map_0 = new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh);
>
>       // Create PUM objects
>       pum_0 = new pum::PUM(surfaces, mesh, *dof_map_0);
>
>
>       // Push PUM objects to a vector
>       pums.push_back(pum_0);
>     }
>
>   public:
>
>     /// Constructor
>     FunctionSpaces(dolfin::Mesh& mesh, std::vector<const pum::GenericSurface*>& surfaces):
>                   mesh(mesh), surfaces(surfaces)
>     {
>       init();
>     }
>
>     /// Destructor
>     ~FunctionSpaces()
>     {
>
>       delete pum_0;
>       delete dof_map_0;
>     }
>
>     /// Return Enriched Function Space
>     pum::FunctionSpace space()
>     {
>        return FunctionSpace(mesh, pums);
>     }
>
>     /// Return Enriched Function Space
>     dolfin::FunctionSpace auxiliary_space()
>     {
>        return FunctionSpace_auxiliary(mesh);
>     }
>
>     /// Update PUM objects while surfaces evolve
>     void update()
>     {
>        for (unsigned int i = 0; i < pums.size(); ++i)
>            const_cast<pum::GenericPUM*>(pums[i])->update();
>     }
>
>
>   };
> }
>
> #endif

> # Copyright (C) 2008-2009 Mehdi Nikbakht and Garth N. Wells.
> # Licensed under the GNU GPL Version 3.0 or any later version.
> #
> # The bilinear form a(v, u) and linear form L(v) for
> # Poisson's equation with discontinuities.
> #
> # Compile this form with FFC: ffc-pum -l dolfin Poisson.ufl
> #
>
> elem_cont = FiniteElement("CG", triangle, 1)
> elem_discont = ElementRestriction(elem_cont, dc) # or ec[dc]
>
> #element = elem_cont * elem_discont
>
> #(vc, vd) = TestFunctions(element)
> #(uc, ud) = TrialFunctions(element)
>
> #v = vc + vd
> #u = uc + ud
>
> element = elem_cont + elem_discont
>
> v = TestFunction(element)
> u = TrialFunction(element)
>
> k  = Constant(triangle)
> f  = Coefficient(elem_cont)
> w  = Coefficient(elem_cont)
> g  = Coefficient(elem_cont)
>
> a = w*dot(grad(v), grad(u))*dx
> L = v*f*dx - v*g*ds
>

---

From: Mehdi Nikbakht <[email-masked]>
Date: Tue, 13 Apr 2010 13:23:12 +0200

> In the directory dolfin/quadrature/ in DOLFIN.

Ok, I got it now. I didn't have the latest dolfin.

>   1. Evaluate basis functions at given quadrature point

Which are not known on the compile time and you need to evaluate them
inside your tabulate_tensors using evaluate_basis* functions.

>   2. Evaluate integrand at each quadrature points

Note that evaluate_basis* functions are evaluated at the real
coordinates. Therefore if you want to use them to evaluate your
integrand, you need to rebuild the integrands without mappings. 

Mehdi

---

From: Anders Logg <[email-masked]>
Date: Tue, 13 Apr 2010 13:31:04 +0200

> inside your tabulate_tensors using evaluate_basis* functions.

Exactly.

> integrand, you need to rebuild the integrands without mappings.

I think we have discussed before the possibility of having two
versions of evalute_basis. One version that evaluates on the reference
element and one that evaluates on the physical element. It would be
practical to have both. Perhaps a third function should also be added
that maps a given point back to the reference element. Then
evaluate_basis on the physical element could call the other two.

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 13 Apr 2010 20:14:32 +0800

>        AK += w * evaluate_integrand(x)

OK.

> the sandbox):.

OK. We have a member function of the object with which we construct 
forms that determines the intersection

> the need for building complex abstractions at this point.

OK

Garth

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 13 Apr 2010 20:19:33 +0800

> libraries like CGAL. See BarycenterQuadrature.cpp.

What we plan is that the generated code would provide a 1D scheme (or 
possible just a polynomial order) and the C++ quadrature object could 
apply it on a sub-triangulations.

> need for writing a special-purpose code generator).

It's not a problem - just good to scrutinise heavily additions to UFC to 
avoid bloat. I'm happy for it to be added. We should make clear that 
it's for testing/special cases since the performance for computing 
element tensors will be poor compared to tabulate_tensor.

Garth

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 13 Apr 2010 20:22:00 +0800

On 13/04/10 19:31, Anders Logg wrote:
> practical to have both.

We should definitely have both.

Garth

> Perhaps a third function should also be added
>>>> #include<cmath>
>>>> #include<algorithm>
>>>> #include<stdexcept>
>>>> #include<fstream>
>>>> #include<boost/assign/list_of.hpp>
>>>> #include<ufc.h>
>>>> #include<pum/GenericPUM.h>
>>>>
>>>>      for (unsigned int r = 0; r<  1; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  2; r++)
>>>>      {
>>>>        for (unsigned int s = 0; s<  2 - r; s++)
>>>>        {
>>>>      for (unsigned int r = 0; r<  3; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  3; r++)
>>>>      {
>>>>        evaluate_basis(r,&dof_values, coordinates, c);
>>>>        values[r] = dof_values;
>>>>      }// end loop over 'r'
>>>>    }
>>>>
>>>>      for (unsigned int r = 0; r<  n; r++)
>>>>      {
>>>>        num_derivatives *= 2;
>>>>      }// end loop over 'r'
>>>>
>>>>      for (unsigned int row = 0; row<  num_derivatives; row++)
>>>>      {
>>>>        for (unsigned int col = 0; col<  n; col++)
>>>>      for (unsigned int row = 1; row<  num_derivatives; row++)
>>>>      {
>>>>        for (unsigned int num = 0; num<  row; num++)
>>>>        {
>>>>          for (unsigned int col = n-1; col+1>  0; col--)
>>>>          {
>>>>            if (combinations[row][col] + 1>  1)
>>>>      for (unsigned int j = 0; j<  num_derivatives; j++)
>>>>      {
>>>>        for (unsigned int k = 0; k<  num_derivatives; k++)
>>>>          transform[j][k] = 1;
>>>>      }
>>>>
>>>>      for (unsigned int row = 0; row<  num_derivatives; row++)
>>>>      {
>>>>        for (unsigned int col = 0; col<  num_derivatives; col++)
>>>>        {
>>>>          for (unsigned int k = 0; k<  n; k++)
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  1; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  2; r++)
>>>>      {
>>>>        for (unsigned int s = 0; s<  2 - r; s++)
>>>>        {
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>        for (unsigned int t = 0; t<  3; t++)
>>>>        {
>>>>          for (unsigned int u = 0; u<  3; u++)
>>>>          {
>>>>        for (unsigned int s = 0; s<  n; s++)
>>>>        {
>>>>          for (unsigned int t = 0; t<  3; t++)
>>>>          {
>>>>            for (unsigned int u = 0; u<  3; u++)
>>>>            {
>>>>          for (unsigned int t = 0; t<  3; t++)
>>>>          {
>>>>            for (unsigned int u = 0; u<  3; u++)
>>>>            {
>>>>              for (unsigned int tu = 0; tu<  3; tu++)
>>>>              {
>>>>          for (unsigned int t = 0; t<  3; t++)
>>>>          {
>>>>            for (unsigned int u = 0; u<  3; u++)
>>>>            {
>>>>              for (unsigned int tu = 0; tu<  3; tu++)
>>>>              {
>>>>        for (unsigned int s = 0; s<  3; s++)
>>>>        {
>>>>          for (unsigned int t = 0; t<  3; t++)
>>>>          {
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>        for (unsigned int s = 0; s<  num_derivatives; s++)
>>>>        {
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>        delete [] transform[r];
>>>>      }// end loop over 'r'
>>>>      delete [] transform;
>>>>    }
>>>>
>>>>      for (unsigned int r = 0; r<  n; r++)
>>>>      {
>>>>        num_derivatives *= 2;
>>>>      }// end loop over 'r'
>>>>
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  3; r++)
>>>>      {
>>>>        for (unsigned int s = 0; s<  num_derivatives; s++)
>>>>        {
>>>>      if (0<= i&&  i<= 2)
>>>>        for (unsigned int r = 0; r<  1; r++)
>>>>        {
>>>>        for (unsigned int r = 0; r<  2; r++)
>>>>        {
>>>>          for (unsigned int s = 0; s<  2 - r; s++)
>>>>          {
>>>>        for (unsigned int r = 0; r<  3; r++)
>>>>        {
>>>>      if (3<= i&&  i<= 5)
>>>>        for (unsigned int r = 0; r<  1; r++)
>>>>        {
>>>>        for (unsigned int r = 0; r<  2; r++)
>>>>        {
>>>>          for (unsigned int s = 0; s<  2 - r; s++)
>>>>          {
>>>>        for (unsigned int r = 0; r<  3; r++)
>>>>        {
>>>>      for (unsigned int r = 0; r<  6; r++)
>>>>      {
>>>>        for (unsigned int s = 0; s<  2; s++)
>>>>        {
>>>>      for (unsigned int r = 0; r<  n; r++)
>>>>      {
>>>>        num_derivatives *= 2;
>>>>      }// end loop over 'r'
>>>>
>>>>      for (unsigned int row = 0; row<  num_derivatives; row++)
>>>>      {
>>>>        for (unsigned int col = 0; col<  n; col++)
>>>>      for (unsigned int row = 1; row<  num_derivatives; row++)
>>>>      {
>>>>        for (unsigned int num = 0; num<  row; num++)
>>>>        {
>>>>          for (unsigned int col = n-1; col+1>  0; col--)
>>>>          {
>>>>            if (combinations[row][col] + 1>  1)
>>>>      for (unsigned int j = 0; j<  num_derivatives; j++)
>>>>      {
>>>>        for (unsigned int k = 0; k<  num_derivatives; k++)
>>>>          transform[j][k] = 1;
>>>>      }
>>>>
>>>>      for (unsigned int row = 0; row<  num_derivatives; row++)
>>>>      {
>>>>        for (unsigned int col = 0; col<  num_derivatives; col++)
>>>>        {
>>>>          for (unsigned int k = 0; k<  n; k++)
>>>>      for (unsigned int r = 0; r<  2*num_derivatives; r++)
>>>>      {
>>>>      if (0<= i&&  i<= 2)
>>>>        for (unsigned int r = 0; r<  1; r++)
>>>>        {
>>>>        for (unsigned int r = 0; r<  2; r++)
>>>>        {
>>>>          for (unsigned int s = 0; s<  2 - r; s++)
>>>>          {
>>>>        for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>        {
>>>>        for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>        {
>>>>          for (unsigned int t = 0; t<  3; t++)
>>>>          {
>>>>            for (unsigned int u = 0; u<  3; u++)
>>>>            {
>>>>          for (unsigned int s = 0; s<  n; s++)
>>>>          {
>>>>            for (unsigned int t = 0; t<  3; t++)
>>>>            {
>>>>              for (unsigned int u = 0; u<  3; u++)
>>>>              {
>>>>            for (unsigned int t = 0; t<  3; t++)
>>>>            {
>>>>              for (unsigned int u = 0; u<  3; u++)
>>>>              {
>>>>                for (unsigned int tu = 0; tu<  3; tu++)
>>>>                {
>>>>            for (unsigned int t = 0; t<  3; t++)
>>>>            {
>>>>              for (unsigned int u = 0; u<  3; u++)
>>>>              {
>>>>                for (unsigned int tu = 0; tu<  3; tu++)
>>>>                {
>>>>          for (unsigned int s = 0; s<  3; s++)
>>>>          {
>>>>            for (unsigned int t = 0; t<  3; t++)
>>>>            {
>>>>        for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>        {
>>>>          for (unsigned int s = 0; s<  num_derivatives; s++)
>>>>          {
>>>>        for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>        {
>>>>        for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>        {
>>>>          delete [] transform[r];
>>>>        }// end loop over 'r'
>>>>        delete [] transform;
>>>>      }
>>>>
>>>>      if (3<= i&&  i<= 5)
>>>>      {
>>>>        for (unsigned int r = 0; r<  1; r++)
>>>>        {
>>>>        for (unsigned int r = 0; r<  2; r++)
>>>>        {
>>>>          for (unsigned int s = 0; s<  2 - r; s++)
>>>>          {
>>>>        for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>        {
>>>>        for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>        {
>>>>          for (unsigned int t = 0; t<  3; t++)
>>>>          {
>>>>            for (unsigned int u = 0; u<  3; u++)
>>>>            {
>>>>          for (unsigned int s = 0; s<  n; s++)
>>>>          {
>>>>            for (unsigned int t = 0; t<  3; t++)
>>>>            {
>>>>              for (unsigned int u = 0; u<  3; u++)
>>>>              {
>>>>            for (unsigned int t = 0; t<  3; t++)
>>>>            {
>>>>              for (unsigned int u = 0; u<  3; u++)
>>>>              {
>>>>                for (unsigned int tu = 0; tu<  3; tu++)
>>>>                {
>>>>            for (unsigned int t = 0; t<  3; t++)
>>>>            {
>>>>              for (unsigned int u = 0; u<  3; u++)
>>>>              {
>>>>                for (unsigned int tu = 0; tu<  3; tu++)
>>>>                {
>>>>          for (unsigned int s = 0; s<  3; s++)
>>>>          {
>>>>            for (unsigned int t = 0; t<  3; t++)
>>>>            {
>>>>        for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>        {
>>>>          for (unsigned int s = 0; s<  num_derivatives; s++)
>>>>          {
>>>>        for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>        {
>>>>        for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>        {
>>>>          delete [] transform[r];
>>>>        }// end loop over 'r'
>>>>        delete [] transform;
>>>>      }
>>>>
>>>>    }
>>>>
>>>>      for (unsigned int r = 0; r<  n; r++)
>>>>      {
>>>>        num_derivatives *= 2;
>>>>      }// end loop over 'r'
>>>>
>>>>      for (unsigned int r = 0; r<  2*num_derivatives; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  6; r++)
>>>>      {
>>>>        for (unsigned int s = 0; s<  2*num_derivatives; s++)
>>>>        {
>>>>      if (d>  2)
>>>>          if (i>  2)
>>>>      if (d>  2)
>>>>          if (i>  2)
>>>>    mutable std::vector<double>  Aa;
>>>>
>>>>      for (unsigned int r = 0; r<  num_entries; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  3; r++)
>>>>      {
>>>>        F0 += FE1[0][r]*w[0][r];
>>>>      }// end loop over 'r'
>>>>      unsigned int m = 0;
>>>>
>>>>      for (unsigned int j = 0; j<  6; j++)
>>>>      {
>>>>        for (unsigned int k = 0; k<  6; k++)
>>>>        {
>>>>          if (((((0<= j&&  j<  3))&&  ((0<= k&&  k<  3)))))
>>>>          {
>>>>        if  ((((0<= j&&  j<  3))))
>>>>        for (unsigned int i = 0; i<  1; i++)
>>>>        {
>>>>          Wn1[i] = W1[i];
>>>>          for (unsigned int j = 0; j<  2; j++)
>>>>      for (unsigned int ip = 0; ip<  Wn1.size(); ip++)
>>>>      {
>>>>          // Creating a table to save the values of derivatives order 1 at the current guass point for<<CG1 on a<triangle of degree 1>>  +<<CG1 on a<triangle of degree 1>>>|_{dc0}>
>>>>          double value_0[2];
>>>>          double table_1_D1[6][2];
>>>>          for (unsigned int j = 0; j<  6; j++)
>>>>          {
>>>>            for (unsigned int k = 0; k<  2; k++)
>>>>          // Creating a table to save the values of shape functions at the current guass point for<CG1 on a<triangle of degree 1>>
>>>>          double value_1[1];
>>>>          double table_0_D0[3][1];
>>>>          for (unsigned int j = 0; j<  3; j++)
>>>>          {
>>>>            for (unsigned int k = 0; k<  1; k++)
>>>>          // Creating a table to save the values of shape functions at the current guass point for<<CG1 on a<triangle of degree 1>>  +<<CG1 on a<triangle of degree 1>>>|_{dc0}>
>>>>          double value_2[1];
>>>>          double table_1_D0[6][1];
>>>>          for (unsigned int j = 0; j<  6; j++)
>>>>          {
>>>>            for (unsigned int k = 0; k<  1; k++)
>>>>        for (unsigned int r = 0; r<  3; r++)
>>>>        {
>>>>          for (unsigned int j = 0; j<  6; j++)
>>>>          {
>>>>            for (unsigned int k = 0; k<  6; k++)
>>>>            {
>>>>              if (!(((0<= j&&  j<  3))&&  ((0<= k&&  k<  3))))
>>>>              {
>>>>                if ((0<= j&&  j<  3)&&  (3<= k&&  k<  6))
>>>>                {
>>>>                  m = j;
>>>>                  n = k;
>>>>                }
>>>>                else if ((3<= j&&  j<  6)&&  (0<= k&&  k<  3))
>>>>                {
>>>>                  m = j;
>>>>                  n = k;
>>>>                }
>>>>                else if ((3<= j&&  j<  6)&&  (3<= k&&  k<  6))
>>>>                {
>>>>                  m = j;
>>>>                  n = k;
>>>>                }
>>>>
>>>>      for (unsigned int j = 0; j<  6; j++)
>>>>        for (unsigned int k = 0; k<  6; k++)
>>>>          if ((0<= j&&  j<  3)&&  (0<= k&&  k<  3))
>>>>            ++m;
>>>>          else
>>>>          {
>>>>              if (((0<= j&&  j<  3)) || ((0<= k&&  k<  3)) || (it_0_0 != active_dofs_0.end()&&  it_0_1 != active_dofs_0.end()))
>>>>              {
>>>>                A[m] = Aa[j*6 + k];
>>>>                ++m;
>>>>              }
>>>>          }
>>>>    }
>>>>
>>>> };
>>>>
>>>>    mutable std::vector<double>  Aa;
>>>>
>>>>      for (unsigned int r = 0; r<  num_entries; r++)
>>>>      {
>>>>      for (unsigned int ip = 0; ip<  4; ip++)
>>>>      {
>>>>        for (unsigned int r = 0; r<  3; r++)
>>>>        {
>>>>        for (unsigned int j = 0; j<  6; j++)
>>>>        {
>>>>          if (((((0<= j&&  j<  3)))))
>>>>          {
>>>>        for (unsigned int i = 0; i<  4; i++)
>>>>        {
>>>>          Wn4[i] = W4[i];
>>>>          for (unsigned int j = 0; j<  2; j++)
>>>>      for (unsigned int ip = 0; ip<  Wn4.size(); ip++)
>>>>      {
>>>>          // Creating a table to save the values of shape functions at the current guass point for<CG1 on a<triangle of degree 1>>
>>>>          double value_0[1];
>>>>          double table_0_D0[3][1];
>>>>          for (unsigned int j = 0; j<  3; j++)
>>>>          {
>>>>            for (unsigned int k = 0; k<  1; k++)
>>>>          // Creating a table to save the values of shape functions at the current guass point for<<CG1 on a<triangle of degree 1>>  +<<CG1 on a<triangle of degree 1>>>|_{dc0}>
>>>>          double value_1[1];
>>>>          double table_1_D0[6][1];
>>>>          for (unsigned int j = 0; j<  6; j++)
>>>>          {
>>>>            for (unsigned int k = 0; k<  1; k++)
>>>>        for (unsigned int r = 0; r<  3; r++)
>>>>        {
>>>>          for (unsigned int j = 0; j<  6; j++)
>>>>          {
>>>>            if (!(((0<= j&&  j<  3))))
>>>>            {
>>>>              if ((3<= j&&  j<  6))
>>>>      for (unsigned int j = 0; j<  6; j++)
>>>>        if ((0<= j&&  j<  3))
>>>>          ++m;
>>>>        else
>>>>        {
>>>>    mutable std::vector<double>  Aa;
>>>>
>>>>      for (unsigned int r = 0; r<  num_entries; r++)
>>>>      {
>>>>        for (unsigned int ip = 0; ip<  2; ip++)
>>>>        {
>>>>          for (unsigned int r = 0; r<  3; r++)
>>>>          {
>>>>          for (unsigned int j = 0; j<  6; j++)
>>>>          {
>>>>            if (((((0<= j&&  j<  3)))))
>>>>            {
>>>>        for (unsigned int ip = 0; ip<  2; ip++)
>>>>        {
>>>>          for (unsigned int r = 0; r<  3; r++)
>>>>          {
>>>>          for (unsigned int j = 0; j<  6; j++)
>>>>          {
>>>>            if (((((0<= j&&  j<  3)))))
>>>>            {
>>>>        for (unsigned int ip = 0; ip<  2; ip++)
>>>>        {
>>>>          for (unsigned int r = 0; r<  3; r++)
>>>>          {
>>>>          for (unsigned int j = 0; j<  6; j++)
>>>>          {
>>>>            if (((((0<= j&&  j<  3)))))
>>>>            {
>>>>        for (unsigned int i = 0; i<  2; i++)
>>>>        {
>>>>          Wn2[i] = W2[i];
>>>>          for (unsigned int j = 0; j<  2; j++)
>>>>      for (unsigned int ip = 0; ip<  Wn2.size(); ip++)
>>>>      {
>>>>          // Creating a table to save the values of shape functions at the current guass point for<CG1 on a<triangle of degree 1>>
>>>>          double value_0[1];
>>>>          double table_0_D0[3][1];
>>>>          for (unsigned int j = 0; j<  3; j++)
>>>>          {
>>>>            for (unsigned int k = 0; k<  1; k++)
>>>>          // Creating a table to save the values of shape functions at the current guass point for<<CG1 on a<triangle of degree 1>>  +<<CG1 on a<triangle of degree 1>>>|_{dc0}>
>>>>          double value_1[1];
>>>>          double table_1_D0[6][1];
>>>>          for (unsigned int j = 0; j<  6; j++)
>>>>          {
>>>>            for (unsigned int k = 0; k<  1; k++)
>>>>        for (unsigned int r = 0; r<  3; r++)
>>>>        {
>>>>          for (unsigned int j = 0; j<  6; j++)
>>>>          {
>>>>            if (!(((0<= j&&  j<  3))))
>>>>            {
>>>>              if ((3<= j&&  j<  6))
>>>>      for (unsigned int j = 0; j<  6; j++)
>>>>        if ((0<= j&&  j<  3))
>>>>          ++m;
>>>>        else
>>>>        {
>>>>      for (unsigned int r = 0; r<  1; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  2; r++)
>>>>      {
>>>>        for (unsigned int s = 0; s<  2 - r; s++)
>>>>        {
>>>>      for (unsigned int r = 0; r<  3; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  3; r++)
>>>>      {
>>>>        evaluate_basis(r,&dof_values, coordinates, c);
>>>>        values[r] = dof_values;
>>>>      }// end loop over 'r'
>>>>    }
>>>>
>>>>      for (unsigned int r = 0; r<  n; r++)
>>>>      {
>>>>        num_derivatives *= 2;
>>>>      }// end loop over 'r'
>>>>
>>>>      for (unsigned int row = 0; row<  num_derivatives; row++)
>>>>      {
>>>>        for (unsigned int col = 0; col<  n; col++)
>>>>      for (unsigned int row = 1; row<  num_derivatives; row++)
>>>>      {
>>>>        for (unsigned int num = 0; num<  row; num++)
>>>>        {
>>>>          for (unsigned int col = n-1; col+1>  0; col--)
>>>>          {
>>>>            if (combinations[row][col] + 1>  1)
>>>>      for (unsigned int j = 0; j<  num_derivatives; j++)
>>>>      {
>>>>        for (unsigned int k = 0; k<  num_derivatives; k++)
>>>>          transform[j][k] = 1;
>>>>      }
>>>>
>>>>      for (unsigned int row = 0; row<  num_derivatives; row++)
>>>>      {
>>>>        for (unsigned int col = 0; col<  num_derivatives; col++)
>>>>        {
>>>>          for (unsigned int k = 0; k<  n; k++)
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  1; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  2; r++)
>>>>      {
>>>>        for (unsigned int s = 0; s<  2 - r; s++)
>>>>        {
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>        for (unsigned int t = 0; t<  3; t++)
>>>>        {
>>>>          for (unsigned int u = 0; u<  3; u++)
>>>>          {
>>>>        for (unsigned int s = 0; s<  n; s++)
>>>>        {
>>>>          for (unsigned int t = 0; t<  3; t++)
>>>>          {
>>>>            for (unsigned int u = 0; u<  3; u++)
>>>>            {
>>>>          for (unsigned int t = 0; t<  3; t++)
>>>>          {
>>>>            for (unsigned int u = 0; u<  3; u++)
>>>>            {
>>>>              for (unsigned int tu = 0; tu<  3; tu++)
>>>>              {
>>>>          for (unsigned int t = 0; t<  3; t++)
>>>>          {
>>>>            for (unsigned int u = 0; u<  3; u++)
>>>>            {
>>>>              for (unsigned int tu = 0; tu<  3; tu++)
>>>>              {
>>>>        for (unsigned int s = 0; s<  3; s++)
>>>>        {
>>>>          for (unsigned int t = 0; t<  3; t++)
>>>>          {
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>        for (unsigned int s = 0; s<  num_derivatives; s++)
>>>>        {
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>        delete [] transform[r];
>>>>      }// end loop over 'r'
>>>>      delete [] transform;
>>>>    }
>>>>
>>>>      for (unsigned int r = 0; r<  n; r++)
>>>>      {
>>>>        num_derivatives *= 2;
>>>>      }// end loop over 'r'
>>>>
>>>>      for (unsigned int r = 0; r<  num_derivatives; r++)
>>>>      {
>>>>      for (unsigned int r = 0; r<  3; r++)
>>>>      {
>>>>        for (unsigned int s = 0; s<  num_derivatives; s++)
>>>>        {
>>>>      if (d>  2)
>>>>          if (i>  2)
>>>> #include<dolfin/common/NoDeleter.h>
>>>> #include<dolfin/mesh/Mesh.h>
>>>> #include<dolfin/fem/DofMap.h>
>>>> #include<dolfin/la/GenericVector.h>
>>>>
>>>> #include<pum/PostProcess.h>
>>>> #include<pum/FunctionSpace.h>
>>>> #include<pum/PUM.h>
>>>> #include<pum/GenericSurface.h>
>>>>
>>>> namespace Poisson
>>>> {
>>>>        for (unsigned int i = 0; i<  num_standard_dofs_0; ++i)
>>>>          positions_0[i] = i;
>>>>
>>>>        std::vector<std::vector<unsigned int>  >  enhanced_dof_maps_0;
>>>>            x.get(&value, 1,&pos_n);
>>>>            x0.add(&value, 1,&pos);
>>>>        for (unsigned int i = 0; i<  num_standard_dofs_0; ++i)
>>>>          positions_0[i] = i;
>>>>
>>>>        std::vector<std::vector<unsigned int>  >  enhanced_dof_maps_0;
>>>>            x.get(&value, 1,&pos_n);
>>>>            j.set(&value, 1,&pos);
>>>> #include<string>
>>>> #include<dolfin/common/NoDeleter.h>
>>>> #include<dolfin/fem/FiniteElement.h>
>>>> #include<dolfin/fem/DofMap.h>
>>>> #include<dolfin/fem/Form.h>
>>>> #include<dolfin/function/FunctionSpace.h>
>>>> #include<dolfin/function/Function.h>
>>>> #include<dolfin/function/GenericFunction.h>
>>>> #include<dolfin/function/CoefficientAssigner.h>
>>>>
>>>> namespace Poisson
>>>> {
>>>>
>>>>    CoefficientSpace_f(const dolfin::Mesh&  mesh):
>>>>    CoefficientSpace_f(dolfin::Mesh&  mesh):
>>>>    CoefficientSpace_g(const dolfin::Mesh&  mesh):
>>>>    CoefficientSpace_g(dolfin::Mesh&  mesh):
>>>>    CoefficientSpace_w(const dolfin::Mesh&  mesh):
>>>>    CoefficientSpace_w(dolfin::Mesh&  mesh):
>>>>    Form_0_FunctionSpace_0(const dolfin::Mesh&  mesh ,const std::vector<const pum::GenericPUM*>&  pums):
>>>>    Form_0_FunctionSpace_0(dolfin::Mesh&  mesh ,const std::vector<const pum::GenericPUM*>&  pums):
>>>>    Form_0_FunctionSpace_1(const dolfin::Mesh&  mesh ,const std::vector<const pum::GenericPUM*>&  pums):
>>>>    Form_0_FunctionSpace_1(dolfin::Mesh&  mesh ,const std::vector<const pum::GenericPUM*>&  pums):
>>>>    Form_0(const pum::FunctionSpace&  V0, const pum::FunctionSpace&  V1, dolfin::GenericFunction&  w):
>>>>    Form_0(boost::shared_ptr<const pum::FunctionSpace>  V0, boost::shared_ptr<const pum::FunctionSpace>  V1, dolfin::GenericFunction&  w):
>>>>    Form_1_FunctionSpace_0(const dolfin::Mesh&  mesh ,const std::vector<const pum::GenericPUM*>&  pums):
>>>>    Form_1_FunctionSpace_0(dolfin::Mesh&  mesh ,const std::vector<const pum::GenericPUM*>&  pums):
>>>>    Form_1(const pum::FunctionSpace&  V0, dolfin::GenericFunction&  f, dolfin::GenericFunction&  g):
>>>>    Form_1(boost::shared_ptr<const pum::FunctionSpace>  V0, dolfin::GenericFunction&  f, dolfin::GenericFunction&  g):
>>>> #include<string>
>>>> #include<dolfin/common/NoDeleter.h>
>>>> #include<dolfin/fem/FiniteElement.h>
>>>> #include<dolfin/fem/DofMap.h>
>>>> #include<dolfin/fem/Form.h>
>>>> #include<dolfin/function/FunctionSpace.h>
>>>> #include<dolfin/function/Function.h>
>>>> #include<dolfin/function/GenericFunction.h>
>>>> #include<dolfin/function/CoefficientAssigner.h>
>>>>
>>>> namespace Poisson
>>>> {
>>>>
>>>>    Form_auxiliary_0_FunctionSpace_auxiliary_0(const dolfin::Mesh&  mesh):
>>>>    Form_auxiliary_0_FunctionSpace_auxiliary_0(dolfin::Mesh&  mesh):
>>>>    Form_auxiliary_1_FunctionSpace_auxiliary_0(const dolfin::Mesh&  mesh):
>>>>    Form_auxiliary_1_FunctionSpace_auxiliary_0(dolfin::Mesh&  mesh):
>>>> #include<dolfin/common/NoDeleter.h>
>>>> #include<dolfin/mesh/Mesh.h>
>>>> #include<dolfin/fem/DofMap.h>
>>>>
>>>> #include<pum/FunctionSpace.h>
>>>> #include<pum/PUM.h>
>>>> #include<pum/GenericSurface.h>
>>>>
>>>> namespace Poisson
>>>> {
>>>>    class FunctionSpaces
>>>>    {
>>>>      dolfin::Mesh&  mesh;
>>>>         for (unsigned int i = 0; i<  pums.size(); ++i)

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 13 Apr 2010 20:27:11 +0800

>                                    const double* coordinates) const = 0;

What about

     virtual void tabulate_tensor(double* A,
                          const double * const * w,
                          const cell&  c,
                          uint num_quad_points,
                          const double * const * coordinates,
                          const double * weights) const = 0;

instead? This would allow more of the quadrature code optimisations to 
be employed.

Garth

---

From: Anders Logg <[email-masked]>
Date: Tue, 13 Apr 2010 14:53:35 +0200

On Tue, Apr 13, 2010 at 08:27:11PM +0800, Garth N. Wells wrote:
> instead? This would allow more of the quadrature code optimisations
> to be employed.

I thought about this but wasn't sure it was general enough for our
needs. I'll need to think about it.

--
Anders

---

From: Anders Logg <[email-masked]>
Date: Tue, 13 Apr 2010 15:00:54 +0200

On Tue, Apr 13, 2010 at 08:19:33PM +0800, Garth N. Wells wrote:
> What we plan is that the generated code would provide a 1D scheme
> (or possible just a polynomial order) and the C++ quadrature object
> could apply it on a sub-triangulations.

I don't understand. How would the 1D scheme connect to a
triangulation?

Btw, we have considered subtriangulating the polyhedra to get a
quadrature rule, but decided to just to barycentric quadrature for now
to save time, both implementation and run time...

> It's not a problem - just good to scrutinise heavily additions to

Agree.

> UFC to avoid bloat. I'm happy for it to be added. We should make
> clear that it's for testing/special cases since the performance for
> computing element tensors will be poor compared to tabulate_tensor.

ok. So we have the following proposed additions to UFC:

  evaluate_integrand
  tabulate_tensor(x, w)
  evaluate_basis (two versions)

  + quite a few other suggestions as blueprints
  https://blueprints.launchpad.net/ufc

Perhaps it's a good time to think about UFC 2.0 and take care of all
the blueprints. Or perhaps we should wait until after DOLFIN 1.0?

I have limited resources at the moment so we might want to postpone
this.

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 13 Apr 2010 21:06:35 +0800

> I don't understand. How would the 1D scheme connect to a

Construct a rule for triangles from the 1D scheme (like FIAT does) and 
apply it on each sub-triangle.

> to save time, both implementation and run time...

We sub-triangulate. I'm hoping that some of the work can be moved to CGAL.

> the blueprints.

I think that we're settled on evaluate_basis, so this can go into 
UFC-dev. The other needs more discussion (and more work), and I won't 
have time for a while.

Garth

Or perhaps we should wait until after DOLFIN 1.0?

---

From: Anders Logg <[email-masked]>
Date: Tue, 13 Apr 2010 15:13:14 +0200

On Tue, Apr 13, 2010 at 09:06:35PM +0800, Garth N. Wells wrote:
> Construct a rule for triangles from the 1D scheme (like FIAT does)
> and apply it on each sub-triangle.

ok. Then I think I understand. So the input to your specialized
version of tabulate_tensor is a list of triangles?

> We sub-triangulate. I'm hoping that some of the work can be moved to CGAL.

I think CGAL can handle this but we decided to skip it for the
moment. It should be possible to make it fast, especially since it is
completely parallel (no need to match sub triangulations on
neighboring triangles).

> UFC-dev. The other needs more discussion (and more work), and I
> won't have time for a while.

Have we decided on the name for the second version of evaluate_basis?

Should it be evaluate_basis_reference (as hinted in one of the
blueprints).

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 13 Apr 2010 21:22:39 +0800

> version of tabulate_tensor is a list of triangles?

No, tabulate_tensor for a given cell first asks if that cell is 
intersected, and if it is intersected it then calls a function which 
supplies the quadrature scheme. An argument to this function is a 
standard quadrature scheme (which is sufficient to integrate the 
functions on a normal cell). The FFC generated code is unaware of how 
the quadrature scheme for intersected cells is computed - that's decided 
on the the code which we build on top of DOLFIN. It's in the C++ code 
which we have written that the sub-triangulation of a cell on either 
side of the intersecting surface is performed.

We're writing this up for the FEniCS book. Some of it is in

   http://dx.doi.org/10.3390/a2031008

> completely parallel (no need to match sub triangulations on

Not only fast - sub-triangulation gets tricky in 3D.

> Should it be evaluate_basis_reference (as hinted in one of the

Seems OK at first glance.

Have you added 'evaluate_integrand' and 'tabulate_tensor(x, w)' as 
Blueprints?

Garth

---

From: Mehdi Nikbakht <[email-masked]>
Date: Tue, 13 Apr 2010 10:31:51 +0200

Mehdi

---

From: Mehdi Nikbakht <[email-masked]>
Date: Tue, 13 Apr 2010 12:51:29 +0200

Mehdi

---

From: Mehdi Nikbakht <[email-masked]>
Date: Tue, 13 Apr 2010 13:23:12 +0200

Mehdi

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 13 Apr 2010 21:40:05 +0200 (CET)

> evaluate_integrand (if it is allowed to call tabaulate_basis).

That's how I think it should be done. How could we otherwise get the values of the basis functions? Adding code to evaluate
the basis functions inside the evaluate_integrand function would be a bit silly unless efficiency proves this to be necessary. We can let the relevant elements be private members of the integral class.

> need for writing a special-purpose code generator).

Thanks to the new compiler design, we just need to add a new representation and and code generator (UFL.Transformer), both should be relative simple to add.

Kristian

> --
> Anders
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAkvEQGkACgkQTuwUCDsYZdGGmgCgkh0D4rEKk24cirvBJ+NX5I+v
> YqgAoIWYa3wX9iwJ/82346+AsL/gwcRx
> =wW/z
> -----END PGP SIGNATURE-----
>

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 13 Apr 2010 21:52:31 +0200 (CET)

>> integrand, you need to rebuild the integrands without mappings.

Why should the basis functions not be evaluated on the physical element?
If evaluate_integrand() gets cell& c as an argument why not just compute the basis function directly and use it when evaluating the UFL.integrand? Then only very little translation is needed by the transformer.

> evaluate_basis on the physical element could call the other two.

That would be one possibility, but we have avoided doing things this due to fear of overhead.
Either way, it should not be too difficult to have both.

Kristian

>> > >         // Total nu
> ...
> [Message clipped]  
> iEYEARECAAYFAkvEVfYACgkQTuwUCDsYZdF6RwCfZ29ThwiA20YFxxmi6Gs3TOU4
> gJQAni1WqXPaG4o6oXajLOA8hQwpFE+a
> =TdIb

---

From: Mehdi Nikbakht <[email-masked]>
Date: Tue, 13 Apr 2010 10:31:51 +0200

Mehdi

---

From: Mehdi Nikbakht <[email-masked]>
Date: Tue, 13 Apr 2010 12:51:29 +0200

Mehdi

---

From: Mehdi Nikbakht <[email-masked]>
Date: Tue, 13 Apr 2010 13:23:12 +0200

Mehdi

---

From: Anders Logg <[email-masked]>
Date: Thu, 15 Apr 2010 20:44:03 +0200

On Tue, Apr 13, 2010 at 09:22:39PM +0800, Garth N. Wells wrote:
> functions on a normal cell). The FFC generated code is unaware of
> how the quadrature scheme for intersected cells is computed - that's
> decided on the the code which we build on top of DOLFIN. It's in the
> C++ code which we have written that the sub-triangulation of a cell
> on either side of the intersecting surface is performed.
>

I have now:

  https://blueprints.launchpad.net/ufc/+spec/evaluate-integrand

--
Anders