Most of the above described methods for the mesh generation are accompanied by some sort of mesh optimization technique [27,33,41,53,63] in order to improve the shape of the elements, their connectivity and to eliminate poor quality elements. This mesh enhancement yields the improved accuracy and stability of the numerical solution and reduces the number of elements required to capture the underlying physical phenomenon. The Laplacian smoothing technique, changing the position of nodes without modifying the topology of the mesh, is the most commonly used method applied to 2D and 3D meshes. It consists in solving Laplacian equation for the locations of the interior nodes with given positions of boundary nodes. The solution can be accomplished in a computationally inexpensive way by iterative process in which each interior node is repositioned into the centre (appropriate metric should be used if an anisotropic mesh is considered) of a polygon (in 2D) or a polyhedron (in 3D) formed by connected elements. This process is repeated until there is (almost) no movement in the mesh. The concept of the Laplacian smoothing can be also extended to 3D surfaces. In the case of smoothing the anisotropic mesh in the parametric space, the relevant metric must be considered. In the real space, however, regardless whether working with the polygon projected to the tangent plane at the node being smoothed or with the raw ``non-planar'' polygon, the smoothed node has to be projected back to the surface to satisfy the surface constraint. Although the application of the Laplacian smoothing typically improves the quality of the mesh in average, the quality of the worst element may be spoiled even more. It is therefore essential to combine the smoothing with some topology and/or optimization based technique. The topology based approaches for triangular meshes rest usually upon the edge swapping . The optimization based techniques, on the other hand, exploit either variational methods  for the detection of the global optimum with respect to a particular cost function or line search methods  to optimize the elements in the local neighbourhood of the node subjected to the smoothing. However, the penalty paid for the high effectiveness of optimization based methods in the elimination of the most severely distorted elements is their computational expense.