codesnippets.py

finite element library for mathematic majored research

"Code snippets for code generation"

__author__ = "Anders Logg ([e-mail])"
__date__ = "2007-02-28 -- 2007-04-18"
__copyright__ = "Copyright (C) 2007 Anders Logg"
__license__  = "GNU GPL version 3 or any later version"

# Modified by Kristian Oelgaard 2007
# Modified by Marie Rognes 2007

# Code snippet for computing the Jacobian, its inverse and determinant in 2D
jacobian_2D = """\
// Extract vertex coordinates
const double * const * x%(restriction)s = c%(restriction)s.coordinates;

// Compute Jacobian of affine map from reference cell
const double J%(restriction)s_00 = x%(restriction)s[1][0] - x%(restriction)s[0][0];
const double J%(restriction)s_01 = x%(restriction)s[2][0] - x%(restriction)s[0][0];
const double J%(restriction)s_10 = x%(restriction)s[1][1] - x%(restriction)s[0][1];
const double J%(restriction)s_11 = x%(restriction)s[2][1] - x%(restriction)s[0][1];
  
// Compute determinant of Jacobian
double detJ%(restriction)s = J%(restriction)s_00*J%(restriction)s_11 - J%(restriction)s_01*J%(restriction)s_10;
  
// Compute inverse of Jacobian
const double Jinv%(restriction)s_00 =  J%(restriction)s_11 / detJ%(restriction)s;
const double Jinv%(restriction)s_01 = -J%(restriction)s_01 / detJ%(restriction)s;
const double Jinv%(restriction)s_10 = -J%(restriction)s_10 / detJ%(restriction)s;
const double Jinv%(restriction)s_11 =  J%(restriction)s_00 / detJ%(restriction)s;"""

# Code snippet for computing the Jacobian, its inverse and determinant in 3D
jacobian_3D = """\
// Extract vertex coordinates
const double * const * x%(restriction)s = c%(restriction)s.coordinates;

// Compute Jacobian of affine map from reference cell
const double J%(restriction)s_00 = x%(restriction)s[1][0] - x%(restriction)s[0][0];
const double J%(restriction)s_01 = x%(restriction)s[2][0] - x%(restriction)s[0][0];
const double J%(restriction)s_02 = x%(restriction)s[3][0] - x%(restriction)s[0][0];
const double J%(restriction)s_10 = x%(restriction)s[1][1] - x%(restriction)s[0][1];
const double J%(restriction)s_11 = x%(restriction)s[2][1] - x%(restriction)s[0][1];
const double J%(restriction)s_12 = x%(restriction)s[3][1] - x%(restriction)s[0][1];
const double J%(restriction)s_20 = x%(restriction)s[1][2] - x%(restriction)s[0][2];
const double J%(restriction)s_21 = x%(restriction)s[2][2] - x%(restriction)s[0][2];
const double J%(restriction)s_22 = x%(restriction)s[3][2] - x%(restriction)s[0][2];
  
// Compute sub determinants
const double d%(restriction)s_00 = J%(restriction)s_11*J%(restriction)s_22 - J%(restriction)s_12*J%(restriction)s_21;
const double d%(restriction)s_01 = J%(restriction)s_12*J%(restriction)s_20 - J%(restriction)s_10*J%(restriction)s_22;
const double d%(restriction)s_02 = J%(restriction)s_10*J%(restriction)s_21 - J%(restriction)s_11*J%(restriction)s_20;

const double d%(restriction)s_10 = J%(restriction)s_02*J%(restriction)s_21 - J%(restriction)s_01*J%(restriction)s_22;
const double d%(restriction)s_11 = J%(restriction)s_00*J%(restriction)s_22 - J%(restriction)s_02*J%(restriction)s_20;
const double d%(restriction)s_12 = J%(restriction)s_01*J%(restriction)s_20 - J%(restriction)s_00*J%(restriction)s_21;

const double d%(restriction)s_20 = J%(restriction)s_01*J%(restriction)s_12 - J%(restriction)s_02*J%(restriction)s_11;
const double d%(restriction)s_21 = J%(restriction)s_02*J%(restriction)s_10 - J%(restriction)s_00*J%(restriction)s_12;
const double d%(restriction)s_22 = J%(restriction)s_00*J%(restriction)s_11 - J%(restriction)s_01*J%(restriction)s_10;
  
// Compute determinant of Jacobian
double detJ%(restriction)s = J%(restriction)s_00*d%(restriction)s_00 + J%(restriction)s_10*d%(restriction)s_10 + J%(restriction)s_20*d%(restriction)s_20;
  
// Compute inverse of Jacobian
const double Jinv%(restriction)s_00 = d%(restriction)s_00 / detJ%(restriction)s;
const double Jinv%(restriction)s_01 = d%(restriction)s_10 / detJ%(restriction)s;
const double Jinv%(restriction)s_02 = d%(restriction)s_20 / detJ%(restriction)s;
const double Jinv%(restriction)s_10 = d%(restriction)s_01 / detJ%(restriction)s;
const double Jinv%(restriction)s_11 = d%(restriction)s_11 / detJ%(restriction)s;
const double Jinv%(restriction)s_12 = d%(restriction)s_21 / detJ%(restriction)s;
const double Jinv%(restriction)s_20 = d%(restriction)s_02 / detJ%(restriction)s;
const double Jinv%(restriction)s_21 = d%(restriction)s_12 / detJ%(restriction)s;
const double Jinv%(restriction)s_22 = d%(restriction)s_22 / detJ%(restriction)s;"""

# Code snippet for computing the scale factor (determinant)
scale_factor = """\
// Set scale factor
const double det = std::abs(detJ);
"""

# Code snippet for computing the determinant of the facet mapping in 2D
facet_determinant_2D = """\
// Vertices on edges
static unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};

// Get vertices
const unsigned int v0 = edge_vertices[%(facet)s][0];
const unsigned int v1 = edge_vertices[%(facet)s][1];

// Compute scale factor (length of edge scaled by length of reference interval)
const double dx0 = x%(restriction)s[v1][0] - x%(restriction)s[v0][0];
const double dx1 = x%(restriction)s[v1][1] - x%(restriction)s[v0][1];
const double det = std::sqrt(dx0*dx0 + dx1*dx1);
"""

# Code snippet for computing the determinant of the facet mapping in 3D
facet_determinant_3D = """\
// Vertices on faces
static unsigned int face_vertices[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};

// Get vertices
const unsigned int v0 = face_vertices[%(facet)s][0];
const unsigned int v1 = face_vertices[%(facet)s][1];
const unsigned int v2 = face_vertices[%(facet)s][2];

// Compute scale factor (area of face scaled by area of reference triangle)
const double a0 = (x%(restriction)s[v0][1]*x%(restriction)s[v1][2] + x%(restriction)s[v0][2]*x%(restriction)s[v2][1] + x%(restriction)s[v1][1]*x%(restriction)s[v2][2])
              - (x%(restriction)s[v2][1]*x%(restriction)s[v1][2] + x%(restriction)s[v2][2]*x%(restriction)s[v0][1] + x%(restriction)s[v1][1]*x%(restriction)s[v0][2]);
const double a1 = (x%(restriction)s[v0][2]*x%(restriction)s[v1][0] + x%(restriction)s[v0][0]*x%(restriction)s[v2][2] + x%(restriction)s[v1][2]*x%(restriction)s[v2][0])
              - (x%(restriction)s[v2][2]*x%(restriction)s[v1][0] + x%(restriction)s[v2][0]*x%(restriction)s[v0][2] + x%(restriction)s[v1][2]*x%(restriction)s[v0][0]);
const double a2 = (x%(restriction)s[v0][0]*x%(restriction)s[v1][1] + x%(restriction)s[v0][1]*x%(restriction)s[v2][0] + x%(restriction)s[v1][0]*x%(restriction)s[v2][1])
              - (x%(restriction)s[v2][0]*x%(restriction)s[v1][1] + x%(restriction)s[v2][1]*x%(restriction)s[v0][0] + x%(restriction)s[v1][0]*x%(restriction)s[v0][1]);
const double det = std::sqrt(a0*a0 + a1*a1 + a2*a2);
"""

# Code snippet for evaluate_dof in 2D
evaluate_dof_2D = """\
double values[%d];
double coordinates[2];

// Nodal coordinates on reference cell
static double X[%d][2] = %s;

// Components for each dof
static unsigned int components[%d] = %s;

// Extract vertex coordinates
const double * const * x = c.coordinates;

// Evaluate basis functions for affine mapping
const double w0 = 1.0 - X[i][0] - X[i][1];
const double w1 = X[i][0];
const double w2 = X[i][1];

// Compute affine mapping x = F(X)
coordinates[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0];
coordinates[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1];

// Evaluate function at coordinates
f.evaluate(values, coordinates, c);

// Pick component for evaluation
return values[components[i]];"""

# Code snippet for evaluate_dof in 3D
evaluate_dof_3D = """\
double values[%d];
double coordinates[3];

// Nodal coordinates on reference cell
static double X[%d][3] = %s;

// Components for each dof
static unsigned int components[%d] = %s;

// Extract vertex coordinates
const double * const * x = c.coordinates;

// Evaluate basis functions for affine mapping
const double w0 = 1.0 - X[i][0] - X[i][1] - X[i][2];
const double w1 = X[i][0];
const double w2 = X[i][1];
const double w3 = X[i][2];

// Compute affine mapping x = F(X)
coordinates[0] = w0*x[0][0] + w1*x[1][0] + w2*x[2][0] + w3*x[3][0];
coordinates[1] = w0*x[0][1] + w1*x[1][1] + w2*x[2][1] + w3*x[3][1];
coordinates[2] = w0*x[0][2] + w1*x[1][2] + w2*x[2][2] + w3*x[3][2];

// Evaluate function at coordinates
f.evaluate(values, coordinates, c);

// Pick component for evaluation
return values[components[i]];"""

# Code snippet for eta basis on triangles (reproduced from FIAT reference.py)
eta_triangle_snippet = """\
if (std::abs(y - 1.0) < %s)
  x = -1.0;
else
  x = 2.0 *x/(1.0 - y) - 1.0;
y = 2.0*y - 1.0;"""

# Code snippet for eta basis on tetrahedra (reproduced from FIAT reference.py)
eta_tetrahedron_snippet = """\
if (std::abs(y + z - 1.0) < %s)
  x = 1.0;
else
  x = -2.0 * x/(y + z - 1.0) - 1.0;
if (std::abs(z - 1.0) < %s)
  y = -1.0;
else
  y = 2.0 * y/(1.0 - z) - 1.0;
z = 2.0 * z - 1.0;"""

# Map basis function to local basisfunction, used by 'evaluate_basis.py'
evaluate_basis_dof_map = """\
unsigned int element = 0;
unsigned int tmp = 0;
for (unsigned int j = 0; j < %d; j++)
{
  if (tmp +  dofs_per_element[j] > i)
  {
    i -= tmp;
    element = element_types[j];
    break;
  }
  else
    tmp += dofs_per_element[j];


}"""


# Inverse affine map from physical cell to UFC reference cell in 2D
map_coordinates_2D = """\
// 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;

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

# Inverse affine map from physical cell to the UFC reference cell in 3D
map_coordinates_3D = """\
// 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_02 = element_coordinates[3][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];
const double J_12 = element_coordinates[3][1] - element_coordinates[0][1];
const double J_20 = element_coordinates[1][2] - element_coordinates[0][2];
const double J_21 = element_coordinates[2][2] - element_coordinates[0][2];
const double J_22 = element_coordinates[3][2] - element_coordinates[0][2];
  
// Compute sub determinants
const double d00 = J_11*J_22 - J_12*J_21;
const double d01 = J_12*J_20 - J_10*J_22;
const double d02 = J_10*J_21 - J_11*J_20;

const double d10 = J_02*J_21 - J_01*J_22;
const double d11 = J_00*J_22 - J_02*J_20;
const double d12 = J_01*J_20 - J_00*J_21;

const double d20 = J_01*J_12 - J_02*J_11;
const double d21 = J_02*J_10 - J_00*J_12;
const double d22 = J_00*J_11 - J_01*J_10;
  
// Compute determinant of Jacobian
double detJ = J_00*d00 + J_10*d10 + J_20*d20;

// Compute constants