Having the limit surface the goal is now to generate a new triangular mesh over it, respecting a given mesh density distribution (typically prescribed at nodes of the STL mesh serving as the control grid and/or curvature based) that makes the mesh suitable for subsequent computational analysis. The discretization is carried out in a hierarchical manner. Firstly, the model vertices are discretized. Then the (limit) curves are segmented using the mass curve of the required element density along that curve. And finally, the individual (limit) surfaces are triangulated.
In order to control the element size distribution, an octree is built around the domain to be discretized. The size of individual octants corresponds approximately to the required element spacing while the nodes (corners) of the octants are storing the required spacing exactly. To ensure the gradual variation of the element size, the maximum one-level octree difference of octants sharing an edge is enforced. This will guarantee creation of well shaped triangles. During the actual mesh generation the required element size is extracted from the octree for a given location using the interpolation of octree nodal values of the element size.
In the presented implementation, the surfaces are discretized by the advancing front technique constrained directly to the surface and modified to reflect surface curvature . Firstly, the initial front consisting of mesh edges at boundary curves of the surface (including inner loops) is established. Once the initial front has been set up, the mesh generation continues on the basis of an edge removal algorithm according to the following steps until the front becomes empty:
The generated mesh is then subjected to an optimization in order to improve the quality of the final mesh. The Laplacian smoothing technique in the combination with topological transformations (diagonal edge swapping) is adopted. This yields the optimized grid after only a few cycles of smoothing (typically up to six). Note that, contradictory to the smoothing carried out in 2D, the repositioning of a node during the smoothing is likely to shift the node out of the surface. Therefore, the node has to be projected back to the surface to satisfy the surface constraint.
A crucial aspect of the proposed mesh generation strategy is related to the point-to-surface projection. Simple and efficient algorithms available for the projection to parametric surfaces cannot be adopted simply because of missing parameterization of the limit surface. Moreover, the situation is further complicated by the fact, that the normal to the limit surface can be evaluated only at the bg nodes of the original or refined control grid. Therefore in order to make the projection sufficiently accurate (in terms of the distance from the limit surface and match with the exact normal), it is necessary to subdivide the original control grid up to a high level. This results in a huge amount of data to be stored, which is not acceptable. In , an efficient and reliable approach for the projection of a point to the limit surface has been proposed. This approach is based on the localized progressive refinement of the control grid towards the actual projection. A recursive implementation of this algorithm enables virtually an unlimited number of refinements with constant memory requirements to be performed. Note that the refinement is of temporary character and it is discarded after the projection is completed. Since some of the projections during the mesh generation can be accomplished with considerably lower accuracy (without any significant impact on the resulting mesh), an alternative approximate, but more efficient projection technique was suggested. This is based on approximating bg faces of the control grid by quadratic Bezier triangular patches and on the employment of a standard projection technique applicable to parametric surfaces (see  for details).
With respect to the application of the above meshing technique to the limit surfaces reconstructed from the STL meshes, it should be noted that planar surfaces (seemingly the simplest case) must be handled in a special way. The reason is that planar surface, having zero curvature in all directions, prevents the use of the curvature based mechanism to tessellate it into STL bg faces. Therefore the aspect ratio and orientation (in plane) of individual bg faces cannot be related to the curvature. This allows two neighbouring bg faces forming the same planar surface to considerably differ in the aspect ratio measured with respect to the shared bg side. As a direct consequence, the limit surface folds over itself, possibly crossing the boundary of the surface. Even thought such a surface still seems to be planar, it is not any more, since the normal is changing orientation from point to point, which is fatal for the meshing algorithm. Therefore, when triangulating a planar surface, the subdivision is not invoked and bg faces of the control grid are used for localization purposes only.