Unexpected stiffness matrix for 3rd order Lagrange element

bjoni · April 28, 2025, 1:49pm

Hi,

when testing Poisson equation for 3rd order Lagrange triangle, I get an unexpected stiffness matrix (first and second order work fine). I called the generated ffcx file directly. Thanks!

The matrix:

+1.0111111111111080e+00 -7.2222222222222771e-02 -7.2222222222222604e-02 -8.3333333333336562e-02 -8.3333333333336618e-02 -7.4690399499995019e-01 +2.4690399499995494e-01 -7.4690399499995019e-01 +2.4690399499995497e-01 +3.0000000000000093e-01 
-7.2222222222222771e-02 +5.1111111111111340e-01 +1.1111111111110895e-02 -8.3333333333330761e-02 -8.3333333333333759e-02 -7.2723166354163801e-02 -1.0610166979168627e-02 +1.8479099562496087e-01 -6.8479099562496304e-01 +2.9999999999999816e-01 
-7.2222222222222604e-02 +1.1111111111110900e-02 +5.1111111111111318e-01 -8.3333333333333870e-02 -8.3333333333331178e-02 +1.8479099562496096e-01 -6.8479099562496282e-01 -7.2723166354164370e-02 -1.0610166979168702e-02 +2.9999999999999855e-01 
-8.3333333333336562e-02 -8.3333333333330789e-02 -8.3333333333333842e-02 +3.0555555555555616e+00 +2.7777777777777973e-01 +3.7613943506936098e-01 +1.3888888888888892e-01 +1.3888888888889717e-01 -1.4872505461804770e+00 -2.2500000000000044e+00 
-8.3333333333336618e-02 -8.3333333333333759e-02 -8.3333333333331205e-02 +2.7777777777777968e-01 +3.0555555555555598e+00 +1.3888888888889733e-01 -1.4872505461804759e+00 +3.7613943506935965e-01 +1.3888888888888901e-01 -2.2500000000000027e+00 
-7.4690399499995019e-01 -7.2723166354163815e-02 +1.8479099562496099e-01 +3.7613943506936098e-01 +1.3888888888889730e-01 +2.7777777777777763e+00 -6.9444444444445130e-01 -2.5569824035898137e-15 -4.0592529337857286e-16 -1.9635254915624254e+00 
+2.4690399499995494e-01 -1.0610166979168627e-02 -6.8479099562496282e-01 +1.3888888888888892e-01 -1.4872505461804759e+00 -6.9444444444445130e-01 +2.7777777777777888e+00 +1.4571677198205180e-16 -2.6784130469081902e-15 -2.8647450843757577e-01 
-7.4690399499995019e-01 +1.8479099562496087e-01 -7.2723166354164356e-02 +1.3888888888889717e-01 +3.7613943506935965e-01 -2.5535129566378600e-15 +1.1796119636642288e-16 +2.7777777777777746e+00 -6.9444444444445053e-01 -1.9635254915624241e+00 
+2.4690399499995497e-01 -6.8479099562496304e-01 -1.0610166979168702e-02 -1.4872505461804768e+00 +1.3888888888888901e-01 -4.3368086899420177e-16 -2.6714741530042829e-15 -6.9444444444445075e-01 +2.7777777777777883e+00 -2.8647450843757633e-01 
+3.0000000000000093e-01 +2.9999999999999821e-01 +2.9999999999999855e-01 -2.2500000000000044e+00 -2.2500000000000027e+00 -1.9635254915624254e+00 -2.8647450843757594e-01 -1.9635254915624243e+00 -2.8647450843757644e-01 +8.1000000000000103e+00

while I would expect (from comsol):

+8.5000000000000497e-01 -8.7500000000000702e-02 -8.7500000000000702e-02 -7.5000000000000996e-02 -7.5000000000000996e-02 -6.3750000000000295e-01 +3.7500000000000200e-01 -6.3750000000000295e-01 +3.7500000000000200e-01 +0.0000000000000000e+00 
-8.7500000000000702e-02 +4.2500000000000199e-01 +0.0000000000000000e+00 +3.7500000000000498e-02 +3.7500000000000699e-02 -3.7500000000000297e-02 -3.7500000000000699e-02 +3.3750000000000202e-01 -6.7500000000000404e-01 +0.0000000000000000e+00 
-8.7500000000000702e-02 +0.0000000000000000e+00 +4.2500000000000199e-01 +3.7500000000000699e-02 +3.7500000000000498e-02 +3.3750000000000202e-01 -6.7500000000000404e-01 -3.7500000000000297e-02 -3.7500000000000699e-02 +0.0000000000000000e+00 
-7.5000000000000996e-02 +3.7500000000000498e-02 +3.7500000000000699e-02 +3.3750000000000102e+00 -6.7500000000000604e-01 +3.3750000000000002e-01 +3.3750000000000302e-01 +3.3750000000000502e-01 -1.6875000000000100e+00 -2.0249999999999999e+00 
-7.5000000000000996e-02 +3.7500000000000699e-02 +3.7500000000000498e-02 -6.7500000000000604e-01 +3.3750000000000102e+00 +3.3750000000000502e-01 -1.6875000000000100e+00 +3.3750000000000002e-01 +3.3750000000000302e-01 -2.0249999999999999e+00 
-6.3750000000000295e-01 -3.7500000000000297e-02 +3.3750000000000202e-01 +3.3750000000000002e-01 +3.3750000000000502e-01 +3.3750000000000102e+00 -1.6875000000000100e+00 +0.0000000000000000e+00 +0.0000000000000000e+00 -2.0249999999999999e+00 
+3.7500000000000200e-01 -3.7500000000000699e-02 -6.7500000000000404e-01 +3.3750000000000302e-01 -1.6875000000000100e+00 -1.6875000000000100e+00 +3.3750000000000102e+00 +0.0000000000000000e+00 +0.0000000000000000e+00 +0.0000000000000000e+00 
-6.3750000000000295e-01 +3.3750000000000202e-01 -3.7500000000000297e-02 +3.3750000000000502e-01 +3.3750000000000002e-01 +0.0000000000000000e+00 +0.0000000000000000e+00 +3.3750000000000102e+00 -1.6875000000000100e+00 -2.0249999999999999e+00 
+3.7500000000000200e-01 -6.7500000000000404e-01 -3.7500000000000699e-02 -1.6875000000000100e+00 +3.3750000000000302e-01 +0.0000000000000000e+00 +0.0000000000000000e+00 -1.6875000000000100e+00 +3.3750000000000102e+00 +0.0000000000000000e+00 
+0.0000000000000000e+00 +0.0000000000000000e+00 +0.0000000000000000e+00 -2.0249999999999999e+00 -2.0249999999999999e+00 -2.0249999999999999e+00 +0.0000000000000000e+00 -2.0249999999999999e+00 +0.0000000000000000e+00 +8.0999999999999996e+00

C-File:


// This code conforms with the UFC specification version 2018.2.0.dev0
// and was automatically generated by FFCx version 0.9.0.
//
// This code was generated with the following options:
//
//  {'epsilon': 1e-14,
//   'output_directory': '.',
//   'profile': False,
//   'scalar_type': 'float64',
//   'sum_factorization': False,
//   'table_atol': 1e-09,
//   'table_rtol': 1e-06,
//   'ufl_file': ['Poisson.ufl'],
//   'verbosity': 30,
//   'visualise': False}

#include <math.h>
#include <stdalign.h>
#include <stdlib.h>
#include <string.h>
#include <ufcx.h>
#define restrict



// Code for integral integral_8e97e9655e7fd2790f3bb75e995f66871ca8f8d1

void tabulate_tensor_integral_8e97e9655e7fd2790f3bb75e995f66871ca8f8d1(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 FE1_C0_D10_Q39d[1][1][1][3] = {{{{-1.0, 1.0, 0.0}}}};
static const double FE1_C1_D01_Q39d[1][1][1][3] = {{{{-1.0, 0.0, 1.0}}}};
static const double FE3_C0_D01_Q39d[1][1][6][10] = {{{{-0.2100309081141055, 0.0, 0.2100309081140965, 4.93577169254329, -0.8515338276409833, -0.09376059857271754, 0.09376059857272634, 0.851533827641004, -4.935771692543287, 0.0},
  {-0.2764485129815143, -0.06641760486740975, 2.773706027588994, -0.0623791259376049, 1.516524266497604, 0.789154701703399, -3.419247426045698, 0.222311203444784, 0.316071802017506, -1.793275331420061},
  {-2.77370602758898, 0.06641760486740975, 0.2764485129815072, -0.3160718020175055, -0.2223112034447813, 3.419247426045671, -0.7891547017033785, -1.516524266497604, 0.06237912593760606, 1.793275331420057},
  {0.4764340609062518, 0.0, -0.4764340609062518, -0.2034003103928528, 0.7439154012332023, -2.223434825677669, 2.223434825677669, -0.7439154012332014, 0.2034003103928523, 0.0},
  {0.6270957396409664, 0.1506616787347146, 0.244925435109706, 0.5576400002159992, -0.5878227266568932, -0.1862754010172019, -0.384422416264041, -3.356551230639105, -1.133116404961437, 4.067865325837293},
  {-0.2449254351097055, -0.1506616787347144, -0.6270957396409665, 1.133116404961438, 3.356551230639106, 0.3844224162640411, 0.1862754010172023, 0.5878227266568933, -0.5576400002159997, -4.067865325837293}}}};
static const double FE3_C0_D10_Q39d[1][1][6][10] = {{{{-0.2764485129815141, 2.773706027588995, -0.06641760486740997, 1.516524266497606, -0.06237912593760572, 0.2223112034447852, 0.3160718020175064, 0.789154701703398, -3.419247426045698, -1.79327533142006},
  {-0.2100309081141054, 0.2100309081140966, 0.0, -0.8515338276409836, 4.935771692543288, 0.8515338276410047, -4.935771692543287, -0.09376059857271883, 0.09376059857272631, 0.0},
  {-2.773706027588981, 0.2764485129815072, 0.06641760486740919, -0.2223112034447822, -0.3160718020175055, -1.516524266497604, 0.06237912593760604, 3.419247426045669, -0.7891547017033786, 1.793275331420058},
  {0.6270957396409658, 0.2449254351097055, 0.1506616787347145, -0.5878227266568946, 0.5576400002159995, -3.356551230639105, -1.133116404961437, -0.1862754010172011, -0.3844224162640414, 4.067865325837293},
  {0.4764340609062521, -0.4764340609062521, 0.0, 0.7439154012332021, -0.2034003103928527, -0.7439154012332021, 0.2034003103928519, -2.223434825677668, 2.223434825677669, 0.0},
  {-0.2449254351097053, -0.6270957396409661, -0.1506616787347146, 3.356551230639107, 1.133116404961437, 0.5878227266568937, -0.5576400002160002, 0.3844224162640408, 0.1862754010172021, -4.067865325837293}}}};
// ------------------------ 
// Section: Jacobian
// Inputs: coordinate_dofs, FE1_C1_D01_Q39d, FE1_C0_D10_Q39d
// Outputs: J_c1, J_c0, J_c2, J_c3
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] * FE1_C0_D10_Q39d[0][0][0][ic];
    J_c3 += coordinate_dofs[(ic) * 3 + 1] * FE1_C1_D01_Q39d[0][0][0][ic];
    J_c1 += coordinate_dofs[(ic) * 3] * FE1_C1_D01_Q39d[0][0][0][ic];
    J_c2 += coordinate_dofs[(ic) * 3 + 1] * FE1_C0_D10_Q39d[0][0][0][ic];
  }
}
// ------------------------ 
double sp_39d_0 = J_c0 * J_c3;
double sp_39d_1 = J_c1 * J_c2;
double sp_39d_2 = -sp_39d_1;
double sp_39d_3 = sp_39d_0 + sp_39d_2;
double sp_39d_4 = J_c0 / sp_39d_3;
double sp_39d_5 = -J_c1;
double sp_39d_6 = sp_39d_5 / sp_39d_3;
double sp_39d_7 = sp_39d_4 * sp_39d_4;
double sp_39d_8 = sp_39d_4 * sp_39d_6;
double sp_39d_9 = sp_39d_6 * sp_39d_6;
double sp_39d_10 = J_c3 / sp_39d_3;
double sp_39d_11 = -J_c2;
double sp_39d_12 = sp_39d_11 / sp_39d_3;
double sp_39d_13 = sp_39d_12 * sp_39d_12;
double sp_39d_14 = sp_39d_10 * sp_39d_12;
double sp_39d_15 = sp_39d_10 * sp_39d_10;
double sp_39d_16 = sp_39d_7 + sp_39d_13;
double sp_39d_17 = sp_39d_8 + sp_39d_14;
double sp_39d_18 = sp_39d_15 + sp_39d_9;
double sp_39d_19 = fabs(sp_39d_3);
double sp_39d_20 = sp_39d_16 * sp_39d_19;
double sp_39d_21 = sp_39d_17 * sp_39d_19;
double sp_39d_22 = sp_39d_18 * sp_39d_19;
for (int iq = 0; iq < 6; ++iq)
{
  // ------------------------ 
  // Section: Intermediates
  // Inputs: 
  // Outputs: fw0, fw1, fw2
  double fw0 = 0;
  double fw1 = 0;
  double fw2 = 0;
  {
    fw0 = sp_39d_22 * weights_39d[iq];
    fw1 = sp_39d_21 * weights_39d[iq];
    fw2 = sp_39d_20 * weights_39d[iq];
  }
  // ------------------------ 
  // ------------------------ 
  // Section: Tensor Computation
  // Inputs: fw1, FE3_C0_D01_Q39d, fw2, fw0, FE3_C0_D10_Q39d
  // Outputs: A
  {
    double temp_0[10] = {0};
    for (int j = 0; j < 10; ++j)
    {
      temp_0[j] = fw0 * FE3_C0_D10_Q39d[0][0][iq][j];
    }
    double temp_1[10] = {0};
    for (int j = 0; j < 10; ++j)
    {
      temp_1[j] = fw1 * FE3_C0_D01_Q39d[0][0][iq][j];
    }
    double temp_2[10] = {0};
    for (int j = 0; j < 10; ++j)
    {
      temp_2[j] = fw1 * FE3_C0_D10_Q39d[0][0][iq][j];
    }
    double temp_3[10] = {0};
    for (int j = 0; j < 10; ++j)
    {
      temp_3[j] = fw2 * FE3_C0_D01_Q39d[0][0][iq][j];
    }
    for (int j = 0; j < 10; ++j)
    {
      for (int i = 0; i < 10; ++i)
      {
        A[10 * (i) + (j)] += FE3_C0_D10_Q39d[0][0][iq][i] * temp_0[j];
        A[10 * (i) + (j)] += FE3_C0_D10_Q39d[0][0][iq][i] * temp_1[j];
        A[10 * (i) + (j)] += FE3_C0_D01_Q39d[0][0][iq][i] * temp_2[j];
        A[10 * (i) + (j)] += FE3_C0_D01_Q39d[0][0][iq][i] * temp_3[j];
      }
    }
  }
  // ------------------------ 
}

}



ufcx_integral integral_8e97e9655e7fd2790f3bb75e995f66871ca8f8d1 =
{
  .enabled_coefficients = NULL,
  .tabulate_tensor_float32 = NULL,
  .tabulate_tensor_float64 = tabulate_tensor_integral_8e97e9655e7fd2790f3bb75e995f66871ca8f8d1,
  .tabulate_tensor_complex64 = NULL,
  .tabulate_tensor_complex128 = NULL,
  .needs_facet_permutations = false,
  .coordinate_element_hash = UINT64_C(0),
};

// End of code for integral integral_8e97e9655e7fd2790f3bb75e995f66871ca8f8d1

// Code for form form_ebd69fd14e49da296310b21a0ead2d2b25d3c854


uint64_t finite_element_hashes_form_ebd69fd14e49da296310b21a0ead2d2b25d3c854[2] = {UINT64_C(0), UINT64_C(0)};
int form_integral_offsets_form_ebd69fd14e49da296310b21a0ead2d2b25d3c854[4] = {0, 1, 1, 1};
static ufcx_integral* form_integrals_form_ebd69fd14e49da296310b21a0ead2d2b25d3c854[1] = {&integral_8e97e9655e7fd2790f3bb75e995f66871ca8f8d1};
int form_integral_ids_form_ebd69fd14e49da296310b21a0ead2d2b25d3c854[1] = {-1};




ufcx_form form_ebd69fd14e49da296310b21a0ead2d2b25d3c854 =
{

  .signature = "fe750be0aa4c02f6f7840ede9806da7a1e1f17c6333b09502dd096b9942e1d2b6e2edf8506cf5aed083d73c57a4f997d2a87472eafb0a63796c9de353061bc1e",
  .rank = 2,
  .num_coefficients = 0,
  .num_constants = 0,
  .original_coefficient_positions = NULL,

  .coefficient_name_map = NULL,
  .constant_name_map = NULL,

  .finite_element_hashes = finite_element_hashes_form_ebd69fd14e49da296310b21a0ead2d2b25d3c854,

  .form_integrals = form_integrals_form_ebd69fd14e49da296310b21a0ead2d2b25d3c854,
  .form_integral_ids = form_integral_ids_form_ebd69fd14e49da296310b21a0ead2d2b25d3c854,
  .form_integral_offsets = form_integral_offsets_form_ebd69fd14e49da296310b21a0ead2d2b25d3c854
};

// Alias name
ufcx_form* form_Poisson_a = &form_ebd69fd14e49da296310b21a0ead2d2b25d3c854;

// End of code for form form_ebd69fd14e49da296310b21a0ead2d2b25d3c854

UFL:

from dolfinx import mesh, fem
from mpi4py import MPI
from ufl import TestFunction, TrialFunction, grad, dx, inner
domain = mesh.create_unit_square(MPI.COMM_WORLD, 1, 1, mesh.CellType.triangle)
V = fem.functionspace(domain, ("Lagrange", 3,(1,)))
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx(domain=domain)

dokken · May 2, 2025, 12:07pm

Please note that when you use 3rd order basis functions, one has to consider what Lagrange-variants to use for the basis function.
The motivation for this is covered in: Variants of Lagrange elements — DOLFINx 0.9.0 documentation

However, for completeness, lets consider:

from dolfinx import mesh, fem
from mpi4py import MPI
from ufl import TestFunction, TrialFunction, grad, dx, inner
from basix.ufl import element
from basix import LagrangeVariant

domain = mesh.create_unit_square(MPI.COMM_WORLD, 1, 1, mesh.CellType.triangle)
el = element(
    "Lagrange",
    domain.basix_cell(),
    3,
    shape=(1,),
    lagrange_variant=LagrangeVariant.equispaced,
)
V = fem.functionspace(domain, el)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx(domain=domain)
A = fem.assemble_matrix(fem.form(a))
print(A.to_dense())

This gives the same result as comsol.
However, it has less favourable properties for higher order, thus, if DOLFINx defaults to gll_warped as Lagrange variant.

bjoni · May 5, 2025, 9:42am

That makes sense, many thanks!

Topic		Replies	Views
Solution Error with Lagrange 1 elements dolfinx	2	39	August 16, 2024
Problem with stiffness matrices mesh	2	388	April 7, 2020
Ordering of the stiffness matrix? Legacy DOLFIN	7	48	October 31, 2024
Element stiffness matrix variational formulation	2	2072	April 17, 2024
FEniCS : How to interpolate data at vertices of (3D) cells? Linear Algebra	6	1626	October 19, 2020

Unexpected stiffness matrix for 3rd order Lagrange element

Related topics