fenics team mailing list archive
-
fenics team
-
Mailing list archive
-
Message #01223
Re: Patch for changed Python interface of Expressions
Patch applied
--
Anders
On Mon, Sep 06, 2010 at 08:44:59AM -0700, Johan Hake wrote:
> I also draged in some changes done in the C++ source code of la/eigenvalue.
>
> Johan
> # Bazaar merge directive format 2 (Bazaar 0.90)
> # revision_id: hake.dev@xxxxxxxxx-20100906154049-idue489k83pizp0n
> # target_branch: bzr+ssh://bazaar.launchpad.net/~fenics-core/fenics-\
> # doc/main/
> # testament_sha1: 624bff4ed92bd2ceeef59bda34f7782bf777a0e4
> # timestamp: 2010-09-06 08:42:44 -0700
> # base_revision_id: k.b.oelgaard@xxxxxxxxx-20100902125002-\
> # ciqdciv4nc1e145w
> #
> # Begin patch
> === modified file 'source/demos/la/eigenvalue/cpp/StiffnessMatrix.h'
> --- source/demos/la/eigenvalue/cpp/StiffnessMatrix.h 2010-09-01 15:13:40 +0000
> +++ source/demos/la/eigenvalue/cpp/StiffnessMatrix.h 2010-09-06 15:40:49 +0000
> @@ -1,5 +1,5 @@
> -// This code conforms with the UFC specification version 1.4
> -// and was automatically generated by FFC version 0.9.2+.
> +// This code conforms with the UFC specification version 1.4.2
> +// and was automatically generated by FFC version 0.9.4+.
> //
> // This code was generated with the option '-l dolfin' and
> // contains DOLFIN-specific wrappers that depend on DOLFIN.
> @@ -9,12 +9,13 @@
> // cache_dir: ''
> // convert_exceptions_to_warnings: False
> // cpp_optimize: False
> +// cpp_optimize_flags: '-O2'
> // epsilon: 1e-14
> // form_postfix: True
> // format: 'dolfin'
> -// log_level: 10
> +// log_level: 20
> // log_prefix: ''
> -// optimize: True
> +// optimize: False
> // output_dir: '.'
> // precision: 15
> // quadrature_degree: 'auto'
> @@ -51,19 +52,19 @@
> /// Return a string identifying the finite element
> virtual const char* signature() const
> {
> - return "FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1)";
> + return "FiniteElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1)";
> }
>
> /// Return the cell shape
> virtual ufc::shape cell_shape() const
> {
> - return ufc::triangle;
> + return ufc::tetrahedron;
> }
>
> /// Return the dimension of the finite element function space
> virtual unsigned int space_dimension() const
> {
> - return 3;
> + return 4;
> }
>
> /// Return the rank of the value space
> @@ -90,21 +91,40 @@
> // 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_02 = x[3][0] - x[0][0];
> const double J_10 = x[1][1] - x[0][1];
> const double J_11 = x[2][1] - x[0][1];
> + const double J_12 = x[3][1] - x[0][1];
> + const double J_20 = x[1][2] - x[0][2];
> + const double J_21 = x[2][2] - x[0][2];
> + const double J_22 = x[3][2] - x[0][2];
> +
> + // Compute sub determinants
> + const double d_00 = J_11*J_22 - J_12*J_21;
> + const double d_01 = J_12*J_20 - J_10*J_22;
> + const double d_02 = J_10*J_21 - J_11*J_20;
> + const double d_10 = J_02*J_21 - J_01*J_22;
> + const double d_11 = J_00*J_22 - J_02*J_20;
> + const double d_12 = J_01*J_20 - J_00*J_21;
> + const double d_20 = J_01*J_12 - J_02*J_11;
> + const double d_21 = J_02*J_10 - J_00*J_12;
> + const double d_22 = J_00*J_11 - J_01*J_10;
>
> // Compute determinant of Jacobian
> - double detJ = J_00*J_11 - J_01*J_10;
> + double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
>
> // Compute inverse of Jacobian
>
> // Compute constants
> - const double C0 = x[1][0] + x[2][0];
> - const double C1 = x[1][1] + x[2][1];
> + const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0];
> + const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1];
> + const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2];
>
> // 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;
> + double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ;
> + double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ;
> + double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ;
> +
>
> // Reset values.
> *values = 0.000000000000000;
> @@ -114,37 +134,49 @@
> {
>
> // Array of basisvalues.
> - double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
> + double basisvalues[4] = {0.000000000000000, 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;
> + double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
>
> // 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));
> + rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
> + ss = r*(r + 1)*(r + 2)/6;
> + basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
> + }// end loop over 'r'
> + for (unsigned int r = 0; r < 1; r++)
> + {
> + for (unsigned int s = 0; s < 1 - r; s++)
> + {
> + rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
> + ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
> + basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
> + }// end loop over 's'
> }// 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));
> + for (unsigned int t = 0; t < 2 - r - s; t++)
> + {
> + rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
> + basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
> + }// end loop over 't'
> }// end loop over 's'
> }// end loop over 'r'
>
> // Table(s) of coefficients.
> - static const double coefficients0[3] = \
> - {0.471404520791032, -0.288675134594813, -0.166666666666667};
> + static const double coefficients0[4] = \
> + {0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993};
>
> // Compute value(s).
> - for (unsigned int r = 0; r < 3; r++)
> + for (unsigned int r = 0; r < 4; r++)
> {
> *values += coefficients0[r]*basisvalues[r];
> }// end loop over 'r'
> @@ -154,37 +186,49 @@
> {
>
> // Array of basisvalues.
> - double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
> + double basisvalues[4] = {0.000000000000000, 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;
> + double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
>
> // 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));
> + rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
> + ss = r*(r + 1)*(r + 2)/6;
> + basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
> + }// end loop over 'r'
> + for (unsigned int r = 0; r < 1; r++)
> + {
> + for (unsigned int s = 0; s < 1 - r; s++)
> + {
> + rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
> + ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
> + basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
> + }// end loop over 's'
> }// 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));
> + for (unsigned int t = 0; t < 2 - r - s; t++)
> + {
> + rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
> + basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
> + }// end loop over 't'
> }// end loop over 's'
> }// end loop over 'r'
>
> // Table(s) of coefficients.
> - static const double coefficients0[3] = \
> - {0.471404520791032, 0.288675134594813, -0.166666666666667};
> + static const double coefficients0[4] = \
> + {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993};
>
> // Compute value(s).
> - for (unsigned int r = 0; r < 3; r++)
> + for (unsigned int r = 0; r < 4; r++)
> {
> *values += coefficients0[r]*basisvalues[r];
> }// end loop over 'r'
> @@ -194,37 +238,101 @@
> {
>
> // 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] = \
> - {0.471404520791032, 0.000000000000000, 0.333333333333333};
> -
> - // Compute value(s).
> - for (unsigned int r = 0; r < 3; r++)
> + double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000};
> +
> + // Declare helper variables.
> + unsigned int rr = 0;
> + unsigned int ss = 0;
> + double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
> +
> + // Compute basisvalues.
> + basisvalues[0] = 1.000000000000000;
> + basisvalues[1] = tmp0;
> + for (unsigned int r = 0; r < 1; r++)
> + {
> + rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
> + ss = r*(r + 1)*(r + 2)/6;
> + basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
> + }// end loop over 'r'
> + for (unsigned int r = 0; r < 1; r++)
> + {
> + for (unsigned int s = 0; s < 1 - r; s++)
> + {
> + rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
> + ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
> + basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
> + }// end loop over 's'
> + }// end loop over 'r'
> + for (unsigned int r = 0; r < 2; r++)
> + {
> + for (unsigned int s = 0; s < 2 - r; s++)
> + {
> + for (unsigned int t = 0; t < 2 - r - s; t++)
> + {
> + rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
> + basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
> + }// end loop over 't'
> + }// end loop over 's'
> + }// end loop over 'r'
> +
> + // Table(s) of coefficients.
> + static const double coefficients0[4] = \
> + {0.288675134594813, 0.000000000000000, 0.210818510677892, -0.074535599249993};
> +
> + // Compute value(s).
> + for (unsigned int r = 0; r < 4; r++)
> + {
> + *values += coefficients0[r]*basisvalues[r];
> + }// end loop over 'r'
> + break;
> + }
> + case 3:
> + {
> +
> + // Array of basisvalues.
> + double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000};
> +
> + // Declare helper variables.
> + unsigned int rr = 0;
> + unsigned int ss = 0;
> + double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
> +
> + // Compute basisvalues.
> + basisvalues[0] = 1.000000000000000;
> + basisvalues[1] = tmp0;
> + for (unsigned int r = 0; r < 1; r++)
> + {
> + rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
> + ss = r*(r + 1)*(r + 2)/6;
> + basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
> + }// end loop over 'r'
> + for (unsigned int r = 0; r < 1; r++)
> + {
> + for (unsigned int s = 0; s < 1 - r; s++)
> + {
> + rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
> + ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
> + basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
> + }// end loop over 's'
> + }// end loop over 'r'
> + for (unsigned int r = 0; r < 2; r++)
> + {
> + for (unsigned int s = 0; s < 2 - r; s++)
> + {
> + for (unsigned int t = 0; t < 2 - r - s; t++)
> + {
> + rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
> + basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
> + }// end loop over 't'
> + }// end loop over 's'
> + }// end loop over 'r'
> +
> + // Table(s) of coefficients.
> + static const double coefficients0[4] = \
> + {0.288675134594813, 0.000000000000000, 0.000000000000000, 0.223606797749979};
> +
> + // Compute value(s).
> + for (unsigned int r = 0; r < 4; r++)
> {
> *values += coefficients0[r]*basisvalues[r];
> }// end loop over 'r'
> @@ -243,7 +351,7 @@
> double dof_values = 0.000000000000000;
>
> // Loop dofs and call evaluate_basis.
> - for (unsigned int r = 0; r < 3; r++)
> + for (unsigned int r = 0; r < 4; r++)
> {
> evaluate_basis(r, &dof_values, coordinates, c);
> values[r] = dof_values;
> @@ -263,31 +371,55 @@
> // 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_02 = x[3][0] - x[0][0];
> const double J_10 = x[1][1] - x[0][1];
> const double J_11 = x[2][1] - x[0][1];
> + const double J_12 = x[3][1] - x[0][1];
> + const double J_20 = x[1][2] - x[0][2];
> + const double J_21 = x[2][2] - x[0][2];
> + const double J_22 = x[3][2] - x[0][2];
> +
> + // Compute sub determinants
> + const double d_00 = J_11*J_22 - J_12*J_21;
> + const double d_01 = J_12*J_20 - J_10*J_22;
> + const double d_02 = J_10*J_21 - J_11*J_20;
> + const double d_10 = J_02*J_21 - J_01*J_22;
> + const double d_11 = J_00*J_22 - J_02*J_20;
> + const double d_12 = J_01*J_20 - J_00*J_21;
> + const double d_20 = J_01*J_12 - J_02*J_11;
> + const double d_21 = J_02*J_10 - J_00*J_12;
> + const double d_22 = J_00*J_11 - J_01*J_10;
>
> // Compute determinant of Jacobian
> - double detJ = J_00*J_11 - J_01*J_10;
> + double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
>
> // 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;
> + const double K_00 = d_00 / detJ;
> + const double K_01 = d_10 / detJ;
> + const double K_02 = d_20 / detJ;
> + const double K_10 = d_01 / detJ;
> + const double K_11 = d_11 / detJ;
> + const double K_12 = d_21 / detJ;
> + const double K_20 = d_02 / detJ;
> + const double K_21 = d_12 / detJ;
> + const double K_22 = d_22 / detJ;
>
> // Compute constants
> - const double C0 = x[1][0] + x[2][0];
> - const double C1 = x[1][1] + x[2][1];
> + const double C0 = x[3][0] + x[2][0] + x[1][0] - x[0][0];
> + const double C1 = x[3][1] + x[2][1] + x[1][1] - x[0][1];
> + const double C2 = x[3][2] + x[2][2] + x[1][2] - x[0][2];
>
> // 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;
> + double X = (d_00*(2.0*coordinates[0] - C0) + d_10*(2.0*coordinates[1] - C1) + d_20*(2.0*coordinates[2] - C2)) / detJ;
> + double Y = (d_01*(2.0*coordinates[0] - C0) + d_11*(2.0*coordinates[1] - C1) + d_21*(2.0*coordinates[2] - C2)) / detJ;
> + double Z = (d_02*(2.0*coordinates[0] - C0) + d_12*(2.0*coordinates[1] - C1) + d_22*(2.0*coordinates[2] - C2)) / detJ;
> +
>
> // Compute number of derivatives.
> unsigned int num_derivatives = 1;
> for (unsigned int r = 0; r < n; r++)
> {
> - num_derivatives *= 2;
> + num_derivatives *= 3;
> }// end loop over 'r'
>
> // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
> @@ -306,7 +438,7 @@
> {
> for (unsigned int col = n-1; col+1 > 0; col--)
> {
> - if (combinations[row][col] + 1 > 1)
> + if (combinations[row][col] + 1 > 2)
> combinations[row][col] = 0;
> else
> {
> @@ -318,7 +450,7 @@
> }
>
> // Compute inverse of Jacobian
> - const double Jinv[2][2] = {{K_00, K_01}, {K_10, K_11}};
> + const double Jinv[3][3] = {{K_00, K_01, K_02}, {K_10, K_11, K_12}, {K_20, K_21, K_22}};
>
> // Declare transformation matrix
> // Declare pointer to two dimensional array and initialise
> @@ -353,45 +485,65 @@
> {
>
> // Array of basisvalues.
> - double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
> + double basisvalues[4] = {0.000000000000000, 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;
> + double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
>
> // 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));
> + rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
> + ss = r*(r + 1)*(r + 2)/6;
> + basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
> + }// end loop over 'r'
> + for (unsigned int r = 0; r < 1; r++)
> + {
> + for (unsigned int s = 0; s < 1 - r; s++)
> + {
> + rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
> + ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
> + basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
> + }// end loop over 's'
> }// 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));
> + for (unsigned int t = 0; t < 2 - r - s; t++)
> + {
> + rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
> + basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
> + }// end loop over 't'
> }// end loop over 's'
> }// end loop over 'r'
>
> // Table(s) of coefficients.
> - static const double coefficients0[3] = \
> - {0.471404520791032, -0.288675134594813, -0.166666666666667};
> + static const double coefficients0[4] = \
> + {0.288675134594813, -0.182574185835055, -0.105409255338946, -0.074535599249993};
>
> // 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}};
> + static const double dmats0[4][4] = \
> + {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
> +
> + static const double dmats1[4][4] = \
> + {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
> +
> + static const double dmats2[4][4] = \
> + {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
>
> // Compute reference derivatives.
> // Declare pointer to array of derivatives on FIAT element.
> @@ -402,24 +554,26 @@
> }// 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}};
> + double dmats[4][4] = \
> + {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.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}};
> + double dmats_old[4][4] = \
> + {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.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 t = 0; t < 4; t++)
> {
> - for (unsigned int u = 0; u < 3; u++)
> + for (unsigned int u = 0; u < 4; u++)
> {
> dmats[t][u] = 0.000000000000000;
> if (t == u)
> @@ -434,9 +588,9 @@
> 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 t = 0; t < 4; t++)
> {
> - for (unsigned int u = 0; u < 3; u++)
> + for (unsigned int u = 0; u < 4; u++)
> {
> dmats_old[t][u] = dmats[t][u];
> dmats[t][u] = 0.000000000000000;
> @@ -446,11 +600,11 @@
> // Update dmats using an inner product.
> if (combinations[r][s] == 0)
> {
> - for (unsigned int t = 0; t < 3; t++)
> + for (unsigned int t = 0; t < 4; t++)
> {
> - for (unsigned int u = 0; u < 3; u++)
> + for (unsigned int u = 0; u < 4; u++)
> {
> - for (unsigned int tu = 0; tu < 3; tu++)
> + for (unsigned int tu = 0; tu < 4; tu++)
> {
> dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
> }// end loop over 'tu'
> @@ -460,11 +614,11 @@
>
> if (combinations[r][s] == 1)
> {
> - for (unsigned int t = 0; t < 3; t++)
> + for (unsigned int t = 0; t < 4; t++)
> {
> - for (unsigned int u = 0; u < 3; u++)
> + for (unsigned int u = 0; u < 4; u++)
> {
> - for (unsigned int tu = 0; tu < 3; tu++)
> + for (unsigned int tu = 0; tu < 4; tu++)
> {
> dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
> }// end loop over 'tu'
> @@ -472,10 +626,24 @@
> }// end loop over 't'
> }
>
> + if (combinations[r][s] == 2)
> + {
> + for (unsigned int t = 0; t < 4; t++)
> + {
> + for (unsigned int u = 0; u < 4; u++)
> + {
> + for (unsigned int tu = 0; tu < 4; tu++)
> + {
> + dmats[t][u] += dmats2[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 s = 0; s < 4; s++)
> {
> - for (unsigned int t = 0; t < 3; t++)
> + for (unsigned int t = 0; t < 4; t++)
> {
> derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
> }// end loop over 't'
> @@ -511,45 +679,65 @@
> {
>
> // Array of basisvalues.
> - double basisvalues[3] = {0.000000000000000, 0.000000000000000, 0.000000000000000};
> + double basisvalues[4] = {0.000000000000000, 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;
> + double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
>
> // 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));
> + rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
> + ss = r*(r + 1)*(r + 2)/6;
> + basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
> + }// end loop over 'r'
> + for (unsigned int r = 0; r < 1; r++)
> + {
> + for (unsigned int s = 0; s < 1 - r; s++)
> + {
> + rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
> + ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
> + basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
> + }// end loop over 's'
> }// 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));
> + for (unsigned int t = 0; t < 2 - r - s; t++)
> + {
> + rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
> + basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
> + }// end loop over 't'
> }// end loop over 's'
> }// end loop over 'r'
>
> // Table(s) of coefficients.
> - static const double coefficients0[3] = \
> - {0.471404520791032, 0.288675134594813, -0.166666666666667};
> + static const double coefficients0[4] = \
> + {0.288675134594813, 0.182574185835055, -0.105409255338946, -0.074535599249993};
>
> // 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}};
> + static const double dmats0[4][4] = \
> + {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
> +
> + static const double dmats1[4][4] = \
> + {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
> +
> + static const double dmats2[4][4] = \
> + {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {5.163977794943223, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
>
> // Compute reference derivatives.
> // Declare pointer to array of derivatives on FIAT element.
> @@ -560,24 +748,26 @@
> }// 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}};
> + double dmats[4][4] = \
> + {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.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}};
> + double dmats_old[4][4] = \
> + {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.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 t = 0; t < 4; t++)
> {
> - for (unsigned int u = 0; u < 3; u++)
> + for (unsigned int u = 0; u < 4; u++)
> {
> dmats[t][u] = 0.000000000000000;
> if (t == u)
> @@ -592,9 +782,9 @@
> 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 t = 0; t < 4; t++)
> {
> - for (unsigned int u = 0; u < 3; u++)
> + for (unsigned int u = 0; u < 4; u++)
> {
> dmats_old[t][u] = dmats[t][u];
> dmats[t][u] = 0.000000000000000;
> @@ -604,11 +794,11 @@
> // Update dmats using an inner product.
> if (combinations[r][s] == 0)
> {
> - for (unsigned int t = 0; t < 3; t++)
> + for (unsigned int t = 0; t < 4; t++)
> {
> - for (unsigned int u = 0; u < 3; u++)
> + for (unsigned int u = 0; u < 4; u++)
> {
> - for (unsigned int tu = 0; tu < 3; tu++)
> + for (unsigned int tu = 0; tu < 4; tu++)
> {
> dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
> }// end loop over 'tu'
> @@ -618,11 +808,11 @@
>
> if (combinations[r][s] == 1)
> {
> - for (unsigned int t = 0; t < 3; t++)
> + for (unsigned int t = 0; t < 4; t++)
> {
> - for (unsigned int u = 0; u < 3; u++)
> + for (unsigned int u = 0; u < 4; u++)
> {
> - for (unsigned int tu = 0; tu < 3; tu++)
> + for (unsigned int tu = 0; tu < 4; tu++)
> {
> dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
> }// end loop over 'tu'
> @@ -630,10 +820,24 @@
> }// end loop over 't'
> }
>
> + if (combinations[r][s] == 2)
> + {
> + for (unsigned int t = 0; t < 4; t++)
> + {
> + for (unsigned int u = 0; u < 4; u++)
> + {
> + for (unsigned int tu = 0; tu < 4; tu++)
> + {
> + dmats[t][u] += dmats2[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 s = 0; s < 4; s++)
> {
> - for (unsigned int t = 0; t < 3; t++)
> + for (unsigned int t = 0; t < 4; t++)
> {
> derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
> }// end loop over 't'
> @@ -669,129 +873,359 @@
> {
>
> // 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] = \
> - {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++)
> + double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000};
> +
> + // Declare helper variables.
> + unsigned int rr = 0;
> + unsigned int ss = 0;
> + double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
> +
> + // Compute basisvalues.
> + basisvalues[0] = 1.000000000000000;
> + basisvalues[1] = tmp0;
> + for (unsigned int r = 0; r < 1; r++)
> + {
> + rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
> + ss = r*(r + 1)*(r + 2)/6;
> + basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
> + }// end loop over 'r'
> + for (unsigned int r = 0; r < 1; r++)
> + {
> + for (unsigned int s = 0; s < 1 - r; s++)
> + {
> + rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
> + ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
> + basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
> + }// end loop over 's'
> + }// end loop over 'r'
> + for (unsigned int r = 0; r < 2; r++)
> + {
> + for (unsigned int s = 0; s < 2 - r; s++)
> + {
> + for (unsigned int t = 0; t < 2 - r - s; t++)
> + {
> + rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
> + basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
> + }// end loop over 't'
> + }// end loop over 's'
> + }// end loop over 'r'
> +
> + // Table(s) of coefficients.
> + static const double coefficients0[4] = \
> + {0.288675134594813, 0.000000000000000, 0.210818510677892, -0.074535599249993};
> +
> + // Tables of derivatives of the polynomial base (transpose).
> + static const double dmats0[4][4] = \
> + {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
> +
> + static const double dmats1[4][4] = \
> + {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
> +
> + static const double dmats2[4][4] = \
> + {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {5.163977794943223, 0.000000000000000, 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[4][4] = \
> + {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}};
> +
> + // Declare (auxiliary) derivative matrix (of polynomial basis).
> + double dmats_old[4][4] = \
> + {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.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 < 4; t++)
> + {
> + for (unsigned int u = 0; u < 4; 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 < 4; t++)
> + {
> + for (unsigned int u = 0; u < 4; 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 < 4; t++)
> + {
> + for (unsigned int u = 0; u < 4; u++)
> + {
> + for (unsigned int tu = 0; tu < 4; 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 < 4; t++)
> + {
> + for (unsigned int u = 0; u < 4; u++)
> + {
> + for (unsigned int tu = 0; tu < 4; tu++)
> + {
> + dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
> + }// end loop over 'tu'
> + }// end loop over 'u'
> + }// end loop over 't'
> + }
> +
> + if (combinations[r][s] == 2)
> + {
> + for (unsigned int t = 0; t < 4; t++)
> + {
> + for (unsigned int u = 0; u < 4; u++)
> + {
> + for (unsigned int tu = 0; tu < 4; tu++)
> + {
> + dmats[t][u] += dmats2[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 < 4; s++)
> + {
> + for (unsigned int t = 0; t < 4; t++)
> + {
> + derivatives[r] += coefficients0[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;
> + break;
> + }
> + case 3:
> + {
> +
> + // Array of basisvalues.
> + double basisvalues[4] = {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000};
> +
> + // Declare helper variables.
> + unsigned int rr = 0;
> + unsigned int ss = 0;
> + double tmp0 = 0.500000000000000*(2.000000000000000 + Y + Z + 2.000000000000000*X);
> +
> + // Compute basisvalues.
> + basisvalues[0] = 1.000000000000000;
> + basisvalues[1] = tmp0;
> + for (unsigned int r = 0; r < 1; r++)
> + {
> + rr = (r + 1)*(r + 1 + 1)*(r + 1 + 2)/6 + 1*(1 + 1)/2;
> + ss = r*(r + 1)*(r + 2)/6;
> + basisvalues[rr] = basisvalues[ss]*(r*(1.000000000000000 + Y) + (2.000000000000000 + Z + 3.000000000000000*Y)/2.000000000000000);
> + }// end loop over 'r'
> + for (unsigned int r = 0; r < 1; r++)
> + {
> + for (unsigned int s = 0; s < 1 - r; s++)
> + {
> + rr = (r + s + 1)*(r + s + 1 + 1)*(r + s + 1 + 2)/6 + (s + 1)*(s + 1 + 1)/2 + 1;
> + ss = (r + s)*(r + s + 1)*(r + s + 2)/6 + s*(s + 1)/2;
> + basisvalues[rr] = basisvalues[ss]*(1.000000000000000 + r + s + Z*(2.000000000000000 + r + s));
> + }// end loop over 's'
> + }// end loop over 'r'
> + for (unsigned int r = 0; r < 2; r++)
> + {
> + for (unsigned int s = 0; s < 2 - r; s++)
> + {
> + for (unsigned int t = 0; t < 2 - r - s; t++)
> + {
> + rr = (r + s + t)*(r + s + t + 1)*(r + s + t + 2)/6 + (s + t)*(s + t + 1)/2 + t;
> + basisvalues[rr] *= std::sqrt((0.500000000000000 + r)*(1.000000000000000 + r + s)*(1.500000000000000 + r + s + t));
> + }// end loop over 't'
> + }// end loop over 's'
> + }// end loop over 'r'
> +
> + // Table(s) of coefficients.
> + static const double coefficients0[4] = \
> + {0.288675134594813, 0.000000000000000, 0.000000000000000, 0.223606797749979};
> +
> + // Tables of derivatives of the polynomial base (transpose).
> + static const double dmats0[4][4] = \
> + {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {6.324555320336760, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
> +
> + static const double dmats1[4][4] = \
> + {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {5.477225575051662, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000}};
> +
> + static const double dmats2[4][4] = \
> + {{0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {3.162277660168380, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {1.825741858350554, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {5.163977794943223, 0.000000000000000, 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[4][4] = \
> + {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 0.000000000000000, 1.000000000000000}};
> +
> + // Declare (auxiliary) derivative matrix (of polynomial basis).
> + double dmats_old[4][4] = \
> + {{1.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 1.000000000000000, 0.000000000000000, 0.000000000000000},
> + {0.000000000000000, 0.000000000000000, 1.000000000000000, 0.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 < 4; t++)
> + {
> + for (unsigned int u = 0; u < 4; 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 < 4; t++)
> + {
> + for (unsigned int u = 0; u < 4; 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 < 4; t++)
> + {
> + for (unsigned int u = 0; u < 4; u++)
> + {
> + for (unsigned int tu = 0; tu < 4; 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 < 4; t++)
> + {
> + for (unsigned int u = 0; u < 4; u++)
> + {
> + for (unsigned int tu = 0; tu < 4; tu++)
> + {
> + dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
> + }// end loop over 'tu'
> + }// end loop over 'u'
> + }// end loop over 't'
> + }
> +
> + if (combinations[r][s] == 2)
> + {
> + for (unsigned int t = 0; t < 4; t++)
> + {
> + for (unsigned int u = 0; u < 4; u++)
> + {
> + for (unsigned int tu = 0; tu < 4; tu++)
> + {
> + dmats[t][u] += dmats2[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 < 4; s++)
> + {
> + for (unsigned int t = 0; t < 4; t++)
> {
> derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
> }// end loop over 't'
> @@ -837,7 +1271,7 @@
> unsigned int num_derivatives = 1;
> for (unsigned int r = 0; r < n; r++)
> {
> - num_derivatives *= 2;
> + num_derivatives *= 3;
> }// end loop over 'r'
>
> // Helper variable to hold values of a single dof.
> @@ -848,7 +1282,7 @@
> }// end loop over 'r'
>
> // Loop dofs and call evaluate_basis_derivatives.
> - for (unsigned int r = 0; r < 3; r++)
> + for (unsigned int r = 0; r < 4; r++)
> {
> evaluate_basis_derivatives(r, n, dof_values, coordinates, c);
> for (unsigned int s = 0; s < num_derivatives; s++)
> @@ -870,7 +1304,7 @@
> double vals[1];
>
> // Declare variable for physical coordinates.
> - double y[2];
> + double y[3];
> const double * const * x = c.coordinates;
> switch (i)
> {
> @@ -878,6 +1312,7 @@
> {
> y[0] = x[0][0];
> y[1] = x[0][1];
> + y[2] = x[0][2];
> f.evaluate(vals, y, c);
> return vals[0];
> break;
> @@ -886,6 +1321,7 @@
> {
> y[0] = x[1][0];
> y[1] = x[1][1];
> + y[2] = x[1][2];
> f.evaluate(vals, y, c);
> return vals[0];
> break;
> @@ -894,6 +1330,16 @@
> {
> y[0] = x[2][0];
> y[1] = x[2][1];
> + y[2] = x[2][2];
> + f.evaluate(vals, y, c);
> + return vals[0];
> + break;
> + }
> + case 3:
> + {
> + y[0] = x[3][0];
> + y[1] = x[3][1];
> + y[2] = x[3][2];
> f.evaluate(vals, y, c);
> return vals[0];
> break;
> @@ -912,20 +1358,28 @@
> double vals[1];
>
> // Declare variable for physical coordinates.
> - double y[2];
> + double y[3];
> const double * const * x = c.coordinates;
> y[0] = x[0][0];
> y[1] = x[0][1];
> + y[2] = x[0][2];
> f.evaluate(vals, y, c);
> values[0] = vals[0];
> y[0] = x[1][0];
> y[1] = x[1][1];
> + y[2] = x[1][2];
> f.evaluate(vals, y, c);
> values[1] = vals[0];
> y[0] = x[2][0];
> y[1] = x[2][1];
> + y[2] = x[2][2];
> f.evaluate(vals, y, c);
> values[2] = vals[0];
> + y[0] = x[3][0];
> + y[1] = x[3][1];
> + y[2] = x[3][2];
> + f.evaluate(vals, y, c);
> + values[3] = vals[0];
> }
>
> /// Interpolate vertex values from dof values
> @@ -937,6 +1391,7 @@
> vertex_values[0] = dof_values[0];
> vertex_values[1] = dof_values[1];
> vertex_values[2] = dof_values[2];
> + vertex_values[3] = dof_values[3];
> }
>
> /// Return the number of sub elements (for a mixed element)
> @@ -978,7 +1433,7 @@
> /// 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 "FFC dofmap for FiniteElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1)";
> }
>
> /// Return true iff mesh entities of topological dimension d are needed
> @@ -1001,6 +1456,11 @@
> return false;
> break;
> }
> + case 3:
> + {
> + return false;
> + break;
> + }
> }
>
> return false;
> @@ -1035,25 +1495,25 @@
> /// 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 4;
> }
>
> /// Return the maximum dimension of the local finite element function space
> virtual unsigned int max_local_dimension() const
> {
> - return 3;
> + return 4;
> }
>
> // Return the geometric dimension of the coordinates this dof map provides
> virtual unsigned int geometric_dimension() const
> {
> - return 2;
> + return 3;
> }
>
> /// Return the number of dofs on each cell facet
> virtual unsigned int num_facet_dofs() const
> {
> - return 2;
> + return 3;
> }
>
> /// Return the number of dofs associated with each cell entity of dimension d
> @@ -1076,6 +1536,11 @@
> return 0;
> break;
> }
> + case 3:
> + {
> + return 0;
> + break;
> + }
> }
>
> return 0;
> @@ -1089,6 +1554,7 @@
> dofs[0] = c.entity_indices[0][0];
> dofs[1] = c.entity_indices[0][1];
> dofs[2] = c.entity_indices[0][2];
> + dofs[3] = c.entity_indices[0][3];
> }
>
> /// Tabulate the local-to-local mapping from facet dofs to cell dofs
> @@ -1101,18 +1567,28 @@
> {
> dofs[0] = 1;
> dofs[1] = 2;
> + dofs[2] = 3;
> break;
> }
> case 1:
> {
> dofs[0] = 0;
> dofs[1] = 2;
> + dofs[2] = 3;
> break;
> }
> case 2:
> {
> dofs[0] = 0;
> dofs[1] = 1;
> + dofs[2] = 3;
> + break;
> + }
> + case 3:
> + {
> + dofs[0] = 0;
> + dofs[1] = 1;
> + dofs[2] = 2;
> break;
> }
> }
> @@ -1123,18 +1599,18 @@
> virtual void tabulate_entity_dofs(unsigned int* dofs,
> unsigned int d, unsigned int i) const
> {
> - if (d > 2)
> + if (d > 3)
> {
> - throw std::runtime_error("d is larger than dimension (2)");
> + throw std::runtime_error("d is larger than dimension (3)");
> }
>
> switch (d)
> {
> case 0:
> {
> - if (i > 2)
> + if (i > 3)
> {
> - throw std::runtime_error("i is larger than number of entities (2)");
> + throw std::runtime_error("i is larger than number of entities (3)");
> }
>
> switch (i)
> @@ -1154,6 +1630,11 @@
> dofs[0] = 2;
> break;
> }
> + case 3:
> + {
> + dofs[0] = 3;
> + break;
> + }
> }
>
> break;
> @@ -1168,6 +1649,11 @@
>
> break;
> }
> + case 3:
> + {
> +
> + break;
> + }
> }
>
> }
> @@ -1180,10 +1666,16 @@
>
> coordinates[0][0] = x[0][0];
> coordinates[0][1] = x[0][1];
> + coordinates[0][2] = x[0][2];
> coordinates[1][0] = x[1][0];
> coordinates[1][1] = x[1][1];
> + coordinates[1][2] = x[1][2];
> coordinates[2][0] = x[2][0];
> coordinates[2][1] = x[2][1];
> + coordinates[2][2] = x[2][2];
> + coordinates[3][0] = x[3][0];
> + coordinates[3][1] = x[3][1];
> + coordinates[3][2] = x[3][2];
> }
>
> /// Return the number of sub dof maps (for a mixed element)
> @@ -1225,10 +1717,10 @@
> const double * const * w,
> const ufc::cell& c) const
> {
> - // Number of operations (multiply-add pairs) for Jacobian data: 11
> - // Number of operations (multiply-add pairs) for geometry tensor: 8
> - // Number of operations (multiply-add pairs) for tensor contraction: 11
> - // Total number of operations (multiply-add pairs): 30
> + // Number of operations (multiply-add pairs) for Jacobian data: 32
> + // Number of operations (multiply-add pairs) for geometry tensor: 27
> + // Number of operations (multiply-add pairs) for tensor contraction: 28
> + // Total number of operations (multiply-add pairs): 87
>
> // Extract vertex coordinates
> const double * const * x = c.coordinates;
> @@ -1236,37 +1728,70 @@
> // 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_02 = x[3][0] - x[0][0];
> const double J_10 = x[1][1] - x[0][1];
> const double J_11 = x[2][1] - x[0][1];
> + const double J_12 = x[3][1] - x[0][1];
> + const double J_20 = x[1][2] - x[0][2];
> + const double J_21 = x[2][2] - x[0][2];
> + const double J_22 = x[3][2] - x[0][2];
> +
> + // Compute sub determinants
> + const double d_00 = J_11*J_22 - J_12*J_21;
> + const double d_01 = J_12*J_20 - J_10*J_22;
> + const double d_02 = J_10*J_21 - J_11*J_20;
> + const double d_10 = J_02*J_21 - J_01*J_22;
> + const double d_11 = J_00*J_22 - J_02*J_20;
> + const double d_12 = J_01*J_20 - J_00*J_21;
> + const double d_20 = J_01*J_12 - J_02*J_11;
> + const double d_21 = J_02*J_10 - J_00*J_12;
> + const double d_22 = J_00*J_11 - J_01*J_10;
>
> // Compute determinant of Jacobian
> - double detJ = J_00*J_11 - J_01*J_10;
> + double detJ = J_00*d_00 + J_10*d_10 + J_20*d_20;
>
> // 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;
> + const double K_00 = d_00 / detJ;
> + const double K_01 = d_10 / detJ;
> + const double K_02 = d_20 / detJ;
> + const double K_10 = d_01 / detJ;
> + const double K_11 = d_11 / detJ;
> + const double K_12 = d_21 / detJ;
> + const double K_20 = d_02 / detJ;
> + const double K_21 = d_12 / detJ;
> + const double K_22 = d_22 / detJ;
>
> // Set scale factor
> const double det = std::abs(detJ);
>
> // Compute geometry tensor
> - const double G0_0_0 = det*(K_00*K_00 + K_01*K_01);
> - const double G0_0_1 = det*(K_00*K_10 + K_01*K_11);
> - const double G0_1_0 = det*(K_10*K_00 + K_11*K_01);
> - const double G0_1_1 = det*(K_10*K_10 + K_11*K_11);
> + const double G0_0_0 = det*(K_00*K_00 + K_01*K_01 + K_02*K_02);
> + const double G0_0_1 = det*(K_00*K_10 + K_01*K_11 + K_02*K_12);
> + const double G0_0_2 = det*(K_00*K_20 + K_01*K_21 + K_02*K_22);
> + const double G0_1_0 = det*(K_10*K_00 + K_11*K_01 + K_12*K_02);
> + const double G0_1_1 = det*(K_10*K_10 + K_11*K_11 + K_12*K_12);
> + const double G0_1_2 = det*(K_10*K_20 + K_11*K_21 + K_12*K_22);
> + const double G0_2_0 = det*(K_20*K_00 + K_21*K_01 + K_22*K_02);
> + const double G0_2_1 = det*(K_20*K_10 + K_21*K_11 + K_22*K_12);
> + const double G0_2_2 = det*(K_20*K_20 + K_21*K_21 + K_22*K_22);
>
> // Compute element tensor
> - A[0] = 0.500000000000000*G0_0_0 + 0.500000000000000*G0_0_1 + 0.500000000000000*G0_1_0 + 0.500000000000000*G0_1_1;
> - A[1] = -0.500000000000000*G0_0_0 - 0.500000000000000*G0_1_0;
> - A[2] = -0.500000000000000*G0_0_1 - 0.500000000000000*G0_1_1;
> - A[3] = -0.500000000000000*G0_0_0 - 0.500000000000000*G0_0_1;
> - A[4] = 0.500000000000000*G0_0_0;
> - A[5] = 0.500000000000000*G0_0_1;
> - A[6] = -0.500000000000000*G0_1_0 - 0.500000000000000*G0_1_1;
> - A[7] = 0.500000000000000*G0_1_0;
> - A[8] = 0.500000000000000*G0_1_1;
> + A[0] = 0.166666666666666*G0_0_0 + 0.166666666666666*G0_0_1 + 0.166666666666666*G0_0_2 + 0.166666666666666*G0_1_0 + 0.166666666666666*G0_1_1 + 0.166666666666666*G0_1_2 + 0.166666666666666*G0_2_0 + 0.166666666666666*G0_2_1 + 0.166666666666666*G0_2_2;
> + A[1] = -0.166666666666666*G0_0_0 - 0.166666666666666*G0_1_0 - 0.166666666666666*G0_2_0;
> + A[2] = -0.166666666666666*G0_0_1 - 0.166666666666666*G0_1_1 - 0.166666666666666*G0_2_1;
> + A[3] = -0.166666666666666*G0_0_2 - 0.166666666666666*G0_1_2 - 0.166666666666666*G0_2_2;
> + A[4] = -0.166666666666666*G0_0_0 - 0.166666666666666*G0_0_1 - 0.166666666666666*G0_0_2;
> + A[5] = 0.166666666666666*G0_0_0;
> + A[6] = 0.166666666666666*G0_0_1;
> + A[7] = 0.166666666666666*G0_0_2;
> + A[8] = -0.166666666666666*G0_1_0 - 0.166666666666666*G0_1_1 - 0.166666666666666*G0_1_2;
> + A[9] = 0.166666666666666*G0_1_0;
> + A[10] = 0.166666666666667*G0_1_1;
> + A[11] = 0.166666666666666*G0_1_2;
> + A[12] = -0.166666666666666*G0_2_0 - 0.166666666666666*G0_2_1 - 0.166666666666666*G0_2_2;
> + A[13] = 0.166666666666666*G0_2_0;
> + A[14] = 0.166666666666666*G0_2_1;
> + A[15] = 0.166666666666666*G0_2_2;
> }
>
> };
> @@ -1305,7 +1830,7 @@
> /// Return a string identifying the form
> virtual const char* signature() const
> {
> - return "Form([Integral(IndexSum(Product(Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(0),), {Index(0): 2})), MultiIndex((Index(1),), {Index(1): 2})), Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 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 "Form([Integral(IndexSum(Product(Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1), 0), MultiIndex((Index(0),), {Index(0): 3})), MultiIndex((Index(0),), {Index(0): 3})), MultiIndex((Index(1),), {Index(1): 3})), Indexed(ComponentTensor(SpatialDerivative(Argument(FiniteElement('Lagrange', Cell('tetrahedron', 1, Space(3)), 1), 1), MultiIndex((Index(2),), {Index(2): 3})), MultiIndex((Index(2),), {Index(2): 3})), MultiIndex((Index(1),), {Index(1): 3}))), MultiIndex((Index(1),), {Index(1): 3})), Measure('cell', 0, None))])";
> }
>
> /// Return the rank of the global tensor (r)
>
> === modified file 'source/demos/la/eigenvalue/cpp/StiffnessMatrix.ufl'
> --- source/demos/la/eigenvalue/cpp/StiffnessMatrix.ufl 2010-09-01 15:13:40 +0000
> +++ source/demos/la/eigenvalue/cpp/StiffnessMatrix.ufl 2010-09-06 15:40:49 +0000
> @@ -2,11 +2,11 @@
> # Licensed under the GNU LGPL Version 2.1
> #
> # First added: 2005-06-05
> -# Last changed: 2006-03-28
> +# Last changed: 2010-09-05
> #
> # The bilinear form for a stiffness matrix (Poisson).
>
> -element = FiniteElement("Lagrange", "triangle", 1)
> +element = FiniteElement("Lagrange", "tetrahedron", 1)
>
> v = TestFunction(element)
> u = TrialFunction(element)
>
> === modified file 'source/demos/la/eigenvalue/cpp/main.cpp'
> --- source/demos/la/eigenvalue/cpp/main.cpp 2010-09-01 15:17:46 +0000
> +++ source/demos/la/eigenvalue/cpp/main.cpp 2010-09-06 15:40:49 +0000
> @@ -2,9 +2,10 @@
> // Licensed under the GNU LGPL Version 2.1.
> //
> // Modified by Anders Logg, 2008.
> +// Modified by Marie E. Rognes, 2010.
> //
> // First added: 2007-03-08
> -// Last changed: 2010-01-02
> +// Last changed: 2010-09-05
> //
> // This simple program illustrates the use of the SLEPc eigenvalue solver.
>
> @@ -21,34 +22,36 @@
> parameters["linear_algebra_backend"] = "PETSc";
>
> // Create mesh
> - UnitSquare mesh(64, 64);
> + Mesh mesh("box_with_dent.xml.gz");
>
> - // Build stiftness matrix
> + // Build stiffness matrix
> Matrix A;
> StiffnessMatrix::FunctionSpace V(mesh);
> StiffnessMatrix::BilinearForm a(V, V);
> assemble(A, a);
>
> - cout << A << endl;
> -
> // Get PETSc matrix
> PETScMatrix& AA(A.down_cast<PETScMatrix>());
>
> - // Compute the first n eigenvalues
> - unsigned int n = 10;
> + // Create eigensolver
> SLEPcEigenSolver esolver;
> - esolver.parameters["spectrum"] = "largest magnitude";
> - esolver.solve(AA, n);
> -
> - cout << "Solver converted in " << esolver.get_iteration_number() << " iterations" << endl;
> -
> - // Display eigenvalues
> - for (unsigned int i = 0; i < n; i++)
> - {
> - double lr, lc;
> - esolver.get_eigenvalue(lr, lc, i);
> - cout << "Eigenvalue " << i << ": " << lr << endl;
> - }
> +
> + // Compute all eigenvalues of A x = \lambda x
> + esolver.solve(AA);
> +
> + // Extract largest (first) eigenpair
> + double r, c;
> + PETScVector rx(A.size(1));
> + PETScVector cx(A.size(1));
> + esolver.get_eigenpair(r, c, rx, cx, 0);
> +
> + std::cout << "Largest eigenvalue: " << r << std::endl;
> +
> + // Initialize function with eigenvector
> + Function u(V, rx);
> +
> + // Plot eigenfunction
> + plot(u);
>
> #else
>
>
> === modified file 'source/demos/pde/biharmonic/python/demo.py'
> --- source/demos/pde/biharmonic/python/demo.py 2010-09-01 20:32:22 +0000
> +++ source/demos/pde/biharmonic/python/demo.py 2010-09-06 15:40:49 +0000
> @@ -52,7 +52,7 @@
> h = CellSize(mesh)
> h_avg = (h('+') + h('-'))/2.0
> n = FacetNormal(mesh)
> -f = Source(V)
> +f = Source()
>
> # Penalty parameter
> alpha = Constant(8.0)
>
> === modified file 'source/demos/pde/biharmonic/python/documentation.rst'
> --- source/demos/pde/biharmonic/python/documentation.rst 2010-09-01 22:59:53 +0000
> +++ source/demos/pde/biharmonic/python/documentation.rst 2010-09-06 15:40:49 +0000
> @@ -85,7 +85,7 @@
> h = CellSize(mesh)
> h_avg = (h('+') + h('-'))/2.0
> n = FacetNormal(mesh)
> - f = Source(V)
> + f = Source()
>
> # Penalty parameter
> alpha = Constant(8.0)
>
> === modified file 'source/demos/pde/cahn-hilliard/python/demo.py'
> --- source/demos/pde/cahn-hilliard/python/demo.py 2010-09-01 20:32:22 +0000
> +++ source/demos/pde/cahn-hilliard/python/demo.py 2010-09-06 15:40:49 +0000
> @@ -18,8 +18,8 @@
> def eval(self, values, x):
> values[0] = 0.63 + 0.02*(0.5 - random.random())
> values[1] = 0.0
> - def dim(self):
> - return 2
> + def value_shape(self):
> + return (2,)
>
> # Class for interfacing with the Newton solver
> class CahnHilliardEquation(NonlinearProblem):
>
> === modified file 'source/demos/pde/cahn-hilliard/python/documentation.rst'
> --- source/demos/pde/cahn-hilliard/python/documentation.rst 2010-09-01 20:10:55 +0000
> +++ source/demos/pde/cahn-hilliard/python/documentation.rst 2010-09-06 15:40:49 +0000
> @@ -34,8 +34,8 @@
> def eval(self, values, x):
> values[0] = 0.63 + 0.02*(0.5 - random.random())
> values[1] = 0.0
> - def dim(self):
> - return 2
> + def value_shape(self):
> + return (2,)
>
> It is a subclass of ``Expression``. In the constructor (``__init__``),
> the random number generator is seeded. If the program is run in
> @@ -43,8 +43,8 @@
> number to ensure a different sequence of numbers on each process. The
> function ``eval`` returns values for a function of dimension two. For
> the first component of the function, a randomized value is returned.
> -The function ``dim`` declares that the ``Expression`` is of dimension
> -two.
> +The method ``value_shape`` declares that the ``Expression`` is vector
> +valued with dimension two.
>
> .. index:: NonlinearProblem
>
>
> === modified file 'source/demos/pde/mixed-poisson/python/demo.py'
> --- source/demos/pde/mixed-poisson/python/demo.py 2010-09-01 15:17:46 +0000
> +++ source/demos/pde/mixed-poisson/python/demo.py 2010-09-06 15:40:49 +0000
> @@ -53,10 +53,9 @@
> g = sin(5*data.x()[0])
> values[0] = g*data.normal()[0]
> values[1] = g*data.normal()[1]
> - def rank(self):
> - return 1
> - def dim(self):
> - return 2
> + def value_shape(self):
> + return (2,)
> +
> G = BoundarySource()
>
> # Define essential boundary
>
> === modified file 'source/demos/pde/mixed-poisson/python/documentation.rst'
> --- source/demos/pde/mixed-poisson/python/documentation.rst 2010-08-31 21:08:44 +0000
> +++ source/demos/pde/mixed-poisson/python/documentation.rst 2010-09-06 15:40:49 +0000
> @@ -104,7 +104,7 @@
> ``Expression`` class. Overloading the ``eval_data`` method (instead of
> the usual ``eval``) allows us to extract more geometry information
> such as the facet normals. Since this is a vector-valued expression,
> -the methods ``rank`` and ``dim`` must also be specified.
> +we need to overload the ``value_shape`` method.
>
> .. index:: Expression
>
> @@ -116,10 +116,9 @@
> g = sin(5*data.x()[0])
> values[0] = g*data.normal()[0]
> values[1] = g*data.normal()[1]
> - def rank(self):
> - return 1
> - def dim(self):
> - return 2
> + def value_shape(self):
> + return (2,)
> +
> G = BoundarySource()
>
> Specifying the relevant part of the boundary can be done as for the
>
> === modified file 'source/demos/undocumented/dielectric/python/demo.py'
> --- source/demos/undocumented/dielectric/python/demo.py 2010-09-01 15:13:40 +0000
> +++ source/demos/undocumented/dielectric/python/demo.py 2010-09-06 15:40:49 +0000
> @@ -12,7 +12,7 @@
> # u(x,y) = V for y = 0
>
> __author__ = "Kristen Kaasbjerg (cosby@xxxxxxxxx)"
> -__date__ = "2008-02-14 -- 2009-10-07"
> +__date__ = "2008-02-14 -- 2010-09-05"
> __copyright__ = ""
> __license__ = "GNU LGPL Version 2.1"
>
> @@ -104,12 +104,12 @@
> v = TestFunction(V2)
> u = TrialFunction(V2)
> f = Constant(0.0)
> -b = Coefficient(V0)
> +b = Coefficient()
> a = b*dot(grad(v), grad(u))*dx
> L = v*f*dx
>
> # Define boundary condition
> -u0 = DirichletFunction(V2)
> +u0 = DirichletFunction()
> bc = DirichletBC(V2, u0, DirichletBoundary())
>
> # Solve PDE and plot solution
> @@ -123,7 +123,7 @@
> # because it is an interpolation of the exact solution
> # in the finite element space!
> Pk = FunctionSpace(mesh, "CG", 5)
> -exact = Exact(Pk)
> +exact = Exact()
>
> e = u - exact
> L2_norm = e*e*dx
>
> === modified file 'source/demos/undocumented/elasticity/python/demo.py'
> --- source/demos/undocumented/elasticity/python/demo.py 2010-09-01 22:59:53 +0000
> +++ source/demos/undocumented/elasticity/python/demo.py 2010-09-06 15:40:49 +0000
> @@ -3,7 +3,7 @@
> for a gear clamped at two of its ends and twisted 30 degrees."""
>
> __author__ = "Kristian B. Oelgaard (k.b.oelgaard@xxxxxxxxxx)"
> -__date__ = "2007-11-14 -- 2009-10-07"
> +__date__ = "2007-11-14 -- 2010-09-05"
> __copyright__ = "Copyright (C) 2007 Kristian B. Oelgaard"
> __license__ = "GNU LGPL Version 2.1"
>
> @@ -45,8 +45,8 @@
> values[1] = y - x[1]
> values[2] = z - x[2]
>
> - def dim(self):
> - return 3
> + def value_shape(self):
> + return (3,)
>
> # Sub domain for rotation at right end
> def right(x, on_boundary):
>
> === modified file 'source/demos/undocumented/elastodynamics/python/demo.py'
> --- source/demos/undocumented/elastodynamics/python/demo.py 2010-09-01 22:59:53 +0000
> +++ source/demos/undocumented/elastodynamics/python/demo.py 2010-09-06 15:40:49 +0000
> @@ -41,23 +41,22 @@
> self.t = t
> self.dt = dt
> self.old = old
> +
> def eval(self, values, x):
>
> # 'Shift' time for n-1 values
> - t_tmp = float(self.t)
> + t_tmp = self.t
> if self.old and t > 0.0:
> t_tmp -= self.dt
>
> cutoff_t = 10.0*1.0/32.0;
> - if t_tmp < cutoff_t:
> - values[0] = 1.0*t_tmp/cutoff_t
> - values[1] = 0.0
> - else:
> - values[0] = 1.0;
> - values[1] = 0.0;
> -
> - def dim(self):
> - return 2
> + weight = t_tmp/cutoff_t if t_tmp < cutoff_t else 1.0
> +
> + values[0] = 1.0*weight
> + values[1] = 0.0
> +
> + def value_shape(self):
> + return (2,)
>
> # Sub domain for clamp at left end
> def left(x, on_boundary):
> @@ -93,7 +92,7 @@
> beta = 0.36
> gamma = 0.7
> dt = 1.0/32.0
> -t = Time(0.0)
> +t = 0.0
> T = 10*dt
>
> # Some useful factors
> @@ -144,11 +143,13 @@
> problem = VariationalProblem(a, L, bcs=bc, exterior_facet_domains=boundary_subdomains)
>
> vtk_file = File("elasticity.pvd")
> -while t < T:
> +while t <= T:
>
> t += dt
> print "Time: ", t
>
> + p.t = t
> + p0.t = t
> # Solve and update functions
> u = problem.solve()
> update(u, u0, v0, a0, beta, gamma, dt)
>
> === modified file 'source/demos/undocumented/periodic/python/demo.py'
> --- source/demos/undocumented/periodic/python/demo.py 2010-09-01 15:13:40 +0000
> +++ source/demos/undocumented/periodic/python/demo.py 2010-09-06 15:40:49 +0000
> @@ -8,7 +8,7 @@
> # Original implementation: ../cpp/main.cpp by Anders Logg
> #
> __author__ = "Kristian B. Oelgaard (k.b.oelgaard@xxxxxxxxxx)"
> -__date__ = "2007-11-15 -- 2009-10-07"
> +__date__ = "2007-11-15 -- 2010-09-05"
> __copyright__ = "Copyright (C) 2007 Kristian B. Oelgaard"
> __license__ = "GNU LGPL Version 2.1"
>
> @@ -55,7 +55,7 @@
> # Define variational problem
> v = TestFunction(V)
> u = TrialFunction(V)
> -f = Source(V)
> +f = Source()
>
> a = dot(grad(v), grad(u))*dx
> L = v*f*dx
>
> === modified file 'source/demos/undocumented/sym-dirichlet-bc/python/demo.py'
> --- source/demos/undocumented/sym-dirichlet-bc/python/demo.py 2010-09-01 15:13:40 +0000
> +++ source/demos/undocumented/sym-dirichlet-bc/python/demo.py 2010-09-06 15:40:49 +0000
> @@ -4,7 +4,7 @@
> # Modified by Kristian Oelgaard, 2008
>
> __author__ = "Kent-Andre Mardal (kent-and@xxxxxxxxx)"
> -__date__ = "2008-08-13 -- 2009-10-07"
> +__date__ = "2008-08-13 -- 2010-09-05"
> __copyright__ = "Copyright (C) 2008 Kent-Andre Mardal"
> __license__ = "GNU LGPL Version 2.1"
>
> @@ -37,8 +37,8 @@
> # Define variational problem
> v = TestFunction(V)
> u = TrialFunction(V)
> -f = Source(V)
> -g = Flux(V)
> +f = Source()
> +g = Flux()
>
> a = inner(grad(v), grad(u))*dx
> L = v*f*dx + v*g*ds
>
> === modified file 'source/programmers-reference/cpp/mesh/index.rst'
> --- source/programmers-reference/cpp/mesh/index.rst 2010-08-30 13:52:07 +0000
> +++ source/programmers-reference/cpp/mesh/index.rst 2010-09-06 15:40:49 +0000
> @@ -34,7 +34,6 @@
> MeshPrimitive
> MeshTopology
> Point
> - Point
> PrimitiveIntersector
> PrimitiveTraits
> Rectangle
>
> # Begin bundle
> IyBCYXphYXIgcmV2aXNpb24gYnVuZGxlIHY0CiMKQlpoOTFBWSZTWR7dD5kANcN/gF19GABY////
> f//ffv////pgKG4+oB94WB2jXzw1rnvSyqA33veFA9seQHbmus8Tru2vc9Sj733eCSS7BnsCtS0y
> UoXXd3F5Smlbu50Nr600KyW7jvb1EemvdtJUQmtOtUABVBQVRttbuEkiNVP1NpoyGp6NJtJ6I0ZG
> U8qeakGQaPSAAABieoJQEASE1MEIynlMhp6jTTQ9QAAAAAAAaaAhSTyn6U9QDQAAAAAAAAAAACTS
> SFMSaNAETQ0yaep6gGjR6g09QAGgAyaBoESU0CaTKYyVN4mjSYyaTap5NBGMp6h6INqDQ0ZANAVJ
> EBACAEaBNA0m1GE001TJ5T1NNqeUw1Gg9TNIyYSQgWfDlbKc/QGiFrN5CDabZIHV43EaLt6KLLmc
> oxXSs8jzYaTYfczju/HsfBwXfrTuZaf8cinY6ZufHmHDVOx4QWa5+tk9HE5CYxOHqTsdJpDhM1xt
> rVtWKB2J8X0vL5btZ4fJcy6+IfL184dc2bNfHPKSeDSqVqueNDdMZaFrrB/nKEw6Wb9PC7bGwUiz
> PHSGPK2Mh4xgogqoep5h4KkUXMsc6uzcslimW6llKXMYyFFFlz+eE0wW2sSnRZtGNjGLOy3nO+Sq
> iUor6iH8RQlEYQ3pSKK/2WU6KXRa9pz4YwbWLcvdZWwuxKMxh18uWPBpUGabIo81neUosyYQrrKN
> opvULKWYu1TDAZCtYk1F4Uj5GxcYxklIs2rDMlzMUuxKGAuLC/CqnrtgtXs88e4MqHpkFobjpjSV
> CVKtyd6r8DjtfaS7fty8N9JoApfDOAIJIakKXLqthRXZ1KAl3785vYlUb3aLsej5MqkFnLUuPRQg
> YpezDlbVsSzDJRUQiZl4tDLE3vO8RWk3HeWacMFjAWUlqITe9LUs9VvCqyOqYwMSEYZyXbEVNKKA
> WEAQYgqD/SCIo8UBFRwB3wERkQDxqEItQip+iWRhUpS/2nVPt9eoRePsFfJXzv8kaJvtyyPIxBoG
> LNqEmT/AXeQ7Q69vba1rW0GglDShcULiwobxQsLChYZjMXGAoYDEf3hpDCVKlSp9Yf3ihkNpLDsF
> BQ0ETaMQm0biVIc+aOmiXCBhMJhLLLJEDcDUE5ec73IJuGGRsNQ6huC4OUe1X7RcBZIRSAKSLCCI
> CJFgCkUWQWEFgHUdURJAgkQCDPJ6bbbbaLW21qtfW0TQFYCgEnwwgTp1vg7vzl1y9vSr29rC2eda
> nlxEeodNl+XFSwFShCHsZQBIgdvLpx1dEh3VTO3MMTZNIaRJAkNOiZFFEkyZmYV0ybYrUvD5YCu8
> kxN3e6irIovhRi0lqxVpCzTML2c4GGsj0bKuJw73e+FxFKDVtTfIDmCbRE1jZptvkjTxltuyrTAJ
> E+vYbAEj5B5A7HzkNj44bdOGcZVVoekHRaqOrKZVN8g6BrNTXM5gX3thmN9Iy7HJVoyXXKkLCnUD
> qcxOiJgL0rm/KCNNziu6zRXZ6RNMzLOpHsbp0u9Dkwbtmu5O0gB1wm3HPi9s4fUizLUaFUK6YKLK
> yzFVZSqvdhLMxdZMreREJN1gqUW8sEqrsrWxar3u1wrPjBeIl3GvecRcxFybYhgwoXSlL4uxiRYr
> dQIth55DyJ2ixZFFhP6A/0uwkrDxvbs/yP/D2vcyE0HskjcKsXgVIag3Asj92Fj7+F3P1t7c3NNM
> aesMgXyRAsNhvJbbbbcnQh4OcwDeM3jKbyzCyxggli7y8mryu88Ywgwg+MPDZBlpX4eFu5J5MELY
> 43uQwYwY3tNDBj1MabmNNNzBykLN6WWhZaWWhZaWX1cQ55zygMo86fbtMSrw1opqKUJShUXiSuFO
> 92zKvxKf0TE5qGgaNktV+tnkXGK3Y1NfTCHg7fRenI3iMRiM7xSylLKUspS0sstLOn87z63U5OP6
> 6jE5bIWnWJ65t/zwbFEQ3EEQQQSlhf1dAAnf+Q8ed3i+Xn8VV+ARy2OGkk9uyf6+nDyW8Lxnt9C9
> 1oST1+HRlvxCL5nSsvQ+Ekv3MVmA29mnBzKtKnGMbx1HDV1BHTPkcZJfLJMzNm9APYPwiTl1xJgO
> oYRJrkSm9vDM19TUI0kMJCwLSGg13eWRw3+TIb+Es2ibQX0BcS4KGzHf1cUks6duu3q1bhsNDRTL
> wxt5suBNhKJ20RYi11yLar2zah1dwwKCq48NCQu/37uRTSQiG4aGDBX8RMOAcgbeHHud1v7tne4g
> /v28bXckbSM/yDnwc5rOjTm9SqqqqEl1NryU0ikeKIm28d+y3Z20vmmZt7x/599yB4gQhtMyfsdF
> QUkBZAv3uXy/Cqr7JIsiy1lo+5SKQ79huP+m6NpJKKIz8dixloWH3SVrVsCjQfxF10sCjWLIlChS
> 5V1p4MmUmls4fbsb7bIkwChQKFAoKGFCwoKDkEQ3QCSnzTkeIfZMTv3MCW1ZUhfNaRHisxLxG9SH
> ll+oZpmngnqcO+Te+nJvT9mBUaVkMxk5dQwFxFChgKChWRVNtRbnJdabGSaLhvFBgI+lMbzLCq+b
> wv9b9eNWrtq49tiv6xxkNg5jrFw2MGiyyxZZ7B+Q64ZBoKEmgochQ8kZyHzUG2bJtlpaWlNMYNMY
> NMYxjTTTTTTTrAMyyk4KQ/LmQJAD2ACTrDudqqKdlKqirWXhLCyoVdSFIfTlI8SSdfwOD6T0uvv/
> byOozKFihuLFiZlClBFKQSlGzSMy4nGxfLrqqqwtY9TAXfXXyEiVBfC8z0sALcDAhLTtE18Wubg5
> naNi8TG+YtGuNyjE62HePDObWKwYjnDgji4jIx0ZPiwRijWaRD9Wxa/KrWub1FuMh7nu9m0uZI88
> nU5c/FBgdbxdfFuvGPGd/vw7OCblPNKSlJ+ZS7+WW6uiFqFChQoUKDiLlxQShSEUKEiiilSQUqlK
> KQsOwtCN9E6zpliUMAvYt2zuPCSHbOMg9vofJj6p98j+o7hL6zwWWPi014ezKbp3/lIk94lBH1z9
> 66j+1hUBhdR5yrNriw2tX7Au8wAnvgCYqqJmEdONJOR6TPodhLGAUVrxmtg4jEAheCaWLlttuQsp
> MKU6ySbpICn0XzsnUlTTN02QhiQ0qLMT22QNnmmyE6kwhaNoNRS6Km0D04OLokKNbH1mq80w7nP0
> OCP3QeDc5FJYkYtttt7SBkyZM7ocRjjjlZOUlWf2iKQKFERQpIJ7pBEYosUgsgqwoHMCwO/x6OQi
> o8SN11Ilwio7gHYICMdHXiUwgKnGU4psiVBDkKQR3CBoGO7ddxgWFUC6gh1uckwnB74Utqi1MnAw
> j6MGLbqXK8hANC7KtRwkewZHBUYGUF4MnqeFcQBtDmC27wg0G2libDJ0Er8BWsm9lBzYcBShsLX7
> SpM9lZjiqBWoMCIFZwNkknCEHcm4IiLT5gD/quEVTBPtViSSVFSTscTk/D9/Bi7nVTmspv6ThfHY
> 7DogU8WPFq3ZWlF6im6IaY7LrooZTGWvvq5ONgE9YBAOe29ttu0LKTQGjDB4dXjDl0Y6iw1TM24v
> BLKVcWgzxDYs6OjoMMZRU85jHDZrIuUavcBhsIoBGABGAYGBkUhSINttttuBTDmYbAT5QFkkoFCh
> 0U04KTfFclviAO+tUNk5gvr+czIOxqkki37zdMubjBPAqyPLjTRRoiYFKRpmRjVQUD0owgjI6Cpr
> dk3o5klESULjQ4iF06QYUlGCEcUFB0UooL1akpy9NLdeWJ3HlanEcJyg+o/KdNOqhebFcZDe4r91
> jpIxjC6wvLdMFqJwoGR/CsfRQCw11NnRGB+1vc79um+9uF6dna4QJVKy3oLOGFpz0ts0WOeNp5ql
> leQqimBYhadYwZF5t2WvIdi921XgrYpxmZhP+ih0Xcl0O3LQE9Ha5U4AiLC/F5hjnxbyR1e5xb9z
> yLpaZRHlJPanfBhUSZSTNzZOflxYNris2HA7pZhZi2PO3OdJD3u+Tbhn1DaV30a7sLs7WKqGGViy
> ow6PiFMpimIhGuGBmzOSNdnbO9MzYqohitaaamtOeEsmEBKNG5ufEWC7a58H+tExHxlAJDLLOqqq
> qEKKSxIxbbS3odYbQ882QBO9GcGiEOY1H/euUjOoOTwOfruykqBOPU3LLGo6TW6B9dKwOiHa8rbo
> SIncRJWqdE7vObPRkyJo85XAQBrvW8c2aQLKBBsokyJyOskh6RgZKbAVrojRGMBsrRQ94gx6Iud1
> U4WKorg7gp5GNT2nQTg9AxcvVqo/fIXDyASyS0KyKqC9C3X4uLVqti4xchpA8BqRBWOuyXXR3AXX
> wx4OSzTr92aGfcN6XnkGB7h7kjCiJlJNdex4P4PJZzfLO7orF07zzeRyXZJbWuvtnj2oq6IyDUhB
> UqqvwwzyeNoz5XjA65TfLZ6RNFAURrpskD4ALmkLDWVZoDRyWgBjMp0Ot9y8Rg9HRzfbwdWboX3O
> vm1MvM5LTWlqK08ehbuam4+LSTFrYNxrdeYSsMq2Cr7sWSgjdA6mx0upXFj1YlchBRKRKurdgrLZ
> V5PWEgB26FwrZKSUW4cetwD0SoIEEi4KFXxkPBZtY73NxODVzYMmiD2Tmp8H4PlOtdg47O9YN7t0
> AuJAuIQ46WpC3FjMy7OY2sRJL0WOyz2mHs9r0tjyuHammR62vZ1wYQ8bCAa34cDQGUAUZNqNliy+
> JBsDkMZM4I00pbBAxHaWO4hK1bXr5hrHqSgcC8al+xTQoUE0GyhyAwHgekhN3wu8S/AaMhEEgvlN
> IjwQyi8DHcCykEnXwDFMUyBxxPI5iAcFYO7HVB0p3Mm5nRyzYJ1uZg4G8NIg9GFCTsaOSnFy3sHR
> i5U8JuZPBzdPq+ueydtbXDqnUhGYhBoQoCcii6iqB4msmAaRMZ7QXzCml7tM2hr2awOs1Q9YF48w
> CUSxbIXAp2ud9a1i/DwPdwEkB2uCjMactSanuMOXOEsBDiM54xKkBlMCAZCSYOZCTYnjzg1VSiC3
> DgOIDYYgNCNqWe1yDuHKGKsC95gkwBQGlUdWVMgNJYhHXJADUTo1SwmDRUPLjYo5s0UHJyfAcGCx
> YPMfJPhOwLFTKLy3uF4lEAGS62ksGdZlQIL8YAANJk3s2zXDqdSsKpiqRhq0mJwVve4h4hNr3LhA
> Duq3FS7jriLqXqpuQWemQvrqrUxdzEW3117XBfbVy0QqWTuITPbrjODGRwyHBoECfRqzaAgOaZEL
> RDucmnDduZcxxOK+Cz9LfxrRWzMVys28HBWTcpdg7erZJQdg8Ojo0aq6mcwcihs1wUG83q72bR4h
> U+YaQ6p7iIHpSaz4moPAhYKNpZDpVtpmJKvazjTd6Y2gPh6WZrCrSybEZeFjSCGgbL20yzdSChjs
> 7JPAvrkVrdHWsXSpBoG0KaReC1Je0AeQuJVM4O0vUGJ0AqF9Hk9GeOvDhqa465C+rsb+nW6OHLjr
> wMm/iLJxcE543y0yc2rZTV3GurC2Ba6s2s6+uDw4NHTs4KFkOwUMHQ4O85Ea5SdGZVae556zzasQ
> PjMUtao0w04xa9LUapaa3shMVUsBh9FS+TB+NGz7m3G+QrONysVX22y7nsuUZuajKo4KSczdeYGa
> LkCaFT0iAWStrEwSNxWaqotDBlgPOLysoFzsPFD9+K89tgi5kllzDKYcrsDtY7MM3Aw6miDjuz8f
> GzrbXWzpwe18jv3OrmZckNKAqrroats3iYlbTD2nFHmZsojCAeC5pCCBogYQBih2lJW+SuInyDFq
> kHxoNwLyEguD5FfPIUftKiNC7ZTlssGQuB3FQ5nIx7q0DgqWiPfSpeI4HOAckg4M7gtX0ZMYEwTY
> 1F3e7+bwZOpwU0aN3q6Ox8sz7HTkt0QbQ4j3YL8oSGenwXb+c3HbZe5ycFA3g7dvIW4eOiQ7ceqX
> vrFSdphYMLQVpmtkFhtQ0trbW0Ed+hqCkWcdV0U4HqaNvIYmZErUCCASZUDreaIhYpYdH4wMzMzD
> Mw4YAVAkAhISWMaGmzYJfCBqE3wH2hllAwgX0YXBZ8smSzAV1Duy/Rwx3jXQvRcYbRcYe8bUsi/f
> WyJOacE0GwtYwkMHR5FlDHAFxLCxLEoiiP4+uCSlHqVkgA/WT+UQEKIBAj7iIJfBU/MMhEBJMAZL
> osAnX3/GViqiCsRFT2QtAoDG92SXFQDYh6/uySzzdxvD0KB6VL/b87976wCixrDerBJGT6rvuUJP
> O7jAnuCb/7GF/jyvV8j4cD4WvMa+X15UFiwUFigoKKDvPTJueoYxESx8eKjQiYsQgwYMQiEGAxgL
> Eid6ym03yYAxBiZMyYAxBiZMyYAxDl8ANk9oST1M4Q94QKbH4vNNeQU3uqwz+lZvFIUKRPgb0qNb
> SW/eqVg9CCkjomBSMLSXQgagYyhSMhqhdCxVGSGoGM1ZhouoKAoqp54+sZEYj4HvMC/tqqUyjftj
> T0shw+0wFgTUDKQK8ww5mEnMTYCgIjki/AB4MTEyKKKLewCWGw2GkpKU9jN8P0mehUKnt+P2+2uM
> 0qJYnAPIImACtBaBWUMpeWmQw5D3R9tg19NmwaZ1NirVKqyNspzhZLbbdjM4lzjghucQ4HeDMEwQ
> yDMEwYyzLcwy0ttywHFvMIfRmhnmI5gfUGOrVRRVVW9/yxIQ+cGmMCwRKEoWCWMpQlDWjaAW0pbb
> tNTaWFEpLGLGMsKJSUElhSxlKSxjrRtAKq3Jom0KWUssbKSxJYUspQoSwpYyjKRzTNisk8D2Gh/d
> 5E3YlSofdS1VIp7KQ+3kH7KywVvAiv5yPiqekLap4G43eBTMAMjA6TmOd2rrO1wz7WTzblM7tTPl
> Ifu6zexU3qatx2C3h8q0mMGG9Oe9azeq7iN4OwIoKDEbsZZKIP9j0WTF7ER0TcP9POQdyZejS568
> Pifp9NpEhc/SocgbfErpr0mhNdBTSQi3vcuAzwIGhnScoFHPef0xM3w83IUQZuvfvrynJh8Fy/MT
> Hk6EWxW21GRkTkIyjlSqiAJ+BAXz9sy6A3nAvnYW9hgMpmkIhjYVkhrOQpo6jR7nrYLMphIXfIpw
> avQ9/fE/gKH3/zH+I/kMAw2uDY7F25wc6dHJZxcmTKMV3N2u19WRPEhPmUg0KUnwcJ+v+rtk7urx
> ne9SEjcR17MDdLPQrDue+bHg72ayl3cMMVlRnKwmK7pnLhISjDSz/8VeGSk9a8EtL5S+Yu5zfF5z
> FdzfN4s7HaWp4M3e9fg9v7OfNwZt8w8JAC7pe87C1lKvV9fkVCnyodFPMgMAijEfINQZBgZWw6GG
> HmTTR3O0KzDjHHv8KzDjjYkk3aBHMp6NmN1XdtndNPF7/T8jPxZrbmr3qeS7hwOzL3qXg1er1PW5
> QU4Po9rTYbTtB4FXcHIL4HmTmIkIkKaBYwXLCUD9x5wYnMRzaMJUeKlDYPM7qSbpDtNbMNdA9tVV
> VVXbt8h14Bh4ihKB9XqNYakgeeQeiRPSruDni2iYLA0/1G161MBZgPmvxnKrOP19I6qYQSyWFpBa
> USwqhZZ1vRKT18zCyD0bn4PyfT2Ke5UP7rJpFE+SD1JREoCoB2RR+LI1lyYBATUzIiZlCPhC/tqO
> KSQVdBUD6DDF/ZdCrhgU1mRo9aavXzZOcgzknwc53LOHHfz+ExXWfSeUkKyLWQMIgIa03ZE9257n
> 3Ya5iQuDG3aXyUccWTtaOTx3YWknpknSY1DxotaO5RgxvJO0gkAykxiJxxVmzAeZozlMBAkKDES+
> 0RzWL+Cem2um0BJZIegsLSssIisVwykW2Wq1oYRZDmNSS5wp0ouajgwXTuXOR46Gi9jjNIYWKEWG
> pKgOgp/ceBp6kawncuEfP6T308VeP7FXHuppTXvXkl02ff2TNJjBsObCCPONWDabXymS304ZTE4u
> 3JO5QNM56lpn8hxu+WlKV2rRzudzCkTQMDkTIPmSWCAxiBijAlVBTP/Zo+Hc7m2WwWOtflDfUiXt
> JaXbkroZoYSDuZXQ0MITP02XXkk03xMW5rAE9jOK4KYh5XybSB09dhUe5orO2HOaBxKMdRgwfqfn
> 0b2rJ1+DXB6deT5rMVOe1fFsDi5dHtXmw5QCysAxVn2K86u1X5tj1Z5L3+6odAge9IHvAKCwRkOg
> R9MWIkOsPnwXVqRtHUvfDza3R7X/p2Wkqa+4/BlJ65Cmx95Z1B13HOZrR61Q+QKAVu4XpNmsdRq7
> DfJaymV6uD8X0BUk19m3Ny+zWajCPsppvTM/H9PYz17mbVJp0jQu732TyQR8vXnk7/QtJ3jzjTXd
> 2hvLAknj7cCld53+tVigWtyK4niQ6+rqOXWAYLjXGvYQ2gK86odnAr5oK/qWbGqjORtmk7ZiJ2k+
> KHj8td3p2PNuT9scAivUHT7fJ16stTJAkSSIECQA5XjPlg7u3fMz9KHBYD2x8yVt8K3T0T2iLiKk
> YPk0WqSVbtQepm+x+dg9k2vtWlO4vg+1B/Y7iR+XqQ/JI+pgnKnlaLWfbJDq6HvmybxDtDNvTyaj
> u71Nf4blO+6G3EA7EVaVkgsVhFdWj7khYffpd8R2mJtbk8XqLHDf07e6sAE/j0agwSTTJTlLlFPH
> 5ek3hjmFDwGcaFCGTzEUKuhYH4lQ83bOi9w4U5w54DZGyMYlpJ2Owp1jFSk02443Z2++ybZ2cNMn
> x+hn6G+EvJOG8MZD6DdrI+J2OD8XBuAu9qwZyxMtmTBNSRWH0BCbRkuiqj51nUqCpJhJDqnF9LZ2
> FwixyhgIfYiRjLWzvAVzrSHaYN/OuA5YfjUZXRrDudWGEbqfjI5pg1g0jiJJeSy4JaZCSVIM7GK8
> YXxl5JPybh+5TGJMr64E+VIwJGDGAtheDT26Eqh3xZg9DxLIGFQ9RU04qhsH/H8PnTX4h9BdLyQ0
> ol/DsG7YsIvxvNxx2sTiSXKQMyBdED/3o4DVaCoE+46BH7y4mnhPBCN5ey3V7et0CC9yqkzyLNhD
> EyAcD0Rdt70UGw9mPOBB/ab9sEb3bXPObXvshG5P68IVR3IK0jDabPJR8aklQXfyMw2js6G+Pztb
> P3IPvexh5lHapNrltXdsnl3Rtn2sXVsmDg2pNo+UT7pHCYTfGW2XOA+DljnVgigk+cI4igoUKChS
> Q5ElRJeMI7KWUKo/fJFjPEZiohVYx5T6vsmeJNg+soKie1ghSb/ZUW9tRtkzTb4nJNUSpFtR9O5o
> xQWd73T9n7VVLRVfo9X5Zhof4KNnw8L3EdUiuc+gpLSuiGH/Dm+RNlu1q8PfHWP3DtHQxE6OZzdC
> nvLUPiLuvtWNZ+aQw7NVlpj+YpkPaKEoKEoUJQtaQ+QlSFlhUwztu0pjfGkaGHv94ds2m0TEsveR
> M/g+eyNQgQE9cF95+Y679kAzA2lDxRJoUhUVRVUk1KMqqFQlqgoO4ZBS5C5Au8x0N+W8kA3mJdxD
> r5MQ1pGlds7yre/HIST+yR7lu5Xehk5EHoayTlB7epkY6cXLBctTiPWLnsGU7grtUnFXFqCW5RK9
> 5GsT+v41W9C8rTUDyYcq2vf83TATSgd1QPZno09DR8VAnVXi00EjbQMD3RFo7BWNKzx3CdkKBYYz
> EnqPzHwmAGmQLgcoLVK7BW/mmOOvWiYjCaR1RjFLKks17YtGWqjzzmAnh0PDq0ST27EtDR7CfV57
> eeklmlRFa8/uUcm97HM0eFVVKiqDYZLCPzvKdJI7ZmbUM5211pkehskYYVLFFFCilRE5fQvB0xMx
> O4UTyC6HzL49CWIolEZDQFOkakjJIDAPB2KG41nqD9IeVNgBoKRQydU8hB9EXVIKkiPvb/mn8P2S
> nzqfB7FoPDzI+pnwnSOP0ScyL1hFCBAB63Pr6yuxDqX+VJte3OlP1exMHukHPe2LqXpP+yvb/CzC
> nFPjRVH0fVVh91IjE3iOcUnKi+tb7TDoEUrBsySTJD3e/3PNg0MNUSrWZA61fUL1PSuXlfvCNFXq
> SC9MGopgwNQntgvD0STMwbMe0rEBCcOChX0wwMzsseBvgW+9dn5+5Ld4puWDAu8Xoerwr5FVSlSF
> UZ7xyCm891IxLXpU/mRsMaEaipwgOK6uopSU7YfOKkVPyUhSPlGB/iTwFhFggDS1CgSQwiYAYHjX
> +oK/WXqmABuCMk5AJwQ34X+9tlth9u0wEzKByiBH/i7kinChID26HzI=
> _______________________________________________
> Mailing list: https://launchpad.net/~fenics
> Post to : fenics@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~fenics
> More help : https://help.launchpad.net/ListHelp
--
Anders
References
