[cig-commits] commit: MADDs-1 compiles but pressure field is wrong
Mercurial
hg at geodynamics.org
Wed Dec 9 20:51:07 PST 2009
changeset: 94:204cfefa3897
user: Marc Spiegelman <mspieg at ldeo.columbia.edu>
date: Sat Dec 05 15:29:38 2009 -0500
files: MADDs-1/cpp/Div.form MADDs-1/cpp/Div.h MADDs-1/cpp/Div.ufl MADDs-1/cpp/Stokes_AQP2_P1.form MADDs-1/cpp/Stokes_P2P1.form MADDs-1/cpp/main.cpp
description:
MADDs-1 compiles but pressure field is wrong
diff -r d80d53c357e0 -r 204cfefa3897 MADDs-1/cpp/Div.form
--- a/MADDs-1/cpp/Div.form Sat Dec 05 08:23:03 2009 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,19 +0,0 @@
-# Copyright (c)2008 Marc Spiegelman
-# Licensed under the GNU LGPL Version 2.1
-#
-# First added: 2008
-# Last changed: 2008
-#
-# The Functional Div(u) for calculating the global Divergence
-#
-# Compile this form with FFC: ffc -l dolfin Div.form
-
-
-P2= VectorElement("Lagrange", "triangle", 2)
-
-u = Function(P2)
-n = FacetNormal("triangle")
-
-a = dot(u,n)*ds
-
-
diff -r d80d53c357e0 -r 204cfefa3897 MADDs-1/cpp/Div.h
--- a/MADDs-1/cpp/Div.h Sat Dec 05 08:23:03 2009 -0500
+++ b/MADDs-1/cpp/Div.h Sat Dec 05 15:29:38 2009 -0500
@@ -1635,1304 +1635,6 @@ public:
break;
case 1:
return new div_0_finite_element_0_1();
- break;
- }
- return 0;
- }
-
-};
-
-/// This class defines the interface for a finite element.
-
-class div_0_finite_element_1_0: public ufc::finite_element
-{
-public:
-
- /// Constructor
- div_0_finite_element_1_0() : ufc::finite_element()
- {
- // Do nothing
- }
-
- /// Destructor
- virtual ~div_0_finite_element_1_0()
- {
- // Do nothing
- }
-
- /// Return a string identifying the finite element
- virtual const char* signature() const
- {
- return "FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)";
- }
-
- /// Return the cell shape
- virtual ufc::shape cell_shape() const
- {
- return ufc::triangle;
- }
-
- /// Return the dimension of the finite element function space
- virtual unsigned int space_dimension() const
- {
- return 1;
- }
-
- /// Return the rank of the value space
- virtual unsigned int value_rank() const
- {
- return 0;
- }
-
- /// Return the dimension of the value space for axis i
- virtual unsigned int value_dimension(unsigned int i) const
- {
- return 1;
- }
-
- /// Evaluate basis function i at given point in cell
- virtual void evaluate_basis(unsigned int i,
- double* values,
- const double* coordinates,
- const ufc::cell& c) const
- {
- // Extract vertex coordinates
- const double * const * element_coordinates = c.coordinates;
-
- // Compute Jacobian of affine map from reference cell
- const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
- const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
- const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
- const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
-
- // Compute determinant of Jacobian
- const double detJ = J_00*J_11 - J_01*J_10;
-
- // Compute inverse of Jacobian
-
- // Get coordinates and map to the reference (UFC) element
- double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
- element_coordinates[0][0]*element_coordinates[2][1] +\
- J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
- double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
- element_coordinates[1][0]*element_coordinates[0][1] -\
- J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
-
- // Map coordinates to the reference square
- if (std::abs(y - 1.0) < 1e-14)
- x = -1.0;
- else
- x = 2.0 *x/(1.0 - y) - 1.0;
- y = 2.0*y - 1.0;
-
- // Reset values
- *values = 0;
-
- // Map degree of freedom to element degree of freedom
- const unsigned int dof = i;
-
- // Generate scalings
- const double scalings_y_0 = 1;
-
- // Compute psitilde_a
- const double psitilde_a_0 = 1;
-
- // Compute psitilde_bs
- const double psitilde_bs_0_0 = 1;
-
- // Compute basisvalues
- const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
-
- // Table(s) of coefficients
- static const double coefficients0[1][1] = \
- {{1.41421356237309}};
-
- // Extract relevant coefficients
- const double coeff0_0 = coefficients0[dof][0];
-
- // Compute value(s)
- *values = coeff0_0*basisvalue0;
- }
-
- /// Evaluate all basis functions at given point in cell
- virtual void evaluate_basis_all(double* values,
- const double* coordinates,
- const ufc::cell& c) const
- {
- throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
- }
-
- /// Evaluate order n derivatives of basis function i at given point in cell
- virtual void evaluate_basis_derivatives(unsigned int i,
- unsigned int n,
- double* values,
- const double* coordinates,
- const ufc::cell& c) const
- {
- // Extract vertex coordinates
- const double * const * element_coordinates = c.coordinates;
-
- // Compute Jacobian of affine map from reference cell
- const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
- const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
- const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
- const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
-
- // Compute determinant of Jacobian
- const double detJ = J_00*J_11 - J_01*J_10;
-
- // Compute inverse of Jacobian
-
- // Get coordinates and map to the reference (UFC) element
- double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
- element_coordinates[0][0]*element_coordinates[2][1] +\
- J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
- double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
- element_coordinates[1][0]*element_coordinates[0][1] -\
- J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
-
- // Map coordinates to the reference square
- if (std::abs(y - 1.0) < 1e-14)
- x = -1.0;
- else
- x = 2.0 *x/(1.0 - y) - 1.0;
- y = 2.0*y - 1.0;
-
- // Compute number of derivatives
- unsigned int num_derivatives = 1;
-
- for (unsigned int j = 0; j < n; j++)
- num_derivatives *= 2;
-
-
- // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
- unsigned int **combinations = new unsigned int *[num_derivatives];
-
- for (unsigned int j = 0; j < num_derivatives; j++)
- {
- combinations[j] = new unsigned int [n];
- for (unsigned int k = 0; k < n; k++)
- combinations[j][k] = 0;
- }
-
- // Generate combinations of derivatives
- for (unsigned int row = 1; row < num_derivatives; row++)
- {
- for (unsigned int num = 0; num < row; num++)
- {
- for (unsigned int col = n-1; col+1 > 0; col--)
- {
- if (combinations[row][col] + 1 > 1)
- combinations[row][col] = 0;
- else
- {
- combinations[row][col] += 1;
- break;
- }
- }
- }
- }
-
- // Compute inverse of Jacobian
- const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
-
- // Declare transformation matrix
- // Declare pointer to two dimensional array and initialise
- double **transform = new double *[num_derivatives];
-
- for (unsigned int j = 0; j < num_derivatives; j++)
- {
- transform[j] = new double [num_derivatives];
- for (unsigned int k = 0; k < num_derivatives; k++)
- transform[j][k] = 1;
- }
-
- // Construct transformation matrix
- for (unsigned int row = 0; row < num_derivatives; row++)
- {
- for (unsigned int col = 0; col < num_derivatives; col++)
- {
- for (unsigned int k = 0; k < n; k++)
- transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
- }
- }
-
- // Reset values
- for (unsigned int j = 0; j < 1*num_derivatives; j++)
- values[j] = 0;
-
- // Map degree of freedom to element degree of freedom
- const unsigned int dof = i;
-
- // Generate scalings
- const double scalings_y_0 = 1;
-
- // Compute psitilde_a
- const double psitilde_a_0 = 1;
-
- // Compute psitilde_bs
- const double psitilde_bs_0_0 = 1;
-
- // Compute basisvalues
- const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
-
- // Table(s) of coefficients
- static const double coefficients0[1][1] = \
- {{1.41421356237309}};
-
- // Interesting (new) part
- // Tables of derivatives of the polynomial base (transpose)
- static const double dmats0[1][1] = \
- {{0}};
-
- static const double dmats1[1][1] = \
- {{0}};
-
- // Compute reference derivatives
- // Declare pointer to array of derivatives on FIAT element
- double *derivatives = new double [num_derivatives];
-
- // Declare coefficients
- double coeff0_0 = 0;
-
- // Declare new coefficients
- double new_coeff0_0 = 0;
-
- // Loop possible derivatives
- for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
- {
- // Get values from coefficients array
- new_coeff0_0 = coefficients0[dof][0];
-
- // Loop derivative order
- for (unsigned int j = 0; j < n; j++)
- {
- // Update old coefficients
- coeff0_0 = new_coeff0_0;
-
- if(combinations[deriv_num][j] == 0)
- {
- new_coeff0_0 = coeff0_0*dmats0[0][0];
- }
- if(combinations[deriv_num][j] == 1)
- {
- new_coeff0_0 = coeff0_0*dmats1[0][0];
- }
-
- }
- // Compute derivatives on reference element as dot product of coefficients and basisvalues
- derivatives[deriv_num] = new_coeff0_0*basisvalue0;
- }
-
- // Transform derivatives back to physical element
- for (unsigned int row = 0; row < num_derivatives; row++)
- {
- for (unsigned int col = 0; col < num_derivatives; col++)
- {
- values[row] += transform[row][col]*derivatives[col];
- }
- }
- // Delete pointer to array of derivatives on FIAT element
- delete [] derivatives;
-
- // Delete pointer to array of combinations of derivatives and transform
- for (unsigned int row = 0; row < num_derivatives; row++)
- {
- delete [] combinations[row];
- delete [] transform[row];
- }
-
- delete [] combinations;
- delete [] transform;
- }
-
- /// Evaluate order n derivatives of all basis functions at given point in cell
- virtual void evaluate_basis_derivatives_all(unsigned int n,
- double* values,
- const double* coordinates,
- const ufc::cell& c) const
- {
- throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
- }
-
- /// Evaluate linear functional for dof i on the function f
- virtual double evaluate_dof(unsigned int i,
- const ufc::function& f,
- const ufc::cell& c) const
- {
- // The reference points, direction and weights:
- static const double X[1][1][2] = {{{0.333333333333333, 0.333333333333333}}};
- static const double W[1][1] = {{1}};
- static const double D[1][1][1] = {{{1}}};
-
- const double * const * x = c.coordinates;
- double result = 0.0;
- // Iterate over the points:
- // Evaluate basis functions for affine mapping
- const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
- const double w1 = X[i][0][0];
- const double w2 = X[i][0][1];
-
- // Compute affine mapping y = F(X)
- double y[2];
- y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
- y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
-
- // Evaluate function at physical points
- double values[1];
- f.evaluate(values, y, c);
-
- // Map function values using appropriate mapping
- // Affine map: Do nothing
-
- // Note that we do not map the weights (yet).
-
- // Take directional components
- for(int k = 0; k < 1; k++)
- result += values[k]*D[i][0][k];
- // Multiply by weights
- result *= W[i][0];
-
- return result;
- }
-
- /// Evaluate linear functionals for all dofs on the function f
- virtual void evaluate_dofs(double* values,
- const ufc::function& f,
- const ufc::cell& c) const
- {
- throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
- }
-
- /// Interpolate vertex values from dof values
- virtual void interpolate_vertex_values(double* vertex_values,
- const double* dof_values,
- const ufc::cell& c) const
- {
- // Evaluate at vertices and use affine mapping
- vertex_values[0] = dof_values[0];
- vertex_values[1] = dof_values[0];
- vertex_values[2] = dof_values[0];
- }
-
- /// Return the number of sub elements (for a mixed element)
- virtual unsigned int num_sub_elements() const
- {
- return 1;
- }
-
- /// Create a new finite element for sub element i (for a mixed element)
- virtual ufc::finite_element* create_sub_element(unsigned int i) const
- {
- return new div_0_finite_element_1_0();
- }
-
-};
-
-/// This class defines the interface for a finite element.
-
-class div_0_finite_element_1_1: public ufc::finite_element
-{
-public:
-
- /// Constructor
- div_0_finite_element_1_1() : ufc::finite_element()
- {
- // Do nothing
- }
-
- /// Destructor
- virtual ~div_0_finite_element_1_1()
- {
- // Do nothing
- }
-
- /// Return a string identifying the finite element
- virtual const char* signature() const
- {
- return "FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)";
- }
-
- /// Return the cell shape
- virtual ufc::shape cell_shape() const
- {
- return ufc::triangle;
- }
-
- /// Return the dimension of the finite element function space
- virtual unsigned int space_dimension() const
- {
- return 1;
- }
-
- /// Return the rank of the value space
- virtual unsigned int value_rank() const
- {
- return 0;
- }
-
- /// Return the dimension of the value space for axis i
- virtual unsigned int value_dimension(unsigned int i) const
- {
- return 1;
- }
-
- /// Evaluate basis function i at given point in cell
- virtual void evaluate_basis(unsigned int i,
- double* values,
- const double* coordinates,
- const ufc::cell& c) const
- {
- // Extract vertex coordinates
- const double * const * element_coordinates = c.coordinates;
-
- // Compute Jacobian of affine map from reference cell
- const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
- const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
- const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
- const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
-
- // Compute determinant of Jacobian
- const double detJ = J_00*J_11 - J_01*J_10;
-
- // Compute inverse of Jacobian
-
- // Get coordinates and map to the reference (UFC) element
- double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
- element_coordinates[0][0]*element_coordinates[2][1] +\
- J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
- double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
- element_coordinates[1][0]*element_coordinates[0][1] -\
- J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
-
- // Map coordinates to the reference square
- if (std::abs(y - 1.0) < 1e-14)
- x = -1.0;
- else
- x = 2.0 *x/(1.0 - y) - 1.0;
- y = 2.0*y - 1.0;
-
- // Reset values
- *values = 0;
-
- // Map degree of freedom to element degree of freedom
- const unsigned int dof = i;
-
- // Generate scalings
- const double scalings_y_0 = 1;
-
- // Compute psitilde_a
- const double psitilde_a_0 = 1;
-
- // Compute psitilde_bs
- const double psitilde_bs_0_0 = 1;
-
- // Compute basisvalues
- const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
-
- // Table(s) of coefficients
- static const double coefficients0[1][1] = \
- {{1.41421356237309}};
-
- // Extract relevant coefficients
- const double coeff0_0 = coefficients0[dof][0];
-
- // Compute value(s)
- *values = coeff0_0*basisvalue0;
- }
-
- /// Evaluate all basis functions at given point in cell
- virtual void evaluate_basis_all(double* values,
- const double* coordinates,
- const ufc::cell& c) const
- {
- throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
- }
-
- /// Evaluate order n derivatives of basis function i at given point in cell
- virtual void evaluate_basis_derivatives(unsigned int i,
- unsigned int n,
- double* values,
- const double* coordinates,
- const ufc::cell& c) const
- {
- // Extract vertex coordinates
- const double * const * element_coordinates = c.coordinates;
-
- // Compute Jacobian of affine map from reference cell
- const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
- const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
- const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
- const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
-
- // Compute determinant of Jacobian
- const double detJ = J_00*J_11 - J_01*J_10;
-
- // Compute inverse of Jacobian
-
- // Get coordinates and map to the reference (UFC) element
- double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
- element_coordinates[0][0]*element_coordinates[2][1] +\
- J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
- double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
- element_coordinates[1][0]*element_coordinates[0][1] -\
- J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
-
- // Map coordinates to the reference square
- if (std::abs(y - 1.0) < 1e-14)
- x = -1.0;
- else
- x = 2.0 *x/(1.0 - y) - 1.0;
- y = 2.0*y - 1.0;
-
- // Compute number of derivatives
- unsigned int num_derivatives = 1;
-
- for (unsigned int j = 0; j < n; j++)
- num_derivatives *= 2;
-
-
- // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
- unsigned int **combinations = new unsigned int *[num_derivatives];
-
- for (unsigned int j = 0; j < num_derivatives; j++)
- {
- combinations[j] = new unsigned int [n];
- for (unsigned int k = 0; k < n; k++)
- combinations[j][k] = 0;
- }
-
- // Generate combinations of derivatives
- for (unsigned int row = 1; row < num_derivatives; row++)
- {
- for (unsigned int num = 0; num < row; num++)
- {
- for (unsigned int col = n-1; col+1 > 0; col--)
- {
- if (combinations[row][col] + 1 > 1)
- combinations[row][col] = 0;
- else
- {
- combinations[row][col] += 1;
- break;
- }
- }
- }
- }
-
- // Compute inverse of Jacobian
- const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
-
- // Declare transformation matrix
- // Declare pointer to two dimensional array and initialise
- double **transform = new double *[num_derivatives];
-
- for (unsigned int j = 0; j < num_derivatives; j++)
- {
- transform[j] = new double [num_derivatives];
- for (unsigned int k = 0; k < num_derivatives; k++)
- transform[j][k] = 1;
- }
-
- // Construct transformation matrix
- for (unsigned int row = 0; row < num_derivatives; row++)
- {
- for (unsigned int col = 0; col < num_derivatives; col++)
- {
- for (unsigned int k = 0; k < n; k++)
- transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
- }
- }
-
- // Reset values
- for (unsigned int j = 0; j < 1*num_derivatives; j++)
- values[j] = 0;
-
- // Map degree of freedom to element degree of freedom
- const unsigned int dof = i;
-
- // Generate scalings
- const double scalings_y_0 = 1;
-
- // Compute psitilde_a
- const double psitilde_a_0 = 1;
-
- // Compute psitilde_bs
- const double psitilde_bs_0_0 = 1;
-
- // Compute basisvalues
- const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
-
- // Table(s) of coefficients
- static const double coefficients0[1][1] = \
- {{1.41421356237309}};
-
- // Interesting (new) part
- // Tables of derivatives of the polynomial base (transpose)
- static const double dmats0[1][1] = \
- {{0}};
-
- static const double dmats1[1][1] = \
- {{0}};
-
- // Compute reference derivatives
- // Declare pointer to array of derivatives on FIAT element
- double *derivatives = new double [num_derivatives];
-
- // Declare coefficients
- double coeff0_0 = 0;
-
- // Declare new coefficients
- double new_coeff0_0 = 0;
-
- // Loop possible derivatives
- for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
- {
- // Get values from coefficients array
- new_coeff0_0 = coefficients0[dof][0];
-
- // Loop derivative order
- for (unsigned int j = 0; j < n; j++)
- {
- // Update old coefficients
- coeff0_0 = new_coeff0_0;
-
- if(combinations[deriv_num][j] == 0)
- {
- new_coeff0_0 = coeff0_0*dmats0[0][0];
- }
- if(combinations[deriv_num][j] == 1)
- {
- new_coeff0_0 = coeff0_0*dmats1[0][0];
- }
-
- }
- // Compute derivatives on reference element as dot product of coefficients and basisvalues
- derivatives[deriv_num] = new_coeff0_0*basisvalue0;
- }
-
- // Transform derivatives back to physical element
- for (unsigned int row = 0; row < num_derivatives; row++)
- {
- for (unsigned int col = 0; col < num_derivatives; col++)
- {
- values[row] += transform[row][col]*derivatives[col];
- }
- }
- // Delete pointer to array of derivatives on FIAT element
- delete [] derivatives;
-
- // Delete pointer to array of combinations of derivatives and transform
- for (unsigned int row = 0; row < num_derivatives; row++)
- {
- delete [] combinations[row];
- delete [] transform[row];
- }
-
- delete [] combinations;
- delete [] transform;
- }
-
- /// Evaluate order n derivatives of all basis functions at given point in cell
- virtual void evaluate_basis_derivatives_all(unsigned int n,
- double* values,
- const double* coordinates,
- const ufc::cell& c) const
- {
- throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
- }
-
- /// Evaluate linear functional for dof i on the function f
- virtual double evaluate_dof(unsigned int i,
- const ufc::function& f,
- const ufc::cell& c) const
- {
- // The reference points, direction and weights:
- static const double X[1][1][2] = {{{0.333333333333333, 0.333333333333333}}};
- static const double W[1][1] = {{1}};
- static const double D[1][1][1] = {{{1}}};
-
- const double * const * x = c.coordinates;
- double result = 0.0;
- // Iterate over the points:
- // Evaluate basis functions for affine mapping
- const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
- const double w1 = X[i][0][0];
- const double w2 = X[i][0][1];
-
- // Compute affine mapping y = F(X)
- double y[2];
- y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
- y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
-
- // Evaluate function at physical points
- double values[1];
- f.evaluate(values, y, c);
-
- // Map function values using appropriate mapping
- // Affine map: Do nothing
-
- // Note that we do not map the weights (yet).
-
- // Take directional components
- for(int k = 0; k < 1; k++)
- result += values[k]*D[i][0][k];
- // Multiply by weights
- result *= W[i][0];
-
- return result;
- }
-
- /// Evaluate linear functionals for all dofs on the function f
- virtual void evaluate_dofs(double* values,
- const ufc::function& f,
- const ufc::cell& c) const
- {
- throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
- }
-
- /// Interpolate vertex values from dof values
- virtual void interpolate_vertex_values(double* vertex_values,
- const double* dof_values,
- const ufc::cell& c) const
- {
- // Evaluate at vertices and use affine mapping
- vertex_values[0] = dof_values[0];
- vertex_values[1] = dof_values[0];
- vertex_values[2] = dof_values[0];
- }
-
- /// Return the number of sub elements (for a mixed element)
- virtual unsigned int num_sub_elements() const
- {
- return 1;
- }
-
- /// Create a new finite element for sub element i (for a mixed element)
- virtual ufc::finite_element* create_sub_element(unsigned int i) const
- {
- return new div_0_finite_element_1_1();
- }
-
-};
-
-/// This class defines the interface for a finite element.
-
-class div_0_finite_element_1: public ufc::finite_element
-{
-public:
-
- /// Constructor
- div_0_finite_element_1() : ufc::finite_element()
- {
- // Do nothing
- }
-
- /// Destructor
- virtual ~div_0_finite_element_1()
- {
- // Do nothing
- }
-
- /// Return a string identifying the finite element
- virtual const char* signature() const
- {
- return "VectorElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0, 2)";
- }
-
- /// Return the cell shape
- virtual ufc::shape cell_shape() const
- {
- return ufc::triangle;
- }
-
- /// Return the dimension of the finite element function space
- virtual unsigned int space_dimension() const
- {
- return 2;
- }
-
- /// Return the rank of the value space
- virtual unsigned int value_rank() const
- {
- return 1;
- }
-
- /// Return the dimension of the value space for axis i
- virtual unsigned int value_dimension(unsigned int i) const
- {
- return 2;
- }
-
- /// Evaluate basis function i at given point in cell
- virtual void evaluate_basis(unsigned int i,
- double* values,
- const double* coordinates,
- const ufc::cell& c) const
- {
- // Extract vertex coordinates
- const double * const * element_coordinates = c.coordinates;
-
- // Compute Jacobian of affine map from reference cell
- const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
- const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
- const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
- const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
-
- // Compute determinant of Jacobian
- const double detJ = J_00*J_11 - J_01*J_10;
-
- // Compute inverse of Jacobian
-
- // Get coordinates and map to the reference (UFC) element
- double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
- element_coordinates[0][0]*element_coordinates[2][1] +\
- J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
- double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
- element_coordinates[1][0]*element_coordinates[0][1] -\
- J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
-
- // Map coordinates to the reference square
- if (std::abs(y - 1.0) < 1e-14)
- x = -1.0;
- else
- x = 2.0 *x/(1.0 - y) - 1.0;
- y = 2.0*y - 1.0;
-
- // Reset values
- values[0] = 0;
- values[1] = 0;
-
- if (0 <= i && i <= 0)
- {
- // Map degree of freedom to element degree of freedom
- const unsigned int dof = i;
-
- // Generate scalings
- const double scalings_y_0 = 1;
-
- // Compute psitilde_a
- const double psitilde_a_0 = 1;
-
- // Compute psitilde_bs
- const double psitilde_bs_0_0 = 1;
-
- // Compute basisvalues
- const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
-
- // Table(s) of coefficients
- static const double coefficients0[1][1] = \
- {{1.41421356237309}};
-
- // Extract relevant coefficients
- const double coeff0_0 = coefficients0[dof][0];
-
- // Compute value(s)
- values[0] = coeff0_0*basisvalue0;
- }
-
- if (1 <= i && i <= 1)
- {
- // Map degree of freedom to element degree of freedom
- const unsigned int dof = i - 1;
-
- // Generate scalings
- const double scalings_y_0 = 1;
-
- // Compute psitilde_a
- const double psitilde_a_0 = 1;
-
- // Compute psitilde_bs
- const double psitilde_bs_0_0 = 1;
-
- // Compute basisvalues
- const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
-
- // Table(s) of coefficients
- static const double coefficients0[1][1] = \
- {{1.41421356237309}};
-
- // Extract relevant coefficients
- const double coeff0_0 = coefficients0[dof][0];
-
- // Compute value(s)
- values[1] = coeff0_0*basisvalue0;
- }
-
- }
-
- /// Evaluate all basis functions at given point in cell
- virtual void evaluate_basis_all(double* values,
- const double* coordinates,
- const ufc::cell& c) const
- {
- throw std::runtime_error("The vectorised version of evaluate_basis() is not yet implemented.");
- }
-
- /// Evaluate order n derivatives of basis function i at given point in cell
- virtual void evaluate_basis_derivatives(unsigned int i,
- unsigned int n,
- double* values,
- const double* coordinates,
- const ufc::cell& c) const
- {
- // Extract vertex coordinates
- const double * const * element_coordinates = c.coordinates;
-
- // Compute Jacobian of affine map from reference cell
- const double J_00 = element_coordinates[1][0] - element_coordinates[0][0];
- const double J_01 = element_coordinates[2][0] - element_coordinates[0][0];
- const double J_10 = element_coordinates[1][1] - element_coordinates[0][1];
- const double J_11 = element_coordinates[2][1] - element_coordinates[0][1];
-
- // Compute determinant of Jacobian
- const double detJ = J_00*J_11 - J_01*J_10;
-
- // Compute inverse of Jacobian
-
- // Get coordinates and map to the reference (UFC) element
- double x = (element_coordinates[0][1]*element_coordinates[2][0] -\
- element_coordinates[0][0]*element_coordinates[2][1] +\
- J_11*coordinates[0] - J_01*coordinates[1]) / detJ;
- double y = (element_coordinates[1][1]*element_coordinates[0][0] -\
- element_coordinates[1][0]*element_coordinates[0][1] -\
- J_10*coordinates[0] + J_00*coordinates[1]) / detJ;
-
- // Map coordinates to the reference square
- if (std::abs(y - 1.0) < 1e-14)
- x = -1.0;
- else
- x = 2.0 *x/(1.0 - y) - 1.0;
- y = 2.0*y - 1.0;
-
- // Compute number of derivatives
- unsigned int num_derivatives = 1;
-
- for (unsigned int j = 0; j < n; j++)
- num_derivatives *= 2;
-
-
- // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
- unsigned int **combinations = new unsigned int *[num_derivatives];
-
- for (unsigned int j = 0; j < num_derivatives; j++)
- {
- combinations[j] = new unsigned int [n];
- for (unsigned int k = 0; k < n; k++)
- combinations[j][k] = 0;
- }
-
- // Generate combinations of derivatives
- for (unsigned int row = 1; row < num_derivatives; row++)
- {
- for (unsigned int num = 0; num < row; num++)
- {
- for (unsigned int col = n-1; col+1 > 0; col--)
- {
- if (combinations[row][col] + 1 > 1)
- combinations[row][col] = 0;
- else
- {
- combinations[row][col] += 1;
- break;
- }
- }
- }
- }
-
- // Compute inverse of Jacobian
- const double Jinv[2][2] = {{J_11 / detJ, -J_01 / detJ}, {-J_10 / detJ, J_00 / detJ}};
-
- // Declare transformation matrix
- // Declare pointer to two dimensional array and initialise
- double **transform = new double *[num_derivatives];
-
- for (unsigned int j = 0; j < num_derivatives; j++)
- {
- transform[j] = new double [num_derivatives];
- for (unsigned int k = 0; k < num_derivatives; k++)
- transform[j][k] = 1;
- }
-
- // Construct transformation matrix
- for (unsigned int row = 0; row < num_derivatives; row++)
- {
- for (unsigned int col = 0; col < num_derivatives; col++)
- {
- for (unsigned int k = 0; k < n; k++)
- transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
- }
- }
-
- // Reset values
- for (unsigned int j = 0; j < 2*num_derivatives; j++)
- values[j] = 0;
-
- if (0 <= i && i <= 0)
- {
- // Map degree of freedom to element degree of freedom
- const unsigned int dof = i;
-
- // Generate scalings
- const double scalings_y_0 = 1;
-
- // Compute psitilde_a
- const double psitilde_a_0 = 1;
-
- // Compute psitilde_bs
- const double psitilde_bs_0_0 = 1;
-
- // Compute basisvalues
- const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
-
- // Table(s) of coefficients
- static const double coefficients0[1][1] = \
- {{1.41421356237309}};
-
- // Interesting (new) part
- // Tables of derivatives of the polynomial base (transpose)
- static const double dmats0[1][1] = \
- {{0}};
-
- static const double dmats1[1][1] = \
- {{0}};
-
- // Compute reference derivatives
- // Declare pointer to array of derivatives on FIAT element
- double *derivatives = new double [num_derivatives];
-
- // Declare coefficients
- double coeff0_0 = 0;
-
- // Declare new coefficients
- double new_coeff0_0 = 0;
-
- // Loop possible derivatives
- for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
- {
- // Get values from coefficients array
- new_coeff0_0 = coefficients0[dof][0];
-
- // Loop derivative order
- for (unsigned int j = 0; j < n; j++)
- {
- // Update old coefficients
- coeff0_0 = new_coeff0_0;
-
- if(combinations[deriv_num][j] == 0)
- {
- new_coeff0_0 = coeff0_0*dmats0[0][0];
- }
- if(combinations[deriv_num][j] == 1)
- {
- new_coeff0_0 = coeff0_0*dmats1[0][0];
- }
-
- }
- // Compute derivatives on reference element as dot product of coefficients and basisvalues
- derivatives[deriv_num] = new_coeff0_0*basisvalue0;
- }
-
- // Transform derivatives back to physical element
- for (unsigned int row = 0; row < num_derivatives; row++)
- {
- for (unsigned int col = 0; col < num_derivatives; col++)
- {
- values[row] += transform[row][col]*derivatives[col];
- }
- }
- // Delete pointer to array of derivatives on FIAT element
- delete [] derivatives;
-
- // Delete pointer to array of combinations of derivatives and transform
- for (unsigned int row = 0; row < num_derivatives; row++)
- {
- delete [] combinations[row];
- delete [] transform[row];
- }
-
- delete [] combinations;
- delete [] transform;
- }
-
- if (1 <= i && i <= 1)
- {
- // Map degree of freedom to element degree of freedom
- const unsigned int dof = i - 1;
-
- // Generate scalings
- const double scalings_y_0 = 1;
-
- // Compute psitilde_a
- const double psitilde_a_0 = 1;
-
- // Compute psitilde_bs
- const double psitilde_bs_0_0 = 1;
-
- // Compute basisvalues
- const double basisvalue0 = 0.707106781186548*psitilde_a_0*scalings_y_0*psitilde_bs_0_0;
-
- // Table(s) of coefficients
- static const double coefficients0[1][1] = \
- {{1.41421356237309}};
-
- // Interesting (new) part
- // Tables of derivatives of the polynomial base (transpose)
- static const double dmats0[1][1] = \
- {{0}};
-
- static const double dmats1[1][1] = \
- {{0}};
-
- // Compute reference derivatives
- // Declare pointer to array of derivatives on FIAT element
- double *derivatives = new double [num_derivatives];
-
- // Declare coefficients
- double coeff0_0 = 0;
-
- // Declare new coefficients
- double new_coeff0_0 = 0;
-
- // Loop possible derivatives
- for (unsigned int deriv_num = 0; deriv_num < num_derivatives; deriv_num++)
- {
- // Get values from coefficients array
- new_coeff0_0 = coefficients0[dof][0];
-
- // Loop derivative order
- for (unsigned int j = 0; j < n; j++)
- {
- // Update old coefficients
- coeff0_0 = new_coeff0_0;
-
- if(combinations[deriv_num][j] == 0)
- {
- new_coeff0_0 = coeff0_0*dmats0[0][0];
- }
- if(combinations[deriv_num][j] == 1)
- {
- new_coeff0_0 = coeff0_0*dmats1[0][0];
- }
-
- }
- // Compute derivatives on reference element as dot product of coefficients and basisvalues
- derivatives[deriv_num] = new_coeff0_0*basisvalue0;
- }
-
- // Transform derivatives back to physical element
- for (unsigned int row = 0; row < num_derivatives; row++)
- {
- for (unsigned int col = 0; col < num_derivatives; col++)
- {
- values[num_derivatives + row] += transform[row][col]*derivatives[col];
- }
- }
- // Delete pointer to array of derivatives on FIAT element
- delete [] derivatives;
-
- // Delete pointer to array of combinations of derivatives and transform
- for (unsigned int row = 0; row < num_derivatives; row++)
- {
- delete [] combinations[row];
- delete [] transform[row];
- }
-
- delete [] combinations;
- delete [] transform;
- }
-
- }
-
- /// Evaluate order n derivatives of all basis functions at given point in cell
- virtual void evaluate_basis_derivatives_all(unsigned int n,
- double* values,
- const double* coordinates,
- const ufc::cell& c) const
- {
- throw std::runtime_error("The vectorised version of evaluate_basis_derivatives() is not yet implemented.");
- }
-
- /// Evaluate linear functional for dof i on the function f
- virtual double evaluate_dof(unsigned int i,
- const ufc::function& f,
- const ufc::cell& c) const
- {
- // The reference points, direction and weights:
- static const double X[2][1][2] = {{{0.333333333333333, 0.333333333333333}}, {{0.333333333333333, 0.333333333333333}}};
- static const double W[2][1] = {{1}, {1}};
- static const double D[2][1][2] = {{{1, 0}}, {{0, 1}}};
-
- const double * const * x = c.coordinates;
- double result = 0.0;
- // Iterate over the points:
- // Evaluate basis functions for affine mapping
- const double w0 = 1.0 - X[i][0][0] - X[i][0][1];
- const double w1 = X[i][0][0];
- const double w2 = X[i][0][1];
-
- // Compute affine mapping y = F(X)
- double y[2];
- y[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
- y[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];
-
- // Evaluate function at physical points
- double values[2];
- f.evaluate(values, y, c);
-
- // Map function values using appropriate mapping
- // Affine map: Do nothing
-
- // Note that we do not map the weights (yet).
-
- // Take directional components
- for(int k = 0; k < 2; k++)
- result += values[k]*D[i][0][k];
- // Multiply by weights
- result *= W[i][0];
-
- return result;
- }
-
- /// Evaluate linear functionals for all dofs on the function f
- virtual void evaluate_dofs(double* values,
- const ufc::function& f,
- const ufc::cell& c) const
- {
- throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
- }
-
- /// Interpolate vertex values from dof values
- virtual void interpolate_vertex_values(double* vertex_values,
- const double* dof_values,
- const ufc::cell& c) const
- {
- // Evaluate at vertices and use affine mapping
- vertex_values[0] = dof_values[0];
- vertex_values[2] = dof_values[0];
- vertex_values[4] = dof_values[0];
- // Evaluate at vertices and use affine mapping
- vertex_values[1] = dof_values[1];
- vertex_values[3] = dof_values[1];
- vertex_values[5] = dof_values[1];
- }
-
- /// Return the number of sub elements (for a mixed element)
- virtual unsigned int num_sub_elements() const
- {
- return 2;
- }
-
- /// Create a new finite element for sub element i (for a mixed element)
- virtual ufc::finite_element* create_sub_element(unsigned int i) const
- {
- switch ( i )
- {
- case 0:
- return new div_0_finite_element_1_0();
- break;
- case 1:
- return new div_0_finite_element_1_1();
break;
}
return 0;
@@ -3521,496 +2223,6 @@ public:
};
-/// This class defines the interface for a local-to-global mapping of
-/// degrees of freedom (dofs).
-
-class div_0_dof_map_1_0: public ufc::dof_map
-{
-private:
-
- unsigned int __global_dimension;
-
-public:
-
- /// Constructor
- div_0_dof_map_1_0() : ufc::dof_map()
- {
- __global_dimension = 0;
- }
-
- /// Destructor
- virtual ~div_0_dof_map_1_0()
- {
- // Do nothing
- }
-
- /// Return a string identifying the dof map
- virtual const char* signature() const
- {
- return "FFC dof map for FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)";
- }
-
- /// Return true iff mesh entities of topological dimension d are needed
- virtual bool needs_mesh_entities(unsigned int d) const
- {
- switch ( d )
- {
- case 0:
- return false;
- break;
- case 1:
- return false;
- break;
- case 2:
- return true;
- break;
- }
- return false;
- }
-
- /// Initialize dof map for mesh (return true iff init_cell() is needed)
- virtual bool init_mesh(const ufc::mesh& m)
- {
- __global_dimension = m.num_entities[2];
- return false;
- }
-
- /// Initialize dof map for given cell
- virtual void init_cell(const ufc::mesh& m,
- const ufc::cell& c)
- {
- // Do nothing
- }
-
- /// Finish initialization of dof map for cells
- virtual void init_cell_finalize()
- {
- // Do nothing
- }
-
- /// Return the dimension of the global finite element function space
- virtual unsigned int global_dimension() const
- {
- return __global_dimension;
- }
-
- /// Return the dimension of the local finite element function space for a cell
- virtual unsigned int local_dimension(const ufc::cell& c) const
- {
- return 1;
- }
-
- /// Return the maximum dimension of the local finite element function space
- virtual unsigned int max_local_dimension() const
- {
- return 1;
- }
-
- // Return the geometric dimension of the coordinates this dof map provides
- virtual unsigned int geometric_dimension() const
- {
- return 2;
- }
-
- /// Return the number of dofs on each cell facet
- virtual unsigned int num_facet_dofs() const
- {
- return 0;
- }
-
- /// Return the number of dofs associated with each cell entity of dimension d
- virtual unsigned int num_entity_dofs(unsigned int d) const
- {
- throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
- }
-
- /// Tabulate the local-to-global mapping of dofs on a cell
- virtual void tabulate_dofs(unsigned int* dofs,
- const ufc::mesh& m,
- const ufc::cell& c) const
- {
- dofs[0] = c.entity_indices[2][0];
- }
-
- /// Tabulate the local-to-local mapping from facet dofs to cell dofs
- virtual void tabulate_facet_dofs(unsigned int* dofs,
- unsigned int facet) const
- {
- switch ( facet )
- {
- case 0:
-
- break;
- case 1:
-
- break;
- case 2:
-
- break;
- }
- }
-
- /// Tabulate the local-to-local mapping of dofs on entity (d, i)
- virtual void tabulate_entity_dofs(unsigned int* dofs,
- unsigned int d, unsigned int i) const
- {
- throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
- }
-
- /// Tabulate the coordinates of all dofs on a cell
- virtual void tabulate_coordinates(double** coordinates,
- const ufc::cell& c) const
- {
- const double * const * x = c.coordinates;
- coordinates[0][0] = 0.333333333333333*x[0][0] + 0.333333333333333*x[1][0] + 0.333333333333333*x[2][0];
- coordinates[0][1] = 0.333333333333333*x[0][1] + 0.333333333333333*x[1][1] + 0.333333333333333*x[2][1];
- }
-
- /// Return the number of sub dof maps (for a mixed element)
- virtual unsigned int num_sub_dof_maps() const
- {
- return 1;
- }
-
- /// Create a new dof_map for sub dof map i (for a mixed element)
- virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
- {
- return new div_0_dof_map_1_0();
- }
-
-};
-
-/// This class defines the interface for a local-to-global mapping of
-/// degrees of freedom (dofs).
-
-class div_0_dof_map_1_1: public ufc::dof_map
-{
-private:
-
- unsigned int __global_dimension;
-
-public:
-
- /// Constructor
- div_0_dof_map_1_1() : ufc::dof_map()
- {
- __global_dimension = 0;
- }
-
- /// Destructor
- virtual ~div_0_dof_map_1_1()
- {
- // Do nothing
- }
-
- /// Return a string identifying the dof map
- virtual const char* signature() const
- {
- return "FFC dof map for FiniteElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0)";
- }
-
- /// Return true iff mesh entities of topological dimension d are needed
- virtual bool needs_mesh_entities(unsigned int d) const
- {
- switch ( d )
- {
- case 0:
- return false;
- break;
- case 1:
- return false;
- break;
- case 2:
- return true;
- break;
- }
- return false;
- }
-
- /// Initialize dof map for mesh (return true iff init_cell() is needed)
- virtual bool init_mesh(const ufc::mesh& m)
- {
- __global_dimension = m.num_entities[2];
- return false;
- }
-
- /// Initialize dof map for given cell
- virtual void init_cell(const ufc::mesh& m,
- const ufc::cell& c)
- {
- // Do nothing
- }
-
- /// Finish initialization of dof map for cells
- virtual void init_cell_finalize()
- {
- // Do nothing
- }
-
- /// Return the dimension of the global finite element function space
- virtual unsigned int global_dimension() const
- {
- return __global_dimension;
- }
-
- /// Return the dimension of the local finite element function space for a cell
- virtual unsigned int local_dimension(const ufc::cell& c) const
- {
- return 1;
- }
-
- /// Return the maximum dimension of the local finite element function space
- virtual unsigned int max_local_dimension() const
- {
- return 1;
- }
-
- // Return the geometric dimension of the coordinates this dof map provides
- virtual unsigned int geometric_dimension() const
- {
- return 2;
- }
-
- /// Return the number of dofs on each cell facet
- virtual unsigned int num_facet_dofs() const
- {
- return 0;
- }
-
- /// Return the number of dofs associated with each cell entity of dimension d
- virtual unsigned int num_entity_dofs(unsigned int d) const
- {
- throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
- }
-
- /// Tabulate the local-to-global mapping of dofs on a cell
- virtual void tabulate_dofs(unsigned int* dofs,
- const ufc::mesh& m,
- const ufc::cell& c) const
- {
- dofs[0] = c.entity_indices[2][0];
- }
-
- /// Tabulate the local-to-local mapping from facet dofs to cell dofs
- virtual void tabulate_facet_dofs(unsigned int* dofs,
- unsigned int facet) const
- {
- switch ( facet )
- {
- case 0:
-
- break;
- case 1:
-
- break;
- case 2:
-
- break;
- }
- }
-
- /// Tabulate the local-to-local mapping of dofs on entity (d, i)
- virtual void tabulate_entity_dofs(unsigned int* dofs,
- unsigned int d, unsigned int i) const
- {
- throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
- }
-
- /// Tabulate the coordinates of all dofs on a cell
- virtual void tabulate_coordinates(double** coordinates,
- const ufc::cell& c) const
- {
- const double * const * x = c.coordinates;
- coordinates[0][0] = 0.333333333333333*x[0][0] + 0.333333333333333*x[1][0] + 0.333333333333333*x[2][0];
- coordinates[0][1] = 0.333333333333333*x[0][1] + 0.333333333333333*x[1][1] + 0.333333333333333*x[2][1];
- }
-
- /// Return the number of sub dof maps (for a mixed element)
- virtual unsigned int num_sub_dof_maps() const
- {
- return 1;
- }
-
- /// Create a new dof_map for sub dof map i (for a mixed element)
- virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
- {
- return new div_0_dof_map_1_1();
- }
-
-};
-
-/// This class defines the interface for a local-to-global mapping of
-/// degrees of freedom (dofs).
-
-class div_0_dof_map_1: public ufc::dof_map
-{
-private:
-
- unsigned int __global_dimension;
-
-public:
-
- /// Constructor
- div_0_dof_map_1() : ufc::dof_map()
- {
- __global_dimension = 0;
- }
-
- /// Destructor
- virtual ~div_0_dof_map_1()
- {
- // Do nothing
- }
-
- /// Return a string identifying the dof map
- virtual const char* signature() const
- {
- return "FFC dof map for VectorElement('Discontinuous Lagrange', Cell('triangle', 1, Space(2)), 0, 2)";
- }
-
- /// Return true iff mesh entities of topological dimension d are needed
- virtual bool needs_mesh_entities(unsigned int d) const
- {
- switch ( d )
- {
- case 0:
- return false;
- break;
- case 1:
- return false;
- break;
- case 2:
- return true;
- break;
- }
- return false;
- }
-
- /// Initialize dof map for mesh (return true iff init_cell() is needed)
- virtual bool init_mesh(const ufc::mesh& m)
- {
- __global_dimension = 2*m.num_entities[2];
- return false;
- }
-
- /// Initialize dof map for given cell
- virtual void init_cell(const ufc::mesh& m,
- const ufc::cell& c)
- {
- // Do nothing
- }
-
- /// Finish initialization of dof map for cells
- virtual void init_cell_finalize()
- {
- // Do nothing
- }
-
- /// Return the dimension of the global finite element function space
- virtual unsigned int global_dimension() const
- {
- return __global_dimension;
- }
-
- /// Return the dimension of the local finite element function space for a cell
- virtual unsigned int local_dimension(const ufc::cell& c) const
- {
- return 2;
- }
-
- /// Return the maximum dimension of the local finite element function space
- virtual unsigned int max_local_dimension() const
- {
- return 2;
- }
-
- // Return the geometric dimension of the coordinates this dof map provides
- virtual unsigned int geometric_dimension() const
- {
- return 2;
- }
-
- /// Return the number of dofs on each cell facet
- virtual unsigned int num_facet_dofs() const
- {
- return 0;
- }
-
- /// Return the number of dofs associated with each cell entity of dimension d
- virtual unsigned int num_entity_dofs(unsigned int d) const
- {
- throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
- }
-
- /// Tabulate the local-to-global mapping of dofs on a cell
- virtual void tabulate_dofs(unsigned int* dofs,
- const ufc::mesh& m,
- const ufc::cell& c) const
- {
- dofs[0] = c.entity_indices[2][0];
- unsigned int offset = m.num_entities[2];
- dofs[1] = offset + c.entity_indices[2][0];
- }
-
- /// Tabulate the local-to-local mapping from facet dofs to cell dofs
- virtual void tabulate_facet_dofs(unsigned int* dofs,
- unsigned int facet) const
- {
- switch ( facet )
- {
- case 0:
-
- break;
- case 1:
-
- break;
- case 2:
-
- break;
- }
- }
-
- /// Tabulate the local-to-local mapping of dofs on entity (d, i)
- virtual void tabulate_entity_dofs(unsigned int* dofs,
- unsigned int d, unsigned int i) const
- {
- throw std::runtime_error("Not implemented (introduced in UFC v1.1).");
- }
-
- /// Tabulate the coordinates of all dofs on a cell
- virtual void tabulate_coordinates(double** coordinates,
- const ufc::cell& c) const
- {
- const double * const * x = c.coordinates;
- coordinates[0][0] = 0.333333333333333*x[0][0] + 0.333333333333333*x[1][0] + 0.333333333333333*x[2][0];
- coordinates[0][1] = 0.333333333333333*x[0][1] + 0.333333333333333*x[1][1] + 0.333333333333333*x[2][1];
- coordinates[1][0] = 0.333333333333333*x[0][0] + 0.333333333333333*x[1][0] + 0.333333333333333*x[2][0];
- coordinates[1][1] = 0.333333333333333*x[0][1] + 0.333333333333333*x[1][1] + 0.333333333333333*x[2][1];
- }
-
- /// Return the number of sub dof maps (for a mixed element)
- virtual unsigned int num_sub_dof_maps() const
- {
- return 2;
- }
-
- /// Create a new dof_map for sub dof map i (for a mixed element)
- virtual ufc::dof_map* create_sub_dof_map(unsigned int i) const
- {
- switch ( i )
- {
- case 0:
- return new div_0_dof_map_1_0();
- break;
- case 1:
- return new div_0_dof_map_1_1();
- break;
- }
- return 0;
- }
-
-};
-
/// This class defines the interface for the tabulation of the
/// exterior facet tensor corresponding to the local contribution to
/// a form from the integral over an exterior facet.
@@ -4058,8 +2270,11 @@ public:
const double dx1 = x[v1][1] - x[v0][1];
const double det = std::sqrt(dx0*dx0 + dx1*dx1);
+ const bool direction = dx1*(x[facet][0] - x[v0][0]) - dx0*(x[facet][1] - x[v0][1]) < 0;
// Compute facet normals from the facet scale factor constants
+ const double n0 = direction ? dx1 / det : -dx1 / det;
+ const double n1 = direction ? -dx0 / det : dx0 / det;
// Array of quadrature weights
@@ -4067,27 +2282,27 @@ public:
// Quadrature points on the UFC reference element: (0.211324865405187), (0.788675134594813)
// Value of basis functions at quadrature points.
- static const double FE1_f0_C0[2][3] = \
+ static const double FE0_f0_C0[2][3] = \
{{0.45534180126148, -0.122008467928146, 0.666666666666667},
{-0.122008467928146, 0.45534180126148, 0.666666666666667}};
// Array of non-zero columns
- static const unsigned int nzc2[3] = {1, 2, 3};
+ static const unsigned int nzc1[3] = {7, 8, 9};
// Array of non-zero columns
- static const unsigned int nzc3[3] = {7, 8, 9};
+ static const unsigned int nzc0[3] = {1, 2, 3};
// Array of non-zero columns
- static const unsigned int nzc6[3] = {0, 1, 5};
+ static const unsigned int nzc5[3] = {6, 7, 11};
// Array of non-zero columns
- static const unsigned int nzc7[3] = {6, 7, 11};
+ static const unsigned int nzc4[3] = {0, 1, 5};
// Array of non-zero columns
- static const unsigned int nzc5[3] = {6, 8, 10};
+ static const unsigned int nzc2[3] = {0, 2, 4};
// Array of non-zero columns
- static const unsigned int nzc4[3] = {0, 2, 4};
+ static const unsigned int nzc3[3] = {6, 8, 10};
// Number of operations to compute geometry constants: 2
// Should be added to total operation count.
- const double G0 = det*w[1][1];
- const double G1 = det*w[1][0];
+ const double G0 = det*n1;
+ const double G1 = det*n0;
// Compute element tensor using UFL quadrature representation
// Optimisations: ('simplify expressions', True), ('ignore zero tables', True), ('non zero columns', True), ('remove zero terms', True), ('ignore ones', True)
@@ -4109,8 +2324,8 @@ public:
// Total number of operations to compute function values = 12
for (unsigned int r = 0; r < 3; r++)
{
- F0 += FE1_f0_C0[ip][r]*w[0][nzc2[r]];
- F1 += FE1_f0_C0[ip][r]*w[0][nzc3[r]];
+ F0 += FE0_f0_C0[ip][r]*w[0][nzc0[r]];
+ F1 += FE0_f0_C0[ip][r]*w[0][nzc1[r]];
}// end loop over 'r'
// Number of operations to compute ip constants: 4
@@ -4140,8 +2355,8 @@ public:
// Total number of operations to compute function values = 12
for (unsigned int r = 0; r < 3; r++)
{
- F0 += FE1_f0_C0[ip][r]*w[0][nzc4[r]];
- F1 += FE1_f0_C0[ip][r]*w[0][nzc5[r]];
+ F0 += FE0_f0_C0[ip][r]*w[0][nzc2[r]];
+ F1 += FE0_f0_C0[ip][r]*w[0][nzc3[r]];
}// end loop over 'r'
// Number of operations to compute ip constants: 4
@@ -4171,8 +2386,8 @@ public:
// Total number of operations to compute function values = 12
for (unsigned int r = 0; r < 3; r++)
{
- F0 += FE1_f0_C0[ip][r]*w[0][nzc6[r]];
- F1 += FE1_f0_C0[ip][r]*w[0][nzc7[r]];
+ F0 += FE0_f0_C0[ip][r]*w[0][nzc4[r]];
+ F1 += FE0_f0_C0[ip][r]*w[0][nzc5[r]];
}// end loop over 'r'
// Number of operations to compute ip constants: 4
@@ -4264,7 +2479,7 @@ public:
/// Return a string identifying the form
virtual const char* signature() const
{
- return "Form([Integral(IndexSum(Product(Indexed(Function(VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 2, 2), 0), MultiIndex((Index(0),), {Index(0): 2})), Indexed(VectorConstant(Cell('triangle', 1, Space(2)), 2, 1), MultiIndex((Index(0),), {Index(0): 2}))), MultiIndex((Index(0),), {Index(0): 2})), Measure('exterior_facet', 0, None))])";
+ return "Form([Integral(IndexSum(Product(Indexed(FacetNormal(Cell('triangle', 1, Space(2))), MultiIndex((Index(0),), {Index(0): 2})), Indexed(Function(VectorElement('Lagrange', Cell('triangle', 1, Space(2)), 2, 2), 0), MultiIndex((Index(0),), {Index(0): 2}))), MultiIndex((Index(0),), {Index(0): 2})), Measure('exterior_facet', 0, None))])";
}
/// Return the rank of the global tensor (r)
@@ -4276,7 +2491,7 @@ public:
/// Return the number of coefficients (n)
virtual unsigned int num_coefficients() const
{
- return 2;
+ return 1;
}
/// Return the number of cell integrals
@@ -4300,31 +2515,13 @@ public:
/// Create a new finite element for argument function i
virtual ufc::finite_element* create_finite_element(unsigned int i) const
{
- switch ( i )
- {
- case 0:
- return new div_0_finite_element_0();
- break;
- case 1:
- return new div_0_finite_element_1();
- break;
- }
- return 0;
+ return new div_0_finite_element_0();
}
/// Create a new dof map for argument function i
virtual ufc::dof_map* create_dof_map(unsigned int i) const
{
- switch ( i )
- {
- case 0:
- return new div_0_dof_map_0();
- break;
- case 1:
- return new div_0_dof_map_1();
- break;
- }
- return 0;
+ return new div_0_dof_map_0();
}
/// Create a new cell integral on sub domain i
@@ -4364,49 +2561,6 @@ namespace Div
namespace Div
{
-class CoefficientSpace_n: public dolfin::FunctionSpace
-{
-public:
-
- CoefficientSpace_n(const dolfin::Mesh& mesh):
- dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
- boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new div_0_finite_element_1()))),
- boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new div_0_dof_map_1()), mesh)))
- {
- // Do nothing
- }
-
- CoefficientSpace_n(dolfin::Mesh& mesh):
- dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
- boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new div_0_finite_element_1()))),
- boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new div_0_dof_map_1()), mesh)))
- {
- // Do nothing
- }
-
- CoefficientSpace_n(boost::shared_ptr<dolfin::Mesh> mesh):
- dolfin::FunctionSpace(mesh,
- boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new div_0_finite_element_1()))),
- boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new div_0_dof_map_1()), *mesh)))
- {
- // Do nothing
- }
-
- CoefficientSpace_n(boost::shared_ptr<const dolfin::Mesh> mesh):
- dolfin::FunctionSpace(mesh,
- boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new div_0_finite_element_1()))),
- boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dof_map>(new div_0_dof_map_1()), *mesh)))
- {
- // Do nothing
- }
-
-
- ~CoefficientSpace_n()
- {
- }
-
-};
-
class CoefficientSpace_u: public dolfin::FunctionSpace
{
public:
@@ -4452,68 +2606,62 @@ public:
typedef CoefficientSpace_u Form_0_FunctionSpace_0;
-typedef CoefficientSpace_n Form_0_FunctionSpace_1;
-
class Form_0: public dolfin::Form
{
public:
// Constructor
Form_0(const dolfin::Mesh& mesh):
- dolfin::Form(0, 2), u(*this, 0), n(*this, 1)
+ dolfin::Form(0, 1), u(*this, 0)
{
_mesh = reference_to_no_delete_pointer(mesh);
_ufc_form = boost::shared_ptr<const ufc::form>(new div_form_0());
}
// Constructor
- Form_0(const dolfin::Mesh& mesh, const dolfin::GenericFunction& u, const dolfin::GenericFunction& n):
- dolfin::Form(0, 2), u(*this, 0), n(*this, 1)
+ Form_0(const dolfin::Mesh& mesh, const dolfin::GenericFunction& u):
+ dolfin::Form(0, 1), u(*this, 0)
{
_mesh = reference_to_no_delete_pointer(mesh);
this->u = u;
- this->n = n;
_ufc_form = boost::shared_ptr<const ufc::form>(new div_form_0());
}
// Constructor
- Form_0(const dolfin::Mesh& mesh, boost::shared_ptr<const dolfin::GenericFunction> u, boost::shared_ptr<const dolfin::GenericFunction> n):
- dolfin::Form(0, 2), u(*this, 0), n(*this, 1)
+ Form_0(const dolfin::Mesh& mesh, boost::shared_ptr<const dolfin::GenericFunction> u):
+ dolfin::Form(0, 1), u(*this, 0)
{
_mesh = reference_to_no_delete_pointer(mesh);
this->u = *u;
- this->n = *n;
_ufc_form = boost::shared_ptr<const ufc::form>(new div_form_0());
}
// Constructor
Form_0(boost::shared_ptr<const dolfin::Mesh> mesh):
- dolfin::Form(0, 2), u(*this, 0), n(*this, 1)
+ dolfin::Form(0, 1), u(*this, 0)
{
_mesh = mesh;
_ufc_form = boost::shared_ptr<const ufc::form>(new div_form_0());
}
// Constructor
- Form_0(boost::shared_ptr<const dolfin::Mesh> mesh, const dolfin::GenericFunction& u, const dolfin::GenericFunction& n):
- dolfin::Form(0, 2), u(*this, 0), n(*this, 1)
+ Form_0(boost::shared_ptr<const dolfin::Mesh> mesh, const dolfin::GenericFunction& u):
+ dolfin::Form(0, 1), u(*this, 0)
{
_mesh = mesh;
this->u = u;
- this->n = n;
_ufc_form = boost::shared_ptr<const ufc::form>(new div_form_0());
}
// Constructor
- Form_0(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::GenericFunction> u, boost::shared_ptr<const dolfin::GenericFunction> n):
- dolfin::Form(0, 2), u(*this, 0), n(*this, 1)
+ Form_0(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::GenericFunction> u):
+ dolfin::Form(0, 1), u(*this, 0)
{
_mesh = mesh;
this->u = *u;
- this->n = *n;
_ufc_form = boost::shared_ptr<const ufc::form>(new div_form_0());
}
@@ -4527,8 +2675,6 @@ public:
{
if (name == "u")
return 0;
- else if (name == "n")
- return 1;
dolfin::error("Invalid coefficient.");
return 0;
@@ -4541,8 +2687,6 @@ public:
{
case 0:
return "u";
- case 1:
- return "n";
}
dolfin::error("Invalid coefficient.");
@@ -4551,11 +2695,9 @@ public:
// Typedefs
typedef Form_0_FunctionSpace_0 CoefficientSpace_u;
- typedef Form_0_FunctionSpace_1 CoefficientSpace_n;
// Coefficients
dolfin::CoefficientAssigner u;
- dolfin::CoefficientAssigner n;
};
// Class typedefs
diff -r d80d53c357e0 -r 204cfefa3897 MADDs-1/cpp/Div.ufl
--- a/MADDs-1/cpp/Div.ufl Sat Dec 05 08:23:03 2009 -0500
+++ b/MADDs-1/cpp/Div.ufl Sat Dec 05 15:29:38 2009 -0500
@@ -12,7 +12,7 @@ P2= VectorElement("Lagrange", triangle,
P2= VectorElement("Lagrange", triangle, 2)
u = Function(P2)
-n = VectorConstant(triangle)
+n = P2.cell().n # facet normal for divergence
a = inner(u,n)*ds
diff -r d80d53c357e0 -r 204cfefa3897 MADDs-1/cpp/Stokes_AQP2_P1.form
--- a/MADDs-1/cpp/Stokes_AQP2_P1.form Sat Dec 05 08:23:03 2009 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-# Copyright (c) 2005-2007 Anders Logg
-# Licensed under the GNU LGPL Version 2.1
-#
-# First added: 2005
-# Last changed: 2007-04-23
-#
-# The bilinear form a(v, u) for the block preconditioning matrix [ A Q ]
-# equations using a mixed formulation (Taylor-Hood elements).
-#
-# Compile this form with FFC: ffc -l dolfin Stokes.form
-
-P2 = VectorElement("Lagrange", "triangle", 2)
-P1 = FiniteElement("Lagrange", "triangle", 1)
-TH = P2 + P1
-
-(v, q) = TestFunctions(TH)
-(u, p) = TrialFunctions(TH)
-
-# alpha is scaling constant (L/h)^2 where L^2=\eta\U_0/\delta\rho g
-
-alpha = Constant("triangle")
-
-#a = (alpha*dot(grad(v), grad(u)) + q*p)*dx
-a = (alpha*dot(grad(v), grad(u)) - div(v)*p + q*p)*dx
diff -r d80d53c357e0 -r 204cfefa3897 MADDs-1/cpp/Stokes_P2P1.form
--- a/MADDs-1/cpp/Stokes_P2P1.form Sat Dec 05 08:23:03 2009 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,25 +0,0 @@
-# Copyright (c) 2005-2007 Anders Logg
-# Licensed under the GNU LGPL Version 2.1
-#
-# First added: 2005
-# Last changed: 2007-04-23
-#
-# The bilinear form a(v, u) and Linear form L(v) for the Stokes
-# equations using a mixed formulation (Taylor-Hood elements).
-#
-# Compile this form with FFC: ffc -l dolfin Stokes.form
-
-P2 = VectorElement("Lagrange", "triangle", 2)
-P1 = FiniteElement("Lagrange", "triangle", 1)
-TH = P2 + P1
-
-(v, q) = TestFunctions(TH)
-(u, p) = TrialFunctions(TH)
-
-# alpha is scaling constant (L/h)^2 where L^2=\eta\U_0/\delta\rho g
-
-alpha = Constant("triangle")
-f = Function(P2)
-
-a = (alpha*dot(grad(v), grad(u)) - div(v)*p + q*div(u))*dx
-L = dot(v, f)*dx
diff -r d80d53c357e0 -r 204cfefa3897 MADDs-1/cpp/main.cpp
--- a/MADDs-1/cpp/main.cpp Sat Dec 05 08:23:03 2009 -0500
+++ b/MADDs-1/cpp/main.cpp Sat Dec 05 15:29:38 2009 -0500
@@ -39,13 +39,13 @@ using namespace dolfin;
};
// Function for Top velocity boundary condition for velocity
- class TopVelocity : public Function
+ class TopVelocity : public Expression
{
public:
- TopVelocity(double U0, double lambda) : U0(U0), lambda(lambda) {}
+ TopVelocity(double U0, double lambda) : Expression(2), U0(U0), lambda(lambda) {}
- void eval(double* values, const double* x) const
+ void eval(double* values, const std::vector<double>& x) const
{
//values[0] = U0*tanh(x[0]/lambda);
values[0] = U0*erf(x[0]/lambda);
@@ -59,13 +59,13 @@ using namespace dolfin;
};
// Function for Evaluating analytic corner flow
- class CornerFlowVelocity : public Function
+ class CornerFlowVelocity : public Expression
{
public:
- CornerFlowVelocity(double U0) : U0(U0) {}
+ CornerFlowVelocity(double U0) : Expression(2), U0(U0) {}
- void eval(double* values, const double* x) const
+ void eval(double* values, const std::vector<double>& x) const
{
double X = x[0];
double Z = 1. - x[1]; //transformed z coordinate
@@ -113,7 +113,7 @@ int main(int argc, char* argv[])
PetscOptionsGetReal("-lambda",&lambda,&flg);
//Create Function space and subspaces
- Stokes_P2P1FunctionSpace W(mesh);
+ Stokes_P2P1::FunctionSpace W(mesh);
SubSpace W0(W, 0); //velocity P2
SubSpace W1(W, 1); //pressure P1
@@ -150,13 +150,13 @@ int main(int argc, char* argv[])
// set forms
printf("Using P2_P1 mixed elements\n");
- Stokes_P2P1BilinearForm a(W, W);
+ Stokes_P2P1::BilinearForm a(W, W);
a.alpha = alpha;
- Stokes_P2P1LinearForm L(W);
+ Stokes_P2P1::LinearForm L(W);
L.f = f;
//preconditioner matrix
- Stokes_AQP2_P1BilinearForm pca(W,W);
+ Stokes_AQP2_P1::BilinearForm pca(W,W);
pca.alpha = alpha;
// Solve PDE and time
@@ -174,15 +174,13 @@ int main(int argc, char* argv[])
//calculate global divergence
- FacetNormal n;
- DivFunctional Div;
- Div.u = u;
- Div.n = n;
+ Div::Functional Div(mesh,u);
double divU = assemble(Div);
//get corner pressures
- const double x_corner[] = { 1., 1.}, x_center[]={ 0., 1.};
+ std::vector<double> x_corner (2,1.), x_center(2,1.);
+ x_center[0] = 0.;
double p_center, p_corner;
p.eval(&p_center,x_center);
p.eval(&p_corner,x_corner);
@@ -206,9 +204,9 @@ int main(int argc, char* argv[])
pfile << p.vector();
// Save solution in VTK format
- File ufile_pvd(output_dir+mesh_filename+"_velocity.pvd");
+ File ufile_pvd(output_dir+mesh_filename+"_velocity.pvd","compressed");
ufile_pvd << u;
- File pfile_pvd(output_dir+mesh_filename+"_pressure.pvd");
+ File pfile_pvd(output_dir+mesh_filename+"_pressure.pvd","compressed");
pfile_pvd << p;
//plot fields
More information about the CIG-COMMITS
mailing list