Write a program which estimates the error in the potential field at each node using a simple energy norm estimator. Utilize a simple averaging procedure to identify which elements need to be refined (can you think of a better way to do this?). Add an extra node at the centroid of the element, which has been identified for refinement (Be careful about how you treat the elements which border a surface. Why?).
Use the sweep2 program to re-triangulate the mesh. After removing the extraneous elements, calculate the FE solution and compare it to the original solution. As you continue the process, you are bound to notice a couple of things. One, after a few iterations, the amount of change in your solution from the previous solution will be minimal. Secondly, the number of degrees of freedom in your system will increase rapidly, such that you may soon not be able to calculate the solution on your workstation (especially if it has limited memory or you have limited time). This may mean that you are converging to the true solution (if you have done the refinements correctly).