Form with domain markers creates duplicated kernels

[IgorBaratta]

MWE:

coord_element = VectorElement("Lagrange", interval, 1)
mesh = Mesh(coord_element)

element = FiniteElement("Lagrange", interval, 1)
V = FunctionSpace(mesh, element)

u = TrialFunction(V)
v = TestFunction(V)

a = inner(u, v) * dx((1,2))

Produces 2 identitical kernels:

// Code for integral integral_8799b8ef6f5269ca9a672cec402eefb7d5222475

void tabulate_tensor_integral_8799b8ef6f5269ca9a672cec402eefb7d5222475(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_4a8[2] = { 0.5, 0.5 };
    // Precomputed values of basis functions and precomputations
    // FE* dimensions: [permutation][entities][points][dofs]
    static const double FE0_C0_Q4a8[1][1][2][2] =
        { { { { 0.7886751345948129, 0.2113248654051871 },
              { 0.2113248654051871, 0.7886751345948129 } } } };
    static const double FE1_C0_D1_Q4a8[1][1][1][2] = { { { { -1.0, 1.0 } } } };
    // Quadrature loop independent computations for quadrature rule 4a8
    const double J_c0 = coordinate_dofs[0] * FE1_C0_D1_Q4a8[0][0][0][0] + coordinate_dofs[3] * FE1_C0_D1_Q4a8[0][0][0][1];
    double sp_4a8[1];
    sp_4a8[0] = fabs(J_c0);
    for (int iq = 0; iq < 2; ++iq)
    {
        const double fw0 = sp_4a8[0] * weights_4a8[iq];
        double t0[2];
        for (int i = 0; i < 2; ++i)
            t0[i] = fw0 * FE0_C0_Q4a8[0][0][iq][i];
        for (int i = 0; i < 2; ++i)
            for (int j = 0; j < 2; ++j)
                A[2 * i + j] += FE0_C0_Q4a8[0][0][iq][j] * t0[i];
    }
}


void tabulate_tensor_integral_b1ad34778d813d876edcdf491953bf64df0b3699(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_4a8[2] = { 0.5, 0.5 };
    // Precomputed values of basis functions and precomputations
    // FE* dimensions: [permutation][entities][points][dofs]
    static const double FE0_C0_Q4a8[1][1][2][2] =
        { { { { 0.7886751345948129, 0.2113248654051871 },
              { 0.2113248654051871, 0.7886751345948129 } } } };
    static const double FE1_C0_D1_Q4a8[1][1][1][2] = { { { { -1.0, 1.0 } } } };
    // Quadrature loop independent computations for quadrature rule 4a8
    const double J_c0 = coordinate_dofs[0] * FE1_C0_D1_Q4a8[0][0][0][0] + coordinate_dofs[3] * FE1_C0_D1_Q4a8[0][0][0][1];
    double sp_4a8[1];
    sp_4a8[0] = fabs(J_c0);
    for (int iq = 0; iq < 2; ++iq)
    {
        const double fw0 = sp_4a8[0] * weights_4a8[iq];
        double t0[2];
        for (int i = 0; i < 2; ++i)
            t0[i] = fw0 * FE0_C0_Q4a8[0][0][iq][i];
        for (int i = 0; i < 2; ++i)
            for (int j = 0; j < 2; ++j)
                A[2 * i + j] += FE0_C0_Q4a8[0][0][iq][j] * t0[i];
    }
}
 

[jorgensd]

The problem seems to be that ufl splits integrals by domain marker.
I.e.

from ufl import *
coord_element = VectorElement("Lagrange", interval, 1)
mesh = Mesh(coord_element)

element = FiniteElement("Lagrange", interval, 1)
V = FunctionSpace(mesh, element)

u = TrialFunction(V)
v = TestFunction(V)

a = inner(u, v) * dx((1,2))
print(a.integrals())

returns two integrals.

[jorgensd]

This is mentioned in FEniCS/ufl#70
The big question is, why do we split integrals by domain id?

[michalhabera]

The problem seems to be that ufl splits integrals by domain marker. I.e.

from ufl import *
coord_element = VectorElement("Lagrange", interval, 1)
mesh = Mesh(coord_element)

element = FiniteElement("Lagrange", interval, 1)
V = FunctionSpace(mesh, element)

u = TrialFunction(V)
v = TestFunction(V)

a = inner(u, v) * dx((1,2))
print(a.integrals())

returns two integrals.

Yep, It is expected as we have one kernel for each integral_data. It comes from compute_form_data this way. Not sure whether UFL's group_form_integrals is effective here.

[jorgensd]

@michalhabera
The issue is already in the initialization of the form: https://github.com/FEniCS/ufl/blob/main/ufl/form.py#L36-L37
As the integrals are split there, I think it is really hard to reverse engineer them back together.

[wence-]

at the cfd stage you can merge integral data where the only difference is in the subdomain id. we can definitely handle that integration in firedrake (it might need minor adaptations)