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)