The spatial localization for the intersection check and the front management, in terms of the edge selection and searching, are considered to be general bottlenecks of the advancing front procedure. Both aspects are addressed in the presented approach. The spatial localization is implemented using the octree data structure. This exhibits the computational complexity which becomes nearly constant for a reasonable octree depth. When searching for the most suitable candidate in the neighbourhood, the local Delaunay property is used as a criterion. This reduces the probability of the consequent occurrence of an intersection. Since gradual variation of the mesh size is ensured no special algorithm for the edge selection strategy is required. Therefore, the simplest approach has been chosen which always uses the first available edge in the front. This obviously leads to the constant computational complexity of the front management.
Another crucial aspect is related to the point-to-surface projection algorithm. To make the projection accurate (in terms of the distance from the limit surface) it is necessary to subdivide the control grid up to a sufficiently high level. This results in a huge amount of data to be stored, which is not acceptable. In the presented approach, the grid is refined (temporarily) only on a single triangle localized from the original control grid (Fig. 3a). The refinement of this triangle is accomplished progressively towards the current projection of the point being projected (Fig. 3b). In order to enable fast setup of the connectivity associated with that triangle and required during the local subdivision, one level of subdivision is computed globally (and stored) prior to the actual mesh generation. A robust implementation of the projection requires a reliable algorithm for triangle localization (including its sub-triangles selection) and a tool for selection correction. Currently, the selection algorithm is enhanced by taking into account the relevant midside nodes on appropriate level of the subdivision which determine approximately the bow of sides of triangles on that level. Simultaneously, the sub-triangle selection is traced during the refinement to enable backtracking if an invalid selection has been detected.
Since the stopping criterion of the progressive refinement is typically related to the target element size, it is evident that the point-to-surface projection is computationally sensitive to the mesh size. With the increasing mesh density the costs of the projection algorithm become strongly prohibitive. The remedy to this problem consists in the implementation of an ``approximate'' projection with much lower computational demands. This may be accomplished by using Bezier triangular patches to approximate the limit surface over individual elements of the control grid. The Bezier triangular patch of degree can be described by the following parametric formula 
The projection of a point to the Bezier triangular patch is accomplished in an iterative manner starting from the barycenter as the first approximation. In each iteration, the point is projected to the tangential plane (constructed at the most recent approximation) and then the increments of barycentric coordinates and , calculated using the patch gradients, are used to get better approximation of the projection. The whole process is repeated until the desired accuracy is achieved.
In the current implementation, the quadratic Bezier triangles have been employed. Their control polygon can be easily calculated using only nodes from the the first level of global subdivision without the necessity to proceed to further levels of subdivision. Note that keeping the order of the patch low increases the computational performance. On the other hand, quadratic patch is not capable to properly capture the change of the curvature from convex to concave. In such a case, the projection gives very unsatisfactory results. The criterion whether to use a patch is based on the match between the patch normal and the normal calculated on the control grid using normal evaluation masks. Should the normals deviate by more than 30 degrees at least at one of the patch vertices, the recursive subdivision based projection algorithm is used. However it is important to realize that for the final location of a node (typically in the last cycle of the mesh smoothing) the ``exact'' projection technique based on recursive subdivision has to be used in order to satisfy the constraint to the limit surface (unless it is planar).