code: def quasi_newton_uniform(points, cells, *args, **kwargs):
      
  
      def get_new_points(mesh):
          # do one Newton step
          # TODO need copy?
          x = mesh.node_coords.copy()
          cells = mesh.cells["nodes"]
          jac_x = jac_uniform(x, cells)
          x -= solve_hessian_approx_uniform(x, cells, jac_x)
          return x
  
      mesh = MeshTri(points, cells)
      runner(get_new_points, mesh, *args, **kwargs)
      return mesh.node_coords, mesh.cells["nodes"]
docs: Like linear_solve above, but assuming rho==1. Note that the energy gradient
  
          \\partial E_i = 2/(d+1) sum_{tau_j in omega_i} (x_i - b_j) \\int_{tau_j} rho
  
      becomes
  
          \\partial E_i = 2/(d+1) sum_{tau_j in omega_i} (x_i - b_j) |tau_j|.
  
      Because of the dependence of |tau_j| on the point coordinates, this is a nonlinear
      problem.
  
      This method makes the simplifying assumption that |tau_j| does in fact _not_ depend
      on the point coordinates. With this, one still only needs to solve a linear system.