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)