I have been trying to set up an adaptive mesh for the finite element method for the laplace equation (and thus AU = F where A is the stiffness matrix) however have been running into some trouble with the error estimates and thus keep ending up with meshes like the following:
[ 0. 0.11111111 0.22222222 0.27777778 0.33333333 0.44444444 0.47222222 0.47916667 0.48090278 0.48111979 0.48133681 0.48177083 0.48263889 0.48611111 0.5 0.55555556 0.66666667 0.77777778 0.88888889 1. ]
As you can see a lot of the nodes are grouping around the same place due to the error estimate not decreasing at this point.The very basic idea for the error estimate is that we have two meshes, our adaptive mesh and a rich/bigger mesh. We find A, Z and our approximation of F over the rich mesh and U over our adaptive mesh and then interpolate it over the fine mesh and calculate the following:
eta_i = (f_i - sum_j (A[i,j]U[j]))*z[i]
And a large eta_i will mean element_i needs refining. Here is my code in which i am certain Poisson_Stiffness (finds A), Nodal_Quadrature (finds F), Solver(finds U) and DualSolution (finds Z) are all correct. I cannot tell if the problem is in the SubdivisionIndicators function or the Mesh refinement, however I'm assuming the former as usually the problem stems from no reduction in error and thus a build up of elements around the same place.
import numpy as np
import math
import scipy
from scipy.sparse import diags
import scipy.sparse.linalg
from scipy.interpolate import interp1d
import matplotlib
import matplotlib.pyplot as plt
def Poisson_Stiffness(x0):
"""Finds the Poisson equation stiffness matrix with any non uniform mesh x0"""
x0 = np.array(x0)
N = len(x0) - 1 # The amount of elements; x0, x1, ..., xN
h = x0[1:] - x0[:-1]
a = np.zeros(N+1)
a[0] = 1 #BOUNDARY CONDITIONS
a[1:-1] = 1/h[1:] + 1/h[:-1]
a[-1] = 1/h[-1]
a[N] = 1 #BOUNDARY CONDITIONS
b = -1/h
b[0] = 0 #BOUNDARY CONDITIONS
c = -1/h
c[N-1] = 0 #BOUNDARY CONDITIONS: DIRICHLET
data = [a.tolist(), b.tolist(), c.tolist()]
Positions = [0, 1, -1]
Stiffness_Matrix = diags(data, Positions, (N+1,N+1))
return Stiffness_Matrix
def NodalQuadrature(x0):
"""Finds the Nodal Quadrature Approximation of sin(pi x)"""
x0 = np.array(x0)
h = x0[1:] - x0[:-1]
N = len(x0) - 1
approx = np.zeros(len(x0))
approx[0] = 0 #BOUNDARY CONDITIONS
for i in range(1,N):
approx[i] = math.sin(math.pi*x0[i])
approx[i] = (approx[i]*h[i-1] + approx[i]*h[i])/2
approx[N] = 0 #BOUNDARY CONDITIONS
return approx
def Solver(x0):
Stiff_Matrix = Poisson_Stiffness(x0)
NodalApproximation = NodalQuadrature(x0)
NodalApproximation[0] = 0
U = scipy.sparse.linalg.spsolve(Stiff_Matrix, NodalApproximation)
return U
def Dualsolution(rich_mesh,qoi):
"""Find Z from stiffness matrix Z = K^-1 Q over richer mesh"""
K = Poisson_Stiffness(rich_mesh)
Q = np.zeros(len(rich_mesh))
if qoi in rich_mesh:
qoi_rich_node = rich_mesh.tolist().index(qoi)
Q[qoi_rich_node] = 1.0
else:
nearest = find_nearest(rich_mesh,qoi)
if nearest < qoi:
qoi_lower_node = rich_mesh.tolist().index(nearest)
qoi_upper_node = qoi_lower_node+1
else:
qoi_upper_node = rich_mesh.tolist().index(nearest)
qoi_lower_node = qoi_upper_node-1
Q[qoi_upper_node] = (qoi - rich_mesh[qoi_lower_node])/(rich_mesh[qoi_upper_node]-rich_mesh[qoi_lower_node])
Q[qoi_lower_node] = (rich_mesh[qoi_upper_node] - qoi)/(rich_mesh[qoi_upper_node]-rich_mesh[qoi_lower_node])
Z = scipy.sparse.linalg.spsolve(K,Q)
return Z
def SubdivisionIndicators(richx0,basex0,qoi):
"""interpolate U, eta_i = r_i z_i = (f_i - sum(AijUj)Zi"""
U = Solver(basex0) #finds U over base mesh
U_inter = interp1d(basex0,U,kind='cubic') #interpolates U over rich mesh
U_rich = U_inter(richx0)
f = NodalQuadrature(richx0) #finds f over rich mesh
A = Poisson_Stiffness(richx0).tocsr() #finds A over rich mesh
Z = Dualsolution(richx0,qoi) #finds Z over rich mesh
eta = np.empty(len(richx0))
for i in range(len(eta)):
eta[i] = abs((f[i] - A[i,:].dot(U_rich))*Z[i])
return eta
def Mesh_refine(base_mesh,amount_of_refinements,z_mesh,QOI):
refinedmesh = base_mesh #initialise
for i in range(amount_of_refinements):
eta = SubdivisionIndicators(z_mesh,refinedmesh,QOI)
maxpositionseta = np.nonzero(eta == eta.max())[0]
maxpositionseta = np.array(maxpositionseta)
ratio = float(len(refinedmesh))/float(len(z_mesh))
for i in range(len(maxpositionseta)):
print eta[maxpositionseta[i]]
maxposition = maxpositionseta[i]*ratio
if maxposition == int(maxposition):
refinedmesh = np.insert(refinedmesh,maxposition,(refinedmesh[maxposition] + refinedmesh[maxposition-1])/2)
else:
refinedmesh = np.insert(refinedmesh,math.ceil(maxposition),(refinedmesh[math.ceil(maxposition)]+refinedmesh[math.floor(maxposition)])/2)
print refinedmesh
return refinedmesh
Mesh_refine(np.linspace(0,1,10),10,np.linspace(0,1,100),0.5)
