mscroggs / symfem

A symbolic finite element definition library
MIT License
53 stars 9 forks source link

Fix rHCT element #207

Closed avieira closed 1 year ago

avieira commented 1 year ago

Describe the bug The rHCT element seems to be wrong once integrated. I noticed it when I tried to build the stiffness matrix for the Poisson equation using these elements. It is clear when one tries to integrate it using symfem or sympy.

To Reproduce Minimal code to reproduce the behavior:

import symfem
from symfem.symbols import x
from sympy import Symbol, integrate, Rational

#Compute integral using symfem
element = symfem.create_element("triangle", "rHCT", 3)
ref = element.reference
f1 = element.get_basis_function(1).directional_derivative((1,0))
f2 = element.get_basis_function(6).directional_derivative((1,0))
integrand = f1*f2
#f1*f2 is 0 everywhere, except on triangle ((0, 1), (0, 0), (1/3, 1/3))
#where it equals (42*x**2 - 24*x*y - 6*x)*(18*x**2 + 6*x*y - 10*x - y + 1)
print(integrand.integral(ref, x)) #=0

#Compute it using sympy
z = Symbol("x")
y = Symbol("y")
expr = (42*z**2 - 24*z*y - 6*z)*(18*z**2 + 6*z*y - 10*z - y + 1) #Copied from what is given by symfem
int_expr = integrate(expr, (y, z, 1-2*z)) #parametrization of triangle ((0, 1), (0, 0), (1/3, 1/3)) : x<y<1-2*x
print(integrate(int_expr, (z,0, Rational(1,3))))#=1/108

Expected behavior The integral should be correct in this example

Screenshots Not applicable.

Additional context No other context needed (I guess?)

mscroggs commented 1 year ago

Looks like this TODO is the problem. I'll work on a fix

https://github.com/mscroggs/symfem/blob/43fbadd94cbb4698058997fd8cccb1f76ee42bd7/symfem/piecewise_functions.py#L435-L436

mscroggs commented 1 year ago

This should be fixed by #208