0026106: BRepMesh - revision of data model
[occt.git] / src / BRepMesh / BRepMesh_FastDiscretFace.cxx
diff --git a/src/BRepMesh/BRepMesh_FastDiscretFace.cxx b/src/BRepMesh/BRepMesh_FastDiscretFace.cxx
deleted file mode 100644 (file)
index f50bc58..0000000
+++ /dev/null
@@ -1,1692 +0,0 @@
-// Created by: Ekaterina SMIRNOVA
-// Copyright (c) 2008-2014 OPEN CASCADE SAS
-//
-// This file is part of Open CASCADE Technology software library.
-//
-// This library is free software; you can redistribute it and/or modify it under
-// the terms of the GNU Lesser General Public License version 2.1 as published
-// by the Free Software Foundation, with special exception defined in the file
-// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
-// distribution for complete text of the license and disclaimer of any warranty.
-//
-// Alternatively, this file may be used under the terms of Open CASCADE
-// commercial license or contractual agreement.
-
-#include <BRepMesh_FastDiscretFace.hxx>
-
-#include <BRepMesh_PairOfPolygon.hxx>
-#include <BRepMesh_ShapeTool.hxx>
-#include <Poly_PolygonOnTriangulation.hxx>
-#include <Poly_Triangulation.hxx>
-
-#include <BRepAdaptor_Surface.hxx>
-#include <BRepAdaptor_HSurface.hxx>
-#include <BRepAdaptor_Curve.hxx>
-#include <Adaptor3d_IsoCurve.hxx>
-#include <Adaptor3d_HCurve.hxx>
-
-#include <BRep_ListIteratorOfListOfPointRepresentation.hxx>
-#include <BRep_PointRepresentation.hxx>
-#include <BRep_TVertex.hxx>
-#include <BRep_Tool.hxx>
-
-#include <GeomLib.hxx>
-#include <Geom_Surface.hxx>
-#include <Geom_BSplineSurface.hxx>
-#include <Geom_BezierSurface.hxx>
-#include <GCPnts_TangentialDeflection.hxx>
-#include <GCPnts_AbscissaPoint.hxx>
-
-#include <Standard_ErrorHandler.hxx>
-#include <Standard_Failure.hxx>
-#include <TColStd_Array1OfReal.hxx>
-#include <TColStd_ListOfInteger.hxx>
-#include <TColStd_SequenceOfReal.hxx>
-#include <TColStd_Array1OfInteger.hxx>
-#include <TColStd_HArray1OfReal.hxx>
-#include <TColgp_Array1OfPnt2d.hxx>
-#include <TopTools_DataMapOfShapeReal.hxx>
-
-#include <TopExp_Explorer.hxx>
-#include <TopoDS.hxx>
-#include <TopoDS_Vertex.hxx>
-#include <TopExp.hxx>
-
-#include <NCollection_Map.hxx>
-#include <Bnd_Box2d.hxx>
-
-#include <algorithm>
-
-//#define DEBUG_MESH "mesh.tcl"
-#ifdef DEBUG_MESH
-#include <fstream>
-#endif
-
-
-IMPLEMENT_STANDARD_RTTIEXT(BRepMesh_FastDiscretFace,Standard_Transient)
-
-static Standard_Real FUN_CalcAverageDUV(TColStd_Array1OfReal& P, const Standard_Integer PLen)
-{
-  Standard_Integer i, j, n = 0;
-  Standard_Real p, result = 0.;
-
-  for(i = 1; i <= PLen; i++)
-  {
-    // Sort
-    for(j = i + 1; j <= PLen; j++)
-    {
-      if(P(i) > P(j))
-      {
-        p = P(i);
-        P(i) = P(j);
-        P(j) = p;
-      }
-    }
-    // Accumulate
-    if (i != 1)
-    {
-      p = Abs(P(i) - P(i-1));
-      if(p > 1.e-7)
-      {
-        result += p;
-        n++;
-      }
-    }
-  }
-  return (n? (result / (Standard_Real) n) : -1.);
-}
-
-namespace
-{
-  Standard_Real deflectionOfSegment (
-    const gp_Pnt& theFirstPoint,
-    const gp_Pnt& theLastPoint,
-    const gp_Pnt& theMidPoint)
-  {
-    // 23.03.2010 skl for OCC21645 - change precision for comparison
-    if (theFirstPoint.SquareDistance (theLastPoint) > Precision::SquareConfusion ())
-    {
-      gp_Lin aLin (theFirstPoint, gp_Dir (gp_Vec (theFirstPoint, theLastPoint)));
-      return aLin.Distance (theMidPoint);
-    }
-
-    return theFirstPoint.Distance (theMidPoint);
-  }
-
-  Standard_Boolean IsCompexSurface (const GeomAbs_SurfaceType theType)
-  {
-    return (
-      theType != GeomAbs_Sphere   &&
-      theType != GeomAbs_Cylinder &&
-      theType != GeomAbs_Cone     &&
-      theType != GeomAbs_Torus);
-  }
-
-  //! Auxiliary class used to extract geometrical parameters of fixed TopoDS_Vertex.
-  class FixedVExplorer
-  {
-  public:
-
-    DEFINE_STANDARD_ALLOC
-
-    FixedVExplorer(const TopoDS_Vertex& theVertex)
-      : myVertex(theVertex)
-    {
-    }
-
-    const TopoDS_Vertex& Vertex() const
-    {
-      return myVertex;
-    }
-
-    Standard_Boolean IsSameUV() const
-    {
-      return Standard_False;
-    }
-
-    TopoDS_Vertex SameVertex() const
-    {
-      return TopoDS_Vertex();
-    }
-
-    gp_Pnt Point() const
-    {
-      return BRep_Tool::Pnt(myVertex);
-    }
-
-  private:
-
-    void operator =(const FixedVExplorer& /*theOther*/)
-    {
-    }
-
-  private:
-    const TopoDS_Vertex& myVertex;
-  };
-
-  void ComputeErrFactors (const Standard_Real               theDeflection, 
-                          const Handle(Adaptor3d_HSurface)& theFace,
-                          Standard_Real&                    theErrFactorU,
-                          Standard_Real&                    theErrFactorV)
-  {
-    theErrFactorU = theDeflection * 10.;
-    theErrFactorV = theDeflection * 10.;
-
-    switch (theFace->GetType ())
-    {
-    case GeomAbs_Cylinder:
-    case GeomAbs_Cone:
-    case GeomAbs_Sphere:
-    case GeomAbs_Torus:
-      break;
-
-    case GeomAbs_SurfaceOfExtrusion:
-    case GeomAbs_SurfaceOfRevolution:
-      {
-        Handle (Adaptor3d_HCurve) aCurve = theFace->BasisCurve ();
-        if (aCurve->GetType () == GeomAbs_BSplineCurve && aCurve->Degree () > 2)
-        {
-          theErrFactorV /= (aCurve->Degree () * aCurve->NbKnots ());
-        }
-        break;
-      }
-    case GeomAbs_BezierSurface:
-      {
-        if (theFace->UDegree () > 2)
-        {
-          theErrFactorU /= (theFace->UDegree ());
-        }
-        if (theFace->VDegree () > 2)
-        {
-          theErrFactorV /= (theFace->VDegree ());
-        }
-        break;
-      }
-    case GeomAbs_BSplineSurface:
-      {
-        if (theFace->UDegree () > 2)
-        {
-          theErrFactorU /= (theFace->UDegree () * theFace->NbUKnots ());
-        }
-        if (theFace->VDegree () > 2)
-        {
-          theErrFactorV /= (theFace->VDegree () *  theFace->NbVKnots ());
-        }
-        break;
-      }
-
-    case GeomAbs_Plane:
-    default:
-      theErrFactorU = theErrFactorV = 1.;
-    }
-  }
-
-  void AdjustCellsCounts (const Handle(Adaptor3d_HSurface)& theFace,
-                          const Standard_Integer            theNbVertices,
-                          Standard_Integer&                 theCellsCountU,
-                          Standard_Integer&                 theCellsCountV)
-  {
-    const GeomAbs_SurfaceType aType = theFace->GetType ();
-    if (aType == GeomAbs_OtherSurface)
-    {
-      // fallback to the default behavior
-      theCellsCountU = theCellsCountV = -1;
-      return;
-    }
-
-    Standard_Real aSqNbVert = theNbVertices;
-    if (aType == GeomAbs_Plane)
-    {
-      theCellsCountU = theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
-    }
-    else if (aType == GeomAbs_Cylinder || aType == GeomAbs_Cone)
-    {
-      theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
-    }
-    else if (aType == GeomAbs_SurfaceOfExtrusion || aType == GeomAbs_SurfaceOfRevolution)
-    {
-      Handle (Adaptor3d_HCurve) aCurve = theFace->BasisCurve ();
-      if (aCurve->GetType () == GeomAbs_Line ||
-          (aCurve->GetType () == GeomAbs_BSplineCurve && aCurve->Degree () < 2))
-      {
-        // planar, cylindrical, conical cases
-        if (aType == GeomAbs_SurfaceOfExtrusion)
-          theCellsCountU = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
-        else
-          theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
-      }
-      if (aType == GeomAbs_SurfaceOfExtrusion)
-      {
-        // V is always a line
-        theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
-      }
-    }
-    else if (aType == GeomAbs_BezierSurface || aType == GeomAbs_BSplineSurface)
-    {
-      if (theFace->UDegree () < 2)
-      {
-        theCellsCountU = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
-      }
-      if (theFace->VDegree () < 2)
-      {
-        theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
-      }
-    }
-
-    theCellsCountU = Max (theCellsCountU, 2);
-    theCellsCountV = Max (theCellsCountV, 2);
-  }
-}
-
-
-//=======================================================================
-//function : BRepMesh_FastDiscretFace
-//purpose  :
-//=======================================================================
-BRepMesh_FastDiscretFace::BRepMesh_FastDiscretFace(
-  const Standard_Real    theAngle,
-  const Standard_Real    theMinSize,
-  const Standard_Boolean isInternalVerticesMode,
-  const Standard_Boolean isControlSurfaceDeflection)
-: myAngle(theAngle),
-  myInternalVerticesMode(isInternalVerticesMode),
-  myMinSize(theMinSize),
-  myIsControlSurfaceDeflection(isControlSurfaceDeflection)
-{
-}
-
-//=======================================================================
-//function : Perform
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscretFace::Perform(const Handle(BRepMesh_FaceAttribute)& theAttribute)
-{
-  add(theAttribute);
-  commitSurfaceTriangulation();
-}
-
-//=======================================================================
-//function : initDataStructure
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscretFace::initDataStructure()
-{
-  const Standard_Real aTolU = myAttribute->ToleranceU();
-  const Standard_Real aTolV = myAttribute->ToleranceV();
-  const Standard_Real uCellSize = 14.0 * aTolU;
-  const Standard_Real vCellSize = 14.0 * aTolV;
-
-  const Standard_Real deltaX = myAttribute->GetDeltaX();
-  const Standard_Real deltaY = myAttribute->GetDeltaY();
-
-  Handle(NCollection_IncAllocator) aAllocator = 
-    new NCollection_IncAllocator(BRepMesh::MEMORY_BLOCK_SIZE_HUGE);
-  myStructure = new BRepMesh_DataStructureOfDelaun(aAllocator);
-  myStructure->Data()->SetCellSize ( uCellSize / deltaX, vCellSize / deltaY);
-  myStructure->Data()->SetTolerance( aTolU     / deltaX, aTolV     / deltaY);
-
-  myAttribute->ChangeStructure() = myStructure;
-  myAttribute->ChangeSurfacePoints() = new BRepMesh::DMapOfIntegerPnt(1, aAllocator);
-  myAttribute->ChangeSurfaceVertices()= new BRepMesh::DMapOfVertexInteger(1, aAllocator);
-
-  // Check the necessity to fill the map of parameters
-  const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
-  GeomAbs_SurfaceType thetype = gFace->GetType();
-  const Standard_Boolean isBSpline = (thetype == GeomAbs_BezierSurface ||
-                                      thetype == GeomAbs_BSplineSurface);
-  const Standard_Boolean useUVParam = (thetype == GeomAbs_Torus || IsCompexSurface (thetype));
-
-  Handle(NCollection_BaseAllocator) aBaseAlloc = aAllocator;
-  myUParam.Clear (aBaseAlloc);
-  myVParam.Clear (aBaseAlloc);
-
-  // essai de determination de la longueur vraie:
-  // akm (bug OCC16) : We must calculate these measures in non-singular
-  //     parts of face. Let`s try to compute average value of three
-  //     (umin, (umin+umax)/2, umax), and respectively for v.
-  //                 vvvvv
-  //Standard_Real longu = 0.0, longv = 0.0; //, last , first;
-  //gp_Pnt P11, P12, P21, P22, P31, P32;
-  BRepMesh::HVectorOfVertex& aBoundaryNodes = myAttribute->ChangeMeshNodes();
-  BRepMesh::VectorOfVertex::Iterator aNodesIt(*aBoundaryNodes);
-  for (; aNodesIt.More(); aNodesIt.Next())
-  {
-    BRepMesh_Vertex& aNode = aNodesIt.ChangeValue();
-    gp_XY aPnt2d = aNode.Coord();
-
-    if (useUVParam)
-    {
-      myUParam.Add(aPnt2d.X());
-      myVParam.Add(aPnt2d.Y());
-    }
-
-    aNode.ChangeCoord() = myAttribute->Scale(aPnt2d, Standard_True);
-    myStructure->AddNode(aNode, Standard_True);
-  }
-  aBoundaryNodes.Nullify();
-
-  if (isBSpline)
-  {
-    const Standard_Real aRange[2][2] = {
-      {myAttribute->GetUMin(), myAttribute->GetUMax()},
-      {myAttribute->GetVMin(), myAttribute->GetVMax()}
-    };
-
-    const GeomAbs_Shape aContinuity = GeomAbs_CN;
-    for (Standard_Integer i = 0; i < 2; ++i)
-    {
-      const Standard_Boolean isU = (i == 0);
-      const Standard_Integer aIntervalsNb = isU ?
-        gFace->NbUIntervals(aContinuity) :
-        gFace->NbVIntervals(aContinuity);
-
-      BRepMesh::IMapOfReal& aParams = isU ? myUParam : myVParam;
-      if (aIntervalsNb < aParams.Size())
-        continue;
-
-      TColStd_Array1OfReal aIntervals(1, aIntervalsNb + 1);
-      if (isU)
-        gFace->UIntervals(aIntervals, aContinuity);
-      else
-        gFace->VIntervals(aIntervals, aContinuity);
-
-      for (Standard_Integer j = 1; j <= aIntervals.Upper(); ++j)
-      {
-        const Standard_Real aParam = aIntervals(j);
-        if (aParam > aRange[i][0] && aParam < aRange[i][1])
-          aParams.Add(aParam);
-      }
-    }
-  }
-
-  //////////////////////////////////////////////////////////// 
-  //add internal vertices after self-intersection check
-  if ( myInternalVerticesMode )
-  {
-    TopExp_Explorer anExplorer(myAttribute->Face(), TopAbs_VERTEX, TopAbs_EDGE);
-    for ( ; anExplorer.More(); anExplorer.Next() )
-      add(TopoDS::Vertex(anExplorer.Current()));
-  }
-
-  const BRepMesh::HDMapOfShapePairOfPolygon& aEdges = myAttribute->ChangeInternalEdges();
-  TopExp_Explorer aWireIt(myAttribute->Face(), TopAbs_WIRE);
-  for (; aWireIt.More(); aWireIt.Next())
-  {
-    TopExp_Explorer aEdgeIt(aWireIt.Current(), TopAbs_EDGE);
-    for (; aEdgeIt.More(); aEdgeIt.Next())
-    {
-      const TopoDS_Edge& aEdge = TopoDS::Edge(aEdgeIt.Current());
-      BRepMesh_PairOfPolygon aPair;
-      if (!aEdges->Find(aEdge, aPair))
-        continue;
-
-      TopAbs_Orientation aOri = aEdge.Orientation();
-      const Handle(Poly_PolygonOnTriangulation)& aPolygon = 
-        aOri == TopAbs_REVERSED ? aPair.Last() : aPair.First();
-
-      const TColStd_Array1OfInteger& aIndices = aPolygon->Nodes();
-      const Standard_Integer aNodesNb = aPolygon->NbNodes();
-
-      Standard_Integer aPrevId = aIndices(1);
-      for (Standard_Integer i = 2; i <= aNodesNb; ++i)
-      {
-        const Standard_Integer aCurId = aIndices(i);
-        addLinkToMesh(aPrevId, aCurId, aOri);
-        aPrevId = aCurId;
-      }
-    }
-  }
-}
-
-//=======================================================================
-//function : addLinkToMesh
-//purpose  :
-//=======================================================================
-void BRepMesh_FastDiscretFace::addLinkToMesh(
-  const Standard_Integer   theFirstNodeId,
-  const Standard_Integer   theLastNodeId,
-  const TopAbs_Orientation theOrientation)
-{
-  if (theOrientation == TopAbs_FORWARD)
-    myStructure->AddLink(BRepMesh_Edge(theFirstNodeId, theLastNodeId, BRepMesh_Frontier));
-  else if (theOrientation == TopAbs_REVERSED)
-    myStructure->AddLink(BRepMesh_Edge(theLastNodeId, theFirstNodeId, BRepMesh_Frontier));
-  else if (theOrientation == TopAbs_INTERNAL)
-    myStructure->AddLink(BRepMesh_Edge(theFirstNodeId, theLastNodeId, BRepMesh_Fixed));
-}
-
-//=======================================================================
-//function : Add
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscretFace::add(const Handle(BRepMesh_FaceAttribute)& theAttribute)
-{
-  if (!theAttribute->IsValid() || theAttribute->ChangeMeshNodes()->IsEmpty())
-    return;
-
-  myAttribute = theAttribute;
-  initDataStructure();
-
-  BRepMesh::HIMapOfInteger& aVertexEdgeMap = myAttribute->ChangeVertexEdgeMap();
-  Standard_Integer nbVertices = aVertexEdgeMap->Extent();
-  BRepMesh::Array1OfInteger tabvert_corr(1, nbVertices);
-  for ( Standard_Integer i = 1; i <= nbVertices; ++i )
-    tabvert_corr(i) = i;
-
-  Handle (Adaptor3d_HSurface) aSurface (myAttribute->Surface ());
-  GeomAbs_SurfaceType aType = aSurface->GetType ();
-
-  while (aType == GeomAbs_OffsetSurface)
-  {
-    aSurface = aSurface->BasisSurface ();
-    aType = aSurface->GetType ();
-  }
-
-  // For better meshing performance we try to estimate the acceleration circles grid structure sizes:
-  // For each parametric direction (U, V) we estimate firstly an approximate distance between the future points -
-  // this estimation takes into account the required face deflection and the complexity of the face.
-  // Particularly, the complexity of the faces based on BSpline curves and surfaces requires much more points.
-  // At the same time, for planar faces and linear parts of the arbitrary surfaces usually no intermediate points
-  // are necessary.
-  // The general idea for each parametric direction:
-  // cells_count = 2 ^ log10 ( estimated_points_count )
-  // For linear parametric direction we fall back to the initial vertex count:
-  // cells_count = 2 ^ log10 ( initial_vertex_count )
-
-  Standard_Real anErrFactorU, anErrFactorV;
-  ComputeErrFactors(myAttribute->GetDefFace (), aSurface, anErrFactorU, anErrFactorV);
-
-  Standard_Integer aCellsCountU, aCellsCountV;
-  if (aType == GeomAbs_Torus)
-  {
-    aCellsCountU = (Standard_Integer)Ceiling(Pow(2, Log10(
-      (myAttribute->GetUMax() - myAttribute->GetUMin()) / myAttribute->GetDeltaX())));
-    aCellsCountV = (Standard_Integer)Ceiling(Pow(2, Log10(
-      (myAttribute->GetVMax() - myAttribute->GetVMin()) / myAttribute->GetDeltaY())));
-  }
-  else if (aType == GeomAbs_Cylinder)
-  {
-    aCellsCountU = (Standard_Integer)Ceiling(Pow(2, Log10(
-      (myAttribute->GetUMax() - myAttribute->GetUMin()) / myAttribute->GetDeltaX() /
-      (myAttribute->GetVMax() - myAttribute->GetVMin()))));
-    aCellsCountV = (Standard_Integer)Ceiling(Pow(2, Log10(
-      (myAttribute->GetVMax() - myAttribute->GetVMin()) / anErrFactorV)));
-  }
-  else
-  {
-    aCellsCountU = (Standard_Integer)Ceiling(Pow(2, Log10(
-      (myAttribute->GetUMax() - myAttribute->GetUMin()) / myAttribute->GetDeltaX() / anErrFactorU)));
-    aCellsCountV = (Standard_Integer)Ceiling(Pow(2, Log10(
-      (myAttribute->GetVMax() - myAttribute->GetVMin()) / myAttribute->GetDeltaY() / anErrFactorV)));
-  }
-
-  AdjustCellsCounts(aSurface, nbVertices, aCellsCountU, aCellsCountV);
-
-  BRepMesh_Delaun trigu(myStructure, tabvert_corr, aCellsCountU, aCellsCountV);
-
-  //removed all free edges from triangulation
-  const Standard_Integer nbLinks = myStructure->NbLinks();
-  for( Standard_Integer i = 1; i <= nbLinks; i++ ) 
-  {
-    if( myStructure->ElementsConnectedTo(i).Extent() < 1 )
-    {
-      BRepMesh_Edge& anEdge = (BRepMesh_Edge&)trigu.GetEdge(i);
-      if ( anEdge.Movability() == BRepMesh_Deleted )
-        continue;
-
-      anEdge.SetMovability(BRepMesh_Free);
-      myStructure->RemoveLink(i);
-    }
-  }
-
-  Standard_Boolean rajout = 
-    (aType == GeomAbs_Sphere || aType == GeomAbs_Torus);
-
-  // Check the necessity to fill the map of parameters
-  const Standard_Boolean useUVParam = (aType == GeomAbs_Torus         ||
-                                       aType == GeomAbs_BezierSurface ||
-                                       aType == GeomAbs_BSplineSurface);
-
-  const Standard_Real umax = myAttribute->GetUMax();
-  const Standard_Real umin = myAttribute->GetUMin();
-  const Standard_Real vmax = myAttribute->GetVMax();
-  const Standard_Real vmin = myAttribute->GetVMin();
-
-  Standard_Boolean isaline = 
-    ((umax - umin) < Precision::PConfusion() || 
-     (vmax - vmin) < Precision::PConfusion());
-
-  Standard_Real aDef = -1;
-  if ( !isaline && myStructure->ElementsOfDomain().Extent() > 0 )
-  {
-    if (!rajout)
-    {
-      // compute maximal deflection
-      aDef = control(trigu, Standard_True);
-      rajout = (aDef > myAttribute->GetDefFace() || aDef < 0.);
-    }
-    if (aType != GeomAbs_Plane)
-    {
-      if (!rajout && useUVParam)
-      {
-        rajout = (myVParam.Extent() > 2 &&
-          (aSurface->IsUClosed() || aSurface->IsVClosed()));
-      }
-
-      if (rajout)
-      {
-        insertInternalVertices(trigu);
-
-        //control internal points
-        if (myIsControlSurfaceDeflection)
-          aDef = control(trigu, Standard_False);
-      }
-    }
-  }
-
-  //modify myStructure back
-  BRepMesh::HVectorOfVertex& aMeshNodes = myStructure->Data()->ChangeVertices();
-  for ( Standard_Integer i = 1; i <= myStructure->NbNodes(); ++i )
-  {
-    BRepMesh_Vertex& aNode = aMeshNodes->ChangeValue(i - 1);
-    aNode.ChangeCoord() = myAttribute->Scale(aNode.Coord(), Standard_False);
-
-    const BRepMesh::ListOfInteger& alist = myStructure->LinksConnectedTo(i);
-    // Register internal nodes used in triangulation
-    if (!alist.IsEmpty() && aVertexEdgeMap->FindIndex(i) == 0)
-      aVertexEdgeMap->Add(i);
-  }
-
-  if (!(aDef < 0.))
-    myAttribute->SetDefFace(aDef);
-}
-
-//=======================================================================
-//function : addVerticesToMesh
-//purpose  : 
-//=======================================================================
-Standard_Boolean BRepMesh_FastDiscretFace::addVerticesToMesh(
-  const BRepMesh::ListOfVertex& theVertices,
-  BRepMesh_Delaun&              theMeshBuilder)
-{
-  if (theVertices.IsEmpty())
-    return Standard_False;
-
-  BRepMesh::Array1OfVertexOfDelaun aArrayOfNewVertices(1, theVertices.Extent());
-  BRepMesh::ListOfVertex::Iterator aVertexIt(theVertices);
-  for (Standard_Integer aVertexId = 0; aVertexIt.More(); aVertexIt.Next())
-    aArrayOfNewVertices(++aVertexId) = aVertexIt.Value();
-
-  theMeshBuilder.AddVertices(aArrayOfNewVertices);
-  return Standard_True;
-}
-
-//=======================================================================
-//function : filterParameters
-//purpose  : 
-//=======================================================================
-static void filterParameters(const BRepMesh::IMapOfReal& theParams,
-                             const Standard_Real         theMinDist,
-                             const Standard_Real         theFilterDist,
-                             BRepMesh::SequenceOfReal&   theResult)
-{
-  // Sort sequence of parameters
-  const Standard_Integer anInitLen = theParams.Extent();
-    
-  TColStd_Array1OfReal aParamArray(1, anInitLen);
-  Standard_Integer j;
-  for (j = 1; j <= anInitLen; j++)
-    aParamArray(j) = theParams(j);
-
-  std::sort (aParamArray.begin(), aParamArray.end());
-
-  // mandatory pre-filtering using the first (minimal) filter value
-  Standard_Integer aParamLength = 1;
-  for (j = 2; j <= anInitLen; j++) 
-  {
-    if ((aParamArray(j)-aParamArray(aParamLength)) > theMinDist)
-    {
-      if (++aParamLength < j)
-        aParamArray(aParamLength) = aParamArray(j);
-    }
-  }
-  
-  //perform filtering on series
-  Standard_Real aLastAdded, aLastCandidate;
-  Standard_Boolean isCandidateDefined = Standard_False;
-  aLastAdded = aParamArray(1);
-  aLastCandidate = aLastAdded;
-  theResult.Append(aLastAdded);
-  
-  for(j=2; j < aParamLength; j++)
-  {
-    Standard_Real aVal = aParamArray(j);
-    if(aVal-aLastAdded > theFilterDist) 
-    {
-      //adds the parameter
-      if(isCandidateDefined) {
-        aLastAdded = aLastCandidate;
-        isCandidateDefined = Standard_False;
-        j--;
-      }
-      else 
-      {
-        aLastAdded = aVal;
-      }
-      theResult.Append(aLastAdded);
-      continue;
-    }
-    
-    aLastCandidate = aVal;
-    isCandidateDefined = Standard_True;
-  }
-  theResult.Append(aParamArray(aParamLength));
-}
-
-//=======================================================================
-//function : insertInternalVertices
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscretFace::insertInternalVertices(BRepMesh_Delaun& theMeshBuilder)
-{
-  Handle(NCollection_IncAllocator) anAlloc = new NCollection_IncAllocator;
-  BRepMesh::ListOfVertex aNewVertices(anAlloc);
-  const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
-  switch (gFace->GetType())
-  {
-  case GeomAbs_Sphere:
-    insertInternalVerticesSphere(aNewVertices);
-    break;
-
-  case GeomAbs_Cylinder:
-    insertInternalVerticesCylinder(aNewVertices);
-    break;
-
-  case GeomAbs_Cone:
-    insertInternalVerticesCone(aNewVertices);
-    break;
-
-  case GeomAbs_Torus:
-    insertInternalVerticesTorus(aNewVertices);
-    break;
-
-  default:
-    insertInternalVerticesOther(aNewVertices);
-    break;
-  }
-  
-  addVerticesToMesh(aNewVertices, theMeshBuilder);
-}
-
-//=======================================================================
-//function : insertInternalVerticesSphere
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscretFace::insertInternalVerticesSphere(
-  BRepMesh::ListOfVertex& theNewVertices)
-{
-  Standard_Real aRange[] = {
-    myAttribute->GetVMin(), myAttribute->GetVMax(),
-    myAttribute->GetUMin(), myAttribute->GetUMax()
-  };
-
-  gp_Sphere aSphere = myAttribute->Surface()->Sphere();
-
-  // Calculate parameters for iteration in V direction
-  Standard_Real aStep = 0.7 * GCPnts_TangentialDeflection::ArcAngularStep(
-    aSphere.Radius(), myAttribute->GetDefFace(), myAngle, myMinSize);
-
-  Standard_Real aDd[2] = {aStep, aStep};
-  Standard_Real aPasMax[2] = {0., 0.};
-  for (Standard_Integer i = 0; i < 2; ++i)
-  {
-    const Standard_Real aMax = aRange[2 * i + 1];
-    const Standard_Real aDiff = aMax - aRange[2 * i + 0];
-    aDd[i] = aDiff / ((Standard_Integer)(aDiff / aDd[i]) + 1);
-    aPasMax[i] = aMax - Precision::PConfusion();
-  }
-
-  const Standard_Real aHalfDu = aDd[1] * 0.5;
-  Standard_Boolean Shift = Standard_False;
-  Standard_Real aPasV = aRange[0] + aDd[0];
-  for (; aPasV < aPasMax[0]; aPasV += aDd[0])
-  {
-    Shift = !Shift;
-    const Standard_Real d = (Shift) ? aHalfDu : 0.;
-    Standard_Real aPasU = aRange[2] + d;
-    for (; aPasU < aPasMax[1]; aPasU += aDd[1])
-    {
-      tryToInsertAnalyticVertex(gp_Pnt2d(aPasU, aPasV), aSphere, theNewVertices);
-    }
-  }
-}
-
-//=======================================================================
-//function : insertInternalVerticesCylinder
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscretFace::insertInternalVerticesCylinder(
-  BRepMesh::ListOfVertex& theNewVertices)
-{
-  const Standard_Real umax = myAttribute->GetUMax();
-  const Standard_Real umin = myAttribute->GetUMin();
-  const Standard_Real vmax = myAttribute->GetVMax();
-  const Standard_Real vmin = myAttribute->GetVMin();
-
-  gp_Cylinder aCylinder = myAttribute->Surface()->Cylinder();
-  const Standard_Real aRadius = aCylinder.Radius();
-
-  Standard_Integer nbU = 0;
-  Standard_Integer nbV = 0;
-  const Standard_Real su = umax - umin;
-  const Standard_Real sv = vmax - vmin;
-  const Standard_Real aArcLen = su * aRadius;
-  if (aArcLen > myAttribute->GetDefFace ())
-  {
-    // Calculate parameters for iteration in U direction
-    const Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep (
-      aRadius, myAttribute->GetDefFace (), myAngle, myMinSize);
-    nbU = (Standard_Integer)(su / Du);
-
-    // Calculate parameters for iteration in V direction
-    const Standard_Real aDv = nbU*sv / aArcLen;
-    // Protection against overflow during casting to int in case 
-    // of long cylinder with small radius.
-    nbV = aDv > static_cast<Standard_Real> (IntegerLast ()) ? 
-      0 : (Standard_Integer)(aDv);
-    nbV = Min (nbV, 100 * nbU);
-  }
-
-  const Standard_Real Du = su / (nbU + 1);
-  const Standard_Real Dv = sv / (nbV + 1);
-
-  Standard_Real pasu, pasv, pasvmax = vmax - Dv*0.5, pasumax = umax - Du*0.5;
-  for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv)
-  {
-    for (pasu = umin + Du; pasu < pasumax; pasu += Du)
-    {
-      tryToInsertAnalyticVertex(gp_Pnt2d(pasu, pasv), aCylinder, theNewVertices);
-    }
-  }
-}
-
-//=======================================================================
-//function : insertInternalVerticesCone
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscretFace::insertInternalVerticesCone(
-  BRepMesh::ListOfVertex& theNewVertices)
-{
-  const Standard_Real umax = myAttribute->GetUMax();
-  const Standard_Real umin = myAttribute->GetUMin();
-  const Standard_Real vmax = myAttribute->GetVMax();
-  const Standard_Real vmin = myAttribute->GetVMin();
-
-  gp_Cone aCone = myAttribute->Surface()->Cone();
-  Standard_Real RefR = aCone.RefRadius();
-  Standard_Real SAng = aCone.SemiAngle();
-  Standard_Real aRadius = Max(Abs(RefR + vmin*Sin(SAng)), Abs(RefR + vmax*Sin(SAng)));
-
-  Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep(
-    aRadius, myAttribute->GetDefFace(), myAngle, myMinSize);
-
-  Standard_Real Dv, pasu, pasv;
-  Standard_Integer nbU = (Standard_Integer)((umax - umin) / Du);
-  Standard_Integer nbV = (Standard_Integer)(nbU*(vmax - vmin) / ((umax - umin)*aRadius));
-  Du = (umax - umin) / (nbU + 1);
-  Dv = (vmax - vmin) / (nbV + 1);
-
-  Standard_Real pasvmax = vmax - Dv*0.5, pasumax = umax - Du*0.5;
-  for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv)
-  {
-    for (pasu = umin + Du; pasu < pasumax; pasu += Du)
-    {
-      tryToInsertAnalyticVertex(gp_Pnt2d(pasu, pasv), aCone, theNewVertices);
-    }
-  }
-}
-
-//=======================================================================
-//function : insertInternalVerticesTorus
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscretFace::insertInternalVerticesTorus(
-  BRepMesh::ListOfVertex& theNewVertices)
-{
-  const Standard_Real umax     = myAttribute->GetUMax();
-  const Standard_Real umin     = myAttribute->GetUMin();
-  const Standard_Real vmax     = myAttribute->GetVMax();
-  const Standard_Real vmin     = myAttribute->GetVMin();
-  const Standard_Real deltaX   = myAttribute->GetDeltaX();
-  const Standard_Real deltaY   = myAttribute->GetDeltaY();
-  const Standard_Real aDefFace = myAttribute->GetDefFace();
-
-  gp_Torus T = myAttribute->Surface()->Torus();
-
-  Standard_Boolean insert;
-  Standard_Integer i, j, ParamULength, ParamVLength;
-  Standard_Real pp, pasu, pasv;
-  Standard_Real r = T.MinorRadius(), R = T.MajorRadius();
-
-  BRepMesh::SequenceOfReal ParamU, ParamV;
-
-  Standard_Real oldDv = GCPnts_TangentialDeflection::ArcAngularStep(
-    r, aDefFace, myAngle, myMinSize);
-
-  Standard_Real Dv = 0.9*oldDv; //TWOTHIRD * oldDv;
-  Dv = oldDv;
-
-  Standard_Real Du;
-  Standard_Integer nbV = Max((Standard_Integer)((vmax - vmin) / Dv), 2);
-  Dv = (vmax - vmin) / (nbV + 1);
-  Standard_Real ru = R + r;
-  if (ru > 1.e-16)
-  {
-    Du = GCPnts_TangentialDeflection::ArcAngularStep(
-      ru, aDefFace, myAngle, myMinSize);
-
-    Standard_Real aa = sqrt(Du*Du + oldDv*oldDv);
-    if (aa < gp::Resolution())
-      return;
-    Du *= Min(oldDv, Du) / aa;
-  }
-  else Du = Dv;
-
-  Standard_Integer nbU = Max((Standard_Integer)((umax - umin) / Du), 2);
-  nbU = Max(nbU, (int)(nbV*(umax - umin)*R / ((vmax - vmin)*r) / 5.));
-  Du = (umax - umin) / (nbU + 1);
-
-  if (R < r)
-  {
-    // As the points of edges are returned.
-    // in this case, the points are not representative.
-
-    //-- Choose DeltaX and DeltaY so that to avoid skipping points on the grid
-    for (i = 0; i <= nbU; i++) ParamU.Append(umin + i* Du);
-  }//R<r
-  else //U if R > r
-  {
-    //--ofv: U
-    // Number of mapped U parameters
-    const Standard_Integer LenU = myUParam.Extent();
-    // Fill array of U parameters
-    TColStd_Array1OfReal Up(1, LenU);
-    for (j = 1; j <= LenU; j++) Up(j) = myUParam(j);
-
-    // Calculate DU, leave array of parameters
-    Standard_Real aDU = FUN_CalcAverageDUV(Up, LenU);
-    aDU = Max(aDU, Abs(umax - umin) / (Standard_Real)nbU / 2.);
-    Standard_Real dUstd = Abs(umax - umin) / (Standard_Real)LenU;
-    if (aDU > dUstd) dUstd = aDU;
-    // Add U parameters
-    for (j = 1; j <= LenU; j++)
-    {
-      pp = Up(j);
-      insert = Standard_True;
-      ParamULength = ParamU.Length();
-      for (i = 1; i <= ParamULength && insert; i++)
-      {
-        insert = (Abs(ParamU.Value(i) - pp) > (0.5*dUstd));
-      }
-      if (insert) ParamU.Append(pp);
-    }
-  }
-
-  //--ofv: V
-  // Number of mapped V parameters
-  const Standard_Integer LenV = myVParam.Extent();
-  // Fill array of V parameters
-  TColStd_Array1OfReal Vp(1, LenV);
-  for (j = 1; j <= LenV; j++) Vp(j) = myVParam(j);
-  // Calculate DV, sort array of parameters
-  Standard_Real aDV = FUN_CalcAverageDUV(Vp, LenV);
-  aDV = Max(aDV, Abs(vmax - vmin) / (Standard_Real)nbV / 2.);
-
-  Standard_Real dVstd = Abs(vmax - vmin) / (Standard_Real)LenV;
-  if (aDV > dVstd) dVstd = aDV;
-  // Add V parameters
-  for (j = 1; j <= LenV; j++)
-  {
-    pp = Vp(j);
-
-    insert = Standard_True;
-    ParamVLength = ParamV.Length();
-    for (i = 1; i <= ParamVLength && insert; i++)
-    {
-      insert = (Abs(ParamV.Value(i) - pp) > (dVstd*2. / 3.));
-    }
-    if (insert) ParamV.Append(pp);
-  }
-
-  Standard_Integer Lu = ParamU.Length(), Lv = ParamV.Length();
-  Standard_Real uminnew = umin + deltaY*0.1;
-  Standard_Real vminnew = vmin + deltaX*0.1;
-  Standard_Real umaxnew = umax - deltaY*0.1;
-  Standard_Real vmaxnew = vmax - deltaX*0.1;
-
-  for (i = 1; i <= Lu; i++)
-  {
-    pasu = ParamU.Value(i);
-    if (pasu >= uminnew && pasu < umaxnew)
-    {
-      for (j = 1; j <= Lv; j++)
-      {
-        pasv = ParamV.Value(j);
-        if (pasv >= vminnew && pasv < vmaxnew)
-        {
-          tryToInsertAnalyticVertex(gp_Pnt2d(pasu, pasv), T, theNewVertices);
-        }
-      }
-    }
-  }
-}
-
-//=======================================================================
-//function : insertInternalVerticesOther
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscretFace::insertInternalVerticesOther(
-  BRepMesh::ListOfVertex& theNewVertices)
-{
-  const Standard_Real aRange[2][2] = {
-      { myAttribute->GetUMax(), myAttribute->GetUMin() },
-      { myAttribute->GetVMax(), myAttribute->GetVMin() }
-  };
-
-  const Standard_Real aDelta[2] = { 
-    myAttribute->GetDeltaX(),
-    myAttribute->GetDeltaY()
-  };
-
-  const Standard_Real                 aDefFace = myAttribute->GetDefFace();
-  const Handle(BRepAdaptor_HSurface)& gFace    = myAttribute->Surface();
-
-  Handle(NCollection_IncAllocator) anAlloc = new NCollection_IncAllocator;
-  BRepMesh::SequenceOfReal aParams[2] = { BRepMesh::SequenceOfReal(anAlloc), 
-                                          BRepMesh::SequenceOfReal(anAlloc) };
-  for (Standard_Integer i = 0; i < 2; ++i)
-  {
-    Standard_Boolean isU = (i == 0);
-    Standard_Real aRes = isU ?
-      gFace->UResolution(aDefFace) :
-      gFace->VResolution(aDefFace);
-
-    // Sort and filter sequence of parameters
-    Standard_Real aMinDiff = Precision::PConfusion();
-    if (aDelta[i] < 1.)
-      aMinDiff /= aDelta[i];
-
-    aMinDiff = Max(myMinSize, aMinDiff);
-
-    Standard_Real aRangeDiff = aRange[i][0] - aRange[i][1];
-    Standard_Real aDiffMaxLim = 0.1 * aRangeDiff;
-    Standard_Real aDiffMinLim = Max(0.005 * aRangeDiff, 2. * aRes);
-    Standard_Real aDiff = Max(myMinSize, Min(aDiffMaxLim, aDiffMinLim));
-    filterParameters(isU ? myUParam : myVParam, aMinDiff, aDiff, aParams[i]);
-  }
-
-  // check intermediate isolines
-  Handle (Geom_Surface) aSurface = gFace->ChangeSurface ().Surface ().Surface ();
-
-  BRepMesh::MapOfReal aParamsToRemove[2] = { BRepMesh::MapOfReal(1, anAlloc),
-                                             BRepMesh::MapOfReal(1, anAlloc) };
-  BRepMesh::MapOfReal aParamsForbiddenToRemove[2] = { BRepMesh::MapOfReal(1, anAlloc),
-                                                      BRepMesh::MapOfReal(1, anAlloc) };
-
-  // insert additional points where it is needed to conform criteria.
-  // precision for normals calculation
-  const Standard_Real aNormPrec = Precision::Confusion();
-  for (Standard_Integer k = 0; k < 2; ++k)
-  {
-    const Standard_Integer aOtherIndex = (k + 1) % 2;
-    BRepMesh::SequenceOfReal& aParams1 = aParams[k];
-    BRepMesh::SequenceOfReal& aParams2 = aParams[aOtherIndex];
-    const Standard_Boolean isU = (k == 0);
-    for (Standard_Integer i = 2; i < aParams1.Length(); ++i)
-    {
-      const Standard_Real aParam1 = aParams1(i);
-      GeomAdaptor_Curve aIso(isU ?
-        aSurface->UIso(aParam1) : aSurface->VIso(aParam1));
-
-      Standard_Real aPrevParam2 = aParams2(1);
-      gp_Pnt aPrevPnt2;
-      gp_Vec aPrevVec2;
-      aIso.D1(aPrevParam2, aPrevPnt2, aPrevVec2);
-      for (Standard_Integer j = 2; j <= aParams2.Length();)
-      {
-        Standard_Real aParam2 = aParams2(j);
-        gp_Pnt aPnt2;
-        gp_Vec aVec2;
-        aIso.D1(aParam2, aPnt2, aVec2);
-
-        Standard_Real aMidParam = 0.5 * (aPrevParam2 + aParam2);
-        gp_Pnt aMidPnt = aIso.Value(aMidParam);
-
-        Standard_Real aDist = deflectionOfSegment(aPrevPnt2, aPnt2, aMidPnt);
-        if (aDist > aDefFace && aDist > myMinSize)
-        {
-          // insertion 
-          aParams2.InsertBefore(j, aMidParam);
-          continue;
-        }
-        //put regular grig for normals
-        gp_Pnt2d aStPnt1, aStPnt2;
-        if (isU)
-        {
-          aStPnt1 = gp_Pnt2d(aParam1, aPrevParam2);
-          aStPnt2 = gp_Pnt2d(aParam1, aMidParam);
-        }
-        else
-        {
-          aStPnt1 = gp_Pnt2d(aPrevParam2, aParam1);
-          aStPnt2 = gp_Pnt2d(aMidParam, aParam1);
-        }
-
-        gp_Dir N1(0, 0, 1), N2(0, 0, 1);
-        Standard_Integer aSt1 = GeomLib::NormEstim(aSurface, aStPnt1, aNormPrec, N1);
-        Standard_Integer aSt2 = GeomLib::NormEstim(aSurface, aStPnt2, aNormPrec, N2);
-
-        const Standard_Real aAngle = N2.Angle(N1);
-        if (aSt1 < 1 && aSt2 < 1 && aAngle > myAngle)
-        {
-          const Standard_Real aLen = GCPnts_AbscissaPoint::Length(
-            aIso, aPrevParam2, aMidParam, aDefFace);
-
-          if (aLen > myMinSize)
-          {
-            // insertion 
-            aParams2.InsertBefore(j, aMidParam);
-            continue;
-          }
-        }
-        aPrevParam2 = aParam2;
-        aPrevPnt2 = aPnt2;
-        aPrevVec2 = aVec2;
-
-        ++j;
-      }
-    }
-  }
-#ifdef DEBUG_InsertInternal
-  // output numbers of parameters along U and V
-  cout << "aParams: " << aParams[0].Length() << " " << aParams[1].Length() << endl;
-#endif
-  // try to reduce number of points evaluating of isolines sampling
-  for (Standard_Integer k = 0; k < 2; ++k)
-  {
-    const Standard_Integer aOtherIndex = (k + 1) % 2;
-    BRepMesh::SequenceOfReal& aParams1 = aParams[k];
-    BRepMesh::SequenceOfReal& aParams2 = aParams[aOtherIndex];
-    const Standard_Boolean isU = (k == 0);
-    BRepMesh::MapOfReal& aToRemove2          = aParamsToRemove[aOtherIndex];
-    BRepMesh::MapOfReal& aForbiddenToRemove1 = aParamsForbiddenToRemove[k];
-    BRepMesh::MapOfReal& aForbiddenToRemove2 = aParamsForbiddenToRemove[aOtherIndex];
-    for (Standard_Integer i = 2; i < aParams1.Length(); ++i)
-    {
-      const Standard_Real aParam1 = aParams1(i);
-      GeomAdaptor_Curve aIso(isU ?
-        aSurface->UIso (aParam1) : aSurface->VIso (aParam1));
-#ifdef DEBUG_InsertInternal
-      // write polyline containing initial parameters to the file
-      {
-        ofstream ff(DEBUG_InsertInternal, std::ios_base::app);
-        ff << "polyline " << (k == 0 ? "u" : "v") << i << " ";
-        for (Standard_Integer j = 1; j <= aParams2.Length(); j++)
-        {
-          gp_Pnt aP;
-          aIso.D0(aParams2(j), aP);
-          ff << aP.X() << " " << aP.Y() << " " << aP.Z() << " ";
-        }
-        ff << endl;
-      }
-#endif
-
-      Standard_Real aPrevParam2 = aParams2(1);
-      gp_Pnt aPrevPnt2;
-      gp_Vec aPrevVec2;
-      aIso.D1 (aPrevParam2, aPrevPnt2, aPrevVec2);
-      for (Standard_Integer j = 2; j <= aParams2.Length();)
-      {
-        Standard_Real aParam2 = aParams2(j);
-        gp_Pnt aPnt2;
-        gp_Vec aVec2;
-        aIso.D1 (aParam2, aPnt2, aVec2);
-
-        // Here we should leave at least 3 parameters as far as
-        // we must have at least one parameter related to surface
-        // internals in order to prevent movement of triangle body
-        // outside the surface in case of highly curved ones, e.g.
-        // BSpline springs.
-        if (aParams2.Length() > 3 && j < aParams2.Length())
-        {
-          // Remove too dense points
-          const Standard_Real aNextParam = aParams2(j + 1);
-          gp_Pnt aNextPnt;
-          gp_Vec aNextVec;
-          aIso.D1(aNextParam, aNextPnt, aNextVec);
-
-          // Lets check current parameter.
-          // If it fits deflection, we can remove it.
-          Standard_Real aDist = deflectionOfSegment(aPrevPnt2, aNextPnt, aPnt2);
-          if (aDist < aDefFace)
-          {
-            // Lets check parameters for angular deflection.
-            if (aPrevVec2.SquareMagnitude() > gp::Resolution() &&
-                aNextVec.SquareMagnitude() > gp::Resolution() &&
-                aPrevVec2.Angle(aNextVec) < myAngle)
-            {
-              // For current Iso line we can remove this parameter.
-#ifdef DEBUG_InsertInternal
-              // write point of removed parameter
-              {
-                ofstream ff(DEBUG_InsertInternal, std::ios_base::app);
-                ff << "point " << (k == 0 ? "u" : "v") << i << "r_" << j << " "
-                  << aPnt2.X() << " " << aPnt2.Y() << " " << aPnt2.Z() << endl;
-              }
-#endif
-              aToRemove2.Add(aParam2);
-              aPrevParam2 = aNextParam;
-              aPrevPnt2 = aNextPnt;
-              aPrevVec2 = aNextVec;
-              j += 2;
-              continue;
-            }
-            else {
-              // We have found a place on the surface refusing 
-              // removement of this parameter.
-#ifdef DEBUG_InsertInternal
-              // write point of forbidden to remove parameter
-              {
-                ofstream ff(DEBUG_InsertInternal, std::ios_base::app);
-                ff << "vertex " << (k == 0 ? "u" : "v") << i << "f_" << j << " "
-                  << aPnt2.X() << " " << aPnt2.Y() << " " << aPnt2.Z() << endl;
-              }
-#endif
-              aForbiddenToRemove1.Add(aParam1);
-              aForbiddenToRemove2.Add(aParam2);
-            }
-          }
-        }
-        aPrevParam2 = aParam2;
-        aPrevPnt2 = aPnt2;
-        aPrevVec2 = aVec2;
-        ++j;
-      }
-    }
-  }
-  // remove filtered out parameters
-  for (Standard_Integer k = 0; k < 2; ++k)
-  {
-    BRepMesh::SequenceOfReal& aParamsk = aParams[k];
-    for (Standard_Integer i = 1; i <= aParamsk.Length();)
-    {
-      const Standard_Real aParam = aParamsk.Value(i);
-      if (aParamsToRemove[k].Contains(aParam) &&
-        !aParamsForbiddenToRemove[k].Contains(aParam))
-      {
-        aParamsk.Remove(i);
-      }
-      else
-        i++;
-    }
-  }
-#ifdef DEBUG_InsertInternal
-  // write polylines containing remaining parameters
-  {
-    ofstream ff(DEBUG_InsertInternal, std::ios_base::app);
-    for (Standard_Integer k = 0; k < 2; ++k)
-    {
-      for (Standard_Integer i = 1; i <= aParams[k].Length(); i++)
-      {
-        ff << "polyline " << (k == 0 ? "uo" : "vo") << i << " ";
-        for (Standard_Integer j = 1; j <= aParams[1 - k].Length(); j++)
-        {
-          gp_Pnt aP;
-          if (k == 0)
-            gFace->D0(aParams[k](i), aParams[1 - k](j), aP);
-          else
-            gFace->D0(aParams[1 - k](j), aParams[k](i), aP);
-          ff << aP.X() << " " << aP.Y() << " " << aP.Z() << " ";
-        }
-        ff << endl;
-      }
-    }
-  }
-#endif
-
-  // insert nodes of the regular grid
-  const BRepMesh::HClassifier& aClassifier = myAttribute->ChangeClassifier();
-  for (Standard_Integer i = 1; i <= aParams[0].Length(); ++i)
-  {
-    const Standard_Real aParam1 = aParams[0].Value (i);
-    for (Standard_Integer j = 1; j <= aParams[1].Length(); ++j)
-    {
-      const Standard_Real aParam2 = aParams[1].Value (j);
-      gp_Pnt2d aPnt2d(aParam1, aParam2);
-
-      // Classify intersection point
-      if (aClassifier->Perform(aPnt2d) == TopAbs_IN)
-      {
-        gp_Pnt aPnt;
-        gFace->D0(aPnt2d.X(), aPnt2d.Y(), aPnt);
-        insertVertex(aPnt, aPnt2d.Coord(), theNewVertices);
-      }
-    }
-  }
-}
-
-//=======================================================================
-//function : checkDeflectionAndInsert
-//purpose  : 
-//=======================================================================
-Standard_Boolean BRepMesh_FastDiscretFace::checkDeflectionAndInsert(
-  const gp_Pnt&              thePnt3d,
-  const gp_XY&               theUV,
-  const Standard_Boolean     isDeflectionCheckOnly,
-  const Standard_Real        theTriangleDeflection,
-  const Standard_Real        theFaceDeflection,
-  const BRepMesh_CircleTool& theCircleTool,
-  BRepMesh::ListOfVertex&    theVertices,
-  Standard_Real&             theMaxTriangleDeflection,
-  const Handle(NCollection_IncAllocator)& theTempAlloc)
-{
-  if (theTriangleDeflection > theMaxTriangleDeflection)
-    theMaxTriangleDeflection = theTriangleDeflection;
-
-  if (theTriangleDeflection < theFaceDeflection)
-    return Standard_True;
-
-  if (myMinSize > Precision::Confusion())
-  {
-    // Iterator in the list of indexes of circles containing the node
-    BRepMesh::ListOfInteger& aCirclesList = 
-      const_cast<BRepMesh_CircleTool&>(theCircleTool).Select(
-      myAttribute->Scale(theUV, Standard_True));
-    
-    BRepMesh::MapOfInteger aUsedNodes(10, theTempAlloc);
-    BRepMesh::ListOfInteger::Iterator aCircleIt(aCirclesList);
-    for (; aCircleIt.More(); aCircleIt.Next())
-    {
-      const BRepMesh_Triangle& aTriangle = 
-        myStructure->GetElement(aCircleIt.Value());
-
-      Standard_Integer aNodes[3];
-      myStructure->ElementNodes(aTriangle, aNodes);
-
-      for (Standard_Integer i = 0; i < 3; ++i)
-      {
-        const Standard_Integer aNodeId = aNodes[i];
-        if (aUsedNodes.Contains(aNodeId))
-          continue;
-
-        aUsedNodes.Add(aNodeId);
-        const BRepMesh_Vertex& aNode = myStructure->GetNode(aNodeId);
-        const gp_Pnt& aPoint = myAttribute->GetPoint(aNode);
-
-        if (thePnt3d.SquareDistance(aPoint) < myMinSize * myMinSize)
-          return Standard_True;
-      }
-    }
-  }
-
-  if (isDeflectionCheckOnly)
-    return Standard_False;
-
-  insertVertex(thePnt3d, theUV, theVertices);
-  return Standard_True;
-}
-
-//=======================================================================
-//function : control
-//purpose  : 
-//=======================================================================
-Standard_Real BRepMesh_FastDiscretFace::control(
-  BRepMesh_Delaun&         theTrigu,
-  const Standard_Boolean   theIsFirst)
-{
-  Standard_Integer aTrianglesNb = myStructure->ElementsOfDomain().Extent();
-  if (aTrianglesNb < 1)
-    return -1.0;
-
-  //IMPORTANT: Constants used in calculations
-  const Standard_Real MinimalArea2d     = 1.e-9;
-  const Standard_Real MinimalSqLength3d = 1.e-12;
-  const Standard_Real aSqDefFace = myAttribute->GetDefFace() * myAttribute->GetDefFace();
-
-  const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
-
-  Handle(Geom_Surface) aBSpline;
-  const GeomAbs_SurfaceType aSurfType = gFace->GetType ();
-  if (IsCompexSurface (aSurfType) && aSurfType != GeomAbs_SurfaceOfExtrusion)
-    aBSpline = gFace->ChangeSurface ().Surface().Surface();
-
-  Handle(NCollection_IncAllocator) anAlloc =
-    new NCollection_IncAllocator(BRepMesh::MEMORY_BLOCK_SIZE_HUGE);
-  NCollection_DataMap<Standard_Integer, gp_Dir> aNorMap(1, anAlloc);
-  BRepMesh::MapOfIntegerInteger                 aStatMap(1, anAlloc);
-  NCollection_Map<BRepMesh_Edge>                aCouples(3 * aTrianglesNb, anAlloc);
-  const BRepMesh_CircleTool& aCircles = theTrigu.Circles();
-
-  // Perform refinement passes
-  // Define the number of iterations
-  Standard_Integer       aIterationsNb = 11;
-  const Standard_Integer aPassesNb = (theIsFirst ? 1 : aIterationsNb);
-  // Initialize stop condition
-  Standard_Real aMaxSqDef = -1.;
-  Standard_Integer aPass = 1, aInsertedNb = 1;
-  Standard_Boolean isAllDegenerated = Standard_False;
-  Handle(NCollection_IncAllocator) aTempAlloc =
-    new NCollection_IncAllocator(BRepMesh::MEMORY_BLOCK_SIZE_HUGE);
-  for (; aPass <= aPassesNb && aInsertedNb && !isAllDegenerated; ++aPass)
-  {
-    aTempAlloc->Reset(Standard_False);
-    BRepMesh::ListOfVertex aNewVertices(aTempAlloc);
-
-    // Reset stop condition
-    aInsertedNb      = 0;
-    aMaxSqDef        = -1.;
-    isAllDegenerated = Standard_True;
-
-    aTrianglesNb = myStructure->ElementsOfDomain().Extent();
-    if (aTrianglesNb < 1)
-      break;
-
-    // Iterate on current triangles
-    const BRepMesh::MapOfInteger& aTriangles = myStructure->ElementsOfDomain();
-    BRepMesh::MapOfInteger::Iterator aTriangleIt(aTriangles);
-    for (; aTriangleIt.More(); aTriangleIt.Next())
-    {
-      const Standard_Integer aTriangleId = aTriangleIt.Key();
-      const BRepMesh_Triangle& aCurrentTriangle = myStructure->GetElement(aTriangleId);
-
-      if (aCurrentTriangle.Movability() == BRepMesh_Deleted)
-        continue;
-
-      Standard_Integer v[3];
-      myStructure->ElementNodes(aCurrentTriangle, v);
-
-      Standard_Integer e[3];
-      Standard_Boolean o[3];
-      aCurrentTriangle.Edges(e, o);
-
-      gp_XY xy[3];
-      gp_XYZ p[3];
-      Standard_Boolean m[3];
-      for (Standard_Integer i = 0; i < 3; ++i)
-      {
-        m[i] = (myStructure->GetLink(e[i]).Movability() == BRepMesh_Frontier);
-
-        const BRepMesh_Vertex& aVertex = myStructure->GetNode(v[i]);
-        xy[i] = myAttribute->Scale(aVertex.Coord(), Standard_False);
-        p [i] = myAttribute->GetPoint(aVertex).Coord();
-      }
-
-      gp_XYZ aLinkVec[3];
-      Standard_Boolean isDegeneratedTri = Standard_False;
-      for (Standard_Integer i = 0; i < 3 && !isDegeneratedTri; ++i)
-      {
-        aLinkVec[i] = p[(i + 1) % 3] - p[i];
-        isDegeneratedTri = (aLinkVec[i].SquareModulus() < MinimalSqLength3d);
-      }
-
-      if (isDegeneratedTri) 
-        continue;
-
-      isAllDegenerated = Standard_False;
-
-      // Check triangle area in 2d
-      if (Abs((xy[1]-xy[0])^(xy[2]-xy[1])) < MinimalArea2d)
-        continue;
-
-      // Check triangle normal
-      gp_Pnt pDef;
-      Standard_Real aSqDef = -1.;
-      Standard_Boolean isSkipped = Standard_False;
-      gp_XYZ normal(aLinkVec[0] ^ aLinkVec[1]);
-      if (normal.SquareModulus () < gp::Resolution())
-        continue;
-
-      normal.Normalize();
-
-      // Check deflection at triangle centroid
-      gp_XY aCenter2d = (xy[0] + xy[1] + xy[2]) / 3.0;
-      gFace->D0(aCenter2d.X(), aCenter2d.Y(), pDef);
-      aSqDef = Abs(normal * (pDef.XYZ() - p[0]));
-      aSqDef *= aSqDef;
-
-      isSkipped = !checkDeflectionAndInsert(pDef, aCenter2d, theIsFirst, 
-        aSqDef, aSqDefFace, aCircles, aNewVertices, aMaxSqDef, aTempAlloc);
-
-      if (isSkipped)
-        break;
-
-      // Check deflection at triangle links
-      for (Standard_Integer i = 0; i < 3 && !isSkipped; ++i)
-      {
-        if (m[i]) // is a boundary
-          continue;
-
-        // Check if this link was already processed
-        if (aCouples.Add(myStructure->GetLink(e[i])))
-        {
-          // Check deflection on edge 1
-          Standard_Integer j = (i + 1) % 3;
-          gp_XY mi2d = (xy[i] + xy[j]) * 0.5;
-          gFace->D0(mi2d.X(), mi2d.Y(), pDef);
-          gp_Lin aLin(p[i], gp_Vec(p[i], p[j]));
-          aSqDef = aLin.SquareDistance(pDef);
-
-          isSkipped = !checkDeflectionAndInsert(pDef, mi2d, theIsFirst, 
-            aSqDef, aSqDefFace, aCircles, aNewVertices, aMaxSqDef, aTempAlloc);
-        }
-      }
-
-      if (isSkipped)
-        break;
-
-      //check normal on bsplines
-      if (theIsFirst && !aBSpline.IsNull())
-      {
-        gp_Dir N[3] = { gp::DZ(), gp::DZ(), gp::DZ() };
-        Standard_Integer aSt[3];
-
-        for (Standard_Integer i = 0; i < 3; ++i)
-        {
-          if (aNorMap.IsBound(v[i]))
-          {
-            aSt[i] = aStatMap.Find(v[i]);
-            N[i] = aNorMap.Find(v[i]);
-          }
-          else
-          {
-            aSt[i] = GeomLib::NormEstim(aBSpline, gp_Pnt2d(xy[i]), Precision::Confusion(), N[i]);
-            aStatMap.Bind(v[i], aSt[i]);
-            aNorMap.Bind(v[i], N[i]);
-          }
-        }
-
-        Standard_Real aAngle[3];
-        for (Standard_Integer i = 0; i < 3; ++i)
-          aAngle[i] = N[(i + 1) % 3].Angle(N[i]);
-
-        if (aSt[0] < 1 && aSt[1] < 1 && aSt[2] < 1)
-        {
-          if (aAngle[0] > myAngle || aAngle[1] > myAngle || aAngle[2] > myAngle)
-          {
-            aMaxSqDef = -1.;
-            break;
-          }
-        }
-      }
-    }
-
-    if (theIsFirst)
-      continue;
-
-#ifdef DEBUG_MESH
-    // append to the file triangles in the form of polyline commands;
-    // so the file can be sourced in draw to analyze triangles on each pass of the algorithm.
-    // write triangles
-    ofstream ftt(DEBUG_MESH, std::ios_base::app);
-    Standard_Integer nbElem = myStructure->NbElements();
-    for (Standard_Integer i = 1; i <= nbElem; i++)
-    {
-      const BRepMesh_Triangle& aTri = myStructure->GetElement(i);
-      if (aTri.Movability() == BRepMesh_Deleted)
-        continue;
-      Standard_Integer n[3];
-      myStructure->ElementNodes(aTri, n);
-      const BRepMesh_Vertex& aV1 = myStructure->GetNode(n[0]);
-      const BRepMesh_Vertex& aV2 = myStructure->GetNode(n[1]);
-      const BRepMesh_Vertex& aV3 = myStructure->GetNode(n[2]);
-      const gp_Pnt& aP1 = myAttribute->GetPoint(aV1);
-      const gp_Pnt& aP2 = myAttribute->GetPoint(aV2);
-      const gp_Pnt& aP3 = myAttribute->GetPoint(aV3);
-      ftt << "polyline t" << aPass << "_" << i << " "
-        << aP1.X() << " " << aP1.Y() << " " << aP1.Z() << " "
-        << aP2.X() << " " << aP2.Y() << " " << aP2.Z() << " "
-        << aP3.X() << " " << aP3.Y() << " " << aP3.Z() << " "
-        << aP1.X() << " " << aP1.Y() << " " << aP1.Z() << endl;
-    }
-    // write points to insert on the current pass
-    BRepMesh::ListOfVertex::Iterator it(aNewVertices);
-    for (Standard_Integer i=1; it.More(); it.Next(), i++)
-    {
-      const BRepMesh_Vertex& aVertex = it.Value();
-      const gp_Pnt& aP = myAttribute->GetPoint(aVertex);
-      ftt << "vertex vt" << aPass << "_" << i << " "
-        << aP.X() << " " << aP.Y() << " " << aP.Z() << endl;
-    }
-#endif
-
-    if (addVerticesToMesh(aNewVertices, theTrigu))
-      ++aInsertedNb;
-  }
-
-  return (aMaxSqDef < 0) ? aMaxSqDef : Sqrt(aMaxSqDef);
-}
-
-//=======================================================================
-//function : add
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscretFace::add(const TopoDS_Vertex& theVertex)
-{
-  if (theVertex.Orientation() != TopAbs_INTERNAL)
-    return;
-
-  try
-  {
-    OCC_CATCH_SIGNALS
-
-    gp_Pnt2d aPnt2d = BRep_Tool::Parameters(theVertex, myAttribute->Face());
-    // check UV values for internal vertices
-    if (myAttribute->ChangeClassifier()->Perform(aPnt2d) != TopAbs_IN)
-      return;
-
-    NCollection_Handle<FixedVExplorer> aFixedVExplorer = new FixedVExplorer(theVertex);
-    Standard_Integer aIndex = myAttribute->GetVertexIndex(aFixedVExplorer);
-    gp_XY anUV = BRepMesh_ShapeTool::FindUV(aIndex, aPnt2d,
-      BRep_Tool::Tolerance(theVertex), myAttribute);
-
-    Standard_Integer aTmpId1, aTmpId2;
-    anUV = myAttribute->Scale(anUV, Standard_True);
-    myAttribute->AddNode(aIndex, anUV, BRepMesh_Fixed, aTmpId1, aTmpId2);
-  }
-  catch (Standard_Failure)
-  {
-  }
-}
-
-//=======================================================================
-//function : insertVertex
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscretFace::insertVertex(
-  const gp_Pnt&           thePnt3d,
-  const gp_XY&            theUV,
-  BRepMesh::ListOfVertex& theVertices)
-{
-  Standard_Integer aNbLocat = myAttribute->LastPointId();
-  myAttribute->ChangeSurfacePoints()->Bind(++aNbLocat, thePnt3d);
-
-  gp_XY aPnt2d  = myAttribute->Scale(theUV, Standard_True);
-  BRepMesh_Vertex aVertex(aPnt2d, aNbLocat, BRepMesh_Free);
-  theVertices.Append(aVertex);
-}
-
-//=======================================================================
-//function : commitSurfaceTriangulation
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscretFace::commitSurfaceTriangulation()
-{
-  if (myAttribute.IsNull() || !myAttribute->IsValid())
-    return;
-
-  const TopoDS_Face& aFace = myAttribute->Face();
-  BRepMesh_ShapeTool::NullifyFace(aFace);
-
-  Handle(BRepMesh_DataStructureOfDelaun)& aStructure = myAttribute->ChangeStructure();
-  const BRepMesh::MapOfInteger&           aTriangles = aStructure->ElementsOfDomain();
-
-  if (aTriangles.IsEmpty())
-  {
-    myAttribute->SetStatus(BRepMesh_Failure);
-    return;
-  }
-
-  BRepMesh::HIMapOfInteger& aVetrexEdgeMap = myAttribute->ChangeVertexEdgeMap();
-
-  // Store triangles
-  Standard_Integer aVerticesNb  = aVetrexEdgeMap->Extent();
-  Standard_Integer aTrianglesNb = aTriangles.Extent();
-  Handle(Poly_Triangulation) aNewTriangulation =
-    new Poly_Triangulation(aVerticesNb, aTrianglesNb, Standard_True);
-
-  Poly_Array1OfTriangle& aPolyTrianges = aNewTriangulation->ChangeTriangles();
-
-  Standard_Integer aTriangeId = 1;
-  BRepMesh::MapOfInteger::Iterator aTriIt(aTriangles);
-  for (; aTriIt.More(); aTriIt.Next())
-  {
-    const BRepMesh_Triangle& aCurElem = aStructure->GetElement(aTriIt.Key());
-
-    Standard_Integer aNode[3];
-    aStructure->ElementNodes(aCurElem, aNode);
-
-    Standard_Integer aNodeId[3];
-    for (Standard_Integer i = 0; i < 3; ++i)
-      aNodeId[i] = aVetrexEdgeMap->FindIndex(aNode[i]);
-
-    aPolyTrianges(aTriangeId++).Set(aNodeId[0], aNodeId[1], aNodeId[2]);
-  }
-
-  // Store mesh nodes
-  TColgp_Array1OfPnt&   aNodes   = aNewTriangulation->ChangeNodes();
-  TColgp_Array1OfPnt2d& aNodes2d = aNewTriangulation->ChangeUVNodes();
-
-  for (Standard_Integer i = 1; i <= aVerticesNb; ++i)
-  {
-    Standard_Integer       aVertexId = aVetrexEdgeMap->FindKey(i);
-    const BRepMesh_Vertex& aVertex   = aStructure->GetNode(aVertexId);
-    const gp_Pnt&          aPoint    = myAttribute->GetPoint(aVertex);
-
-    aNodes(i)   = aPoint;
-    aNodes2d(i) = aVertex.Coord();
-  }
-
-  aNewTriangulation->Deflection(myAttribute->GetDefFace());
-  BRepMesh_ShapeTool::AddInFace(aFace, aNewTriangulation);
-
-  // Delete unused data
-  myUParam.Clear(0L);
-  myVParam.Clear(0L);
-  myAttribute->ChangeStructure().Nullify();
-  myAttribute->ChangeSurfacePoints().Nullify();
-  myAttribute->ChangeSurfaceVertices().Nullify();
-
-  myAttribute->ChangeClassifier().Nullify();
-  myAttribute->ChangeLocation2D().Nullify();
-  myAttribute->ChangeVertexEdgeMap().Nullify();
-}