# Import Python libraries
from numpy import zeros, array
from scipy.sparse.linalg.dsolve import spsolve
from scipy import sparse
from pylab import clf
# *************** CONNECTIVITY ARRAYS ***************
# Function to construct element connectivity arrays
def create_connectivity_arrays(mesh):
nnode = len(mesh.nodes)
# create a temporary array 0, 1, 2, ..., nnode-1
temp = range(0, nnode)
# Turn off boundary nodes
for i in range(0, nnode):
if mesh.is_boundary_node(i): temp[i] = -1
# re-enumerate temporary array
counter = 0
for i in range(0, nnode):
if temp[i] >= 0:
temp[i] = counter
counter += 1
ndof = counter
# create elements connectivity arrays
con_array = []
for elem in mesh.elems:
a, b, c = elem
con_elem = []
con_elem.append(temp[a])
con_elem.append(temp[b])
con_elem.append(temp[c])
con_array.append(con_elem)
return con_array, ndof
# *************** UTILS ***************
# Function to calculate the area of a triangle.
# Input: vertices a, b, c (coordinate pairs).
def calc_triangle_area(a, b, c):
ux = b[0] - a[0]
uy = b[1] - a[1]
vx = c[0] - a[0]
vy = c[1] - a[1]
return 0.5*(ux*vy - uy*vx)
# *************** REFERENCE MAP AND TRANSFORMATIONS ***************
# Reference map transforming reference triangle ([-1,-1], [1,-1], [-1,1])
# onto physical mesh triangle (a, b, c). Here a, b, c are coordinate pairs.
def refmap(xi1, xi2, a, b, c):
x = a[0]*fn1(xi1, xi2) + b[0]*fn2(xi1, xi2) + c[0]*fn3(xi1, xi2)
y = a[1]*fn1(xi1, xi2) + b[1]*fn2(xi1, xi2) + c[1]*fn3(xi1, xi2)
return x, y
# Calculate inverse map derivatives, needed for transformation of
# derivatives of shape functions from the reference domain to the
# physical mesh triangle. Input: nodes a, b, c (coordinate pairs).
# This function is taken from Hermes2D.
def calc_inv_map(a, b, c):
m = zeros((2,2), float)
m[0][0] = b[0] - a[0]
m[0][1] = c[0] - a[0]
m[1][0] = b[1] - a[1]
m[1][1] = c[1] - a[1]
jacobian = 0.25 * (m[0][0] * m[1][1] - m[0][1] * m[1][0])
ij = 0.5 / jacobian
inv_map = zeros((2,2), float)
inv_map[0][0] = m[1][1] * ij
inv_map[1][0] = -m[0][1] * ij
inv_map[0][1] = -m[1][0] * ij
inv_map[1][1] = m[0][0] * ij
return inv_map
# Transform derivatives from reference domain to physical mesh element
def transf_der(ref_dx, ref_dy, m):
phys_dx = m[0][0]*ref_dx + m[0][1]*ref_dy
phys_dy = m[1][0]*ref_dx + m[1][1]*ref_dy
return phys_dx, phys_dy
# *************** NUMERICAL QUADRATURE ***************
# Fourth-order Gauss quadrature on the reference triangle.
# List of points and weights in the form x_coord, y_coord, weight.
gauss_quad = [
( -0.108103018168070, -0.108103018168070, 0.446763179356022 ),
( -0.108103018168070, -0.783793963663860, 0.446763179356022 ),
( -0.783793963663860, -0.108103018168070, 0.446763179356022 ),
( -0.816847572980458, -0.816847572980458, 0.219903487310644 ),
( -0.816847572980458, 0.633695145960918, 0.219903487310644 ),
( 0.633695145960918, -0.816847572980458, 0.219903487310644 )
]
# *************** WEAK FORMS (EQUATION-DEPENDENT) ***************
# Bilinear form
# Integrals evaluated using 4th-order Gauss quadrature
def a_form(j, i, q, p, a, b, c, m, gauss_quad):
#print "**** a_form, i, j:", i, j
result = 0
for quad_point in gauss_quad:
xi1, xi2, ref_weight = quad_point
ref_val_j = fn[j](xi1, xi2)
phys_val_j = ref_val_j
ref_dx_j = dfndx[j](xi1, xi2)
ref_dy_j = dfndy[j](xi1, xi2)
phys_dx_j, phys_dy_j = transf_der(ref_dx_j, ref_dy_j, m)
ref_val_i = fn[i](xi1, xi2)
phys_val_i = ref_val_i
ref_dx_i = dfndx[i](xi1, xi2)
ref_dy_i = dfndy[i](xi1, xi2)
phys_dx_i, phys_dy_i = transf_der(ref_dx_i, ref_dy_i, m)
jacobian = 0.5 * calc_triangle_area(a, b, c)
phys_weight = ref_weight * jacobian
x1, x2 = refmap(xi1, xi2, a, b, c)
result += (q(x1, x2) * (phys_dx_j * phys_dx_i + phys_dy_j * phys_dy_i) +
q(x1, x2) * phys_val_j * phys_val_i) * phys_weight
return result
# Linear form
# Integrals evaluated using 4th-order Gauss quadrature
def l_form(i, f, a, b, c, gauss_quad):
result = 0
for quad_point in gauss_quad:
xi1, xi2, ref_weight = quad_point
ref_val_i = fn[i](xi1, xi2)
phys_val_i = ref_val_i
jacobian = 0.5 * calc_triangle_area(a, b, c)
phys_weight = ref_weight * jacobian
x1, x2 = refmap(xi1, xi2, a, b, c)
result += f(x1, x2) * phys_val_i * phys_weight
return result
# *************** ASSEMBLING PROCEDURE ***************
# Calculate element stiffness matrix and load vector
def calc_elem_integrals(q, p, f, a, b, c, m, gauss_quad):
mat = zeros((shape_fn_num, shape_fn_num), float)
vec = zeros([shape_fn_num])
for i in range(0, shape_fn_num):
for j in range(0, shape_fn_num):
mat[i][j] = a_form(j, i, q, p, a, b, c, m, gauss_quad)
vec[i] = l_form(i, f, a, b, c, gauss_quad)
return mat, vec
# Calculate stiffness matrices and load vectors in all elements
def assemble(q, p, f, mesh, ndof, con_array, gauss_quad):
glob_I = []
glob_J = []
glob_V = []
glob_vec = zeros(ndof, float)
elem_counter = 0
shape_fn_num = len(con_array[0])
for elem in mesh.elems:
a, b, c, = elem
m = calc_inv_map(mesh.nodes[a], mesh.nodes[b], mesh.nodes[c])
loc_mat, loc_vec = calc_elem_integrals(q, p, f, mesh.nodes[a], mesh.nodes[b], mesh.nodes[c], m, gauss_quad)
for loc_i in range(shape_fn_num):
glob_i = con_array[elem_counter][loc_i]
if glob_i >= 0: # test function active
for loc_j in range(shape_fn_num):
glob_j = con_array[elem_counter][loc_j]
if glob_j >= 0:
glob_I.append(glob_i)
glob_J.append(glob_j)
glob_V.append(loc_mat[loc_i][loc_j])
glob_vec[glob_i] += loc_vec[loc_i]
elem_counter += 1
return glob_I, glob_J, glob_V, glob_vec
# *************** DEFINITION OF SHAPE FUNCTIONS ***************
# Shape functions and their derivatives:
def fn1(xi1, xi2): return (-xi1 - xi2)/2.
def dfn1dx(xi1, xi2): return -1./2.
def dfn1dy(xi1, xi2): return -1./2.
def fn2(xi1, xi2): return (xi1 + 1.)/2.
def dfn2dx(xi1, xi2): return 1./2.
def dfn2dy(xi1, xi2): return 0.
def fn3(xi1, xi2): return (xi2 + 1.)/2.
def dfn3dx(xi1, xi2): return 0.
def dfn3dy(xi1, xi2): return 1./2.
fn = [fn1, fn2, fn3]
dfndx = [dfn1dx, dfn2dx, dfn3dx]
dfndy = [dfn1dy, dfn2dy, dfn3dy]
shape_fn_num = 3
# *************** EQUATION DATA ***************
# Equation solved: -div(q grad u) + p u = f
# Boundary conditions: Zero Dirichlet
# Constant parameters:
def q(x,y):
return 1.
def p(x,y):
return 0.
def f(x,y):
return 1.
