title: Adaptive mesh refinement - Python
body: <p>i'm currently incredibly stuck on what isn't working in my code and have been staring at it for hours. I have created some functions to approximate the solution to the laplace equation adaptively using the finite element method then estimate it's error using the dual weighted residual. The error function should give a vector of errors (one error for each element), i then choose the biggest errors, add more elements around them, solve again and then recheck the error; however i have no idea why my error estimate isn't changing!</p>
  
  <p>My first 4 functions are correct but i will include them incase someone wants to try the code:</p>
  
  <pre><code>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_rich_node): #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))
      Q[qoi_rich_node] = 1.0
  
      Z = scipy.sparse.linalg.spsolve(K,Q)
      return Z
  </code></pre>
  
  <p>My error indicator function takes in an approximation Uh, with the mesh it is solved over, and finds eta = (f - Bu)z. </p>
  
  <pre><code>def Error_Indicators(Uh,U_mesh,Z,Z_mesh,f):
      """Take in U, Interpolate to same mesh as Z then solve for eta vector"""
  
      u_inter = interp1d(U_mesh,Uh) #Interpolation of old mesh
      U2 = u_inter(Z_mesh) #New function u for the new mesh to use in 
  
      Bz = Poisson_Stiffness(Z_mesh)
      Bz = Bz.tocsr()
  
      eta = np.empty(len(Z_mesh))
      for i in range(len(Z_mesh)):
          for j in range(len(Z_mesh)):
              eta[i] += (f[i] - Bz[i,j]*U2[j])
  
      for i in range(len(Z)):
          eta[i] = eta[i]*Z[i]
  
      return eta
  </code></pre>
  
  <p>My next function seems to adapt the mesh very well to the given error indicator! Just no idea why the indicator seems to stay the same regardless?</p>
  
  <pre><code>def Mesh_Refinement(base_mesh,tolerance,refinement,z_mesh,QOI_z_mesh):
      """Solve for U on a normal mesh, Take in Z, Find error indicators, adapt. OUTPUT NEW MESH"""
      New_mesh = base_mesh
      Z = Dualsolution(z_mesh,QOI_z_mesh) #Solve dual solution only once
  
      f = np.empty(len(z_mesh))
      for i in range(len(z_mesh)):
          f[i] = math.sin(math.pi*z_mesh[i])
  
      U = Solver(New_mesh)
      eta = Error_Indicators(U,base_mesh,Z,z_mesh,f) 
  
      while max(abs(k) for k in eta) &gt; tolerance:
  
          orderedeta = np.sort(eta) #Sort error indicators LENGTH 40
          biggest = np.flipud(orderedeta[int((1-refinement)*len(eta)):len(eta)]) 
          position = np.empty(len(biggest))
          ratio = float(len(New_mesh))/float(len(z_mesh))
  
          for i in range(len(biggest)):
              position[i] = eta.tolist().index(biggest[i])*ratio #GIVES WHAT NUMBER NODE TO REFINE
  
          refine = np.zeros(len(position))
          for i in range(len(position)):
              refine[i] = math.floor(position[i])+0.5 #AT WHAT NODE TO PUT NEW ELEMENT 5.5 ETC
          refine = np.flipud(sorted(set(refine)))
  
          for i in range(len(refine)):
              New_mesh = np.insert(New_mesh,refine[i]+0.5,(New_mesh[refine[i]+0.5]+New_mesh[refine[i]-0.5])/2)
  
          U = Solver(New_mesh)
          eta = Error_Indicators(U,New_mesh,Z,z_mesh,f)
          print eta
  </code></pre>
  
  <p>An example input for this would be:
       Mesh_Refinement(np.linspace(0,1,3),0.1,0.2,np.linspace(0,1,60),20)</p>
  
  <p>I understand there is alot of code here but i am at a loss, i have no idea where to turn!</p>