Code generation

All the code in the previous section is Python code symbolically describing the variational form of the projection.

This is why we in FEniCS rely on code generation. One can interpret the variational form written in UFL as a directed acyclic graph of operations, where we for each simple operation can implement it as C code.

FFCx

We use FFCx to generate the code used for finite element assembly.

The previous section of the tutorial can be found in the file introduction.py. We can use the FFCx main module to generate C code for all the objects in this file

import ffcx.main
from pathlib import Path
import os
cwd = Path.cwd()
infile = cwd / "ufl_formulation.py"

ffcx.main.main(["-o", str(cwd), "--visualise", str(infile)])

IntegralData over domain(cell, ('otherwise',)) with integrals:
{ weight * |(J[0, 0] * J[1, 1] + -1 * J[0, 1] * J[1, 0])| * (reference_value(w_0)) / (reference_value(w_1)) * (reference_value(v_0)) } * dx(<Mesh #0>[('otherwise',)], {'estimated_polynomial_degree': 3})
and metadata:
{metadata}

0

This computes the computational graph of the bi-linear and linear form

Bilinear graph

Linear graph

With this graph, we also get compiled c code:

os.system("ls ufl_formulation.*")

ufl_formulation.c
ufl_formulation.h
ufl_formulation.py

0

We can look at the assembly code for the local matrix. We start by inspecting the signature of the tabulate_tensor function, that computes the local element matrix

os.system("head -344 ufl_formulation.c | tail +272")

void tabulate_tensor_integral_444ae3de485a81d987c5e8c21247536bf812e211(double* restrict A,
                                    const double* restrict w,
                                    const double* restrict c,
                                    const double* restrict coordinate_dofs,
                                    const int* restrict entity_local_index,
                                    const uint8_t* restrict quadrature_permutation)
{
// Quadrature rules
static const double weights_39d[6] = {0.054975871827661, 0.054975871827661, 0.054975871827661, 0.1116907948390055, 0.1116907948390055, 0.1116907948390055};
// Precomputed values of basis functions and precomputations
// FE* dimensions: [permutation][entities][points][dofs]
static const double FE0_C0_Q39d[1][1][6][6] = {{{{-0.07480380774819603, 0.5176323419876736, -0.07480380774819671, 0.2992152309927871, 0.03354481152314834, 0.2992152309927839},
  {-0.07480380774819613, -0.0748038077481966, 0.5176323419876735, 0.2992152309927871, 0.2992152309927838, 0.03354481152314828},
  {0.5176323419876713, -0.0748038077481967, -0.07480380774819674, 0.03354481152314866, 0.2992152309927869, 0.2992152309927868},
  {-0.04820837781551195, -0.08473049309397784, -0.04820837781551192, 0.1928335112620479, 0.7954802262009061, 0.1928335112620478},
  {-0.04820837781551193, -0.048208377815512, -0.08473049309397786, 0.1928335112620479, 0.192833511262048, 0.7954802262009062},
  {-0.08473049309397794, -0.04820837781551188, -0.04820837781551195, 0.7954802262009061, 0.1928335112620479, 0.1928335112620479}}}};
static const double FE2_C0_D10_Q39d[1][1][1][3] = {{{{-1.0, 1.0, 0.0}}}};
static const double FE2_C1_D01_Q39d[1][1][1][3] = {{{{-1.0, 0.0, 1.0}}}};
// ------------------------ 
// Section: Jacobian
// Inputs: coordinate_dofs, FE2_C0_D10_Q39d, FE2_C1_D01_Q39d
// Outputs: J_c2, J_c0, J_c3, J_c1
double J_c0 = 0.0;
double J_c3 = 0.0;
double J_c1 = 0.0;
double J_c2 = 0.0;
{
  for (int ic = 0; ic < 3; ++ic)
  {
    J_c0 += coordinate_dofs[(ic) * 3] * FE2_C0_D10_Q39d[0][0][0][ic];
    J_c3 += coordinate_dofs[(ic) * 3 + 1] * FE2_C1_D01_Q39d[0][0][0][ic];
    J_c1 += coordinate_dofs[(ic) * 3] * FE2_C1_D01_Q39d[0][0][0][ic];
    J_c2 += coordinate_dofs[(ic) * 3 + 1] * FE2_C0_D10_Q39d[0][0][0][ic];
  }
}
// ------------------------ 
double sp_39d_0 = J_c0 * J_c3;
double sp_39d_1 = J_c1 * J_c2;
double sp_39d_2 = -sp_39d_1;
double sp_39d_3 = sp_39d_0 + sp_39d_2;
double sp_39d_4 = fabs(sp_39d_3);
for (int iq = 0; iq < 6; ++iq)
{
  // ------------------------ 
  // Section: Intermediates
  // Inputs: 
  // Outputs: fw0
  double fw0 = 0;
  {
    fw0 = sp_39d_4 * weights_39d[iq];
  }
  // ------------------------ 
  // ------------------------ 
  // Section: Tensor Computation
  // Inputs: FE0_C0_Q39d, fw0
  // Outputs: A
  {
    double temp_0[6] = {0};
    for (int j = 0; j < 6; ++j)
    {
      temp_0[j] = fw0 * FE0_C0_Q39d[0][0][iq][j];
    }
    for (int j = 0; j < 6; ++j)
    {
      for (int i = 0; i < 6; ++i)
      {
        A[6 * (i) + (j)] += FE0_C0_Q39d[0][0][iq][i] * temp_0[j];
      }
    }
  }
  // ------------------------ 
}

0