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:
{}

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 -523 ufl_formulation.c | tail +476")

void tabulate_tensor_integral_d6d1f05c892850147487cbd2aced82d778f6f49d(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 FE4_C0_D10_Q39d[1][1][1][3] = { { { { -1.0, 1.0, 0.0 } } } };
  static const double FE4_C1_D01_Q39d[1][1][1][3] = { { { { -1.0, 0.0, 1.0 } } } };
  static const double FE6_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 } } } };
  // Quadrature loop independent computations for quadrature rule 39d
  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] * FE4_C0_D10_Q39d[0][0][0][ic];
    J_c3 += coordinate_dofs[ic * 3 + 1] * FE4_C1_D01_Q39d[0][0][0][ic];
    J_c1 += coordinate_dofs[ic * 3] * FE4_C1_D01_Q39d[0][0][0][ic];
    J_c2 += coordinate_dofs[ic * 3 + 1] * FE4_C0_D10_Q39d[0][0][0][ic];
  }
  double sp_39d[4];
  sp_39d[0] = J_c0 * J_c3;
  sp_39d[1] = J_c1 * J_c2;
  sp_39d[2] = sp_39d[0] + -1 * sp_39d[1];
  sp_39d[3] = fabs(sp_39d[2]);
  for (int iq = 0; iq < 6; ++iq)
  {
    const double fw0 = sp_39d[3] * weights_39d[iq];
    double t0[6];
    for (int i = 0; i < 6; ++i)
      t0[i] = fw0 * FE6_C0_Q39d[0][0][iq][i];
    for (int i = 0; i < 6; ++i)
      for (int j = 0; j < 6; ++j)
        A[6 * i + j] += FE6_C0_Q39d[0][0][iq][j] * t0[i];
  }
}

0