After the a priori boundary recovery was successfully completed, which means that there are no cyclic dependencies in the dependency graph, the actual domain discretization is quite straightforward. The initial Delaunay triangulation is setuped in the form of 5 tetrahedra filling the parallelepiped constructed during the process of boundary recovery. Then all the points from set (eventually updated due to newly formed tetrahedra) are incrementally inserted into the triangulation using the modified Watson's algorithm. This is accomplished by traversing the list of constraining faces (also eventually updated) and inserting the not yet inserted vertices of each constraining face, following the dependency hierarchy of the dependency graph. If a particular constraining face does not appear in the intermediate constrained Delaunay triangulation, it is recovered using the 3-to-2 topological transformation. After all the points have been inserted, the outer tetrahedra, except those formed during the a priori boundary recovery process, which are now reclassified as interior, are identified in the resulting constrained Delaunay triangulation and removed from it. In the next phase, new points are incrementally inserted in the interior of the domain, while keeping the boundary constraints, to make the tetrahedra of appropriate size and of aspect ratio close to one. The new point is always located at the barycenter of those tetrahedra size of which does not comply with the user prescribed mesh density distribution or whose aspect ratio exceeds a user prescribed limit. After all internal points have been inserted, the perturbed points are returned to their original positions and the resulting mesh is then subjected to the optimization in terms of the combination of Laplacian smoothing and topological transformations, in order to remove the potential slivers and to improve the overall mesh quality.