codesnippets.py

finite element library for mathematic majored research

const double C0 = d00*(element_coordinates[0][0] - element_coordinates[2][0] - element_coordinates[3][0]) \\
                + d10*(element_coordinates[0][1] - element_coordinates[2][1] - element_coordinates[3][1]) \\
                + d20*(element_coordinates[0][2] - element_coordinates[2][2] - element_coordinates[3][2]);

const double C1 = d01*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[3][0]) \\
                + d11*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[3][1]) \\
                + d21*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[3][2]);

const double C2 = d02*(element_coordinates[0][0] - element_coordinates[1][0] - element_coordinates[2][0]) \\
                + d12*(element_coordinates[0][1] - element_coordinates[1][1] - element_coordinates[2][1]) \\
                + d22*(element_coordinates[0][2] - element_coordinates[1][2] - element_coordinates[2][2]);

// Get coordinates and map to the UFC reference element
double x = (C0 + d00*coordinates[0] + d10*coordinates[1] + d20*coordinates[2]) / detJ;
double y = (C1 + d01*coordinates[0] + d11*coordinates[1] + d21*coordinates[2]) / detJ;
double z = (C2 + d02*coordinates[0] + d12*coordinates[1] + d22*coordinates[2]) / detJ;"""


# Inverse affine map from physical cell to the FIAT reference cell in 2D
map_coordinates_FIAT_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;

// Compute constants
const double C0 = element_coordinates[1][0] + element_coordinates[2][0];
const double C1 = element_coordinates[1][1] + element_coordinates[2][1];

// Get coordinates and map to the reference (FIAT) element
double x = (J_01*C1 - J_11*C0 + 2.0*J_11*coordinates[0] - 2.0*J_01*coordinates[1]) / detJ;
double y = (J_10*C0 - J_00*C1 - 2.0*J_10*coordinates[0] + 2.0*J_00*coordinates[1]) / detJ;"""

# Inverse affine map from physical cell to (FIAT) reference cell in 3D
map_coordinates_FIAT_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
const double C0 = element_coordinates[3][0] + element_coordinates[2][0] \\
                + element_coordinates[1][0] - element_coordinates[0][0];
const double C1 = element_coordinates[3][1] + element_coordinates[2][1] \\
                + element_coordinates[1][1] - element_coordinates[0][1];
const double C2 = element_coordinates[3][2] + element_coordinates[2][2] \\
                + element_coordinates[1][2] - element_coordinates[0][2];

// Get coordinates and map to the reference (FIAT) element
double x = coordinates[0];
double y = coordinates[1];
double z = coordinates[2];

x = (2.0*d00*x + 2.0*d10*y + 2.0*d20*z - d00*C0 - d10*C1 - d20*C2) / detJ;
y = (2.0*d01*x + 2.0*d11*y + 2.0*d21*z - d01*C0 - d11*C1 - d21*C2) / detJ;
z = (2.0*d02*x + 2.0*d12*y + 2.0*d22*z - d02*C0 - d12*C1 - d22*C2) / detJ;"""

# Snippet to generate combinations of derivatives of order n
combinations_snippet = """\
// Declare pointer to two dimensional array that holds combinations of derivatives and initialise
unsigned int **%(combinations)s = new unsigned int *[%(num_derivatives)s];
    
for (unsigned int j = 0; j < %(num_derivatives)s; j++)
{
  %(combinations)s[j] = new unsigned int [%(n)s];
  for (unsigned int k = 0; k < %(n)s; k++)
    %(combinations)s[j][k] = 0;
}
    
// Generate combinations of derivatives
for (unsigned int row = 1; row < %(num_derivatives)s; row++)
{
  for (unsigned int num = 0; num < row; num++)
  {
    for (unsigned int col = %(n)s-1; col+1 > 0; col--)
    {
      if (%(combinations)s[row][col] + 1 > %(shape-1)s)
        %(combinations)s[row][col] = 0;
      else
      {
        %(combinations)s[row][col] += 1;
        break;
      }
    }
  }
}"""

# Snippet to transform of derivatives of order n
transform2D_snippet = """\
// Compute inverse of Jacobian
const double %(Jinv)s[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)s = new double *[%(num_derivatives)s];
    
for (unsigned int j = 0; j < %(num_derivatives)s; j++)
{
  %(transform)s[j] = new double [%(num_derivatives)s];
  for (unsigned int k = 0; k < %(num_derivatives)s; k++)
    %(transform)s[j][k] = 1;
}

// Construct transformation matrix
for (unsigned int row = 0; row < %(num_derivatives)s; row++)
{
  for (unsigned int col = 0; col < %(num_derivatives)s; col++)
  {
    for (unsigned int k = 0; k < %(n)s; k++)
      %(transform)s[row][col] *= %(Jinv)s[%(combinations)s[col][k]][%(combinations)s[row][k]];
  }
}"""

transform3D_snippet = """\
// Compute inverse of Jacobian
const double %(Jinv)s[3][3] =\
{{d00 / detJ, d10 / detJ, d20 / detJ},\
 {d01 / detJ, d11 / detJ, d21 / detJ},\
 {d02 / detJ, d12 / detJ, d22 / detJ}};

// Declare transformation matrix
// Declare pointer to two dimensional array and initialise
double **%(transform)s = new double *[%(num_derivatives)s];
    
for (unsigned int j = 0; j < %(num_derivatives)s; j++)
{
  %(transform)s[j] = new double [%(num_derivatives)s];
  for (unsigned int k = 0; k < %(num_derivatives)s; k++)
    %(transform)s[j][k] = 1;
}

// Construct transformation matrix
for (unsigned int row = 0; row < %(num_derivatives)s; row++)
{
  for (unsigned int col = 0; col < %(num_derivatives)s; col++)
  {
    for (unsigned int k = 0; k < %(n)s; k++)
      %(transform)s[row][col] *= %(Jinv)s[%(combinations)s[col][k]][%(combinations)s[row][k]];
  }
}"""

# Snippet to transform of derivatives of order n
transform2D_FIAT_snippet = """\
// Compute inverse of Jacobian, components are scaled because of difference in FFC/FIAT reference elements
const double %(Jinv)s[2][2] =  {{2*J_11 / detJ, -2*J_01 / detJ}, {-2*J_10 / detJ, 2*J_00 / detJ}};

// Declare transformation matrix
// Declare pointer to two dimensional array and initialise
double **%(transform)s = new double *[%(num_derivatives)s];
    
for (unsigned int j = 0; j < %(num_derivatives)s; j++)
{
  %(transform)s[j] = new double [%(num_derivatives)s];
  for (unsigned int k = 0; k < %(num_derivatives)s; k++)
    %(transform)s[j][k] = 1;
}

// Construct transformation matrix
for (unsigned int row = 0; row < %(num_derivatives)s; row++)
{
  for (unsigned int col = 0; col < %(num_derivatives)s; col++)
  {
    for (unsigned int k = 0; k < %(n)s; k++)
      %(transform)s[row][col] *= %(Jinv)s[%(combinations)s[col][k]][%(combinations)s[row][k]];
  }
}"""

transform3D_FIAT_snippet = """\
// Compute inverse of Jacobian, components are scaled because of difference in FFC/FIAT reference elements
const double %(Jinv)s[3][3] =\
{{2*d00 / detJ, 2*d10 / detJ, 2*d20 / detJ},\
 {2*d01 / detJ, 2*d11 / detJ, 2*d21 / detJ},\
 {2*d02 / detJ, 2*d12 / detJ, 2*d22 / detJ}};

// Declare transformation matrix
// Declare pointer to two dimensional array and initialise
double **%(transform)s = new double *[%(num_derivatives)s];
    
for (unsigned int j = 0; j < %(num_derivatives)s; j++)
{
  %(transform)s[j] = new double [%(num_derivatives)s];
  for (unsigned int k = 0; k < %(num_derivatives)s; k++)
    %(transform)s[j][k] = 1;
}

// Construct transformation matrix
for (unsigned int row = 0; row < %(num_derivatives)s; row++)
{
  for (unsigned int col = 0; col < %(num_derivatives)s; col++)
  {
    for (unsigned int k = 0; k < %(n)s; k++)
      %(transform)s[row][col] *= %(Jinv)s[%(combinations)s[col][k]][%(combinations)s[row][k]];
  }
}"""

# Code snippet for computing the inverse of the Jacobian and determinant in 2D
inverse_jacobian_2D = """\
// 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 inverse of the Jacobian and determinant in 3D
inverse_jacobian_3D = """\
// Compute sub determinants
const double d00 = J%(restriction)s_11*J%(restriction)s_22 - J%(restriction)s_12*J%(restriction)s_21;
const double d01 = J%(restriction)s_12*J%(restriction)s_20 - J%(restriction)s_10*J%(restriction)s_22;
const double d02 = J%(restriction)s_10*J%(restriction)s_21 - J%(restriction)s_11*J%(restriction)s_20;

const double d10 = J%(restriction)s_02*J%(restriction)s_21 - J%(restriction)s_01*J%(restriction)s_22;
const double d11 = J%(restriction)s_00*J%(restriction)s_22 - J%(restriction)s_02*J%(restriction)s_20;
const double d12 = J%(restriction)s_01*J%(restriction)s_20 - J%(restriction)s_00*J%(restriction)s_21;

const double d20 = J%(restriction)s_01*J%(restriction)s_12 - J%(restriction)s_02*J%(restriction)s_11;
const double d21 = J%(restriction)s_02*J%(restriction)s_10 - J%(restriction)s_00*J%(restriction)s_12;
const double d22 = 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*d00 + J%(restriction)s_10*d10 + J%(restriction)s_20*d20;
  
// Compute inverse of Jacobian
const double Jinv%(restriction)s_00 = d00 / detJ%(restriction)s;
const double Jinv%(restriction)s_01 = d10 / detJ%(restriction)s;
const double Jinv%(restriction)s_02 = d20 / detJ%(restriction)s;
const double Jinv%(restriction)s_10 = d01 / detJ%(restriction)s;
const double Jinv%(restriction)s_11 = d11 / detJ%(restriction)s;
const double Jinv%(restriction)s_12 = d21 / detJ%(restriction)s;
const double Jinv%(restriction)s_20 = d02 / detJ%(restriction)s;
const double Jinv%(restriction)s_21 = d12 / detJ%(restriction)s;
const double Jinv%(restriction)s_22 = d22 / detJ%(restriction)s;"""