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