fenics team mailing list archive
-
fenics team
-
Mailing list archive
-
Message #01222
Patch for changed Python interface of Expressions
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=
Follow ups
