Approximating function by 2D FEA

I am trying to approximate function by 2D finite element method in Python. For this test case I am using 1 triangle. I am comparing my result with FEniCS, and unfortunately they do not match. My stiffness matrix A and right vector b are completely different from FEniCS matrices. Tried to calculate elements of matrices in two coordinates: 1. Global coordinate (physical). 2. Local coordinate (natural). In my opinion, the problem is in integration of product of basis functions. I do not know how to properly integrate product of basis functions in global coordinate as well as in local coordinate, what are the limits of integration. Might as well I am missing something in integral expression or multiplying basis functions incorrectly. Can you please help me with the realization of this method, thanks.

For using local coordinates just uncomment Method 2, A matrix, b vector block of code and comment another one

FEniCS code:

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

def approximate(f, V):
	
	u = TrialFunction(V)
	v = TestFunction(V)
	a = u*v*dx
	L = f*v*dx
	
	#matrix A
	A = Matrix()
	assemble(a, A)
	print('\nFEniCS matrix A\n')
	print(A.array())

	#vector b
	b = Vector()
	assemble(L, b)
	print('\nFEniCS vector b\n')
	print(b.get_local())

	#solution vector
	u = Function(V)
	solve(a == L, u)
	print('\nFEniCS solution vector\n')
	print(u.vector().get_local())
	
	return u

def calculate():
	
	#function which need to approximate
	f = Expression('15*x[0]*x[0]', degree=1)

	#creating 1 triangle element
	mesh = Mesh()
	editor = MeshEditor()

	editor.open(mesh, 'triangle', 2, 2)
	editor.init_vertices(3)
	editor.init_cells(1)

	editor.add_vertex(0, np.array([0.25, 0.0]))
	editor.add_vertex(1, np.array([1.5, 0.0]))
	editor.add_vertex(2, np.array([0.25, 1.85]))
	editor.add_cell(0, np.array([0, 1, 2], dtype=np.uintp))

	editor.close()
	plot(mesh)
	
	#cell and its coordinates
	mesh.init()
	for i in range(mesh.num_cells()):
		cell = Cell(mesh, i)
		info('Cell number %i' % i)
		vis = cell.entities(0)
		for vi in vis:
			vertex = Vertex(mesh, vi)
			xx = vertex.point().x()
			yy = vertex.point().y()
			info('Vertex %i point(%g, %g)' % (vi, xx, yy))

	#approximating, using linear basis function
	V = FunctionSpace(mesh, 'P', 1)
	print('\nBasis functions:', V.dim())
	u = approximate(f, V)
	#plot(u, title='u')

calculate()
plt.show()

Code:

import sympy as sym
import numpy as np

x, y, n, e = sym.symbols('x y n e')

#function which need to approximate
f = 15*x*x

#x,y coordinates of 3 vertexes of triangle
i = [0.25, 0.0]															#vertex 0	
j = [1.5, 0.0]															#vertex 1
k = [0.25, 1.85]														#vertex 2

#calculating determinant of Jacobian matrix
detJ = i[0]*(j[1]-k[1]) + j[0]*(k[1]-i[1]) + k[0]*(i[1]-j[1])
print('Determinant of Jacobian matrix:', detJ)


#---------------------------- Method 1. Using basis functions in global coordinate ---------------------------------# 

#phi0 basis function in global coordinate
Ai0 = j[0]*k[1] - k[0]*j[1]
Bi0 = j[1] - k[1]
Si0 = k[0] - j[0]
phi_global0 = (1/detJ)*(Ai0+Bi0*x+Si0*y)
print('phi0 basis function:', phi_global0)

#phi1 basis function in global coordinate 
Ai1 = k[0]*i[1] - i[0]*k[1]
Bi1 = k[1] - i[1]
Si1 = i[0] - k[0]
phi_global1 = (1/detJ)*(Ai1+Bi1*x+Si1*y)
print('phi1 basis function:', phi_global1)

#phi2 basis function in global coordinate
Ai2 = i[0]*j[1] - j[0]*i[1]
Bi2 = i[1] - j[1]
Si2 = j[0] - i[0]
phi_global2 = (1/detJ)*(Ai2+Bi2*x+Si2*y)
print('phi2 basis function:', phi_global2)

#phi0+phi1+phi2 has to be equal 1, which proofs correctness of basis functions
print('phi0+phi1+phi2:', phi_global0+phi_global1+phi_global2)

#-------------------------------------------------------------------------------------------------------------------# 


'''
#---------------------------- Method 2. Using basis functions in local coordinate ----------------------------------# 

#phi0 basis function in local coordinate
phi_local0 = 1-n-e
print('phi0 basis function:', phi_local0)

#phi1 basis function in local coordinate
phi_local1 = n
print('phi0 basis function:', phi_local1)

#phi2 basis function in local coordinate
phi_local2 = e
print('phi0 basis function:', phi_local2)

#phi0+phi1+phi2 has to be equal 1, which proofs correctness of basis functions
print('phi0+phi1+phi2:', phi_local0+phi_local1+phi_local2)

#--------------------------------------------------------------------------------------------------------------------# 
'''


#A matrix
A = sym.zeros(3, 3)


#A matrix, integrating in global (physical) coordinates
A[0,0] = sym.integrate(phi_global0*phi_global0, (x, 0.25, 1.5), (y, 0, 1.85))
A[0,1] = sym.integrate(phi_global0*phi_global1, (x, 0.25, 1.5), (y, 0, 1.85))
A[0,2] = sym.integrate(phi_global0*phi_global2, (x, 0.25, 1.5), (y, 0, 1.85))

A[1,0] = sym.integrate(phi_global1*phi_global0, (x, 0.25, 1.5), (y, 0, 1.85))
A[1,1] = sym.integrate(phi_global1*phi_global1, (x, 0.25, 1.5), (y, 0, 1.85))
A[1,2] = sym.integrate(phi_global1*phi_global2, (x, 0.25, 1.5), (y, 0, 1.85))

A[2,0] = sym.integrate(phi_global2*phi_global0, (x, 0.25, 1.5), (y, 0, 1.85))
A[2,1] = sym.integrate(phi_global2*phi_global1, (x, 0.25, 1.5), (y, 0, 1.85))
A[2,2] = sym.integrate(phi_global2*phi_global2, (x, 0.25, 1.5), (y, 0, 1.85))
print('\nMatrix A:', A)

'''
#A matrix, integrating in local (natural) coordinates
A[0,0] = sym.integrate(phi_local0*phi_local0*detJ, (n, 0, 1), (e, 0, 1))
A[0,1] = sym.integrate(phi_local0*phi_local1*detJ, (n, 0, 1), (e, 0, 1))
A[0,2] = sym.integrate(phi_local0*phi_local2*detJ, (n, 0, 1), (e, 0, 1))

A[1,0] = sym.integrate(phi_local1*phi_local0*detJ, (n, 0, 1), (e, 0, 1))
A[1,1] = sym.integrate(phi_local1*phi_local1*detJ, (n, 0, 1), (e, 0, 1))
A[1,2] = sym.integrate(phi_local1*phi_local2*detJ, (n, 0, 1), (e, 0, 1))

A[2,0] = sym.integrate(phi_local2*phi_local0*detJ, (n, 0, 1), (e, 0, 1))
A[2,1] = sym.integrate(phi_local2*phi_local1*detJ, (n, 0, 1), (e, 0, 1))
A[2,2] = sym.integrate(phi_local2*phi_local2*detJ, (n, 0, 1), (e, 0, 1))
print('\nMatrix A:', A)
'''

#b vector
b = sym.zeros(3, 1)

#b vector, integrating in global (physical) coordinates
b[0] = sym.integrate(phi_global0 * f, (x, 0.25, 1.5), (y, 0, 1.85))
b[1] = sym.integrate(phi_global1 * f, (x, 0.25, 1.5), (y, 0, 1.85))
b[2] = sym.integrate(phi_global2 * f, (x, 0.25, 1.5), (y, 0, 1.85))
print('\nVector b:', b)

'''
#b vector, integrating in local (natural) coordinates
b[0] = sym.integrate(phi_local0 * f * detJ, (n, 0, 1), (e, 0, 1))
b[1] = sym.integrate(phi_local1 * f * detJ, (n, 0, 1), (e, 0, 1))
b[2] = sym.integrate(phi_local2 * f * detJ, (n, 0, 1), (e, 0, 1))
print('\nVector b:', b)
'''

#solution vector
c = A.LUsolve(b)
print('\nSolution vector c:', c)

As far as I can tell by just looking at your code (not executing it), you have an issue when integrating.
You are integrating over a quadrilateral, not a triangle.
This can be illustrated with the following MWE:

import sympy
x, y = sympy.symbols("x y")
f = x**3*y + 2 * y
sympy_int = sympy.integrate(f, (x,0,1), (y, 0, 1))

intf = 1/8*x**4*y**2 + y**2
exact_int = intf.subs(x, 1).subs(y, 1) print(exact_int, sympy_int)
print(exact_int-sympy_int)

yielding:

1.12500000000000 9/8
0

To integrate over a triangle, you would have to change the limits of one of the variables to depend on the other:

# Triangle integration (of triangle with vertices (0,0), (0, 1), (1, 0))
sympy_int = sympy.integrate(f, (y,0,1-x), (x, 0, 1))
print(sympy_int)

yielding:

41/120

Now, it is working, thanks a lot. The stiffness matrix A is identical with FEniCS result.

Corrected code, reference coordinate:

#A matrix, integrating in local (natural) coordinates
A[0,0] = sym.integrate(phi_local0*phi_local0*detJ, (n, 0, 1-e), (e, 0, 1))
A[0,1] = sym.integrate(phi_local0*phi_local1*detJ, (n, 0, 1-e), (e, 0, 1))
A[0,2] = sym.integrate(phi_local0*phi_local2*detJ, (n, 0, 1-e), (e, 0, 1))

A[1,0] = sym.integrate(phi_local1*phi_local0*detJ, (n, 0, 1-e), (e, 0, 1))
A[1,1] = sym.integrate(phi_local1*phi_local1*detJ, (n, 0, 1-e), (e, 0, 1))
A[1,2] = sym.integrate(phi_local1*phi_local2*detJ, (n, 0, 1-e), (e, 0, 1))

A[2,0] = sym.integrate(phi_local2*phi_local0*detJ, (n, 0, 1-e), (e, 0, 1))
A[2,1] = sym.integrate(phi_local2*phi_local1*detJ, (n, 0, 1-e), (e, 0, 1))
A[2,2] = sym.integrate(phi_local2*phi_local2*detJ, (n, 0, 1-e), (e, 0, 1))
print('\nMatrix A:', A)

Now, I have a problem with load vector b. It is still different from FEniCS load vector. For instance, if I use function without X and Y variables:

#function which need to approximate
f = 15

#b vector, integrating in local (natural) coordinates
b[0] = sym.integrate(phi_local1 * f * detJ, (n, 0, 1-e), (e, 0, 1))
b[1] = sym.integrate(phi_local1 * f * detJ, (n, 0, 1-e), (e, 0, 1))
b[2] = sym.integrate(phi_local2 * f * detJ, (n, 0, 1-e), (e, 0, 1))
print('\nVector b:', b)

the load vector b will be identical with the FEniCS load vector. But calculating the function with X and Y variables:

#function which need to approximate
f = 3*x + 4*y

will yield the wrong load vector b which is not equal to FEniCS result.

I suppose the issue is related to linear mapping (affine mapping). Namely, I do not know how to transform X and Y variables in function from global coordinate system to local coordinate system. Do you know how to change these variables?

See for instance: Introduction to finite element methods (Section: Affine mapping of the reference cell)

It is an excellent book, working, thanks again

Code:

#affing mapping
x = (phi_local0 * i[0]) + (phi_local1 * j[0]) + (phi_local2 * k[0])

#function which need to approximate
f = 3*x

#b vector, integrating in local (natural) coordinates
b[0] = sym.integrate(phi_local0 * f * detJ, (n, 0, 1-e), (e, 0, 1))
b[1] = sym.integrate(phi_local1 * f * detJ, (n, 0, 1-e), (e, 0, 1))
b[2] = sym.integrate(phi_local2 * f * detJ, (n, 0, 1-e), (e, 0, 1))
print('\nVector b:', b)