I'm currently attempting to turn my code for solving the laplace equation using finite element approximations into an adaptive one using the dual weighted residual as my error estimator: i.e. my indicator for element i: $_i=(f[i]_0^N A_{ij}U_h[j])z[i]$. i.e. I set up f, A and z over a very fine mesh and solve U_h over the normal mesh then interpolate it so i can perform the above operation.I am solving for the dual solution z relatively simply; i.e. B(v,z) = Q(v) where B is the stiffness matrix and Q is a vector of 0s with 1 at the quantity of interest.However the errors do not behave how they should; as an example i got the mesh to just increase from 10 in uniform multiples of 10 until 150. However these are the error estimates (as solved above):0.03354620650380.03354620650380.03354620650380.03354620650380.03354620650380.03336085919160.03354620650380.03346145048040.03331588307980.03354620650380.03343175378620.03333441042010.03354415448210.0334327061457They seem to drop then grow, drop then grow repeatedly and i have no idea why, does anybody have any clues what could be causing this?The following is my codedef rhs(x0):    rhsfunction = np.empty(len(x0))    for i in range(len(x0)):        rhsfunction[i] = math.sin(math.pi*x0[i])    return rhsfunctiondef 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_Matrixdef 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 CONDITIONSreturn approxdef Solver(x0):    Stiff_Matrix = Poisson_Stiffness(x0)    NodalApproximation = NodalQuadrature(x0)    NodalApproximation[0] = 0    U = scipy.sparse.linalg.spsolve(Stiff_Matrix, NodalApproximation)return Udef Dualsolution(rich_mesh,qoi): #BOUNDARY CONDITIONS?    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 Zdef Indicators(richx0,basex0,qoi):interpolate U, eta_i = r_i z_i = (f_i - sum(AijUj)Zi    U = Solver(basex0)    U_inter = interp1d(basex0,U)    U_rich = U_inter(richx0)    f = rhs(richx0)    A = Poisson_Stiffness(richx0).tocsr()    Z = Dualsolution(richx0,qoi)    eta = np.empty(len(richx0))    for i in range(len(eta)):        eta[i] = f[i]        for j in range(len(richx0)):            eta[i] = eta[i] - A[i,j]*U_rich[j]        eta[i] = eta[i] * Z[i]I understand theres alot of code here but i have no idea where else to turn. I'm almost certain my functions for solving the finite element method are correct (i.e. the Stiffness matrix, nodal quadrature and the solver). Any help would be appreciated alot.Thanks, JamesEDIT, this is my plot of the error indicator which must be incorrect but i cannot figure out what is wrong with it.
Adaptive Finite Element Method - Laplace
finite element;python;mesh generation;adaptive mesh refinement