The implementation of the Natural Element Method (NEM) by means of a Galerkin-based procedure is parallel to that adopted in FEM or Element-Free Galerkin (EFG) method [14], with the key distinction that separates the three is in the construction of the shape functions and their derivatives.

A computational procedure to evaluate the shape functions (n-n coordinates) is outlined in [3], which is extended by [9] to compute the derivatives of the interpolating function . In the application of n-n interpolation to PDEs (fluid-structure interaction), [10] adopted Lasserre's recursive formula [15] to compute the area (volume of a convex polyhedral in ) of the second-order Voronoi cells. The shape functions as well as their derivatives are computed. In [10], it is pointed out that Lasserre's formula [15] is more robust than the one due to [3], which breaks-down if the point lies along the edge of a Delaunay triangle. Expressions for the shape functions and their derivatives are provided [3,16,9,10].