ffc team mailing list archive
-
ffc team
-
Mailing list archive
-
Message #03841
Re: evaluate_integrand
On Tue, 2010-04-13 at 11:59 +0200, Anders Logg wrote:
> On Tue, Apr 13, 2010 at 11:31:38AM +0200, Mehdi Nikbakht wrote:
> > On Tue, 2010-04-13 at 10:43 +0200, Anders Logg wrote:
> > > On Tue, Apr 13, 2010 at 10:31:51AM +0200, Mehdi Nikbakht wrote:
> > > > On Tue, 2010-04-13 at 09:45 +0200, Anders Logg wrote:
> > > > > On Tue, Apr 13, 2010 at 07:17:28AM +0800, Garth N. Wells wrote:
> > > > > >
> > > > > >
> > > > > > On 12/04/10 23:35, Anders Logg wrote:
> > > > > > >On Mon, Apr 12, 2010 at 10:20:13PM +0800, Garth N. Wells wrote:
> > > > > > >>
> > > > > > >>
> > > > > > >>On 12/04/10 21:49, Anders Logg wrote:
> > > > > > >>>On Mon, Apr 12, 2010 at 09:34:38PM +0800, Garth N. Wells wrote:
> > > > > > >>>>
> > > > > > >>>>
> > > > > > >>>>On 12/04/10 21:29, Anders Logg wrote:
> > > > > > >>>>>On Mon, Apr 12, 2010 at 09:21:32PM +0800, Garth N. Wells wrote:
> > > > > > >>>>>>
> > > > > > >>>>>>
> > > > > > >>>>>>On 12/04/10 21:19, Garth N. Wells wrote:
> > > > > > >>>>>>>
> > > > > > >>>>>>>
> > > > > > >>>>>>>On 12/04/10 20:47, Anders Logg wrote:
> > > > > > >>>>>>>>We are doing some work where we need to do run-time quadrature over
> > > > > > >>>>>>>>arbitrary polyhedra.
> > > > > > >>>>>>>
> > > > > > >>>>>>>We (Mehdi and I) do this already (using UFC), so I don't see why a new
> > > > > > >>>>>>>function is required. Can you explain why evaluate_tensor is not enough?
> > > > > > >>>>>>>
> > > > > > >>>>>>
> > > > > > >>>>>>I meant 'tabulate_tensor'.
> > > > > > >>>>>
> > > > > > >>>>>Which function do you call for evaluating the integrand?
> > > > > > >>>>>
> > > > > > >>>>
> > > > > > >>>>We evaluate it inside ufc::tabulate_tensor. We construct our forms
> > > > > > >>>>with an extra argument, say an object "CutCellIntegrator", which can
> > > > > > >>>>provide quadrature schemes which depend on the considered cell.
> > > > > > >>>
> > > > > > >>>That would require a special purpose code generator.
> > > > > > >>
> > > > > > >>What's wrong with that? FFC won't (and shouldn't) be able to do
> > > > > > >>everything. Just adding a function to UFC won't make FFC do what we
> > > > > > >>do now. We reuse FFC (import modules) and add special purpose
> > > > > > >>extensions.
> > > > > > >
> > > > > > >Exactly, it won't make FFC do what we need, but we could *use* FFC to
> > > > > > >do what we need (without adding a special-purpose code generator).
> > > > > > >
> > > > > > >>>Having
> > > > > > >>>evaluate_integrand would allow more flexibility for users to implement
> > > > > > >>>their own special quadrature scheme.
> > > > > > >>>
> > > > > > >>
> > > > > > >>We make "CutCellIntegrator" an abstract base class, so the user has
> > > > > > >>*complete* freedom to define the quadrature scheme and the generated
> > > > > > >>code does not depend on the scheme, since the scheme may depend on
> > > > > > >>things like how the cell 'cut' is represented.
> > > > > > >
> > > > > > >Then it sounds to me like that generated code is not at all special,
> > > > > > >but instead general purpose and should be added to UFC/FFC.
> > > > > > >
> > > > > > >And the most general interface would (I think) be an interface for
> > > > > > >evaluating the integrand at a given point. We already have the same
> > > > > > >for basis functions (evaluate_basis_function) so it is a natural
> > > > > > >extension.
> > > > > > >
> > > > > >
> > > > > > I still don't see the need for 'evaluate_integrand' unless you plan
> > > > > > to call it directly from the assembler side (i.e. DOLFIN). Is that
> > > > > > the case? Perhaps you can give me a concrete example of how you plan
> > > > > > to use it.
> > > > >
> > > > > Yes, that's the plan. In pseudo-code, this is what we want to do:
> > > > >
> > > > > for polyhedron in intersection.cut_cells:
> > > > >
> > > > > quadrature_rule = QuadratureRule(polyhedron)
> > > > > AK = 0
> > > > >
> > > > > for (x, w) in quadrature_rule:
> > > > >
> > > > > AK += w * evaluate_integrand(x)
> > > > >
> > > > > A += AK
> > > > >
> > > > > All data representing the geometry, the polyhedra, mappings from those
> > > > > polyhedra to the original cells etc is in the Intersection class (in
> > > > > the sandbox):
> > > > >
> > > > > intersection = Intersection(mesh0, mesh1)
> > > > >
> > > > > Eventually, we might want to move part of the functionality into
> > > > > either DOLFIN or FFC, but having access to and evaluate_integrand
> > > > > function makes it possible to experiment (from the C++ side) without
> > > > > the need for building complex abstractions at this point.
> > > >
> > > > As far as I understood you want to compute the integration points for
> > > > polyhedrons inside Dolfin and evaluate_integrands will just compute that
> > > > integrand in that specific integration points. If it is the case how
> > > > would you determine what is the order of quadrature rule that you want
> > > > to use?
> > >
> > > The quadrature rule would be a simple option for the user to
> > > set. Currently we only have one general rule implemented which is
> > > barycentric quadrature.
> >
> > Then if you use higher order elements, you need to update your
> > quadrature rule manually. I think it would be nice to compute your own
> > quadrature rule inside tabualte_tensor by using the standard quadrature
> > rule.
>
> Yes, but the problem is that the computation of the quadrature rule is
> nontrivial and it might be better to do it from C++ than as part of
> the generated code, especially when the code depends on external
> libraries like CGAL. See BarycenterQuadrature.cpp.
Where can I find this file?
> > >
> > > > Since to evaluate integrands you need to tabulate basis functions and
> > > > their derivatives on arbitrary points. How do you want to tabulate basis
> > > > functions and their derivatives inside evaluate_integrands?
> > >
> > > That's a good point. We would need to evaluate the basis functions and
> > > their derivatives at a given arbitrary point which is not known at
> > > compile-time.
> >
> > Yes, you need to use the tabulate_basis* functions implemented by
> > Kristian. Then I am not the only one who is using them. Good for
> > Kristian. ;)
>
> ok, good. Then there is no principal problem of generating code for
> evaluate_integrand (if it is allowed to call tabaulate_basis).
>
> I don't see what the problem is of adding evaluate_integrand. It is a
> natural extension (we have evaluate_basis already), it would be
> "simple" to implement (Kristian can correct me if I'm wrong) and it
> makes generated code useful for assembly over cut meshes (without the
> need for writing a special-purpose code generator).
I don't think it would be straightforward in the current ffc modules.
You need to tabulate basis functions and their derivatives in the
compile time which differs from what is already in the place. You should
have your own quadrature transformer to generates entries required for
this case. It is similar to what we have done for XFEM.
Please have a look on the tabulate_Tensor function inside generated code
for Poisson demo in the attached files.
Mehdi
>
> --
> Anders
// This code conforms with the UFC specification version 1.4
// and was automatically generated by FFC version 0.9.2+.
//
// This code was generated with the option '-l dolfin' and
// contains DOLFIN-specific wrappers that depend on DOLFIN.
//
// This code was generated with the following parameters:
//
// cache_dir: ''
// convert_exceptions_to_warnings: False
// cpp_optimize: False
// epsilon: 1e-14
// form_postfix: True
// format: 'dolfin'
// log_level: 20
// log_prefix: ''
// optimize: False
// output_dir: '.'
// precision: 15
// quadrature_degree: 'auto'
// quadrature_rule: 'auto'
// representation: 'quadrature'
// split: False
#ifndef __POISSON_H
#define __POISSON_H
#include <cmath>
#include <algorithm>
#include <stdexcept>
#include <fstream>
#include <boost/assign/list_of.hpp>
#include <ufc.h>
#include <pum/GenericPUM.h>
/// This class defines the interface for a finite element.
class poisson_finite_element_0: public ufc::finite_element
{
public:
/// Constructor
poisson_finite_element_0() : ufc::finite_element()
{
// Do nothing
}
/// Destructor
virtual ~poisson_finite_element_0()
{
// Do nothing
}
/// Return a string identifying the finite element
virtual const char* signature() const
{
return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
}
/// Return the cell shape
virtual ufc::shape cell_shape() const
{
return ufc::triangle;
}
/// Return the dimension of the finite element function space
virtual unsigned int space_dimension() const
{
return 3;
}
/// Return the rank of the value space
virtual unsigned int value_rank() const
{
return 0;
}
/// Return the dimension of the value space for axis i
virtual unsigned int value_dimension(unsigned int i) const
{
return 1;
}
/// Evaluate basis function i at given point in cell
virtual void evaluate_basis(unsigned int i,
double* values,
const double* coordinates,
const ufc::cell& c) const
{
// Extract vertex coordinates
const double * const * x = c.coordinates;
// Compute Jacobian of affine map from reference cell
const double J_00 = x[1][0] - x[0][0];
const double J_01 = x[2][0] - x[0][0];
const double J_10 = x[1][1] - x[0][1];
const double J_11 = x[2][1] - x[0][1];
// Compute determinant of Jacobian
double detJ = J_00*J_11 - J_01*J_10;
// Compute inverse of Jacobian
// Compute constants
const double C0 = x[1][0] + x[2][0];
const double C1 = x[1][1] + x[2][1];
// Get coordinates and map to the reference (FIAT) element
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
// Reset values.
*values = 0.000000000000000;
// Map degree of freedom to element degree of freedom
const unsigned int dof = i;
// Array of basisvalues.
double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
// Declare helper variables.
unsigned int rr = 0;
unsigned int ss = 0;
double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
// Compute basisvalues.
basisvalues[0] = 1.000000000000000;
basisvalues[1] = tmp0;
for (unsigned int r = 0; r < 1; r++)
{
rr = (r + 1)*(r + 1 + 1)/2 + 1;
ss = r*(r + 1)/2;
basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
}// end loop over 'r'
for (unsigned int r = 0; r < 2; r++)
{
for (unsigned int s = 0; s < 2 - r; s++)
{
rr = (r + s)*(r + s + 1)/2 + s;
basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
}// end loop over 's'
}// end loop over 'r'
// Table(s) of coefficients.
static const double coefficients0[3][3] = \
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.000000000000000, 0.333333333333333}};
// Compute value(s).
for (unsigned int r = 0; r < 3; r++)
{
*values += coefficients0[dof][r]*basisvalues[r];
}// end loop over 'r'
}
/// Evaluate all basis functions at given point in cell
virtual void evaluate_basis_all(double* values,
const double* coordinates,
const ufc::cell& c) const
{
// Helper variable to hold values of a single dof.
double dof_values = 0.000000000000000;
// Loop dofs and call evaluate_basis.
for (unsigned int r = 0; r < 3; r++)
{
evaluate_basis(r, &dof_values, coordinates, c);
values[r] = dof_values;
}// end loop over 'r'
}
/// Evaluate order n derivatives of basis function i at given point in cell
virtual void evaluate_basis_derivatives(unsigned int i,
unsigned int n,
double* values,
const double* coordinates,
const ufc::cell& c) const
{
// Extract vertex coordinates
const double * const * x = c.coordinates;
// Compute Jacobian of affine map from reference cell
const double J_00 = x[1][0] - x[0][0];
const double J_01 = x[2][0] - x[0][0];
const double J_10 = x[1][1] - x[0][1];
const double J_11 = x[2][1] - x[0][1];
// Compute determinant of Jacobian
double detJ = J_00*J_11 - J_01*J_10;
// Compute inverse of Jacobian
const double K_00 = J_11 / detJ;
const double K_01 = -J_01 / detJ;
const double K_10 = -J_10 / detJ;
const double K_11 = J_00 / detJ;
// Compute constants
const double C0 = x[1][0] + x[2][0];
const double C1 = x[1][1] + x[2][1];
// Get coordinates and map to the reference (FIAT) element
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
// Compute number of derivatives.
unsigned int num_derivatives = 1;
for (unsigned int r = 0; r < n; r++)
{
num_derivatives *= 2;
}// end loop over 'r'
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
unsigned int **combinations = new unsigned int *[num_derivatives];
for (unsigned int row = 0; row < num_derivatives; row++)
{
combinations[row] = new unsigned int [n];
for (unsigned int col = 0; col < n; col++)
combinations[row][col] = 0;
}
// Generate combinations of derivatives
for (unsigned int row = 1; row < num_derivatives; row++)
{
for (unsigned int num = 0; num < row; num++)
{
for (unsigned int col = n-1; col+1 > 0; col--)
{
if (combinations[row][col] + 1 > 1)
combinations[row][col] = 0;
else
{
combinations[row][col] += 1;
break;
}
}
}
}
// Compute inverse of Jacobian
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
// Declare transformation matrix
// Declare pointer to two dimensional array and initialise
double **transform = new double *[num_derivatives];
for (unsigned int j = 0; j < num_derivatives; j++)
{
transform[j] = new double [num_derivatives];
for (unsigned int k = 0; k < num_derivatives; k++)
transform[j][k] = 1;
}
// Construct transformation matrix
for (unsigned int row = 0; row < num_derivatives; row++)
{
for (unsigned int col = 0; col < num_derivatives; col++)
{
for (unsigned int k = 0; k < n; k++)
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
}
}
// Reset values. Assuming that values is always an array.
for (unsigned int r = 0; r < num_derivatives; r++)
{
values[r] = 0.000000000000000;
}// end loop over 'r'
// Map degree of freedom to element degree of freedom
const unsigned int dof = i;
// Array of basisvalues.
double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
// Declare helper variables.
unsigned int rr = 0;
unsigned int ss = 0;
double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
// Compute basisvalues.
basisvalues[0] = 1.000000000000000;
basisvalues[1] = tmp0;
for (unsigned int r = 0; r < 1; r++)
{
rr = (r + 1)*(r + 1 + 1)/2 + 1;
ss = r*(r + 1)/2;
basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
}// end loop over 'r'
for (unsigned int r = 0; r < 2; r++)
{
for (unsigned int s = 0; s < 2 - r; s++)
{
rr = (r + s)*(r + s + 1)/2 + s;
basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
}// end loop over 's'
}// end loop over 'r'
// Table(s) of coefficients.
static const double coefficients0[3][3] = \
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.000000000000000, 0.333333333333333}};
// Tables of derivatives of the polynomial base (transpose).
static const double dmats0[3][3] = \
{{0.000000000000000, 0.000000000000000, 0.000000000000000},
{4.898979485566356, 0.000000000000000, 0.000000000000000},
{0.000000000000000, 0.000000000000000, 0.000000000000000}};
static const double dmats1[3][3] = \
{{0.000000000000000, 0.000000000000000, 0.000000000000000},
{2.449489742783178, 0.000000000000000, 0.000000000000000},
{4.242640687119285, 0.000000000000000, 0.000000000000000}};
// Compute reference derivatives.
// Declare pointer to array of derivatives on FIAT element.
double *derivatives = new double[num_derivatives];
for (unsigned int r = 0; r < num_derivatives; r++)
{
derivatives[r] = 0.000000000000000;
}// end loop over 'r'
// Declare derivative matrix (of polynomial basis).
double dmats[3][3] = \
{{1.000000000000000, 0.000000000000000, 0.000000000000000},
{0.000000000000000, 1.000000000000000, 0.000000000000000},
{0.000000000000000, 0.000000000000000, 1.000000000000000}};
// Declare (auxiliary) derivative matrix (of polynomial basis).
double dmats_old[3][3] = \
{{1.000000000000000, 0.000000000000000, 0.000000000000000},
{0.000000000000000, 1.000000000000000, 0.000000000000000},
{0.000000000000000, 0.000000000000000, 1.000000000000000}};
// Loop possible derivatives.
for (unsigned int r = 0; r < num_derivatives; r++)
{
// Resetting dmats values to compute next derivative.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats[t][u] = 0.000000000000000;
if (t == u)
{
dmats[t][u] = 1.000000000000000;
}
}// end loop over 'u'
}// end loop over 't'
// Looping derivative order to generate dmats.
for (unsigned int s = 0; s < n; s++)
{
// Updating dmats_old with new values and resetting dmats.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats_old[t][u] = dmats[t][u];
dmats[t][u] = 0.000000000000000;
}// end loop over 'u'
}// end loop over 't'
// Update dmats using an inner product.
if (combinations[r][s] == 0)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
}// end loop over 'tu'
}// end loop over 'u'
}// end loop over 't'
}
if (combinations[r][s] == 1)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
}// end loop over 'tu'
}// end loop over 'u'
}// end loop over 't'
}
}// end loop over 's'
for (unsigned int s = 0; s < 3; s++)
{
for (unsigned int t = 0; t < 3; t++)
{
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
}// end loop over 't'
}// end loop over 's'
}// end loop over 'r'
// Transform derivatives back to physical element
for (unsigned int r = 0; r < num_derivatives; r++)
{
for (unsigned int s = 0; s < num_derivatives; s++)
{
values[r] += transform[r][s]*derivatives[s];
}// end loop over 's'
}// end loop over 'r'
// Delete pointer to array of derivatives on FIAT element
delete [] derivatives;
// Delete pointer to array of combinations of derivatives and transform
for (unsigned int r = 0; r < num_derivatives; r++)
{
delete [] combinations[r];
}// end loop over 'r'
delete [] combinations;
for (unsigned int r = 0; r < num_derivatives; r++)
{
delete [] transform[r];
}// end loop over 'r'
delete [] transform;
}
/// Evaluate order n derivatives of all basis functions at given point in cell
virtual void evaluate_basis_derivatives_all(unsigned int n,
double* values,
const double* coordinates,
const ufc::cell& c) const
{
// Compute number of derivatives.
unsigned int num_derivatives = 1;
for (unsigned int r = 0; r < n; r++)
{
num_derivatives *= 2;
}// end loop over 'r'
// Helper variable to hold values of a single dof.
double *dof_values = new double[num_derivatives];
for (unsigned int r = 0; r < num_derivatives; r++)
{
dof_values[r] = 0.000000000000000;
}// end loop over 'r'
// Loop dofs and call evaluate_basis_derivatives.
for (unsigned int r = 0; r < 3; r++)
{
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
for (unsigned int s = 0; s < num_derivatives; s++)
{
values[r*num_derivatives + s] = dof_values[s];
}// end loop over 's'
}// end loop over 'r'
// Delete pointer.
delete [] dof_values;
}
/// Evaluate linear functional for dof i on the function f
virtual double evaluate_dof(unsigned int i,
const ufc::function& f,
const ufc::cell& c) const
{
// Declare variables for result of evaluation.
double vals[1];
// Declare variable for physical coordinates.
double y[2];
const double * const * x = c.coordinates;
switch (i)
{
case 0:
{
y[0] = x[0][0];
y[1] = x[0][1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
case 1:
{
y[0] = x[1][0];
y[1] = x[1][1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
case 2:
{
y[0] = x[2][0];
y[1] = x[2][1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
}
return 0.000000000000000;
}
/// Evaluate linear functionals for all dofs on the function f
virtual void evaluate_dofs(double* values,
const ufc::function& f,
const ufc::cell& c) const
{
// Declare variables for result of evaluation.
double vals[1];
// Declare variable for physical coordinates.
double y[2];
const double * const * x = c.coordinates;
y[0] = x[0][0];
y[1] = x[0][1];
f.evaluate(vals, y, c);
values[0] = vals[0];
y[0] = x[1][0];
y[1] = x[1][1];
f.evaluate(vals, y, c);
values[1] = vals[0];
y[0] = x[2][0];
y[1] = x[2][1];
f.evaluate(vals, y, c);
values[2] = vals[0];
}
/// Interpolate vertex values from dof values
virtual void interpolate_vertex_values(double* vertex_values,
const double* dof_values,
const ufc::cell& c) const
{
// Evaluate function and change variables
vertex_values[0] = dof_values[0];
vertex_values[1] = dof_values[1];
vertex_values[2] = dof_values[2];
}
/// Return the number of sub elements (for a mixed element)
virtual unsigned int num_sub_elements() const
{
return 0;
}
/// Create a new finite element for sub element i (for a mixed element)
virtual ufc::finite_element* create_sub_element(unsigned int i) const
{
return 0;
}
};
/// This class defines the interface for a finite element.
class poisson_finite_element_1: public ufc::finite_element
{
public:
/// Constructor
poisson_finite_element_1() : ufc::finite_element()
{
// Do nothing
}
/// Destructor
virtual ~poisson_finite_element_1()
{
// Do nothing
}
/// Return a string identifying the finite element
virtual const char* signature() const
{
return "EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None)))";
}
/// Return the cell shape
virtual ufc::shape cell_shape() const
{
return ufc::triangle;
}
/// Return the dimension of the finite element function space
virtual unsigned int space_dimension() const
{
return 6;
}
/// Return the rank of the value space
virtual unsigned int value_rank() const
{
return 0;
}
/// Return the dimension of the value space for axis i
virtual unsigned int value_dimension(unsigned int i) const
{
return 1;
}
/// Evaluate basis function i at given point in cell
virtual void evaluate_basis(unsigned int i,
double* values,
const double* coordinates,
const ufc::cell& c) const
{
// Extract vertex coordinates
const double * const * x = c.coordinates;
// Compute Jacobian of affine map from reference cell
const double J_00 = x[1][0] - x[0][0];
const double J_01 = x[2][0] - x[0][0];
const double J_10 = x[1][1] - x[0][1];
const double J_11 = x[2][1] - x[0][1];
// Compute determinant of Jacobian
double detJ = J_00*J_11 - J_01*J_10;
// Compute inverse of Jacobian
// Compute constants
const double C0 = x[1][0] + x[2][0];
const double C1 = x[1][1] + x[2][1];
// Get coordinates and map to the reference (FIAT) element
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
// Reset values.
values[0] = 0.000000000000000;
values[1] = 0.000000000000000;
if (0 <= i && i <= 2)
{
// Map degree of freedom to element degree of freedom
const unsigned int dof = i;
// Array of basisvalues.
double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
// Declare helper variables.
unsigned int rr = 0;
unsigned int ss = 0;
double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
// Compute basisvalues.
basisvalues[0] = 1.000000000000000;
basisvalues[1] = tmp0;
for (unsigned int r = 0; r < 1; r++)
{
rr = (r + 1)*(r + 1 + 1)/2 + 1;
ss = r*(r + 1)/2;
basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
}// end loop over 'r'
for (unsigned int r = 0; r < 2; r++)
{
for (unsigned int s = 0; s < 2 - r; s++)
{
rr = (r + s)*(r + s + 1)/2 + s;
basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
}// end loop over 's'
}// end loop over 'r'
// Table(s) of coefficients.
static const double coefficients0[3][3] = \
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.000000000000000, 0.333333333333333}};
// Compute value(s).
for (unsigned int r = 0; r < 3; r++)
{
values[0] += coefficients0[dof][r]*basisvalues[r];
}// end loop over 'r'
}
if (3 <= i && i <= 5)
{
// Map degree of freedom to element degree of freedom
const unsigned int dof = i - 3;
// Array of basisvalues.
double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
// Declare helper variables.
unsigned int rr = 0;
unsigned int ss = 0;
double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
// Compute basisvalues.
basisvalues[0] = 1.000000000000000;
basisvalues[1] = tmp0;
for (unsigned int r = 0; r < 1; r++)
{
rr = (r + 1)*(r + 1 + 1)/2 + 1;
ss = r*(r + 1)/2;
basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
}// end loop over 'r'
for (unsigned int r = 0; r < 2; r++)
{
for (unsigned int s = 0; s < 2 - r; s++)
{
rr = (r + s)*(r + s + 1)/2 + s;
basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
}// end loop over 's'
}// end loop over 'r'
// Table(s) of coefficients.
static const double coefficients0[3][3] = \
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.000000000000000, 0.333333333333333}};
// Compute value(s).
for (unsigned int r = 0; r < 3; r++)
{
values[1] += coefficients0[dof][r]*basisvalues[r];
}// end loop over 'r'
}
}
/// Evaluate all basis functions at given point in cell
virtual void evaluate_basis_all(double* values,
const double* coordinates,
const ufc::cell& c) const
{
// Helper variable to hold values of a single dof.
double dof_values[2] = {0.000000000000000, 0.000000000000000};
// Loop dofs and call evaluate_basis.
for (unsigned int r = 0; r < 6; r++)
{
evaluate_basis(r, dof_values, coordinates, c);
for (unsigned int s = 0; s < 2; s++)
{
values[r*2 + s] = dof_values[s];
}// end loop over 's'
}// end loop over 'r'
}
/// Evaluate order n derivatives of basis function i at given point in cell
virtual void evaluate_basis_derivatives(unsigned int i,
unsigned int n,
double* values,
const double* coordinates,
const ufc::cell& c) const
{
// Extract vertex coordinates
const double * const * x = c.coordinates;
// Compute Jacobian of affine map from reference cell
const double J_00 = x[1][0] - x[0][0];
const double J_01 = x[2][0] - x[0][0];
const double J_10 = x[1][1] - x[0][1];
const double J_11 = x[2][1] - x[0][1];
// Compute determinant of Jacobian
double detJ = J_00*J_11 - J_01*J_10;
// Compute inverse of Jacobian
const double K_00 = J_11 / detJ;
const double K_01 = -J_01 / detJ;
const double K_10 = -J_10 / detJ;
const double K_11 = J_00 / detJ;
// Compute constants
const double C0 = x[1][0] + x[2][0];
const double C1 = x[1][1] + x[2][1];
// Get coordinates and map to the reference (FIAT) element
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
// Compute number of derivatives.
unsigned int num_derivatives = 1;
for (unsigned int r = 0; r < n; r++)
{
num_derivatives *= 2;
}// end loop over 'r'
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
unsigned int **combinations = new unsigned int *[num_derivatives];
for (unsigned int row = 0; row < num_derivatives; row++)
{
combinations[row] = new unsigned int [n];
for (unsigned int col = 0; col < n; col++)
combinations[row][col] = 0;
}
// Generate combinations of derivatives
for (unsigned int row = 1; row < num_derivatives; row++)
{
for (unsigned int num = 0; num < row; num++)
{
for (unsigned int col = n-1; col+1 > 0; col--)
{
if (combinations[row][col] + 1 > 1)
combinations[row][col] = 0;
else
{
combinations[row][col] += 1;
break;
}
}
}
}
// Compute inverse of Jacobian
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
// Declare transformation matrix
// Declare pointer to two dimensional array and initialise
double **transform = new double *[num_derivatives];
for (unsigned int j = 0; j < num_derivatives; j++)
{
transform[j] = new double [num_derivatives];
for (unsigned int k = 0; k < num_derivatives; k++)
transform[j][k] = 1;
}
// Construct transformation matrix
for (unsigned int row = 0; row < num_derivatives; row++)
{
for (unsigned int col = 0; col < num_derivatives; col++)
{
for (unsigned int k = 0; k < n; k++)
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
}
}
// Reset values. Assuming that values is always an array.
for (unsigned int r = 0; r < 2*num_derivatives; r++)
{
values[r] = 0.000000000000000;
}// end loop over 'r'
if (0 <= i && i <= 2)
{
// Map degree of freedom to element degree of freedom
const unsigned int dof = i;
// Array of basisvalues.
double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
// Declare helper variables.
unsigned int rr = 0;
unsigned int ss = 0;
double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
// Compute basisvalues.
basisvalues[0] = 1.000000000000000;
basisvalues[1] = tmp0;
for (unsigned int r = 0; r < 1; r++)
{
rr = (r + 1)*(r + 1 + 1)/2 + 1;
ss = r*(r + 1)/2;
basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
}// end loop over 'r'
for (unsigned int r = 0; r < 2; r++)
{
for (unsigned int s = 0; s < 2 - r; s++)
{
rr = (r + s)*(r + s + 1)/2 + s;
basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
}// end loop over 's'
}// end loop over 'r'
// Table(s) of coefficients.
static const double coefficients0[3][3] = \
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.000000000000000, 0.333333333333333}};
// Tables of derivatives of the polynomial base (transpose).
static const double dmats0[3][3] = \
{{0.000000000000000, 0.000000000000000, 0.000000000000000},
{4.898979485566356, 0.000000000000000, 0.000000000000000},
{0.000000000000000, 0.000000000000000, 0.000000000000000}};
static const double dmats1[3][3] = \
{{0.000000000000000, 0.000000000000000, 0.000000000000000},
{2.449489742783178, 0.000000000000000, 0.000000000000000},
{4.242640687119285, 0.000000000000000, 0.000000000000000}};
// Compute reference derivatives.
// Declare pointer to array of derivatives on FIAT element.
double *derivatives = new double[num_derivatives];
for (unsigned int r = 0; r < num_derivatives; r++)
{
derivatives[r] = 0.000000000000000;
}// end loop over 'r'
// Declare derivative matrix (of polynomial basis).
double dmats[3][3] = \
{{1.000000000000000, 0.000000000000000, 0.000000000000000},
{0.000000000000000, 1.000000000000000, 0.000000000000000},
{0.000000000000000, 0.000000000000000, 1.000000000000000}};
// Declare (auxiliary) derivative matrix (of polynomial basis).
double dmats_old[3][3] = \
{{1.000000000000000, 0.000000000000000, 0.000000000000000},
{0.000000000000000, 1.000000000000000, 0.000000000000000},
{0.000000000000000, 0.000000000000000, 1.000000000000000}};
// Loop possible derivatives.
for (unsigned int r = 0; r < num_derivatives; r++)
{
// Resetting dmats values to compute next derivative.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats[t][u] = 0.000000000000000;
if (t == u)
{
dmats[t][u] = 1.000000000000000;
}
}// end loop over 'u'
}// end loop over 't'
// Looping derivative order to generate dmats.
for (unsigned int s = 0; s < n; s++)
{
// Updating dmats_old with new values and resetting dmats.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats_old[t][u] = dmats[t][u];
dmats[t][u] = 0.000000000000000;
}// end loop over 'u'
}// end loop over 't'
// Update dmats using an inner product.
if (combinations[r][s] == 0)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
}// end loop over 'tu'
}// end loop over 'u'
}// end loop over 't'
}
if (combinations[r][s] == 1)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
}// end loop over 'tu'
}// end loop over 'u'
}// end loop over 't'
}
}// end loop over 's'
for (unsigned int s = 0; s < 3; s++)
{
for (unsigned int t = 0; t < 3; t++)
{
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
}// end loop over 't'
}// end loop over 's'
}// end loop over 'r'
// Transform derivatives back to physical element
for (unsigned int r = 0; r < num_derivatives; r++)
{
for (unsigned int s = 0; s < num_derivatives; s++)
{
values[r] += transform[r][s]*derivatives[s];
}// end loop over 's'
}// end loop over 'r'
// Delete pointer to array of derivatives on FIAT element
delete [] derivatives;
// Delete pointer to array of combinations of derivatives and transform
for (unsigned int r = 0; r < num_derivatives; r++)
{
delete [] combinations[r];
}// end loop over 'r'
delete [] combinations;
for (unsigned int r = 0; r < num_derivatives; r++)
{
delete [] transform[r];
}// end loop over 'r'
delete [] transform;
}
if (3 <= i && i <= 5)
{
// Map degree of freedom to element degree of freedom
const unsigned int dof = i - 3;
// Array of basisvalues.
double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
// Declare helper variables.
unsigned int rr = 0;
unsigned int ss = 0;
double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
// Compute basisvalues.
basisvalues[0] = 1.000000000000000;
basisvalues[1] = tmp0;
for (unsigned int r = 0; r < 1; r++)
{
rr = (r + 1)*(r + 1 + 1)/2 + 1;
ss = r*(r + 1)/2;
basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
}// end loop over 'r'
for (unsigned int r = 0; r < 2; r++)
{
for (unsigned int s = 0; s < 2 - r; s++)
{
rr = (r + s)*(r + s + 1)/2 + s;
basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
}// end loop over 's'
}// end loop over 'r'
// Table(s) of coefficients.
static const double coefficients0[3][3] = \
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.000000000000000, 0.333333333333333}};
// Tables of derivatives of the polynomial base (transpose).
static const double dmats0[3][3] = \
{{0.000000000000000, 0.000000000000000, 0.000000000000000},
{4.898979485566356, 0.000000000000000, 0.000000000000000},
{0.000000000000000, 0.000000000000000, 0.000000000000000}};
static const double dmats1[3][3] = \
{{0.000000000000000, 0.000000000000000, 0.000000000000000},
{2.449489742783178, 0.000000000000000, 0.000000000000000},
{4.242640687119285, 0.000000000000000, 0.000000000000000}};
// Compute reference derivatives.
// Declare pointer to array of derivatives on FIAT element.
double *derivatives = new double[num_derivatives];
for (unsigned int r = 0; r < num_derivatives; r++)
{
derivatives[r] = 0.000000000000000;
}// end loop over 'r'
// Declare derivative matrix (of polynomial basis).
double dmats[3][3] = \
{{1.000000000000000, 0.000000000000000, 0.000000000000000},
{0.000000000000000, 1.000000000000000, 0.000000000000000},
{0.000000000000000, 0.000000000000000, 1.000000000000000}};
// Declare (auxiliary) derivative matrix (of polynomial basis).
double dmats_old[3][3] = \
{{1.000000000000000, 0.000000000000000, 0.000000000000000},
{0.000000000000000, 1.000000000000000, 0.000000000000000},
{0.000000000000000, 0.000000000000000, 1.000000000000000}};
// Loop possible derivatives.
for (unsigned int r = 0; r < num_derivatives; r++)
{
// Resetting dmats values to compute next derivative.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats[t][u] = 0.000000000000000;
if (t == u)
{
dmats[t][u] = 1.000000000000000;
}
}// end loop over 'u'
}// end loop over 't'
// Looping derivative order to generate dmats.
for (unsigned int s = 0; s < n; s++)
{
// Updating dmats_old with new values and resetting dmats.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats_old[t][u] = dmats[t][u];
dmats[t][u] = 0.000000000000000;
}// end loop over 'u'
}// end loop over 't'
// Update dmats using an inner product.
if (combinations[r][s] == 0)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
}// end loop over 'tu'
}// end loop over 'u'
}// end loop over 't'
}
if (combinations[r][s] == 1)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
}// end loop over 'tu'
}// end loop over 'u'
}// end loop over 't'
}
}// end loop over 's'
for (unsigned int s = 0; s < 3; s++)
{
for (unsigned int t = 0; t < 3; t++)
{
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
}// end loop over 't'
}// end loop over 's'
}// end loop over 'r'
// Transform derivatives back to physical element
for (unsigned int r = 0; r < num_derivatives; r++)
{
for (unsigned int s = 0; s < num_derivatives; s++)
{
values[num_derivatives + r] += transform[r][s]*derivatives[s];
}// end loop over 's'
}// end loop over 'r'
// Delete pointer to array of derivatives on FIAT element
delete [] derivatives;
// Delete pointer to array of combinations of derivatives and transform
for (unsigned int r = 0; r < num_derivatives; r++)
{
delete [] combinations[r];
}// end loop over 'r'
delete [] combinations;
for (unsigned int r = 0; r < num_derivatives; r++)
{
delete [] transform[r];
}// end loop over 'r'
delete [] transform;
}
}
/// Evaluate order n derivatives of all basis functions at given point in cell
virtual void evaluate_basis_derivatives_all(unsigned int n,
double* values,
const double* coordinates,
const ufc::cell& c) const
{
// Compute number of derivatives.
unsigned int num_derivatives = 1;
for (unsigned int r = 0; r < n; r++)
{
num_derivatives *= 2;
}// end loop over 'r'
// Helper variable to hold values of a single dof.
double *dof_values = new double[2*num_derivatives];
for (unsigned int r = 0; r < 2*num_derivatives; r++)
{
dof_values[r] = 0.000000000000000;
}// end loop over 'r'
// Loop dofs and call evaluate_basis_derivatives.
for (unsigned int r = 0; r < 6; r++)
{
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
for (unsigned int s = 0; s < 2*num_derivatives; s++)
{
values[r*2*num_derivatives + s] = dof_values[s];
}// end loop over 's'
}// end loop over 'r'
// Delete pointer.
delete [] dof_values;
}
/// Evaluate linear functional for dof i on the function f
virtual double evaluate_dof(unsigned int i,
const ufc::function& f,
const ufc::cell& c) const
{
// Declare variables for result of evaluation.
double vals[1];
// Declare variable for physical coordinates.
double y[2];
const double * const * x = c.coordinates;
switch (i)
{
case 0:
{
y[0] = x[0][0];
y[1] = x[0][1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
case 1:
{
y[0] = x[1][0];
y[1] = x[1][1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
case 2:
{
y[0] = x[2][0];
y[1] = x[2][1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
case 3:
{
y[0] = x[0][0];
y[1] = x[0][1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
case 4:
{
y[0] = x[1][0];
y[1] = x[1][1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
case 5:
{
y[0] = x[2][0];
y[1] = x[2][1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
}
return 0.000000000000000;
}
/// Evaluate linear functionals for all dofs on the function f
virtual void evaluate_dofs(double* values,
const ufc::function& f,
const ufc::cell& c) const
{
// Declare variables for result of evaluation.
double vals[1];
// Declare variable for physical coordinates.
double y[2];
const double * const * x = c.coordinates;
y[0] = x[0][0];
y[1] = x[0][1];
f.evaluate(vals, y, c);
values[0] = vals[0];
y[0] = x[1][0];
y[1] = x[1][1];
f.evaluate(vals, y, c);
values[1] = vals[0];
y[0] = x[2][0];
y[1] = x[2][1];
f.evaluate(vals, y, c);
values[2] = vals[0];
y[0] = x[0][0];
y[1] = x[0][1];
f.evaluate(vals, y, c);
values[3] = vals[0];
y[0] = x[1][0];
y[1] = x[1][1];
f.evaluate(vals, y, c);
values[4] = vals[0];
y[0] = x[2][0];
y[1] = x[2][1];
f.evaluate(vals, y, c);
values[5] = vals[0];
}
/// Interpolate vertex values from dof values
virtual void interpolate_vertex_values(double* vertex_values,
const double* dof_values,
const ufc::cell& c) const
{
// Evaluate function and change variables
vertex_values[0] = dof_values[0] + dof_values[3];
vertex_values[1] = dof_values[1] + dof_values[4];
vertex_values[2] = dof_values[2] + dof_values[5];
}
/// Return the number of sub elements (for a mixed element)
virtual unsigned int num_sub_elements() const
{
return 0;
}
/// Create a new finite element for sub element i (for a mixed element)
virtual ufc::finite_element* create_sub_element(unsigned int i) const
{
return 0;
}
};
/// This class defines the interface for a local-to-global mapping of
/// degrees of freedom (dofs).
class poisson_dof_map_0: public ufc::dof_map
{
private:
unsigned int _global_dimension;
public:
/// Constructor
poisson_dof_map_0( ) :ufc::dof_map()
{
_global_dimension = 0;
}
/// Destructor
virtual ~poisson_dof_map_0()
{
// Do nothing
}
/// Return a string identifying the dof map
virtual const char* signature() const
{
return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
}
/// Return true iff mesh entities of topological dimension d are needed
virtual bool needs_mesh_entities(unsigned int d) const
{
switch (d)
{
case 0:
{
return true;
break;
}
case 1:
{
return false;
break;
}
case 2:
{
return false;
break;
}
}
return false;
}
/// Initialize dof map for mesh (return true iff init_cell() is needed)
virtual bool init_mesh(const ufc::mesh& m)
{
_global_dimension = m.num_entities[0];
return false;
}
/// Initialize dof map for given cell
virtual void init_cell(const ufc::mesh& m,
const ufc::cell& c)
{
// Do nothing
}
/// Finish initialization of dof map for cells
virtual void init_cell_finalize()
{
// Do nothing
}
/// Return the dimension of the global finite element function space
virtual unsigned int global_dimension() const
{
return _global_dimension;
}
/// Return the dimension of the local finite element function space for a cell
virtual unsigned int local_dimension(const ufc::cell& c) const
{
return 3;
}
/// Return the maximum dimension of the local finite element function space
virtual unsigned int max_local_dimension() const
{
return 3;
}
// Return the geometric dimension of the coordinates this dof map provides
virtual unsigned int geometric_dimension() const
{
return 2;
}
/// Return the number of dofs on each cell facet
virtual unsigned int num_facet_dofs() const
{
return 2;
}
/// Return the number of dofs associated with each cell entity of dimension d
virtual unsigned int num_entity_dofs(unsigned int d) const
{
switch (d)
{
case 0:
{
return 1;
break;
}
case 1:
{
return 0;
break;
}
case 2:
{
return 0;
break;
}
}
return 0;
}
/// Tabulate the local-to-global mapping of dofs on a cell
virtual void tabulate_dofs(unsigned int* dofs,
const ufc::mesh& m,
const ufc::cell& c) const
{
dofs[0] = c.entity_indices[0][0];
dofs[1] = c.entity_indices[0][1];
dofs[2] = c.entity_indices[0][2];
}
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
virtual void tabulate_facet_dofs(unsigned int* dofs,
unsigned int facet) const
{
switch (facet)
{
case 0:
{
dofs[0] = 1;
dofs[1] = 2;
break;
}
case 1:
{
dofs[0] = 0;
dofs[1] = 2;
break;
}
case 2:
{
dofs[0] = 0;
dofs[1] = 1;
break;
}
}
}
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
virtual void tabulate_entity_dofs(unsigned int* dofs,
unsigned int d, unsigned int i) const
{
if (d > 2)
{
throw std::runtime_error("d is larger than dimension (2)");
}
switch (d)
{
case 0:
{
if (i > 2)
{
throw std::runtime_error("i is larger than number of entities (2)");
}
switch (i)
{
case 0:
{
dofs[0] = 0;
break;
}
case 1:
{
dofs[0] = 1;
break;
}
case 2:
{
dofs[0] = 2;
break;
}
}
break;
}
case 1:
{
break;
}
case 2:
{
break;
}
}
}
/// Tabulate the coordinates of all dofs on a cell
virtual void tabulate_coordinates(double** coordinates,
const ufc::cell& c) const
{
const double * const * x = c.coordinates;
coordinates[0][0] = x[0][0];
coordinates[0][1] = x[0][1];
coordinates[1][0] = x[1][0];
coordinates[1][1] = x[1][1];
coordinates[2][0] = x[2][0];
coordinates[2][1] = x[2][1];
}
/// Return the number of sub dof maps (for a mixed element)
virtual unsigned int num_sub_dof_maps() const
{
return 0;
}
/// Create a new dof_map for sub dof map i (for a mixed element)
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
{
return 0;
}
};
/// This class defines the interface for a local-to-global mapping of
/// degrees of freedom (dofs).
class poisson_dof_map_1: public ufc::dof_map
{
private:
unsigned int _global_dimension;
const std::vector<const pum::GenericPUM*>& pums;
public:
/// Constructor
poisson_dof_map_1( const std::vector<const pum::GenericPUM*>& pums) :ufc::dof_map() , pums(pums)
{
_global_dimension = 0;
}
/// Destructor
virtual ~poisson_dof_map_1()
{
// Do nothing
}
/// Return a string identifying the dof map
virtual const char* signature() const
{
return "FFC dofmap for EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None)))";
}
/// Return true iff mesh entities of topological dimension d are needed
virtual bool needs_mesh_entities(unsigned int d) const
{
switch (d)
{
case 0:
{
return true;
break;
}
case 1:
{
return false;
break;
}
case 2:
{
return false;
break;
}
}
return false;
}
/// Initialize dof map for mesh (return true iff init_cell() is needed)
virtual bool init_mesh(const ufc::mesh& m)
{
_global_dimension = m.num_entities[0];
return false;
}
/// Initialize dof map for given cell
virtual void init_cell(const ufc::mesh& m,
const ufc::cell& c)
{
// Do nothing
}
/// Finish initialization of dof map for cells
virtual void init_cell_finalize()
{
// Do nothing
}
/// Return the dimension of the global finite element function space
virtual unsigned int global_dimension() const
{
return _global_dimension + pums[0]->enriched_global_dimension();
}
/// Return the dimension of the local finite element function space for a cell
virtual unsigned int local_dimension(const ufc::cell& c) const
{
return 3 + pums[0]->enriched_local_dimension(c);
}
/// Return the maximum dimension of the local finite element function space
virtual unsigned int max_local_dimension() const
{
return 3 + pums[0]->enriched_max_local_dimension();
}
// Return the geometric dimension of the coordinates this dof map provides
virtual unsigned int geometric_dimension() const
{
return 2;
}
/// Return the number of dofs on each cell facet
virtual unsigned int num_facet_dofs() const
{
return 2;
}
/// Return the number of dofs associated with each cell entity of dimension d
virtual unsigned int num_entity_dofs(unsigned int d) const
{
switch (d)
{
case 0:
{
return 1;
break;
}
case 1:
{
return 0;
break;
}
case 2:
{
return 0;
break;
}
}
return 0;
}
/// Tabulate the local-to-global mapping of dofs on a cell
virtual void tabulate_dofs(unsigned int* dofs,
const ufc::mesh& m,
const ufc::cell& c) const
{
dofs[0] = c.entity_indices[0][0];
dofs[1] = c.entity_indices[0][1];
dofs[2] = c.entity_indices[0][2];
\
// Generate code for tabulating extra degrees of freedom.
unsigned int local_offset = 3;
unsigned int global_offset = m.num_entities[0];
// Calculate local-to-global mapping for the enriched dofs related to the discontinuous field 0
pums[0]->tabulate_enriched_dofs(dofs, c, local_offset, global_offset);
}
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
virtual void tabulate_facet_dofs(unsigned int* dofs,
unsigned int facet) const
{
switch (facet)
{
case 0:
{
dofs[0] = 1;
dofs[1] = 2;
break;
}
case 1:
{
dofs[0] = 0;
dofs[1] = 2;
break;
}
case 2:
{
dofs[0] = 0;
dofs[1] = 1;
break;
}
}
}
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
virtual void tabulate_entity_dofs(unsigned int* dofs,
unsigned int d, unsigned int i) const
{
if (d > 2)
{
throw std::runtime_error("d is larger than dimension (2)");
}
switch (d)
{
case 0:
{
if (i > 2)
{
throw std::runtime_error("i is larger than number of entities (2)");
}
switch (i)
{
case 0:
{
dofs[0] = 0;
break;
}
case 1:
{
dofs[0] = 1;
break;
}
case 2:
{
dofs[0] = 2;
break;
}
}
break;
}
case 1:
{
break;
}
case 2:
{
break;
}
}
}
/// Tabulate the coordinates of all dofs on a cell
virtual void tabulate_coordinates(double** coordinates,
const ufc::cell& c) const
{
const double * const * x = c.coordinates;
coordinates[0][0] = x[0][0];
coordinates[0][1] = x[0][1];
coordinates[1][0] = x[1][0];
coordinates[1][1] = x[1][1];
coordinates[2][0] = x[2][0];
coordinates[2][1] = x[2][1];
}
/// Return the number of sub dof maps (for a mixed element)
virtual unsigned int num_sub_dof_maps() const
{
return 0;
}
/// Create a new dof_map for sub dof map i (for a mixed element)
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
{
return 0;
}
};
/// This class defines the interface for the tabulation of the cell
/// tensor corresponding to the local contribution to a form from
/// the integral over a cell.
class poisson_cell_integral_0_0: public ufc::cell_integral
{
const std::vector<const pum::GenericPUM*>& pums;
mutable std::vector <double> Aa;
/// Tabulate the regular entities of tensor for the contribution from a local cell
virtual void tabulate_regular_tensor(double* A,
const double * const * w,
const ufc::cell& c) const
{
// Extract vertex coordinates
const double * const * x = c.coordinates;
// Compute Jacobian of affine map from reference cell
const double J_00 = x[1][0] - x[0][0];
const double J_01 = x[2][0] - x[0][0];
const double J_10 = x[1][1] - x[0][1];
const double J_11 = x[2][1] - x[0][1];
// Compute determinant of Jacobian
double detJ = J_00*J_11 - J_01*J_10;
// Compute inverse of Jacobian
const double K_00 = J_11 / detJ;
const double K_01 = -J_01 / detJ;
const double K_10 = -J_10 / detJ;
const double K_11 = J_00 / detJ;
// Set scale factor
const double det = std::abs(detJ);
// Array of quadrature weights.
static const double W1 = 0.500000000000000;
// Quadrature points on the UFC reference element: (0.333333333333333, 0.333333333333333)
// Value of basis functions at quadrature points.
static const double FE0_D01[1][6] = \
{{-1.000000000000000, 0.000000000000000, 1.000000000000000, -1.000000000000000, 0.000000000000000, 1.000000000000000}};
static const double FE0_D10[1][6] = \
{{-1.000000000000000, 1.000000000000000, 0.000000000000000, -1.000000000000000, 1.000000000000000, 0.000000000000000}};
static const double FE1[1][3] = \
{{0.333333333333333, 0.333333333333333, 0.333333333333333}};
// enriched local dimension of the current cell
unsigned int offset = pums[0]->enriched_local_dimension(c);
// Reset values in the element tensor.
const unsigned int num_entries = (3 + offset)*(3 + offset);
for (unsigned int r = 0; r < num_entries; r++)
{
A[r] = 0.000000000000000;
}// end loop over 'r'
// Compute element tensor using UFL quadrature representation
// Optimisations: ('optimisation', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False), ('ignore zero tables', False)
// Loop quadrature points for integral.
// Number of operations to compute element tensor for following IP loop = 25
// Only 1 integration point, omitting IP loop.
// Coefficient declarations.
double F0 = 0.000000000000000;
// Total number of operations to compute function values = 6
for (unsigned int r = 0; r < 3; r++)
{
F0 += FE1[0][r]*w[0][r];
}// end loop over 'r'
unsigned int m = 0;
// Number of operations for primary indices = 19
for (unsigned int j = 0; j < 6; j++)
{
for (unsigned int k = 0; k < 6; k++)
{
// Number of operations to compute entry: 19
if (((((0 <= j && j < 3)) && ((0 <= k && k < 3)))))
{
A[m] += ((((K_00*FE0_D10[0][j] + K_10*FE0_D01[0][j]))*((K_00*FE0_D10[0][k] + K_10*FE0_D01[0][k])) + ((K_01*FE0_D10[0][j] + K_11*FE0_D01[0][j]))*((K_01*FE0_D10[0][k] + K_11*FE0_D01[0][k]))))*F0*W1*det;
++m;
}
}// end loop over 'k'
// Offset the entries corresponding to enriched terms
if ((((0 <= j && j < 3))))
m += offset;
}// end loop over 'j'
}
public:
/// Constructor
poisson_cell_integral_0_0( const std::vector<const pum::GenericPUM*>& pums) : ufc::cell_integral() , pums(pums)
{
// Do nothing
}
/// Destructor
virtual ~poisson_cell_integral_0_0()
{
// Do nothing
}
/// Tabulate the tensor for the contribution from a local cell
virtual void tabulate_tensor(double* A,
const double * const * w,
const ufc::cell& c) const
{
// Tabulate regular entires of element tensor
tabulate_regular_tensor(A, w, c);
// enriched local dimension of the current cell(s)
unsigned int offset = pums[0]->enriched_local_dimension(c);
if (offset == 0)
{
return;
}
// Extract vertex coordinates
const double * const * x = c.coordinates;
// Compute Jacobian of affine map from reference cell
const double J_00 = x[1][0] - x[0][0];
const double J_01 = x[2][0] - x[0][0];
const double J_10 = x[1][1] - x[0][1];
const double J_11 = x[2][1] - x[0][1];
// Compute determinant of Jacobian
double detJ = J_00*J_11 - J_01*J_10;
// Set scale factor
const double det = std::abs(detJ);
// FIXME: It will crash for multiple discontinuities, if we don't have at least one cell which all dofs are enriched
const unsigned int min_entries = 36;
const unsigned int num_entries = std::max((offset + 3)*(offset + 3), min_entries);
// Resizing and reseting auxiliary tensors
Aa.resize(num_entries);
std::fill(Aa.begin(), Aa.end(), 0.0);
// Define an array to save current quadrature point
double coordinate[2];
// Define ufc::finite_element object(s) to evalaute shape functions or their derivatives on runtime
poisson_finite_element_0 element_0;
poisson_finite_element_1 element_1;
// Array of quadrature weights.
static const double W1[1] = {0.500000000000000};
// Quadrature points on the UFC reference element: (0.333333333333333, 0.333333333333333)
// Array of quadrature points.
static const double P1[2] = \
{0.333333333333333, 0.333333333333333};
// Define vectors for quadrature points and weights(note that the sizes are determined in compile time)
std::vector<double> Wn1;
std::vector<double> Pn1;
// Check whether there is any need to use modified integration scheme
if (pums[0]->modified_quadrature(c))
{
const std::vector<double> weight1(W1, W1 + 1);
const std::vector<double> point1(P1, P1 + 2);
ConstQuadratureRule standard_gauss = std::make_pair(point1, weight1);
QuadratureRule modified_gauss;
pums[0]->cell_quadrature_rule(modified_gauss, standard_gauss, c);
Pn1 = modified_gauss.first;
Wn1 = modified_gauss.second;
}
else
{
// Map quadrature points from the reference cell to the physical cell
Wn1.resize(1);
Pn1.resize(2);
for (unsigned int i = 0; i < 1; i++)
{
Wn1[i] = W1[i];
for (unsigned int j = 0; j < 2; j++)
Pn1[2*i + j] = x[0][j]*(1.0 - P1[2*i] - P1[2*i + 1]) + x[1][j]*P1[2*i + 1] + x[2][j]*P1[2*i];
}
}
// Return the values of enriched function at the quadrature points
std::vector<double> enriched_values_1;
pums[0]->tabulate_enriched_basis(enriched_values_1, Pn1, c);
// Define auxilary indices: m, n
unsigned int m = 0;
unsigned int n = 0;
// Loop quadrature points for integral.
for (unsigned int ip = 0; ip < Wn1.size(); ip++)
{
// Evalaute tables and entries in the element tensor, if the enhanced value at this quadrature point is non-zero
if (enriched_values_1[ip] != 0)
{
// Pick up the coordinates of the current quadrature point
coordinate[0] = Pn1[2*ip];
coordinate[1] = Pn1[2*ip + 1];
// Creating a table to save the values of derivatives order 1 at the current guass point for <<CG1 on a <triangle of degree 1>> + <<CG1 on a <triangle of degree 1>>>|_{dc0}>
double value_0[2];
double table_1_D1[6][2];
for (unsigned int j = 0; j < 6; j++)
{
element_1.evaluate_basis_derivatives(j, 1, value_0, coordinate, c);
for (unsigned int k = 0; k < 2; k++)
table_1_D1[j][k] = value_0[k];
}
// Creating a table to save the values of shape functions at the current guass point for <CG1 on a <triangle of degree 1>>
double value_1[1];
double table_0_D0[3][1];
for (unsigned int j = 0; j < 3; j++)
{
element_0.evaluate_basis(j, value_1, coordinate, c);
for (unsigned int k = 0; k < 1; k++)
table_0_D0[j][k] = value_1[k];
}
// Creating a table to save the values of shape functions at the current guass point for <<CG1 on a <triangle of degree 1>> + <<CG1 on a <triangle of degree 1>>>|_{dc0}>
double value_2[1];
double table_1_D0[6][1];
for (unsigned int j = 0; j < 6; j++)
{
element_1.evaluate_basis(j, value_2, coordinate, c);
for (unsigned int k = 0; k < 1; k++)
table_1_D0[j][k] = value_2[k];
}
// Coefficient declarations.
double F0 = 0.000000000000000;
// Total number of operations to compute function values = 6
for (unsigned int r = 0; r < 3; r++)
{
F0 += table_0_D0[r][0]*w[0][r];
}// end loop over 'r'
// Number of operations for primary indices: 7
for (unsigned int j = 0; j < 6; j++)
{
for (unsigned int k = 0; k < 6; k++)
{
if (!(((0 <= j && j < 3)) && ((0 <= k && k < 3))))
{
// Move the indices of discontinuous spaces to the end of mixed space
if ((0 <= j && j < 3) && (3 <= k && k < 6))
{
m = j;
n = k;
}
else if ((3 <= j && j < 6) && (0 <= k && k < 3))
{
m = j;
n = k;
}
else if ((3 <= j && j < 6) && (3 <= k && k < 6))
{
m = j;
n = k;
}
// Number of operations to compute entry: 7
Aa[m*6 + n] += ((table_1_D1[j][1]*table_1_D1[k][1] + table_1_D1[j][0]*table_1_D1[k][0]))*F0*Wn1[ip]*det;
}// end check for enriched entiries
}// end loop over 'k'
}// end loop over 'j'
}
}// end loop over 'ip'
// Pick up entries from the total element tensor for the nodes active in the enrichment
// Determine a vector that contains the local numbering of enriched degrees of freedom in ufc::cell c for the field 0
std::vector<unsigned int> active_dofs_0;
pums[0]->tabulate_enriched_local_dofs(active_dofs_0, c);
std::vector<unsigned int>::iterator it_0_0, it_0_1;
m = 0;
for (unsigned int j = 0; j < 6; j++)
for (unsigned int k = 0; k < 6; k++)
if ((0 <= j && j < 3) && (0 <= k && k < 3))
++m;
else
{
it_0_0 = find(active_dofs_0.begin(), active_dofs_0.end(), j - 3);
it_0_1 = find(active_dofs_0.begin(), active_dofs_0.end(), k - 3);
// Check whether the entry is coressponding to the active enriched node
if (it_0_0 != active_dofs_0.end() || it_0_1 != active_dofs_0.end())
if (((0 <= j && j < 3)) || ((0 <= k && k < 3)) || (it_0_0 != active_dofs_0.end() && it_0_1 != active_dofs_0.end()))
{
A[m] = Aa[j*6 + k];
++m;
}
}
}
};
/// This class defines the interface for the tabulation of the cell
/// tensor corresponding to the local contribution to a form from
/// the integral over a cell.
class poisson_cell_integral_1_0: public ufc::cell_integral
{
const std::vector<const pum::GenericPUM*>& pums;
mutable std::vector <double> Aa;
/// Tabulate the regular entities of tensor for the contribution from a local cell
virtual void tabulate_regular_tensor(double* A,
const double * const * w,
const ufc::cell& c) const
{
// Extract vertex coordinates
const double * const * x = c.coordinates;
// Compute Jacobian of affine map from reference cell
const double J_00 = x[1][0] - x[0][0];
const double J_01 = x[2][0] - x[0][0];
const double J_10 = x[1][1] - x[0][1];
const double J_11 = x[2][1] - x[0][1];
// Compute determinant of Jacobian
double detJ = J_00*J_11 - J_01*J_10;
// Compute inverse of Jacobian
// Set scale factor
const double det = std::abs(detJ);
// Array of quadrature weights.
static const double W4[4] = {0.159020690871988, 0.090979309128011, 0.159020690871988, 0.090979309128011};
// Quadrature points on the UFC reference element: (0.178558728263616, 0.155051025721682), (0.075031110222608, 0.644948974278318), (0.666390246014701, 0.155051025721682), (0.280019915499074, 0.644948974278318)
// Value of basis functions at quadrature points.
static const double FE0[4][6] = \
{{0.666390246014701, 0.178558728263616, 0.155051025721682, 0.666390246014701, 0.178558728263616, 0.155051025721682},
{0.280019915499074, 0.075031110222608, 0.644948974278318, 0.280019915499074, 0.075031110222608, 0.644948974278318},
{0.178558728263616, 0.666390246014701, 0.155051025721682, 0.178558728263616, 0.666390246014701, 0.155051025721682},
{0.075031110222608, 0.280019915499074, 0.644948974278318, 0.075031110222608, 0.280019915499074, 0.644948974278318}};
static const double FE1[4][3] = \
{{0.666390246014701, 0.178558728263616, 0.155051025721682},
{0.280019915499074, 0.075031110222608, 0.644948974278318},
{0.178558728263616, 0.666390246014701, 0.155051025721682},
{0.075031110222608, 0.280019915499074, 0.644948974278318}};
// enriched local dimension of the current cell
unsigned int offset = pums[0]->enriched_local_dimension(c);
// Reset values in the element tensor.
const unsigned int num_entries = (3 + offset);
for (unsigned int r = 0; r < num_entries; r++)
{
A[r] = 0.000000000000000;
}// end loop over 'r'
// Compute element tensor using UFL quadrature representation
// Optimisations: ('optimisation', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False), ('ignore zero tables', False)
// Loop quadrature points for integral.
// Number of operations to compute element tensor for following IP loop = 40
for (unsigned int ip = 0; ip < 4; ip++)
{
// Coefficient declarations.
double F0 = 0.000000000000000;
// Total number of operations to compute function values = 6
for (unsigned int r = 0; r < 3; r++)
{
F0 += FE1[ip][r]*w[0][r];
}// end loop over 'r'
unsigned int m = 0;
// Number of operations for primary indices = 4
for (unsigned int j = 0; j < 6; j++)
{
// Number of operations to compute entry: 4
if (((((0 <= j && j < 3)))))
{
A[m] += FE0[ip][j]*F0*W4[ip]*det;
++m;
}
}// end loop over 'j'
}// end loop over 'ip'
}
public:
/// Constructor
poisson_cell_integral_1_0( const std::vector<const pum::GenericPUM*>& pums) : ufc::cell_integral() , pums(pums)
{
// Do nothing
}
/// Destructor
virtual ~poisson_cell_integral_1_0()
{
// Do nothing
}
/// Tabulate the tensor for the contribution from a local cell
virtual void tabulate_tensor(double* A,
const double * const * w,
const ufc::cell& c) const
{
// Tabulate regular entires of element tensor
tabulate_regular_tensor(A, w, c);
// enriched local dimension of the current cell(s)
unsigned int offset = pums[0]->enriched_local_dimension(c);
if (offset == 0)
{
return;
}
// Extract vertex coordinates
const double * const * x = c.coordinates;
// Compute Jacobian of affine map from reference cell
const double J_00 = x[1][0] - x[0][0];
const double J_01 = x[2][0] - x[0][0];
const double J_10 = x[1][1] - x[0][1];
const double J_11 = x[2][1] - x[0][1];
// Compute determinant of Jacobian
double detJ = J_00*J_11 - J_01*J_10;
// Set scale factor
const double det = std::abs(detJ);
// FIXME: It will crash for multiple discontinuities, if we don't have at least one cell which all dofs are enriched
const unsigned int min_entries = 6;
const unsigned int num_entries = std::max((offset + 3), min_entries);
// Resizing and reseting auxiliary tensors
Aa.resize(num_entries);
std::fill(Aa.begin(), Aa.end(), 0.0);
// Define an array to save current quadrature point
double coordinate[2];
// Define ufc::finite_element object(s) to evalaute shape functions or their derivatives on runtime
poisson_finite_element_0 element_0;
poisson_finite_element_1 element_1;
// Array of quadrature weights.
static const double W4[4] = {0.159020690871988, 0.090979309128011, 0.159020690871988, 0.090979309128011};
// Quadrature points on the UFC reference element: (0.178558728263616, 0.155051025721682), (0.075031110222608, 0.644948974278318), (0.666390246014701, 0.155051025721682), (0.280019915499074, 0.644948974278318)
// Array of quadrature points.
static const double P4[8] = \
{0.178558728263616, 0.155051025721682,
0.075031110222608, 0.644948974278318,
0.666390246014701, 0.155051025721682,
0.280019915499074, 0.644948974278318};
// Define vectors for quadrature points and weights(note that the sizes are determined in compile time)
std::vector<double> Wn4;
std::vector<double> Pn4;
// Check whether there is any need to use modified integration scheme
if (pums[0]->modified_quadrature(c))
{
const std::vector<double> weight4(W4, W4 + 4);
const std::vector<double> point4(P4, P4 + 8);
ConstQuadratureRule standard_gauss = std::make_pair(point4, weight4);
QuadratureRule modified_gauss;
pums[0]->cell_quadrature_rule(modified_gauss, standard_gauss, c);
Pn4 = modified_gauss.first;
Wn4 = modified_gauss.second;
}
else
{
// Map quadrature points from the reference cell to the physical cell
Wn4.resize(4);
Pn4.resize(8);
for (unsigned int i = 0; i < 4; i++)
{
Wn4[i] = W4[i];
for (unsigned int j = 0; j < 2; j++)
Pn4[2*i + j] = x[0][j]*(1.0 - P4[2*i] - P4[2*i + 1]) + x[1][j]*P4[2*i + 1] + x[2][j]*P4[2*i];
}
}
// Return the values of enriched function at the quadrature points
std::vector<double> enriched_values_4;
pums[0]->tabulate_enriched_basis(enriched_values_4, Pn4, c);
// Define an auxilary index: m
unsigned int m = 0;
// Loop quadrature points for integral.
for (unsigned int ip = 0; ip < Wn4.size(); ip++)
{
// Evalaute tables and entries in the element tensor, if the enhanced value at this quadrature point is non-zero
if (enriched_values_4[ip] != 0)
{
// Pick up the coordinates of the current quadrature point
coordinate[0] = Pn4[2*ip];
coordinate[1] = Pn4[2*ip + 1];
// Creating a table to save the values of shape functions at the current guass point for <CG1 on a <triangle of degree 1>>
double value_0[1];
double table_0_D0[3][1];
for (unsigned int j = 0; j < 3; j++)
{
element_0.evaluate_basis(j, value_0, coordinate, c);
for (unsigned int k = 0; k < 1; k++)
table_0_D0[j][k] = value_0[k];
}
// Creating a table to save the values of shape functions at the current guass point for <<CG1 on a <triangle of degree 1>> + <<CG1 on a <triangle of degree 1>>>|_{dc0}>
double value_1[1];
double table_1_D0[6][1];
for (unsigned int j = 0; j < 6; j++)
{
element_1.evaluate_basis(j, value_1, coordinate, c);
for (unsigned int k = 0; k < 1; k++)
table_1_D0[j][k] = value_1[k];
}
// Coefficient declarations.
double F0 = 0.000000000000000;
// Total number of operations to compute function values = 6
for (unsigned int r = 0; r < 3; r++)
{
F0 += table_0_D0[r][0]*w[0][r];
}// end loop over 'r'
// Number of operations for primary indices: 4
for (unsigned int j = 0; j < 6; j++)
{
if (!(((0 <= j && j < 3))))
{
// Move the indices of discontinuous spaces to the end of mixed space
if ((3 <= j && j < 6))
{
m = j;
}
// Number of operations to compute entry: 4
Aa[m] += table_1_D0[j][0]*F0*Wn4[ip]*det;
}// end check for enriched entiries
}// end loop over 'j'
}
}// end loop over 'ip'
// Pick up entries from the total element tensor for the nodes active in the enrichment
// Determine a vector that contains the local numbering of enriched degrees of freedom in ufc::cell c for the field 0
std::vector<unsigned int> active_dofs_0;
pums[0]->tabulate_enriched_local_dofs(active_dofs_0, c);
std::vector<unsigned int>::iterator it_0_0;
m = 0;
for (unsigned int j = 0; j < 6; j++)
if ((0 <= j && j < 3))
++m;
else
{
it_0_0 = find(active_dofs_0.begin(), active_dofs_0.end(), j - 3);
// Check whether the entry is coressponding to the active enriched node
if (it_0_0 != active_dofs_0.end())
{
A[m] = Aa[j];
++m;
}
}
}
};
/// This class defines the interface for the tabulation of the
/// exterior facet tensor corresponding to the local contribution to
/// a form from the integral over an exterior facet.
class poisson_exterior_facet_integral_1_0: public ufc::exterior_facet_integral
{
const std::vector<const pum::GenericPUM*>& pums;
mutable std::vector <double> Aa;
/// Tabulate the regular entities of the tensor for the contribution from a local exterior facet
virtual void tabulate_regular_tensor(double* A,
const double * const * w,
const ufc::cell& c,
unsigned int facet) const
{
// Extract vertex coordinates
const double * const * x = c.coordinates;
// Compute Jacobian of affine map from reference cell
// Compute determinant of Jacobian
// Compute inverse of Jacobian
// Get vertices on edge
static unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};
const unsigned int v0 = edge_vertices[facet][0];
const unsigned int v1 = edge_vertices[facet][1];
// Compute scale factor (length of edge scaled by length of reference interval)
const double dx0 = x[v1][0] - x[v0][0];
const double dx1 = x[v1][1] - x[v0][1];
const double det = std::sqrt(dx0*dx0 + dx1*dx1);
// Array of quadrature weights.
static const double W2[2] = {0.500000000000000, 0.500000000000000};
// Quadrature points on the UFC reference element: (0.211324865405187), (0.788675134594813)
// Value of basis functions at quadrature points.
static const double FE0_f0[2][6] = \
{{0.000000000000000, 0.788675134594813, 0.211324865405187, 0.000000000000000, 0.788675134594813, 0.211324865405187},
{0.000000000000000, 0.211324865405187, 0.788675134594813, 0.000000000000000, 0.211324865405187, 0.788675134594813}};
static const double FE0_f1[2][6] = \
{{0.788675134594813, 0.000000000000000, 0.211324865405187, 0.788675134594813, 0.000000000000000, 0.211324865405187},
{0.211324865405187, 0.000000000000000, 0.788675134594813, 0.211324865405187, 0.000000000000000, 0.788675134594813}};
static const double FE0_f2[2][6] = \
{{0.788675134594813, 0.211324865405187, 0.000000000000000, 0.788675134594813, 0.211324865405187, 0.000000000000000},
{0.211324865405187, 0.788675134594813, 0.000000000000000, 0.211324865405187, 0.788675134594813, 0.000000000000000}};
static const double FE1_f0[2][3] = \
{{0.000000000000000, 0.788675134594813, 0.211324865405187},
{0.000000000000000, 0.211324865405187, 0.788675134594813}};
static const double FE1_f1[2][3] = \
{{0.788675134594813, 0.000000000000000, 0.211324865405187},
{0.211324865405187, 0.000000000000000, 0.788675134594813}};
static const double FE1_f2[2][3] = \
{{0.788675134594813, 0.211324865405187, 0.000000000000000},
{0.211324865405187, 0.788675134594813, 0.000000000000000}};
// enriched local dimension of the current cell
unsigned int offset = pums[0]->enriched_local_dimension(c);
// Reset values in the element tensor.
const unsigned int num_entries = (3 + offset);
for (unsigned int r = 0; r < num_entries; r++)
{
A[r] = 0.000000000000000;
}// end loop over 'r'
// Compute element tensor using UFL quadrature representation
// Optimisations: ('optimisation', False), ('non zero columns', False), ('remove zero terms', False), ('ignore ones', False), ('ignore zero tables', False)
switch (facet)
{
case 0:
{
// Total number of operations to compute element tensor (from this point): 22
// Loop quadrature points for integral.
// Number of operations to compute element tensor for following IP loop = 22
for (unsigned int ip = 0; ip < 2; ip++)
{
// Coefficient declarations.
double F0 = 0.000000000000000;
// Total number of operations to compute function values = 6
for (unsigned int r = 0; r < 3; r++)
{
F0 += FE1_f0[ip][r]*w[1][r];
}// end loop over 'r'
unsigned int m = 0;
// Number of operations for primary indices = 5
for (unsigned int j = 0; j < 6; j++)
{
// Number of operations to compute entry: 5
if (((((0 <= j && j < 3)))))
{
A[m] += FE0_f0[ip][j]*F0*(-1.000000000000000)*W2[ip]*det;
++m;
}
}// end loop over 'j'
}// end loop over 'ip'
break;
}
case 1:
{
// Total number of operations to compute element tensor (from this point): 22
// Loop quadrature points for integral.
// Number of operations to compute element tensor for following IP loop = 22
for (unsigned int ip = 0; ip < 2; ip++)
{
// Coefficient declarations.
double F0 = 0.000000000000000;
// Total number of operations to compute function values = 6
for (unsigned int r = 0; r < 3; r++)
{
F0 += FE1_f1[ip][r]*w[1][r];
}// end loop over 'r'
unsigned int m = 0;
// Number of operations for primary indices = 5
for (unsigned int j = 0; j < 6; j++)
{
// Number of operations to compute entry: 5
if (((((0 <= j && j < 3)))))
{
A[m] += FE0_f1[ip][j]*F0*(-1.000000000000000)*W2[ip]*det;
++m;
}
}// end loop over 'j'
}// end loop over 'ip'
break;
}
case 2:
{
// Total number of operations to compute element tensor (from this point): 22
// Loop quadrature points for integral.
// Number of operations to compute element tensor for following IP loop = 22
for (unsigned int ip = 0; ip < 2; ip++)
{
// Coefficient declarations.
double F0 = 0.000000000000000;
// Total number of operations to compute function values = 6
for (unsigned int r = 0; r < 3; r++)
{
F0 += FE1_f2[ip][r]*w[1][r];
}// end loop over 'r'
unsigned int m = 0;
// Number of operations for primary indices = 5
for (unsigned int j = 0; j < 6; j++)
{
// Number of operations to compute entry: 5
if (((((0 <= j && j < 3)))))
{
A[m] += FE0_f2[ip][j]*F0*(-1.000000000000000)*W2[ip]*det;
++m;
}
}// end loop over 'j'
}// end loop over 'ip'
break;
}
}
}
public:
/// Constructor
poisson_exterior_facet_integral_1_0( const std::vector<const pum::GenericPUM*>& pums) : ufc::exterior_facet_integral() , pums(pums)
{
// Do nothing
}
/// Destructor
virtual ~poisson_exterior_facet_integral_1_0()
{
// Do nothing
}
/// Tabulate the tensor for the contribution from a local exterior facet
virtual void tabulate_tensor(double* A,
const double * const * w,
const ufc::cell& c,
unsigned int facet) const
{
// Tabulate regular entires of element tensor
tabulate_regular_tensor(A, w, c, facet);
// enriched local dimension of the current cell(s)
unsigned int offset = pums[0]->enriched_local_dimension(c);
if (offset == 0)
{
return;
}
// Extract vertex coordinates
const double * const * x = c.coordinates;
// Get vertices on edge
static unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};
const unsigned int v0 = edge_vertices[facet][0];
const unsigned int v1 = edge_vertices[facet][1];
// Compute scale factor (length of edge scaled by length of reference interval)
const double dx0 = x[v1][0] - x[v0][0];
const double dx1 = x[v1][1] - x[v0][1];
const double det = std::sqrt(dx0*dx0 + dx1*dx1);
// FIXME: It will crash for multiple discontinuities, if we don't have at least one cell which all dofs are enriched
const unsigned int min_entries = 6;
const unsigned int num_entries = std::max((offset + 3), min_entries);
// Resizing and reseting auxiliary tensors
Aa.resize(num_entries);
std::fill(Aa.begin(), Aa.end(), 0.0);
// Define an array to save current quadrature point
double coordinate[2];
// Define ufc::finite_element object(s) to evalaute shape functions or their derivatives on runtime
poisson_finite_element_0 element_0;
poisson_finite_element_1 element_1;
// Array of quadrature weights.
static const double W2[2] = {0.500000000000000, 0.500000000000000};
// Quadrature points on the UFC reference element: (0.211324865405187), (0.788675134594813)
// Array of quadrature points.
static const double P2[2] = \
{0.211324865405187,
0.788675134594813};
// Define vectors for quadrature points and weights(note that the sizes are determined in compile time)
std::vector<double> Wn2;
std::vector<double> Pn2;
// Check whether there is any need to use modified integration scheme
if (pums[0]->modified_quadrature(c, facet))
{
const std::vector<double> weight2(W2, W2 + 2);
const std::vector<double> point2(P2, P2 + 2);
ConstQuadratureRule standard_gauss = std::make_pair(point2, weight2);
QuadratureRule modified_gauss;
pums[0]->facet_quadrature_rule(modified_gauss, standard_gauss, c, facet);
Pn2 = modified_gauss.first;
Wn2 = modified_gauss.second;
}
else
{
// Map quadrature points from the reference cell to the physical cell
Wn2.resize(2);
Pn2.resize(4);
for (unsigned int i = 0; i < 2; i++)
{
Wn2[i] = W2[i];
for (unsigned int j = 0; j < 2; j++)
Pn2[2*i + j] = x[v0][j]*(1.0 - P2[i]) + x[v1][j]*P2[i];
}
}
// Return the values of enriched function at the quadrature points
std::vector<double> enriched_values_2;
pums[0]->tabulate_enriched_basis(enriched_values_2, Pn2, c);
// Define an auxilary index: m
unsigned int m = 0;
// Loop quadrature points for integral.
for (unsigned int ip = 0; ip < Wn2.size(); ip++)
{
// Evalaute tables and entries in the element tensor, if the enhanced value at this quadrature point is non-zero
if (enriched_values_2[ip] != 0)
{
// Pick up the coordinates of the current quadrature point
coordinate[0] = Pn2[2*ip];
coordinate[1] = Pn2[2*ip + 1];
// Creating a table to save the values of shape functions at the current guass point for <CG1 on a <triangle of degree 1>>
double value_0[1];
double table_0_D0[3][1];
for (unsigned int j = 0; j < 3; j++)
{
element_0.evaluate_basis(j, value_0, coordinate, c);
for (unsigned int k = 0; k < 1; k++)
table_0_D0[j][k] = value_0[k];
}
// Creating a table to save the values of shape functions at the current guass point for <<CG1 on a <triangle of degree 1>> + <<CG1 on a <triangle of degree 1>>>|_{dc0}>
double value_1[1];
double table_1_D0[6][1];
for (unsigned int j = 0; j < 6; j++)
{
element_1.evaluate_basis(j, value_1, coordinate, c);
for (unsigned int k = 0; k < 1; k++)
table_1_D0[j][k] = value_1[k];
}
// Coefficient declarations.
double F0 = 0.000000000000000;
// Total number of operations to compute function values = 6
for (unsigned int r = 0; r < 3; r++)
{
F0 += table_0_D0[r][0]*w[1][r];
}// end loop over 'r'
// Number of operations for primary indices: 5
for (unsigned int j = 0; j < 6; j++)
{
if (!(((0 <= j && j < 3))))
{
// Move the indices of discontinuous spaces to the end of mixed space
if ((3 <= j && j < 6))
{
m = j;
}
// Number of operations to compute entry: 5
Aa[m] += table_1_D0[j][0]*F0*(-1.000000000000000)*Wn2[ip]*det;
}// end check for enriched entiries
}// end loop over 'j'
}
}// end loop over 'ip'
// Pick up entries from the total element tensor for the nodes active in the enrichment
// Determine a vector that contains the local numbering of enriched degrees of freedom in ufc::cell c for the field 0
std::vector<unsigned int> active_dofs_0;
pums[0]->tabulate_enriched_local_dofs(active_dofs_0, c);
std::vector<unsigned int>::iterator it_0_0;
m = 0;
for (unsigned int j = 0; j < 6; j++)
if ((0 <= j && j < 3))
++m;
else
{
it_0_0 = find(active_dofs_0.begin(), active_dofs_0.end(), j - 3);
// Check whether the entry is coressponding to the active enriched node
if (it_0_0 != active_dofs_0.end())
{
A[m] = Aa[j];
++m;
}
}
}
};
/// This class defines the interface for the assembly of the global
/// tensor corresponding to a form with r + n arguments, that is, a
/// mapping
///
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
///
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
/// global tensor A is defined by
///
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
///
/// where each argument Vj represents the application to the
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
/// fixed functions (coefficients).
class poisson_form_0: public ufc::form
{
const std::vector<const pum::GenericPUM*>& pums;
public:
/// Constructor
poisson_form_0( const std::vector<const pum::GenericPUM*>& pums) : ufc::form() , pums(pums)
{
// Do nothing
}
/// Destructor
virtual ~poisson_form_0()
{
// Do nothing
}
/// Return a string identifying the form
virtual const char* signature() const
{
return "Form([Integral(Product(Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), IndexSum(Product(Indexed(ComponentTensor(SpatialDerivative(Argument(EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None))), 0), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ComponentTensor(SpatialDerivative(Argument(EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None))), 1), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(2),), {Index(2): 2})), MultiIndex((Index(1),), {Index(1): 2}))), MultiIndex((Index(1),), {Index(1): 2}))), Measure('cell', 0, None))])";
}
/// Return the rank of the global tensor (r)
virtual unsigned int rank() const
{
return 2;
}
/// Return the number of coefficients (n)
virtual unsigned int num_coefficients() const
{
return 1;
}
/// Return the number of cell integrals
virtual unsigned int num_cell_integrals() const
{
return 1;
}
/// Return the number of exterior facet integrals
virtual unsigned int num_exterior_facet_integrals() const
{
return 0;
}
/// Return the number of interior facet integrals
virtual unsigned int num_interior_facet_integrals() const
{
return 0;
}
/// Create a new finite element for argument function i
virtual ufc::finite_element* create_finite_element(unsigned int i) const
{
switch (i)
{
case 0:
{
return new poisson_finite_element_1();
break;
}
case 1:
{
return new poisson_finite_element_1();
break;
}
case 2:
{
return new poisson_finite_element_0();
break;
}
}
return 0;
}
/// Create a new dof map for argument function i
virtual ufc::dof_map* create_dof_map(unsigned int i) const
{
switch (i)
{
case 0:
{
return new poisson_dof_map_1(pums);
break;
}
case 1:
{
return new poisson_dof_map_1(pums);
break;
}
case 2:
{
return new poisson_dof_map_0();
break;
}
}
return 0;
}
/// Create a new cell integral on sub domain i
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
{
switch (i)
{
case 0:
{
return new poisson_cell_integral_0_0(pums);
break;
}
}
return 0;
}
/// Create a new exterior facet integral on sub domain i
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
{
return 0;
}
/// Create a new interior facet integral on sub domain i
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
{
return 0;
}
};
/// This class defines the interface for the assembly of the global
/// tensor corresponding to a form with r + n arguments, that is, a
/// mapping
///
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
///
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
/// global tensor A is defined by
///
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
///
/// where each argument Vj represents the application to the
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
/// fixed functions (coefficients).
class poisson_form_1: public ufc::form
{
const std::vector<const pum::GenericPUM*>& pums;
public:
/// Constructor
poisson_form_1( const std::vector<const pum::GenericPUM*>& pums) : ufc::form() , pums(pums)
{
// Do nothing
}
/// Destructor
virtual ~poisson_form_1()
{
// Do nothing
}
/// Return a string identifying the form
virtual const char* signature() const
{
return "Form([Integral(Product(Argument(EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None))), 0), Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0)), Measure('cell', 0, None)), Integral(Product(IntValue(-1, (), (), {}), Product(Argument(EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None))), 0), Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 1))), Measure('exterior_facet', 0, None))])";
}
/// Return the rank of the global tensor (r)
virtual unsigned int rank() const
{
return 1;
}
/// Return the number of coefficients (n)
virtual unsigned int num_coefficients() const
{
return 2;
}
/// Return the number of cell integrals
virtual unsigned int num_cell_integrals() const
{
return 1;
}
/// Return the number of exterior facet integrals
virtual unsigned int num_exterior_facet_integrals() const
{
return 1;
}
/// Return the number of interior facet integrals
virtual unsigned int num_interior_facet_integrals() const
{
return 0;
}
/// Create a new finite element for argument function i
virtual ufc::finite_element* create_finite_element(unsigned int i) const
{
switch (i)
{
case 0:
{
return new poisson_finite_element_1();
break;
}
case 1:
{
return new poisson_finite_element_0();
break;
}
case 2:
{
return new poisson_finite_element_0();
break;
}
}
return 0;
}
/// Create a new dof map for argument function i
virtual ufc::dof_map* create_dof_map(unsigned int i) const
{
switch (i)
{
case 0:
{
return new poisson_dof_map_1(pums);
break;
}
case 1:
{
return new poisson_dof_map_0();
break;
}
case 2:
{
return new poisson_dof_map_0();
break;
}
}
return 0;
}
/// Create a new cell integral on sub domain i
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
{
switch (i)
{
case 0:
{
return new poisson_cell_integral_1_0(pums);
break;
}
}
return 0;
}
/// Create a new exterior facet integral on sub domain i
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
{
switch (i)
{
case 0:
{
return new poisson_exterior_facet_integral_1_0(pums);
break;
}
}
return 0;
}
/// Create a new interior facet integral on sub domain i
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
{
return 0;
}
};
/// This class defines the interface for a finite element.
class poisson_auxiliary_finite_element_0: public ufc::finite_element
{
public:
/// Constructor
poisson_auxiliary_finite_element_0() : ufc::finite_element()
{
// Do nothing
}
/// Destructor
virtual ~poisson_auxiliary_finite_element_0()
{
// Do nothing
}
/// Return a string identifying the finite element
virtual const char* signature() const
{
return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
}
/// Return the cell shape
virtual ufc::shape cell_shape() const
{
return ufc::triangle;
}
/// Return the dimension of the finite element function space
virtual unsigned int space_dimension() const
{
return 3;
}
/// Return the rank of the value space
virtual unsigned int value_rank() const
{
return 0;
}
/// Return the dimension of the value space for axis i
virtual unsigned int value_dimension(unsigned int i) const
{
return 1;
}
/// Evaluate basis function i at given point in cell
virtual void evaluate_basis(unsigned int i,
double* values,
const double* coordinates,
const ufc::cell& c) const
{
// Extract vertex coordinates
const double * const * x = c.coordinates;
// Compute Jacobian of affine map from reference cell
const double J_00 = x[1][0] - x[0][0];
const double J_01 = x[2][0] - x[0][0];
const double J_10 = x[1][1] - x[0][1];
const double J_11 = x[2][1] - x[0][1];
// Compute determinant of Jacobian
double detJ = J_00*J_11 - J_01*J_10;
// Compute inverse of Jacobian
// Compute constants
const double C0 = x[1][0] + x[2][0];
const double C1 = x[1][1] + x[2][1];
// Get coordinates and map to the reference (FIAT) element
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
// Reset values.
*values = 0.000000000000000;
// Map degree of freedom to element degree of freedom
const unsigned int dof = i;
// Array of basisvalues.
double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
// Declare helper variables.
unsigned int rr = 0;
unsigned int ss = 0;
double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
// Compute basisvalues.
basisvalues[0] = 1.000000000000000;
basisvalues[1] = tmp0;
for (unsigned int r = 0; r < 1; r++)
{
rr = (r + 1)*(r + 1 + 1)/2 + 1;
ss = r*(r + 1)/2;
basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
}// end loop over 'r'
for (unsigned int r = 0; r < 2; r++)
{
for (unsigned int s = 0; s < 2 - r; s++)
{
rr = (r + s)*(r + s + 1)/2 + s;
basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
}// end loop over 's'
}// end loop over 'r'
// Table(s) of coefficients.
static const double coefficients0[3][3] = \
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.000000000000000, 0.333333333333333}};
// Compute value(s).
for (unsigned int r = 0; r < 3; r++)
{
*values += coefficients0[dof][r]*basisvalues[r];
}// end loop over 'r'
}
/// Evaluate all basis functions at given point in cell
virtual void evaluate_basis_all(double* values,
const double* coordinates,
const ufc::cell& c) const
{
// Helper variable to hold values of a single dof.
double dof_values = 0.000000000000000;
// Loop dofs and call evaluate_basis.
for (unsigned int r = 0; r < 3; r++)
{
evaluate_basis(r, &dof_values, coordinates, c);
values[r] = dof_values;
}// end loop over 'r'
}
/// Evaluate order n derivatives of basis function i at given point in cell
virtual void evaluate_basis_derivatives(unsigned int i,
unsigned int n,
double* values,
const double* coordinates,
const ufc::cell& c) const
{
// Extract vertex coordinates
const double * const * x = c.coordinates;
// Compute Jacobian of affine map from reference cell
const double J_00 = x[1][0] - x[0][0];
const double J_01 = x[2][0] - x[0][0];
const double J_10 = x[1][1] - x[0][1];
const double J_11 = x[2][1] - x[0][1];
// Compute determinant of Jacobian
double detJ = J_00*J_11 - J_01*J_10;
// Compute inverse of Jacobian
const double K_00 = J_11 / detJ;
const double K_01 = -J_01 / detJ;
const double K_10 = -J_10 / detJ;
const double K_11 = J_00 / detJ;
// Compute constants
const double C0 = x[1][0] + x[2][0];
const double C1 = x[1][1] + x[2][1];
// Get coordinates and map to the reference (FIAT) element
double X = (J_01*(C1 - 2.0*coordinates[1]) + J_11*(2.0*coordinates[0] - C0)) / detJ;
double Y = (J_00*(2.0*coordinates[1] - C1) + J_10*(C0 - 2.0*coordinates[0])) / detJ;
// Compute number of derivatives.
unsigned int num_derivatives = 1;
for (unsigned int r = 0; r < n; r++)
{
num_derivatives *= 2;
}// end loop over 'r'
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
unsigned int **combinations = new unsigned int *[num_derivatives];
for (unsigned int row = 0; row < num_derivatives; row++)
{
combinations[row] = new unsigned int [n];
for (unsigned int col = 0; col < n; col++)
combinations[row][col] = 0;
}
// Generate combinations of derivatives
for (unsigned int row = 1; row < num_derivatives; row++)
{
for (unsigned int num = 0; num < row; num++)
{
for (unsigned int col = n-1; col+1 > 0; col--)
{
if (combinations[row][col] + 1 > 1)
combinations[row][col] = 0;
else
{
combinations[row][col] += 1;
break;
}
}
}
}
// Compute inverse of Jacobian
const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
// Declare transformation matrix
// Declare pointer to two dimensional array and initialise
double **transform = new double *[num_derivatives];
for (unsigned int j = 0; j < num_derivatives; j++)
{
transform[j] = new double [num_derivatives];
for (unsigned int k = 0; k < num_derivatives; k++)
transform[j][k] = 1;
}
// Construct transformation matrix
for (unsigned int row = 0; row < num_derivatives; row++)
{
for (unsigned int col = 0; col < num_derivatives; col++)
{
for (unsigned int k = 0; k < n; k++)
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
}
}
// Reset values. Assuming that values is always an array.
for (unsigned int r = 0; r < num_derivatives; r++)
{
values[r] = 0.000000000000000;
}// end loop over 'r'
// Map degree of freedom to element degree of freedom
const unsigned int dof = i;
// Array of basisvalues.
double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
// Declare helper variables.
unsigned int rr = 0;
unsigned int ss = 0;
double tmp0 = (1.000000000000000 + Y + 2.000000000000000*X)/2.000000000000000;
// Compute basisvalues.
basisvalues[0] = 1.000000000000000;
basisvalues[1] = tmp0;
for (unsigned int r = 0; r < 1; r++)
{
rr = (r + 1)*(r + 1 + 1)/2 + 1;
ss = r*(r + 1)/2;
basisvalues[rr] = basisvalues[ss]*(0.500000000000000 + r + Y*(1.500000000000000 + r));
}// end loop over 'r'
for (unsigned int r = 0; r < 2; r++)
{
for (unsigned int s = 0; s < 2 - r; s++)
{
rr = (r + s)*(r + s + 1)/2 + s;
basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s));
}// end loop over 's'
}// end loop over 'r'
// Table(s) of coefficients.
static const double coefficients0[3][3] = \
{{0.471404520791032, -0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.288675134594813, -0.166666666666667},
{0.471404520791032, 0.000000000000000, 0.333333333333333}};
// Tables of derivatives of the polynomial base (transpose).
static const double dmats0[3][3] = \
{{0.000000000000000, 0.000000000000000, 0.000000000000000},
{4.898979485566356, 0.000000000000000, 0.000000000000000},
{0.000000000000000, 0.000000000000000, 0.000000000000000}};
static const double dmats1[3][3] = \
{{0.000000000000000, 0.000000000000000, 0.000000000000000},
{2.449489742783178, 0.000000000000000, 0.000000000000000},
{4.242640687119285, 0.000000000000000, 0.000000000000000}};
// Compute reference derivatives.
// Declare pointer to array of derivatives on FIAT element.
double *derivatives = new double[num_derivatives];
for (unsigned int r = 0; r < num_derivatives; r++)
{
derivatives[r] = 0.000000000000000;
}// end loop over 'r'
// Declare derivative matrix (of polynomial basis).
double dmats[3][3] = \
{{1.000000000000000, 0.000000000000000, 0.000000000000000},
{0.000000000000000, 1.000000000000000, 0.000000000000000},
{0.000000000000000, 0.000000000000000, 1.000000000000000}};
// Declare (auxiliary) derivative matrix (of polynomial basis).
double dmats_old[3][3] = \
{{1.000000000000000, 0.000000000000000, 0.000000000000000},
{0.000000000000000, 1.000000000000000, 0.000000000000000},
{0.000000000000000, 0.000000000000000, 1.000000000000000}};
// Loop possible derivatives.
for (unsigned int r = 0; r < num_derivatives; r++)
{
// Resetting dmats values to compute next derivative.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats[t][u] = 0.000000000000000;
if (t == u)
{
dmats[t][u] = 1.000000000000000;
}
}// end loop over 'u'
}// end loop over 't'
// Looping derivative order to generate dmats.
for (unsigned int s = 0; s < n; s++)
{
// Updating dmats_old with new values and resetting dmats.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats_old[t][u] = dmats[t][u];
dmats[t][u] = 0.000000000000000;
}// end loop over 'u'
}// end loop over 't'
// Update dmats using an inner product.
if (combinations[r][s] == 0)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
}// end loop over 'tu'
}// end loop over 'u'
}// end loop over 't'
}
if (combinations[r][s] == 1)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
}// end loop over 'tu'
}// end loop over 'u'
}// end loop over 't'
}
}// end loop over 's'
for (unsigned int s = 0; s < 3; s++)
{
for (unsigned int t = 0; t < 3; t++)
{
derivatives[r] += coefficients0[dof][s]*dmats[s][t]*basisvalues[t];
}// end loop over 't'
}// end loop over 's'
}// end loop over 'r'
// Transform derivatives back to physical element
for (unsigned int r = 0; r < num_derivatives; r++)
{
for (unsigned int s = 0; s < num_derivatives; s++)
{
values[r] += transform[r][s]*derivatives[s];
}// end loop over 's'
}// end loop over 'r'
// Delete pointer to array of derivatives on FIAT element
delete [] derivatives;
// Delete pointer to array of combinations of derivatives and transform
for (unsigned int r = 0; r < num_derivatives; r++)
{
delete [] combinations[r];
}// end loop over 'r'
delete [] combinations;
for (unsigned int r = 0; r < num_derivatives; r++)
{
delete [] transform[r];
}// end loop over 'r'
delete [] transform;
}
/// Evaluate order n derivatives of all basis functions at given point in cell
virtual void evaluate_basis_derivatives_all(unsigned int n,
double* values,
const double* coordinates,
const ufc::cell& c) const
{
// Compute number of derivatives.
unsigned int num_derivatives = 1;
for (unsigned int r = 0; r < n; r++)
{
num_derivatives *= 2;
}// end loop over 'r'
// Helper variable to hold values of a single dof.
double *dof_values = new double[num_derivatives];
for (unsigned int r = 0; r < num_derivatives; r++)
{
dof_values[r] = 0.000000000000000;
}// end loop over 'r'
// Loop dofs and call evaluate_basis_derivatives.
for (unsigned int r = 0; r < 3; r++)
{
evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
for (unsigned int s = 0; s < num_derivatives; s++)
{
values[r*num_derivatives + s] = dof_values[s];
}// end loop over 's'
}// end loop over 'r'
// Delete pointer.
delete [] dof_values;
}
/// Evaluate linear functional for dof i on the function f
virtual double evaluate_dof(unsigned int i,
const ufc::function& f,
const ufc::cell& c) const
{
// Declare variables for result of evaluation.
double vals[1];
// Declare variable for physical coordinates.
double y[2];
const double * const * x = c.coordinates;
switch (i)
{
case 0:
{
y[0] = x[0][0];
y[1] = x[0][1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
case 1:
{
y[0] = x[1][0];
y[1] = x[1][1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
case 2:
{
y[0] = x[2][0];
y[1] = x[2][1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
}
return 0.000000000000000;
}
/// Evaluate linear functionals for all dofs on the function f
virtual void evaluate_dofs(double* values,
const ufc::function& f,
const ufc::cell& c) const
{
// Declare variables for result of evaluation.
double vals[1];
// Declare variable for physical coordinates.
double y[2];
const double * const * x = c.coordinates;
y[0] = x[0][0];
y[1] = x[0][1];
f.evaluate(vals, y, c);
values[0] = vals[0];
y[0] = x[1][0];
y[1] = x[1][1];
f.evaluate(vals, y, c);
values[1] = vals[0];
y[0] = x[2][0];
y[1] = x[2][1];
f.evaluate(vals, y, c);
values[2] = vals[0];
}
/// Interpolate vertex values from dof values
virtual void interpolate_vertex_values(double* vertex_values,
const double* dof_values,
const ufc::cell& c) const
{
// Evaluate function and change variables
vertex_values[0] = dof_values[0];
vertex_values[1] = dof_values[1];
vertex_values[2] = dof_values[2];
}
/// Return the number of sub elements (for a mixed element)
virtual unsigned int num_sub_elements() const
{
return 0;
}
/// Create a new finite element for sub element i (for a mixed element)
virtual ufc::finite_element* create_sub_element(unsigned int i) const
{
return 0;
}
};
/// This class defines the interface for a local-to-global mapping of
/// degrees of freedom (dofs).
class poisson_auxiliary_dof_map_0: public ufc::dof_map
{
private:
unsigned int _global_dimension;
public:
/// Constructor
poisson_auxiliary_dof_map_0() :ufc::dof_map()
{
_global_dimension = 0;
}
/// Destructor
virtual ~poisson_auxiliary_dof_map_0()
{
// Do nothing
}
/// Return a string identifying the dof map
virtual const char* signature() const
{
return "FFC dofmap for FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
}
/// Return true iff mesh entities of topological dimension d are needed
virtual bool needs_mesh_entities(unsigned int d) const
{
switch (d)
{
case 0:
{
return true;
break;
}
case 1:
{
return false;
break;
}
case 2:
{
return false;
break;
}
}
return false;
}
/// Initialize dof map for mesh (return true iff init_cell() is needed)
virtual bool init_mesh(const ufc::mesh& m)
{
_global_dimension = m.num_entities[0];
return false;
}
/// Initialize dof map for given cell
virtual void init_cell(const ufc::mesh& m,
const ufc::cell& c)
{
// Do nothing
}
/// Finish initialization of dof map for cells
virtual void init_cell_finalize()
{
// Do nothing
}
/// Return the dimension of the global finite element function space
virtual unsigned int global_dimension() const
{
return _global_dimension;
}
/// Return the dimension of the local finite element function space for a cell
virtual unsigned int local_dimension(const ufc::cell& c) const
{
return 3;
}
/// Return the maximum dimension of the local finite element function space
virtual unsigned int max_local_dimension() const
{
return 3;
}
// Return the geometric dimension of the coordinates this dof map provides
virtual unsigned int geometric_dimension() const
{
return 2;
}
/// Return the number of dofs on each cell facet
virtual unsigned int num_facet_dofs() const
{
return 2;
}
/// Return the number of dofs associated with each cell entity of dimension d
virtual unsigned int num_entity_dofs(unsigned int d) const
{
switch (d)
{
case 0:
{
return 1;
break;
}
case 1:
{
return 0;
break;
}
case 2:
{
return 0;
break;
}
}
return 0;
}
/// Tabulate the local-to-global mapping of dofs on a cell
virtual void tabulate_dofs(unsigned int* dofs,
const ufc::mesh& m,
const ufc::cell& c) const
{
dofs[0] = c.entity_indices[0][0];
dofs[1] = c.entity_indices[0][1];
dofs[2] = c.entity_indices[0][2];
}
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
virtual void tabulate_facet_dofs(unsigned int* dofs,
unsigned int facet) const
{
switch (facet)
{
case 0:
{
dofs[0] = 1;
dofs[1] = 2;
break;
}
case 1:
{
dofs[0] = 0;
dofs[1] = 2;
break;
}
case 2:
{
dofs[0] = 0;
dofs[1] = 1;
break;
}
}
}
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
virtual void tabulate_entity_dofs(unsigned int* dofs,
unsigned int d, unsigned int i) const
{
if (d > 2)
{
throw std::runtime_error("d is larger than dimension (2)");
}
switch (d)
{
case 0:
{
if (i > 2)
{
throw std::runtime_error("i is larger than number of entities (2)");
}
switch (i)
{
case 0:
{
dofs[0] = 0;
break;
}
case 1:
{
dofs[0] = 1;
break;
}
case 2:
{
dofs[0] = 2;
break;
}
}
break;
}
case 1:
{
break;
}
case 2:
{
break;
}
}
}
/// Tabulate the coordinates of all dofs on a cell
virtual void tabulate_coordinates(double** coordinates,
const ufc::cell& c) const
{
const double * const * x = c.coordinates;
coordinates[0][0] = x[0][0];
coordinates[0][1] = x[0][1];
coordinates[1][0] = x[1][0];
coordinates[1][1] = x[1][1];
coordinates[2][0] = x[2][0];
coordinates[2][1] = x[2][1];
}
/// Return the number of sub dof maps (for a mixed element)
virtual unsigned int num_sub_dof_maps() const
{
return 0;
}
/// Create a new dof_map for sub dof map i (for a mixed element)
virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
{
return 0;
}
};
/// This class defines the interface for the assembly of the global
/// tensor corresponding to a form with r + n arguments, that is, a
/// mapping
///
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
///
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
/// global tensor A is defined by
///
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
///
/// where each argument Vj represents the application to the
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
/// fixed functions (coefficients).
class poisson_auxiliary_form_0: public ufc::form
{
public:
/// Constructor
poisson_auxiliary_form_0() : ufc::form()
{
// Do nothing
}
/// Destructor
virtual ~poisson_auxiliary_form_0()
{
// Do nothing
}
/// Return a string identifying the form
virtual const char* signature() const
{
return "auxiliary form";
}
/// Return the rank of the global tensor (r)
virtual unsigned int rank() const
{
return 0;
}
/// Return the number of coefficients (n)
virtual unsigned int num_coefficients() const
{
return 0;
}
/// Return the number of cell integrals
virtual unsigned int num_cell_integrals() const
{
return 0;
}
/// Return the number of exterior facet integrals
virtual unsigned int num_exterior_facet_integrals() const
{
return 0;
}
/// Return the number of interior facet integrals
virtual unsigned int num_interior_facet_integrals() const
{
return 0;
}
/// Create a new finite element for argument function i
virtual ufc::finite_element* create_finite_element(unsigned int i) const
{
switch (i)
{
case 0:
{
return new poisson_finite_element_0();
break;
}
}
return 0;
}
/// Create a new dof map for argument function i
virtual ufc::dof_map* create_dof_map(unsigned int i) const
{
switch (i)
{
case 0:
{
return new poisson_auxiliary_dof_map_0();
break;
}
}
return 0;
}
/// Create a new cell integral on sub domain i
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
{
return 0;
}
/// Create a new exterior facet integral on sub domain i
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
{
return 0;
}
/// Create a new interior facet integral on sub domain i
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
{
return 0;
}
};
/// This class defines the interface for the assembly of the global
/// tensor corresponding to a form with r + n arguments, that is, a
/// mapping
///
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
///
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
/// global tensor A is defined by
///
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
///
/// where each argument Vj represents the application to the
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
/// fixed functions (coefficients).
class poisson_auxiliary_form_1: public ufc::form
{
public:
/// Constructor
poisson_auxiliary_form_1() : ufc::form()
{
// Do nothing
}
/// Destructor
virtual ~poisson_auxiliary_form_1()
{
// Do nothing
}
/// Return a string identifying the form
virtual const char* signature() const
{
return "auxiliary form";
}
/// Return the rank of the global tensor (r)
virtual unsigned int rank() const
{
return 0;
}
/// Return the number of coefficients (n)
virtual unsigned int num_coefficients() const
{
return 0;
}
/// Return the number of cell integrals
virtual unsigned int num_cell_integrals() const
{
return 0;
}
/// Return the number of exterior facet integrals
virtual unsigned int num_exterior_facet_integrals() const
{
return 0;
}
/// Return the number of interior facet integrals
virtual unsigned int num_interior_facet_integrals() const
{
return 0;
}
/// Create a new finite element for argument function i
virtual ufc::finite_element* create_finite_element(unsigned int i) const
{
switch (i)
{
case 0:
{
return new poisson_finite_element_0();
break;
}
}
return 0;
}
/// Create a new dof map for argument function i
virtual ufc::dof_map* create_dof_map(unsigned int i) const
{
switch (i)
{
case 0:
{
return new poisson_auxiliary_dof_map_0();
break;
}
}
return 0;
}
/// Create a new cell integral on sub domain i
virtual ufc::cell_integral* create_cell_integral(unsigned int i) const
{
return 0;
}
/// Create a new exterior facet integral on sub domain i
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(unsigned int i) const
{
return 0;
}
/// Create a new interior facet integral on sub domain i
virtual ufc::interior_facet_integral* create_interior_facet_integral(unsigned int i) const
{
return 0;
}
};
/// This class defines the interface for post-processing on vector x
/// to obtain x0, u and j where,
///
/// - x is the solution vector containing standard and enriched degrees of freedom
/// defined on continuous/discontinuous space
/// - u is the standard part of solution vector defined on continuous space
/// - j is the enriched part pf solution vector defined on continuous space
/// - x0 is the result vector, equall to u + j, defined on continuous space
/// by considering enrichement function
//
// Dolfin includes
#include <dolfin/common/NoDeleter.h>
#include <dolfin/mesh/Mesh.h>
#include <dolfin/fem/DofMap.h>
#include <dolfin/la/GenericVector.h>
// PartitionOfUnity includes
#include <pum/PostProcess.h>
#include <pum/FunctionSpace.h>
#include <pum/PUM.h>
#include <pum/GenericSurface.h>
namespace Poisson
{
class PostProcess: public pum::PostProcess
{
const dolfin::Mesh& mesh;
pum::FunctionSpace& function_space;
public:
/// Constructor
PostProcess(pum::FunctionSpace& function_space):
pum::PostProcess(function_space.mesh()), mesh(function_space.mesh()),
function_space(function_space)
{
// Do nothing
}
/// Destructor
~PostProcess()
{
// Do nothing
}
/// Return a string identifying the underling element
const char* signature() const
{
return "Interpolating results to the continuous space of EnrichedElement(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), ElementRestriction(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), Measure('surface', 0, None)))";
}
/// Obtain result vector 'x0' from solution vector 'x'
void interpolate(const dolfin::GenericVector& x, dolfin::GenericVector& x0) const
{
// Compute number of standard dofs for field 0
dolfin::DofMap dof_map_0(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh);
unsigned int num_standard_dofs_0 = dof_map_0.global_dimension();
double value;
double h;
/// Selecting standard degrees of freedom related to field 0 from the solution vector
double* values_0 = new double[num_standard_dofs_0];
unsigned int* positions_0 = new unsigned int [num_standard_dofs_0];
for (unsigned int i = 0; i < num_standard_dofs_0; ++i)
positions_0[i] = i;
x.get(values_0, num_standard_dofs_0, positions_0);
x0.set(values_0, num_standard_dofs_0, positions_0);
x0.apply("insert");
const std::vector<const pum::GenericPUM*>& pums = function_space.pums;
/// Selecting enriched degrees of freedom related to field 0 from the solution vector
std::vector <std::vector<unsigned int> > enhanced_dof_maps_0;
enhanced_dof_maps_0.resize(num_standard_dofs_0);
std::vector<unsigned int> enhanced_dof_values_0;
enhanced_dof_values_0.resize(num_standard_dofs_0);
compute_enhanced_dof_maps(*pums[0], dof_map_0, enhanced_dof_maps_0);
compute_enhanced_dof_values(*pums[0], dof_map_0, enhanced_dof_values_0);
//compute_enhanced_dof_maps(pum_0, dof_map_0, enhanced_dof_maps_0);
//compute_enhanced_dof_values(pum_0, dof_map_0, enhanced_dof_values_0);
for (unsigned int i = 0; i != num_standard_dofs_0; ++i)
{
unsigned int pos = i;
for (std::vector<unsigned int>::const_iterator it = enhanced_dof_maps_0[i].begin();
it != enhanced_dof_maps_0[i].end(); ++it)
{
h = enhanced_dof_values_0[i];
unsigned int pos_n = *it + num_standard_dofs_0;
x.get(&value, 1, &pos_n);
value *= h;
x0.add(&value, 1, &pos);
}
}
// memory clean up
delete[] values_0;
delete[] positions_0;
x0.apply("add");
}
/// Obtain continuous u and discontinuous j parts of solution vector 'x'
void interpolate(const dolfin::GenericVector& x, dolfin::GenericVector& u, dolfin::GenericVector& j) const
{
// Compute number of standard dofs for field 0
dolfin::DofMap dof_map_0(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh);
unsigned int num_standard_dofs_0 = dof_map_0.global_dimension();
double value;
/// Selecting standard degrees of freedom related to field 0 from the solution vector
double* values_0 = new double[num_standard_dofs_0];
unsigned int* positions_0 = new unsigned int [num_standard_dofs_0];
for (unsigned int i = 0; i < num_standard_dofs_0; ++i)
positions_0[i] = i;
x.get(values_0, num_standard_dofs_0, positions_0);
u.set(values_0, num_standard_dofs_0, positions_0);
const std::vector<const pum::GenericPUM*>& pums = function_space.pums;
/// Selecting enriched degrees of freedom related to field 0 from the solution vector
std::vector <std::vector<unsigned int> > enhanced_dof_maps_0;
enhanced_dof_maps_0.resize(num_standard_dofs_0);
// Compute enhanced dof maps
compute_enhanced_dof_maps(*pums[0], dof_map_0, enhanced_dof_maps_0);
for (unsigned int i = 0; i != num_standard_dofs_0; ++i)
{
unsigned int pos = i ;
for (std::vector<unsigned int>::const_iterator it = enhanced_dof_maps_0[i].begin();
it != enhanced_dof_maps_0[i].end(); ++it)
{
unsigned int pos_n = *it + num_standard_dofs_0;
x.get(&value, 1, &pos_n);
j.set(&value, 1, &pos);
}
}
// memory clean up
delete[] values_0;
delete[] positions_0;
u.apply("insert");
j.apply("insert");
}
};
}
// DOLFIN wrappers
// Standard library includes
#include <string>
// DOLFIN includes
#include <dolfin/common/NoDeleter.h>
#include <dolfin/fem/FiniteElement.h>
#include <dolfin/fem/DofMap.h>
#include <dolfin/fem/Form.h>
#include <dolfin/function/FunctionSpace.h>
#include <dolfin/function/Function.h>
#include <dolfin/function/GenericFunction.h>
#include <dolfin/function/CoefficientAssigner.h>
namespace Poisson
{
class CoefficientSpace_f: public dolfin::FunctionSpace
{
public:
CoefficientSpace_f(const dolfin::Mesh & mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
(new poisson_dof_map_0()), mesh)))
{
// Do nothing
}
CoefficientSpace_f(dolfin::Mesh & mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh)))
{
// Do nothing
}
CoefficientSpace_f(boost::shared_ptr<dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), *mesh)))
{
// Do nothing
}
CoefficientSpace_f(boost::shared_ptr<const dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), *mesh)))
{
// Do nothing
}
~CoefficientSpace_f()
{
}
};
class CoefficientSpace_g: public dolfin::FunctionSpace
{
public:
CoefficientSpace_g(const dolfin::Mesh & mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
(new poisson_dof_map_0()), mesh)))
{
// Do nothing
}
CoefficientSpace_g(dolfin::Mesh & mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh)))
{
// Do nothing
}
CoefficientSpace_g(boost::shared_ptr<dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), *mesh)))
{
// Do nothing
}
CoefficientSpace_g(boost::shared_ptr<const dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), *mesh)))
{
// Do nothing
}
~CoefficientSpace_g()
{
}
};
class CoefficientSpace_w: public dolfin::FunctionSpace
{
public:
CoefficientSpace_w(const dolfin::Mesh & mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
(new poisson_dof_map_0()), mesh)))
{
// Do nothing
}
CoefficientSpace_w(dolfin::Mesh & mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh)))
{
// Do nothing
}
CoefficientSpace_w(boost::shared_ptr<dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), *mesh)))
{
// Do nothing
}
CoefficientSpace_w(boost::shared_ptr<const dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), *mesh)))
{
// Do nothing
}
~CoefficientSpace_w()
{
}
};
class Form_0_FunctionSpace_0: public pum::FunctionSpace
{
public:
Form_0_FunctionSpace_0(const dolfin::Mesh & mesh ,const std::vector<const pum::GenericPUM*>& pums):
pum::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
(new poisson_dof_map_1(pums)), mesh)) ,pums)
{
// Do nothing
}
Form_0_FunctionSpace_0(dolfin::Mesh & mesh ,const std::vector<const pum::GenericPUM*>& pums):
pum::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), mesh)) ,pums)
{
// Do nothing
}
Form_0_FunctionSpace_0(boost::shared_ptr<dolfin::Mesh> mesh ,const std::vector<const pum::GenericPUM*>& pums):
pum::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), *mesh)) ,pums)
{
// Do nothing
}
Form_0_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh ,const std::vector<const pum::GenericPUM*>& pums):
pum::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), *mesh)) ,pums)
{
// Do nothing
}
~Form_0_FunctionSpace_0()
{
}
};
class Form_0_FunctionSpace_1: public pum::FunctionSpace
{
public:
Form_0_FunctionSpace_1(const dolfin::Mesh & mesh ,const std::vector<const pum::GenericPUM*>& pums):
pum::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
(new poisson_dof_map_1(pums)), mesh)) ,pums)
{
// Do nothing
}
Form_0_FunctionSpace_1(dolfin::Mesh & mesh ,const std::vector<const pum::GenericPUM*>& pums):
pum::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), mesh)) ,pums)
{
// Do nothing
}
Form_0_FunctionSpace_1(boost::shared_ptr<dolfin::Mesh> mesh ,const std::vector<const pum::GenericPUM*>& pums):
pum::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), *mesh)) ,pums)
{
// Do nothing
}
Form_0_FunctionSpace_1(boost::shared_ptr<const dolfin::Mesh> mesh ,const std::vector<const pum::GenericPUM*>& pums):
pum::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), *mesh)) ,pums)
{
// Do nothing
}
~Form_0_FunctionSpace_1()
{
}
};
typedef CoefficientSpace_w Form_0_FunctionSpace_2;
class Form_0: public dolfin::Form
{
public:
// Constructor
Form_0(const pum::FunctionSpace& V0, const pum::FunctionSpace& V1):
dolfin::Form(2, 1), w(*this, 0)
{
_function_spaces[0] = reference_to_no_delete_pointer(V0);
_function_spaces[1] = reference_to_no_delete_pointer(V1);
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_0(V0.pums));
}
// Constructor
Form_0(const pum::FunctionSpace& V0, const pum::FunctionSpace& V1, dolfin::GenericFunction & w):
dolfin::Form(2, 1), w(*this, 0)
{
_function_spaces[0] = reference_to_no_delete_pointer(V0);
_function_spaces[1] = reference_to_no_delete_pointer(V1);
this->w = w;
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_0(V0.pums));
}
// Constructor
Form_0(const pum::FunctionSpace& V0, const pum::FunctionSpace& V1, boost::shared_ptr<dolfin::GenericFunction> w):
dolfin::Form(2, 1), w(*this, 0)
{
_function_spaces[0] = reference_to_no_delete_pointer(V0);
_function_spaces[1] = reference_to_no_delete_pointer(V1);
this->w = *w;
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_0(V0.pums));
}
// Constructor
Form_0(boost::shared_ptr<const pum::FunctionSpace> V0, boost::shared_ptr<const pum::FunctionSpace> V1):
dolfin::Form(2, 1), w(*this, 0)
{
_function_spaces[0] = V0;
_function_spaces[1] = V1;
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_0(V0->pums));
}
// Constructor
Form_0(boost::shared_ptr<const pum::FunctionSpace> V0, boost::shared_ptr<const pum::FunctionSpace> V1, dolfin::GenericFunction & w):
dolfin::Form(2, 1), w(*this, 0)
{
_function_spaces[0] = V0;
_function_spaces[1] = V1;
this->w = w;
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_0(V0->pums));
}
// Constructor
Form_0(boost::shared_ptr<const pum::FunctionSpace> V0, boost::shared_ptr<const pum::FunctionSpace> V1, boost::shared_ptr<dolfin::GenericFunction> w):
dolfin::Form(2, 1), w(*this, 0)
{
_function_spaces[0] = V0;
_function_spaces[1] = V1;
this->w = *w;
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_0(V0->pums));
}
// Destructor
~Form_0()
{}
/// Return the number of the coefficient with this name
virtual dolfin::uint coefficient_number(const std::string& name) const
{
if(name == "w") return 0;
dolfin::error("Invalid coefficient.");
return 0;
}
/// Return the name of the coefficient with this number
virtual std::string coefficient_name(dolfin::uint i) const
{
switch(i)
{
case 0: return "w";
}
dolfin::error("Invalid coefficient.");
return "unnamed";
}
// Typedefs
typedef Form_0_FunctionSpace_0 TestSpace;
typedef Form_0_FunctionSpace_1 TrialSpace;
typedef Form_0_FunctionSpace_2 CoefficientSpace_w;
// Coefficients
dolfin::CoefficientAssigner w;
};
class Form_1_FunctionSpace_0: public pum::FunctionSpace
{
public:
Form_1_FunctionSpace_0(const dolfin::Mesh & mesh ,const std::vector<const pum::GenericPUM*>& pums):
pum::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
(new poisson_dof_map_1(pums)), mesh)) ,pums)
{
// Do nothing
}
Form_1_FunctionSpace_0(dolfin::Mesh & mesh ,const std::vector<const pum::GenericPUM*>& pums):
pum::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), mesh)) ,pums)
{
// Do nothing
}
Form_1_FunctionSpace_0(boost::shared_ptr<dolfin::Mesh> mesh ,const std::vector<const pum::GenericPUM*>& pums):
pum::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), *mesh)) ,pums)
{
// Do nothing
}
Form_1_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh ,const std::vector<const pum::GenericPUM*>& pums):
pum::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_finite_element_1()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_1(pums)), *mesh)) ,pums)
{
// Do nothing
}
~Form_1_FunctionSpace_0()
{
}
};
typedef CoefficientSpace_f Form_1_FunctionSpace_1;
typedef CoefficientSpace_g Form_1_FunctionSpace_2;
class Form_1: public dolfin::Form
{
public:
// Constructor
Form_1(const pum::FunctionSpace& V0):
dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
{
_function_spaces[0] = reference_to_no_delete_pointer(V0);
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_1(V0.pums));
}
// Constructor
Form_1(const pum::FunctionSpace& V0, dolfin::GenericFunction & f, dolfin::GenericFunction & g):
dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
{
_function_spaces[0] = reference_to_no_delete_pointer(V0);
this->f = f;
this->g = g;
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_1(V0.pums));
}
// Constructor
Form_1(const pum::FunctionSpace& V0, boost::shared_ptr<dolfin::GenericFunction> f, boost::shared_ptr<dolfin::GenericFunction> g):
dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
{
_function_spaces[0] = reference_to_no_delete_pointer(V0);
this->f = *f;
this->g = *g;
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_1(V0.pums));
}
// Constructor
Form_1(boost::shared_ptr<const pum::FunctionSpace> V0):
dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
{
_function_spaces[0] = V0;
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_1(V0->pums));
}
// Constructor
Form_1(boost::shared_ptr<const pum::FunctionSpace> V0, dolfin::GenericFunction & f, dolfin::GenericFunction & g):
dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
{
_function_spaces[0] = V0;
this->f = f;
this->g = g;
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_1(V0->pums));
}
// Constructor
Form_1(boost::shared_ptr<const pum::FunctionSpace> V0, boost::shared_ptr<dolfin::GenericFunction> f, boost::shared_ptr<dolfin::GenericFunction> g):
dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
{
_function_spaces[0] = V0;
this->f = *f;
this->g = *g;
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_form_1(V0->pums));
}
// Destructor
~Form_1()
{}
/// Return the number of the coefficient with this name
virtual dolfin::uint coefficient_number(const std::string& name) const
{
if(name == "f") return 0;
else if(name == "g") return 1;
dolfin::error("Invalid coefficient.");
return 0;
}
/// Return the name of the coefficient with this number
virtual std::string coefficient_name(dolfin::uint i) const
{
switch(i)
{
case 0: return "f";
case 1: return "g";
}
dolfin::error("Invalid coefficient.");
return "unnamed";
}
// Typedefs
typedef Form_1_FunctionSpace_0 TestSpace;
typedef Form_1_FunctionSpace_1 CoefficientSpace_f;
typedef Form_1_FunctionSpace_2 CoefficientSpace_g;
// Coefficients
dolfin::CoefficientAssigner f;
dolfin::CoefficientAssigner g;
};
// Class typedefs
typedef Form_0 BilinearForm;
typedef Form_1 LinearForm;
typedef Form_0::TestSpace FunctionSpace;
} // namespace Poisson
// DOLFIN wrappers
// Standard library includes
#include <string>
// DOLFIN includes
#include <dolfin/common/NoDeleter.h>
#include <dolfin/fem/FiniteElement.h>
#include <dolfin/fem/DofMap.h>
#include <dolfin/fem/Form.h>
#include <dolfin/function/FunctionSpace.h>
#include <dolfin/function/Function.h>
#include <dolfin/function/GenericFunction.h>
#include <dolfin/function/CoefficientAssigner.h>
namespace Poisson
{
class Form_auxiliary_0_FunctionSpace_auxiliary_0: public dolfin::FunctionSpace
{
public:
Form_auxiliary_0_FunctionSpace_auxiliary_0(const dolfin::Mesh & mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
(new poisson_auxiliary_dof_map_0()), mesh)))
{
// Do nothing
}
Form_auxiliary_0_FunctionSpace_auxiliary_0(dolfin::Mesh & mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_auxiliary_dof_map_0()), mesh)))
{
// Do nothing
}
Form_auxiliary_0_FunctionSpace_auxiliary_0(boost::shared_ptr<dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_auxiliary_dof_map_0()), *mesh)))
{
// Do nothing
}
Form_auxiliary_0_FunctionSpace_auxiliary_0(boost::shared_ptr<const dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_auxiliary_dof_map_0()), *mesh)))
{
// Do nothing
}
~Form_auxiliary_0_FunctionSpace_auxiliary_0()
{
}
};
class Form_auxiliary_0: public dolfin::Form
{
public:
// Constructor
Form_auxiliary_0(const dolfin::FunctionSpace& V0):
dolfin::Form(1, 0)
{
_function_spaces[0] = reference_to_no_delete_pointer(V0);
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_auxiliary_form_0());
}
// Constructor
Form_auxiliary_0(boost::shared_ptr<const dolfin::FunctionSpace> V0):
dolfin::Form(1, 0)
{
_function_spaces[0] = V0;
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_auxiliary_form_0());
}
// Destructor
~Form_auxiliary_0()
{}
/// Return the number of the coefficient with this name
virtual dolfin::uint coefficient_number(const std::string& name) const
{
dolfin::error("No coefficients.");
return 0;
}
/// Return the name of the coefficient with this number
virtual std::string coefficient_name(dolfin::uint i) const
{
dolfin::error("No coefficients.");
return "unnamed";
}
// Typedefs
typedef Form_auxiliary_0_FunctionSpace_auxiliary_0 TestSpace;
// Coefficients
};
class Form_auxiliary_1_FunctionSpace_auxiliary_0: public dolfin::FunctionSpace
{
public:
Form_auxiliary_1_FunctionSpace_auxiliary_0(const dolfin::Mesh & mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement
(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>
(new poisson_auxiliary_dof_map_0()), mesh)))
{
// Do nothing
}
Form_auxiliary_1_FunctionSpace_auxiliary_0(dolfin::Mesh & mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_auxiliary_dof_map_0()), mesh)))
{
// Do nothing
}
Form_auxiliary_1_FunctionSpace_auxiliary_0(boost::shared_ptr<dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_auxiliary_dof_map_0()), *mesh)))
{
// Do nothing
}
Form_auxiliary_1_FunctionSpace_auxiliary_0(boost::shared_ptr<const dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new poisson_auxiliary_finite_element_0()))),
boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_auxiliary_dof_map_0()), *mesh)))
{
// Do nothing
}
~Form_auxiliary_1_FunctionSpace_auxiliary_0()
{
}
};
class Form_auxiliary_1: public dolfin::Form
{
public:
// Constructor
Form_auxiliary_1(const dolfin::FunctionSpace& V0):
dolfin::Form(1, 0)
{
_function_spaces[0] = reference_to_no_delete_pointer(V0);
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_auxiliary_form_1());
}
// Constructor
Form_auxiliary_1(boost::shared_ptr<const dolfin::FunctionSpace> V0):
dolfin::Form(1, 0)
{
_function_spaces[0] = V0;
_ufc_form = boost::shared_ptr<const ufc::form>(new poisson_auxiliary_form_1());
}
// Destructor
~Form_auxiliary_1()
{}
/// Return the number of the coefficient with this name
virtual dolfin::uint coefficient_number(const std::string& name) const
{
dolfin::error("No coefficients.");
return 0;
}
/// Return the name of the coefficient with this number
virtual std::string coefficient_name(dolfin::uint i) const
{
dolfin::error("No coefficients.");
return "unnamed";
}
// Typedefs
typedef Form_auxiliary_1_FunctionSpace_auxiliary_0 TestSpace;
// Coefficients
};
// Class typedefs
typedef Form_auxiliary_0::TestSpace FunctionSpace_auxiliary;
} // namespace Poisson
/// This class defines the interface for computing enriched function space
/// and auxiliary function spaces
///
//
// Dolfin includes
#include <dolfin/common/NoDeleter.h>
#include <dolfin/mesh/Mesh.h>
#include <dolfin/fem/DofMap.h>
// PartitionOfUnity includes
#include <pum/FunctionSpace.h>
#include <pum/PUM.h>
#include <pum/GenericSurface.h>
namespace Poisson
{
class FunctionSpaces
{
dolfin::Mesh& mesh;
std::vector<const pum::GenericSurface*>& surfaces;
std::vector<const pum::GenericPUM*> pums;
pum::PUM* pum_0;
dolfin::DofMap* dof_map_0;
/// Build PUM objects to intialize enriched FunctionSpace
void init()
{
// Create DofMap instance for the field 0
//dolfin::DofMap dof_map_0(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh);
dof_map_0 = new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new poisson_dof_map_0()), mesh);
// Create PUM objects
pum_0 = new pum::PUM(surfaces, mesh, *dof_map_0);
// Push PUM objects to a vector
pums.push_back(pum_0);
}
public:
/// Constructor
FunctionSpaces(dolfin::Mesh& mesh, std::vector<const pum::GenericSurface*>& surfaces):
mesh(mesh), surfaces(surfaces)
{
init();
}
/// Destructor
~FunctionSpaces()
{
delete pum_0;
delete dof_map_0;
}
/// Return Enriched Function Space
pum::FunctionSpace space()
{
return FunctionSpace(mesh, pums);
}
/// Return Enriched Function Space
dolfin::FunctionSpace auxiliary_space()
{
return FunctionSpace_auxiliary(mesh);
}
/// Update PUM objects while surfaces evolve
void update()
{
for (unsigned int i = 0; i < pums.size(); ++i)
const_cast<pum::GenericPUM*>(pums[i])->update();
}
};
}
#endif
# Copyright (C) 2008-2009 Mehdi Nikbakht and Garth N. Wells.
# Licensed under the GNU GPL Version 3.0 or any later version.
#
# The bilinear form a(v, u) and linear form L(v) for
# Poisson's equation with discontinuities.
#
# Compile this form with FFC: ffc-pum -l dolfin Poisson.ufl
#
elem_cont = FiniteElement("CG", triangle, 1)
elem_discont = ElementRestriction(elem_cont, dc) # or ec[dc]
#element = elem_cont * elem_discont
#(vc, vd) = TestFunctions(element)
#(uc, ud) = TrialFunctions(element)
#v = vc + vd
#u = uc + ud
element = elem_cont + elem_discont
v = TestFunction(element)
u = TrialFunction(element)
k = Constant(triangle)
f = Coefficient(elem_cont)
w = Coefficient(elem_cont)
g = Coefficient(elem_cont)
a = w*dot(grad(v), grad(u))*dx
L = v*f*dx - v*g*ds
Follow ups
References
-
Re: evaluate_integrand
From: Anders Logg, 2010-04-12
-
Re: evaluate_integrand
From: Garth N. Wells, 2010-04-12
-
Re: evaluate_integrand
From: Anders Logg, 2010-04-12
-
Re: evaluate_integrand
From: Garth N. Wells, 2010-04-12
-
Re: evaluate_integrand
From: Anders Logg, 2010-04-12
-
Re: evaluate_integrand
From: Garth N. Wells, 2010-04-12
-
Re: evaluate_integrand
From: Anders Logg, 2010-04-13
-
Re: evaluate_integrand
From: Mehdi Nikbakht, 2010-04-13
-
Re: evaluate_integrand
From: Anders Logg, 2010-04-13
-
Re: evaluate_integrand
From: Anders Logg, 2010-04-13
