ffc team mailing list archive
-
ffc team
-
Mailing list archive
-
Message #03843
Re: evaluate_integrand
On Tue, 2010-04-13 at 13:10 +0200, Anders Logg wrote:
> On Tue, Apr 13, 2010 at 12:51:29PM +0200, Mehdi Nikbakht wrote:
> > On Tue, 2010-04-13 at 11:59 +0200, Anders Logg wrote:
> > > 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:
> > > > > On Tue, Apr 13, 2010 at 10:31:51AM +0200, Mehdi Nikbakht wrote:
> > > > > > On Tue, 2010-04-13 at 09:45 +0200, Anders Logg wrote:
> > > > > > > On Tue, Apr 13, 2010 at 07:17:28AM +0800, Garth N. Wells wrote:
> > > > > > > >
> > > > > > > >
> > > > > > > > On 12/04/10 23:35, Anders Logg wrote:
> > > > > > > > >On Mon, Apr 12, 2010 at 10:20:13PM +0800, Garth N. Wells wrote:
> > > > > > > > >>
> > > > > > > > >>
> > > > > > > > >>On 12/04/10 21:49, Anders Logg wrote:
> > > > > > > > >>>On Mon, Apr 12, 2010 at 09:34:38PM +0800, Garth N. Wells wrote:
> > > > > > > > >>>>
> > > > > > > > >>>>
> > > > > > > > >>>>On 12/04/10 21:29, Anders Logg wrote:
> > > > > > > > >>>>>On Mon, Apr 12, 2010 at 09:21:32PM +0800, Garth N. Wells wrote:
> > > > > > > > >>>>>>
> > > > > > > > >>>>>>
> > > > > > > > >>>>>>On 12/04/10 21:19, Garth N. Wells wrote:
> > > > > > > > >>>>>>>
> > > > > > > > >>>>>>>
> > > > > > > > >>>>>>>On 12/04/10 20:47, Anders Logg wrote:
> > > > > > > > >>>>>>>>We are doing some work where we need to do run-time quadrature over
> > > > > > > > >>>>>>>>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?
> > > > > > > > >>>>>>>
> > > > > > > > >>>>>>
> > > > > > > > >>>>>>I meant 'tabulate_tensor'.
> > > > > > > > >>>>>
> > > > > > > > >>>>>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.
> > > > > > > > >>>
> > > > > > > > >>>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.
> > > > > > > > >
> > > > > > > > >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).
> > > > > > > > >
> > > > > > > > >>>Having
> > > > > > > > >>>evaluate_integrand would allow more flexibility for users to implement
> > > > > > > > >>>their own special quadrature scheme.
> > > > > > > > >>>
> > > > > > > > >>
> > > > > > > > >>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.
> > > > > > > > >
> > > > > > > > >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.
> > > > > > > > >
> > > > > > > >
> > > > > > > > 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.
> > > > > >
> > > > > > 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?
> > > > >
> > > > > 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.
> > > >
> > > > 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.
> >
> > Where can I find this file?
>
> In the directory dolfin/quadrature/ in DOLFIN.
Ok, I got it now. I didn't have the latest dolfin.
>
> > > > > > 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?
> > > > >
> > > > > 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.
> > > >
> > > > 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).
> >
> > 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.
>
> 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
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
>
> > Please have a look on the tabulate_Tensor function inside generated code
> > 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 double * const * w,
> > 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 double * const * w,
> > 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 double * const * w,
> > 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 double * const * w,
> > 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 double * const * w,
> > 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 double * const * w,
> > 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
> >
>
Follow ups
References
-
Re: evaluate_integrand
From: Anders Logg, 2010-04-12
-
Re: evaluate_integrand
From: Garth N. Wells, 2010-04-12
-
Re: evaluate_integrand
From: Anders Logg, 2010-04-12
-
Re: evaluate_integrand
From: Garth N. Wells, 2010-04-12
-
Re: evaluate_integrand
From: Anders Logg, 2010-04-13
-
Re: evaluate_integrand
From: Mehdi Nikbakht, 2010-04-13
-
Re: evaluate_integrand
From: Anders Logg, 2010-04-13
-
Re: evaluate_integrand
From: Anders Logg, 2010-04-13
-
Re: evaluate_integrand
From: Mehdi Nikbakht, 2010-04-13
-
Re: evaluate_integrand
From: Anders Logg, 2010-04-13
