I am implementing a finite element solver (in 2D) to solve the generic differential equation :
$$-\nabla(a(x) \nabla u) = f$$
Brief explanation
By integrating and multipling by a test function, the weak form allows to integrate the stiffness matrix :
$$ A_{ij} = \int a(x) \nabla \phi_i(x) \nabla \phi_j(x)\, dx\quad i=1,2,3\ .$$
More simply, we can write the local stiffness matrix (of a 3 node linear triangle) of the element $K$ (the assembly will be done in a next step):
$$A_{ij}^{K} = \int a(x) \nabla \phi_i(x) \nabla \phi_j(x)\, dx\quad i=1,2,3\ .$$
In first place, I use the central quadrature formula that allows to compute:
$$A_{ij}^{K} = a_\text{mean}(b_i b_j + c_i c_j) K_\text{area}$$
with $a_\text{mean}$ = central value of $a(x)$ in the center of the triangle. The coefficient are the one form the shape function: $\phi_i = a_i + b_i x1 + c_i x_2$
For one single local stiffness matrix, it represents this code in Python:
loc2glob = t[0:3,K] #local to global map
x = p[0,loc2glob] #get a list of the x-coordinates
y = p[1,loc2glob] #get a list of the y-coordinates
area,b,c = gradient(x,y) # compute the numerical gradient
xc = np.mean(x), yc = np.mean(y)
a_centroid = a(xc,yc) #simplification
b = np.atleast_2d(b)
c = np.atleast_2d(c)
AK = (np.dot(b.T, b) + np.dot(c.T, c))*a_centroid*area
And ...?
After the verification of this method, I end up with implementing a more general integration method (based on quadrature). The coordinates $x,y$ are transformed in $r,s$ (local integration, to have a general method). The Jacobian is computed numerically.
$$A_{ij}^{K} = \int_{el_K} a(r, s) \nabla \phi_i (r, s) \nabla \phi_i(r,s) |J_K(r,s)|\ dr\ ds$$
A little more explanation on $(r,s)$ coordinate
In the domain of finite element method, computer implementation are used to switch the coordinate $(x,y)$ of a (triangle) element to an $(r,s)$ orthogonal base. The principle is very nice, it allows to do a variable change in the integral, and to have an easiest method to apply the quadrature :
$$\int_{\Omega_K} q(x,y) \phi_i(x,y)\phi_j(x,y)\ dx\ dy = \int_{\triangle} q(\xi, \eta) \psi_i(\xi, \eta)\psi_j(\xi, \eta) \left|\frac{\partial (x, y)}{\partial (\xi, \eta)}\right|\ d\xi\ d\eta$$
Graphically, it corresponds to :
So you can use the coordinate $(r,s)$ of a generic triangle for the integral, providing the Jacobian of the coordinate transformation.
In practice, I use a 3 gauss point method. The basics are here:
AK = np.zeros((3,3))
for q in range(len(qwgts)):
r = rspts[q,0] # r coordinate of the q_th quadrature point
s = rspts[q,1] # s coordinate of the q_th quadrature point
a_gauss_point = a(r,s)
S,dSdx,dSdy,detJ = isopmap(x,y,r,s,shapefcn)
dSdx =np.matrix(dSdx),dSdy =np.matrix(dSdy)
wxarea = a_gauss_point*qwgts[q]*detJ/2 #weight * area (corrected by 0.5)
AK = AK + (dSdx.T.dot(dSdx) + dSdy.T.dot(dSdy))* wxarea
I know that the method "isopmap" return a correct Jacobian, since I have tested with data found in the literature.
My problem...
When I take $a(x,y) = 1$ , the two methods give the same result. But when I put $a(x,y) = x + y$ for example, I end with significant difference between the two matrix.
And so?
Do you have any explanation? Of course the problem is about the value $a(r,s)$ in the second code, but I don't see where the thing fails...
SOLUTION : EDIT (17.07.2017)
In fact it comes that I didn't understand properly the logic of the mapping (x,y) domain to (r,s) domain. In fact the integral can be written like this:
Then I just have to change in my python code the evaluation of the function $a(x,y)$:
x_physical = np.dot(x,S)
y_physical = np.dot(y,S)
a_gauss_point = a(x_physical,y_physical)
If you want to reproduce this behaviour (EDIT : correct code uploaded)
You can download the Python coded that contains 100 lines of code in order to compute very easily the $A_K$ matrix.
The following is the code
from __future__ import division # avoid integer problem of division
#Import zone
import numpy as np
#Define shape function
def P1shapes(r,s):
S = np.array([1-r-s,r,s])
dSdr = np.array([-1,1,0])
dSds = np.array([-1,0,1])
return S,dSdr,dSds
#Jacobian function
def isopmap(x,y,r,s,shapefcn):
# x = vector of x coordinate of the element's point
#shapefcn = P1shapes or P2shapes
S,dSdr,dSds = shapefcn(r,s);
j11=np.dot(dSdr,x)
j12=np.dot(dSdr,y)
j21=np.dot(dSds,x)
j22=np.dot(dSds,y)
detJ=j11*j22-j12*j21
dSdx=( j22*dSdr-j12*dSds)/detJ
dSdy=(-j21*dSdr+j11*dSds)/detJ
return S,dSdx,dSdy,detJ
#Gradient Function
def gradient(x,y):
#for linear shape function (1st order)
area = polyArea(x,y)
b = np.array([y[1]-y[2],y[2]-y[0],y[0]-y[1]])/2/area
c = np.array([x[2]-x[1],x[0]-x[2],x[1]-x[0]])/2/area
return area,b,c
def polyArea(x,y):
#Shoelace formula => Compute Area of Polygons, x and y set of pointd given (vertices coordinates)
return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))
#Define a(x,y)
def a(x,y):
return x + y
#Gauss point quadrature
qwgts=np.array([1/3,1/3,1/3])
rspts=np.array([[1/6,1/6],
[2/3,1/6],
[1/6,2/3]])
#Point of the 3 nodes triangles (test case)
x = [-21.68467035,-21.17695462 ,-22.18700401]
y = [-9.94204652 ,-8.91258056 ,-9.12362242]
#First method (simplifed one)
area,b,c = gradient(x,y) # compute the numerical gradient
xc = np.mean(x)
yc = np.mean(y)
a_centroid = a(xc,yc) #simplification
b = np.atleast_2d(b)
c = np.atleast_2d(c)
AK_simplified = (np.dot(b.T, b) + np.dot(c.T, c))*a_centroid*area
#Second Method
AK_quadrature = np.zeros((3,3))
for q in range(len(qwgts)):
r = rspts[q,0] # r coordinate of the q_th quadrature point
s = rspts[q,1] # s coordinate of the q_th quadrature point
S,dSdx,dSdy,detJ = isopmap(x,y,r,s,P1shapes)
#Map the gauss point to the physical domain (r,s)_gauss to (x,y)_gauss
x_physical = np.dot(x,S)
y_physical = np.dot(y,S)
a_gauss_point = a(x_physical,y_physical)
dSdx =np.matrix(dSdx)
dSdy =np.matrix(dSdy)
wxarea = a_gauss_point*qwgts[q]*detJ/2 #weight * area (corrected by 0.5)
AK_quadrature = AK_quadrature + (dSdx.T.dot(dSdx) + dSdy.T.dot(dSdy))* wxarea
#Show the matrix
print AK_simplified
print " "
print AK_quadrature
