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