How dolfinx calculate integral

dolfinx
1

Hello, I want to ask a question about how dolfinx calculate integral.
Say if I have a function f and I create the following form, where v is the test function:

a = f * f * v * dx

when it assemble the term, will the basis function of both f be integralled? or it just regard “f*f” as a single function and use the quadrature node of " f * f ".

Thanks!

2

As long as f is a dolfinx.fem.Function (not a TrialFunction where this wouldn’t make sense), this will generate an evaluation of f at the quadrature points for both instances.

This can be checked by looking at the generated code of the following MWE

import basix.ufl
import ufl


c_el = basix.ufl.element("Lagrange", "interval", 1, shape=(1,))
mesh = ufl.Mesh(c_el)

el = basix.ufl.element("Lagrange", "interval", 1)
V = ufl.FunctionSpace(mesh, el)
u = ufl.Coefficient(V)
v = ufl.TestFunction(V)

F = u * u * v * ufl.dx

executed with python3 -m ffcx name_of_script.py
which yields a name_of_script.c file with the following kernel:

void tabulate_tensor_integral_260ff49c39f0d831f5c0318c5ce971311757d924(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}}}};
// ------------------------ 
// Section: Jacobian
// Inputs: FE1_C0_D1_Q4a8, coordinate_dofs
// Outputs: J_c0
double J_c0 = 0.0;
{
  for (int ic = 0; ic < 2; ++ic)
  {
    J_c0 += coordinate_dofs[(ic) * 3] * FE1_C0_D1_Q4a8[0][0][0][ic];
  }
}
// ------------------------ 
double sp_4a8_0 = fabs(J_c0);
for (int iq = 0; iq < 2; ++iq)
{
  // ------------------------ 
  // Section: Coefficient
  // Inputs: w, FE0_C0_Q4a8
  // Outputs: w0
  double w0 = 0.0;
  {
    for (int ic = 0; ic < 2; ++ic)
    {
      w0 += w[ic] * FE0_C0_Q4a8[0][0][iq][ic];
    }
  }
  // ------------------------ 
  // ------------------------ 
  // Section: Intermediates
  // Inputs: w0
  // Outputs: fw0
  double fw0 = 0;
  {
    double sv_4a8_0 = w0 * w0;
    double sv_4a8_1 = sv_4a8_0 * sp_4a8_0;
    fw0 = sv_4a8_1 * weights_4a8[iq];
  }
  // ------------------------ 
  // ------------------------ 
  // Section: Tensor Computation
  // Inputs: fw0, FE0_C0_Q4a8
  // Outputs: A
  {
    for (int i = 0; i < 2; ++i)
    {
      A[(i)] += fw0 * FE0_C0_Q4a8[0][0][iq][i];
    }
  }
  // ------------------------ 
}

where you see the evaluation of the f at quadrature points inside the Coefficient , and in section “Intermediates” you see the computation of f(x_q)**2*abs(detJ(x_q))*weight(x_q)

1 Like