The integration regions Ω Ω are then chosen as rectangles centered around node and extending to the 4 nearest neighbor nodes. The gradient ∇ ∇ φ φ is then approximated using centered difference normal to the boundary of the integration region Ω Ω , and average ϵ ϵ over the integration surface ∂ ∂ Ω Ω . This approach allows us to approximate the left hand side of the Poisson equation above in first order as