As far as Im aware these features are in DOLFINx.
For instance, you can adapt Figure 3 to work with ffcx:
from ufl import Cell, VectorElement, FunctionSpace, Mesh, Coefficient
# Define triangle cell embedded in Rˆ3
cell = Cell("triangle", 3)
# Define Lagrange vector element
coordinate_element = VectorElement("Lagrange", cell, 1)
mesh = Mesh(coordinate_element)
V = FunctionSpace(mesh, coordinate_element)
# Arguments defined over V will have 3 components:
u = Coefficient(V)
print(u[0], u[1], u[2])
python3 -m ffcx manifold.py
w_0[0] w_0[1] w_0[2]
or Figure 4:
from ufl import Cell, VectorElement, FunctionSpace, Mesh, Coefficient, FiniteElement, TrialFunction, TestFunction, inner, grad, dx
# Define triangle cell embedded in Rˆ3
cell = Cell("triangle", 3)
# Define Lagrange vector element
coordinate_element = VectorElement("Lagrange", cell, 1)
mesh = Mesh(coordinate_element)
element = FiniteElement("Lagrange", cell, 1)
V = FunctionSpace(mesh, element)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx
which yields the kernel in manifold.c
void tabulate_tensor_integral_6a56ba35a29fca58aae87b107061fd2de6b29fd1(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_083[1] = { 0.5 };
// Precomputed values of basis functions and precomputations
// FE* dimensions: [permutation][entities][points][dofs]
static const double FE10_C0_D010_Q083[1][1][1][3] = { { { { -1.0, 0.0, 1.0 } } } };
static const double FE6_C0_D100_Q083[1][1][1][3] = { { { { -1.0, 1.0, 0.0 } } } };
// Quadrature loop independent computations for quadrature rule 083
double J_c1 = 0.0;
double J_c0 = 0.0;
double J_c2 = 0.0;
double J_c4 = 0.0;
double J_c3 = 0.0;
double J_c5 = 0.0;
for (int ic = 0; ic < 3; ++ic)
{
J_c1 += coordinate_dofs[ic * 3] * FE10_C0_D010_Q083[0][0][0][ic];
J_c0 += coordinate_dofs[ic * 3] * FE6_C0_D100_Q083[0][0][0][ic];
J_c2 += coordinate_dofs[ic * 3 + 1] * FE6_C0_D100_Q083[0][0][0][ic];
J_c4 += coordinate_dofs[ic * 3 + 2] * FE6_C0_D100_Q083[0][0][0][ic];
J_c3 += coordinate_dofs[ic * 3 + 1] * FE10_C0_D010_Q083[0][0][0][ic];
J_c5 += coordinate_dofs[ic * 3 + 2] * FE10_C0_D010_Q083[0][0][0][ic];
}
double sp_083[74];
sp_083[0] = J_c0 * J_c0;
sp_083[1] = J_c2 * J_c2;
sp_083[2] = sp_083[0] + sp_083[1];
sp_083[3] = J_c4 * J_c4;
sp_083[4] = sp_083[2] + sp_083[3];
sp_083[5] = J_c1 * J_c1;
sp_083[6] = J_c3 * J_c3;
sp_083[7] = sp_083[5] + sp_083[6];
sp_083[8] = J_c5 * J_c5;
sp_083[9] = sp_083[7] + sp_083[8];
sp_083[10] = sp_083[4] * sp_083[9];
sp_083[11] = J_c0 * J_c1;
sp_083[12] = J_c2 * J_c3;
sp_083[13] = sp_083[11] + sp_083[12];
sp_083[14] = J_c4 * J_c5;
sp_083[15] = sp_083[13] + sp_083[14];
sp_083[16] = sp_083[15] * sp_083[15];
sp_083[17] = sp_083[10] + -1 * sp_083[16];
sp_083[18] = sp_083[4] / sp_083[17];
sp_083[19] = J_c1 * sp_083[18];
sp_083[20] = (-1 * sp_083[15]) / sp_083[17];
sp_083[21] = J_c0 * sp_083[20];
sp_083[22] = sp_083[19] + sp_083[21];
sp_083[23] = sp_083[9] / sp_083[17];
sp_083[24] = J_c0 * sp_083[23];
sp_083[25] = J_c1 * sp_083[20];
sp_083[26] = sp_083[24] + sp_083[25];
sp_083[27] = sp_083[22] * sp_083[22];
sp_083[28] = sp_083[22] * sp_083[26];
sp_083[29] = sp_083[26] * sp_083[26];
sp_083[30] = J_c3 * sp_083[18];
sp_083[31] = J_c2 * sp_083[20];
sp_083[32] = sp_083[30] + sp_083[31];
sp_083[33] = J_c2 * sp_083[23];
sp_083[34] = J_c3 * sp_083[20];
sp_083[35] = sp_083[33] + sp_083[34];
sp_083[36] = sp_083[32] * sp_083[32];
sp_083[37] = sp_083[32] * sp_083[35];
sp_083[38] = sp_083[35] * sp_083[35];
sp_083[39] = sp_083[27] + sp_083[36];
sp_083[40] = sp_083[28] + sp_083[37];
sp_083[41] = sp_083[29] + sp_083[38];
sp_083[42] = J_c5 * sp_083[18];
sp_083[43] = J_c4 * sp_083[20];
sp_083[44] = sp_083[42] + sp_083[43];
sp_083[45] = J_c4 * sp_083[23];
sp_083[46] = J_c5 * sp_083[20];
sp_083[47] = sp_083[45] + sp_083[46];
sp_083[48] = sp_083[44] * sp_083[44];
sp_083[49] = sp_083[44] * sp_083[47];
sp_083[50] = sp_083[47] * sp_083[47];
sp_083[51] = sp_083[39] + sp_083[48];
sp_083[52] = sp_083[40] + sp_083[49];
sp_083[53] = sp_083[41] + sp_083[50];
sp_083[54] = J_c2 * J_c5;
sp_083[55] = J_c3 * J_c4;
sp_083[56] = sp_083[54] + -1 * sp_083[55];
sp_083[57] = sp_083[56] * sp_083[56];
sp_083[58] = J_c1 * J_c4;
sp_083[59] = J_c0 * J_c5;
sp_083[60] = sp_083[58] + -1 * sp_083[59];
sp_083[61] = sp_083[60] * sp_083[60];
sp_083[62] = sp_083[57] + sp_083[61];
sp_083[63] = J_c0 * J_c3;
sp_083[64] = J_c1 * J_c2;
sp_083[65] = sp_083[63] + -1 * sp_083[64];
sp_083[66] = sp_083[65] * sp_083[65];
sp_083[67] = sp_083[62] + sp_083[66];
sp_083[68] = sqrt(sp_083[67]);
sp_083[69] = 1 * sp_083[68];
sp_083[70] = fabs(sp_083[69]);
sp_083[71] = sp_083[51] * sp_083[70];
sp_083[72] = sp_083[52] * sp_083[70];
sp_083[73] = sp_083[53] * sp_083[70];
for (int iq = 0; iq < 1; ++iq)
{
const double fw0 = sp_083[73] * weights_083[iq];
const double fw1 = sp_083[72] * weights_083[iq];
const double fw2 = sp_083[71] * weights_083[iq];
double t0[3];
double t1[3];
for (int i = 0; i < 3; ++i)
{
t0[i] = fw0 * FE6_C0_D100_Q083[0][0][0][i] + fw1 * FE10_C0_D010_Q083[0][0][0][i];
t1[i] = fw1 * FE6_C0_D100_Q083[0][0][0][i] + fw2 * FE10_C0_D010_Q083[0][0][0][i];
}
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
A[3 * i + j] += FE6_C0_D100_Q083[0][0][0][j] * t0[i] + FE10_C0_D010_Q083[0][0][0][j] * t1[i];
}
}
which is different from using
cell = Cell("triangle", 2)
yielding
void tabulate_tensor_integral_27b24fdf56699cf003a569e3622e9eaddc6a64cf(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_083[1] = { 0.5 };
// Precomputed values of basis functions and precomputations
// FE* dimensions: [permutation][entities][points][dofs]
static const double FE4_C0_D10_Q083[1][1][1][3] = { { { { -1.0, 1.0, 0.0 } } } };
static const double FE8_C0_D01_Q083[1][1][1][3] = { { { { -1.0, 0.0, 1.0 } } } };
// Quadrature loop independent computations for quadrature rule 083
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_Q083[0][0][0][ic];
J_c3 += coordinate_dofs[ic * 3 + 1] * FE8_C0_D01_Q083[0][0][0][ic];
J_c1 += coordinate_dofs[ic * 3] * FE8_C0_D01_Q083[0][0][0][ic];
J_c2 += coordinate_dofs[ic * 3 + 1] * FE4_C0_D10_Q083[0][0][0][ic];
}
double sp_083[20];
sp_083[0] = J_c0 * J_c3;
sp_083[1] = J_c1 * J_c2;
sp_083[2] = sp_083[0] + -1 * sp_083[1];
sp_083[3] = J_c0 / sp_083[2];
sp_083[4] = (-1 * J_c1) / sp_083[2];
sp_083[5] = sp_083[3] * sp_083[3];
sp_083[6] = sp_083[3] * sp_083[4];
sp_083[7] = sp_083[4] * sp_083[4];
sp_083[8] = J_c3 / sp_083[2];
sp_083[9] = (-1 * J_c2) / sp_083[2];
sp_083[10] = sp_083[9] * sp_083[9];
sp_083[11] = sp_083[8] * sp_083[9];
sp_083[12] = sp_083[8] * sp_083[8];
sp_083[13] = sp_083[5] + sp_083[10];
sp_083[14] = sp_083[6] + sp_083[11];
sp_083[15] = sp_083[12] + sp_083[7];
sp_083[16] = fabs(sp_083[2]);
sp_083[17] = sp_083[13] * sp_083[16];
sp_083[18] = sp_083[14] * sp_083[16];
sp_083[19] = sp_083[15] * sp_083[16];
for (int iq = 0; iq < 1; ++iq)
{
const double fw0 = sp_083[19] * weights_083[iq];
const double fw1 = sp_083[18] * weights_083[iq];
const double fw2 = sp_083[17] * weights_083[iq];
double t0[3];
double t1[3];
for (int i = 0; i < 3; ++i)
{
t0[i] = fw0 * FE4_C0_D10_Q083[0][0][0][i] + fw1 * FE8_C0_D01_Q083[0][0][0][i];
t1[i] = fw1 * FE4_C0_D10_Q083[0][0][0][i] + fw2 * FE8_C0_D01_Q083[0][0][0][i];
}
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
A[3 * i + j] += FE4_C0_D10_Q083[0][0][0][j] * t0[i] + FE8_C0_D01_Q083[0][0][0][j] * t1[i];
}
}