Logo Search packages:      
Sourcecode: opencascade version File versions  Download package

BRepMesh_FastDiscret.cxx

// File:    BRepMesh_FastDiscret.cxx
// Created: Tue Feb 27 16:39:53 1996
// Author:  Ekaterina SMIRNOVA
// Copyright: Open CASCADE SAS 2008


#include <BRepMesh_FastDiscret.ixx>

#include <BRepMesh_Array1OfVertexOfDelaun.hxx>
#include <BRepMesh_ListIteratorOfListOfVertex.hxx>
#include <BRepMesh_GeomTool.hxx>
#include <MeshShape_ListIteratorOfListOfSurfacePoint.hxx>
#include <MeshShape_CurvatureType.hxx>
#include <MeshShape_Couple.hxx>
#include <GeomAbs_IsoType.hxx>
#include <GeomAbs_SurfaceType.hxx>
#include <TopAbs.hxx>
//#include <TColStd_ListOfReal.hxx>
#include <TColStd_Array1OfBoolean.hxx>
#include <TColStd_HArray1OfReal.hxx>
#include <TColStd_ListIteratorOfListOfReal.hxx>
#include <TColStd_MapIteratorOfMapOfInteger.hxx>
#include <Precision.hxx>
#include <TColgp_Array1OfXY.hxx>

#include <BRep_Builder.hxx>
#include <BRep_Tool.hxx>
#include <Poly_Triangulation.hxx>
#include <Poly_PolygonOnTriangulation.hxx>
#include <Poly_Connect.hxx>
#include <TColStd_SequenceOfInteger.hxx>
#include <TColStd_Array1OfInteger.hxx>
#include <TColStd_HArray1OfInteger.hxx>

#include <TColgp_Array1OfPnt.hxx>
#include <TColgp_Array1OfPnt2d.hxx>
#include <Precision.hxx>

#include <BRepAdaptor_Curve.hxx>
#include <BRepAdaptor_Surface.hxx>
#include <BRepAdaptor_HSurface.hxx>
#include <BRepTools.hxx>
#include <BndLib_Add3dCurve.hxx>
#include <BRepBndLib.hxx>
#include <Bnd_Box.hxx>
#include <TopoDS.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>

#include <Geom2d_Curve.hxx>

#include <MeshShape_DataMapOfShapeListOfTransient.hxx>
#include <TColStd_ListOfTransient.hxx>
#include <MeshShape_DataMapIteratorOfDataMapOfShapeListOfTransient.hxx>
#include <TColStd_DataMapOfIntegerInteger.hxx>
#include <BRepMesh_ShapeTool.hxx>
#include <ElSLib.hxx>
#include <Geom_Surface.hxx>
#include <Adaptor3d_IsoCurve.hxx>
#include <BRepMesh_IndexedMapOfVertex.hxx>
#include <BRepMesh_Delaun.hxx>
#include <Extrema_LocateExtPC.hxx>

#include <BRepMesh_ListOfXY.hxx>
#include <BRepMesh_ListIteratorOfListOfXY.hxx>

#include <TColStd_Array1OfInteger.hxx>
#include <BRepMesh_IDMapOfNodeOfDataStructureOfDelaun.hxx>
#include <Standard_ErrorHandler.hxx>
#include <Standard_Failure.hxx>
#include <MeshDS_DataMapOfInteger.hxx>
#include <TColGeom2d_SequenceOfCurve.hxx>
#include <TopTools_SequenceOfShape.hxx>
#include <NCollection_IncAllocator.hxx>

#include <BRep_ListIteratorOfListOfPointRepresentation.hxx>
#include <BRep_PointRepresentation.hxx>
#include <BRep_TVertex.hxx>
#include <TColStd_MapOfInteger.hxx>
#include <SortTools_ShellSortOfReal.hxx>
#include <TCollection_CompareOfReal.hxx>

static Standard_Integer LIMITE_TRIANGULATION = 1048576;


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.);
}
//--ofv.

inline Standard_Real MaxFaceTol (const TopoDS_Face& theFace)
{
  Standard_Real T, TMax = BRep_Tool::Tolerance(theFace);
  TopExp_Explorer Ex;

  for (Ex.Init(theFace,TopAbs_EDGE); Ex.More(); Ex.Next())
  {
    T = BRep_Tool::Tolerance(TopoDS::Edge(Ex.Current()));
    if (T > TMax) TMax = T;
  }

  for (Ex.Init(theFace,TopAbs_VERTEX); Ex.More(); Ex.Next())
  {
    T = BRep_Tool::Tolerance(TopoDS::Vertex(Ex.Current()));
    if (T > TMax) TMax = T;
  }

  return TMax;
}


//=======================================================================
//function : BRepMesh_FastDiscret
//purpose  : 
//=======================================================================
BRepMesh_FastDiscret::BRepMesh_FastDiscret(const Standard_Real    defle,
                                   const Standard_Real    angl,
                                   const Bnd_Box&         Box,
                           const Standard_Boolean ws,
                           const Standard_Boolean inshape,
                           const Standard_Boolean relative,
                           const Standard_Boolean shapetrigu): 
  angle(angl), deflection(defle),WithShare(ws),nbDomains(0),  nbLocat(0), 
  myrelative(relative), myshapetrigu(shapetrigu), myinshape(inshape),
  myInternalVerticesMode(Standard_True)
{
  myAllocator = new NCollection_IncAllocator(64000);
  if (relative)
  {
    Standard_Real TXmin, TYmin, TZmin, TXmax, TYmax, TZmax;
    Box.Get(TXmin, TYmin, TZmin, TXmax, TYmax, TZmax);
    mydtotale = TXmax-TXmin;
    const Standard_Real dy = TYmax-TYmin;
    const Standard_Real dz = TZmax-TZmin;
    if (dy > mydtotale) mydtotale = dy;
    if (dz > mydtotale) mydtotale = dz;
  }
}


//=======================================================================
//function : BRepMesh_FastDiscret
//purpose  : 
//=======================================================================

00173 BRepMesh_FastDiscret::BRepMesh_FastDiscret(const Standard_Real    defle,
                           const TopoDS_Shape&    shape,
                                   const Bnd_Box&         Box,
                           const Standard_Real    angl,
                           const Standard_Boolean ws,
                           const Standard_Boolean inshape,
                           const Standard_Boolean relative,
                           const Standard_Boolean shapetrigu): 
 angle(angl), deflection(defle),WithShare(ws),nbDomains(0),  nbLocat(0), 
 myrelative(relative), myshapetrigu(shapetrigu), myinshape(inshape),
 myInternalVerticesMode(Standard_True)
{
  myAllocator = new NCollection_IncAllocator(64000);
  if (relative)
  {
    Standard_Real TXmin, TYmin, TZmin, TXmax, TYmax, TZmax;
    Box.Get(TXmin, TYmin, TZmin, TXmax, TYmax, TZmax);
    mydtotale = TXmax-TXmin;
    const Standard_Real dy = TYmax-TYmin;
    const Standard_Real dz = TZmax-TZmin;
    if (dy > mydtotale) mydtotale = dy;
    if (dz > mydtotale) mydtotale = dz;
  }
  Add(shape);
}


//=======================================================================
//function : Add(shape)
//purpose  : 
//=======================================================================

void BRepMesh_FastDiscret::Add(const TopoDS_Shape& shape)
{
  TopExp_Explorer ex;
  for (ex.Init(shape, TopAbs_FACE); ex.More(); ex.Next()) {
    Add(TopoDS::Face(ex.Current()));
  }
}


//=======================================================================
//function : Add(face)
//purpose  : 
//=======================================================================

#define MESH_FAILURE(theface)        \
{                                    \
  myfacestate = BRepMesh_Failure;    \
  mynottriangulated.Append(theface); \
  return;                            \
}

void BRepMesh_FastDiscret::Add(const TopoDS_Face& theface)
{
#ifndef DEB_MESH
  try
  {
    OCC_CATCH_SIGNALS
#endif
  TopoDS_Face face = theface;
  BRepTools::Update(face);
  face.Orientation(TopAbs_FORWARD);
  nbDomains++;
  Domains.Bind(nbDomains, face);
  structure.Nullify();
  Handle(NCollection_IncAllocator) anAlloc = Handle(NCollection_IncAllocator)::DownCast(myAllocator);
  anAlloc->Reset(Standard_False);  
  structure=new BRepMesh_DataStructureOfDelaun(anAlloc);
  structure->NewDomain(nbDomains);
  BRepAdaptor_Surface  BS(face, Standard_False);
  Handle(BRepAdaptor_HSurface) gFace = new BRepAdaptor_HSurface(BS);
  
  GeomAbs_SurfaceType thetype;
  thetype = BS.GetType();

  gp_Pnt2d uvFirst, uvLast;

  TopAbs_Orientation orFace = face.Orientation();
  Handle(Poly_Triangulation) T;
  TopLoc_Location l;

  if (!WithShare) {          
    vertices.Clear();
    edges.Clear();
  }

  mylistver.Clear();
  myvemap.Clear();
  mylocation2d.Clear();
  internaledges.Clear();

  Standard_Integer i, j;//,i1, i2, k;
  i = 1;

  Standard_Real defedge, defface;
  Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, dx, dy, dz;
  Standard_Integer nbEdge = 0;
  Standard_Real savangle = angle;
  Standard_Real cdef;
  Standard_Real maxdef = 2.* MaxFaceTol(theface);
  defface = 0.;

  if (!myrelative) defface = Max(deflection, maxdef);

  TColStd_SequenceOfReal aFSeq, aLSeq;
  TColGeom2d_SequenceOfCurve aCSeq;
  TopTools_SequenceOfShape aShSeq;

  TopoDS_Iterator exW(face);
  for (; exW.More(); exW.Next()) {
    const TopoDS_Shape& aWire = exW.Value();
    if (aWire.ShapeType() != TopAbs_WIRE)
      continue;
    TopoDS_Iterator ex(aWire);
    for(; ex.More(); ex.Next()) {
      const TopoDS_Edge& edge = TopoDS::Edge(ex.Value());
      nbEdge++;
      if (!mapdefle.IsBound(edge)) {
        if (myrelative) {
          if (edges.IsBound(edge)) {
          const TColStd_ListOfTransient& L = edges.Find(edge);
          const Handle(Poly_PolygonOnTriangulation)& P = 
            *(Handle(Poly_PolygonOnTriangulation)*)&(L.First());
          defedge = P->Deflection();
        }
          else {
            Bnd_Box B;
            BRepBndLib::Add(edge, B);
            B.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
            dx = aXmax-aXmin;
            dy = aYmax-aYmin;
            dz = aZmax-aZmin;
            defedge = dx;
            if (defedge < dy) defedge = dy;
            if (defedge < dz) defedge = dz;
            // ajustement par rapport a la taille totale:
            cdef = mydtotale/(2*defedge);
            if (cdef < 0.5) cdef = 0.5;
            if (cdef > 2.) cdef = 2.;
            defedge = cdef * defedge * deflection;
            angle = savangle * cdef;
          }
          defface = defface + defedge;
        defface = Max(maxdef, defface);
        }
        else defedge = deflection;
    
        defedge = Max(maxdef, defedge);
        defedge = Max(1.e-05 , defedge);
        mapdefle.Bind(edge, defedge);
      }
      else{
        defedge = mapdefle(edge);
        if (myrelative) {defface = defface + defedge; defface = Max(maxdef, defface);}
      }
      Standard_Real f,l;
      Handle(Geom2d_Curve) C = BRep_Tool::CurveOnSurface(edge, face, f, l);
      if (C.IsNull()) continue;

      aFSeq.Append(f);
      aLSeq.Append(l);
      aCSeq.Append(C);
      aShSeq.Append(edge);
      Add(edge, face, gFace, C, defedge, f, l);
      angle = savangle;
    }
  }

  if (nbEdge == 0 || myvemap.Extent() < 3)
  {
      MESH_FAILURE(theface);
  }

  if (myrelative ) defface = defface / nbEdge;
  else             defface = deflection;

  if (WithShare) defface = Max(maxdef, defface);
  
  T = BRep_Tool::Triangulation(face, l);

  if (!myshapetrigu || T.IsNull()) {
    
    Standard_Real xCur, yCur;
    Standard_Real maxX, minX, maxY, minY;
    minX=minY=1.e100;
    maxX=maxY=-1.e100;
    
    Standard_Integer ipn = 0;
    Standard_Integer i =1;
    for (i = 1; i <= myvemap.Extent(); i++) {
      const BRepMesh_Vertex& aV = structure->GetNode(myvemap.FindKey(i));
      ipn++;
      xCur=aV.Coord().X();
      yCur=aV.Coord().Y();
      minX=Min(xCur, minX);
      maxX=Max(xCur, maxX);
      minY=Min(yCur, minY);
      maxY=Max(yCur, maxY);
    }
    myumin = minX;
    myumax = maxX;
    myvmin = minY;
    myvmax = maxY;
    
    const Standard_Real umin = BS.FirstUParameter();
    const Standard_Real umax = BS.LastUParameter();
    const Standard_Real vmin = BS.FirstVParameter();
    const Standard_Real vmax = BS.LastVParameter();

    if (myumin < umin || myumax > umax)
    {
      if (BS.IsUPeriodic())
      {
        if ((myumax - myumin) > (umax - umin))
        {
          myumax = myumin + (umax - umin);
        }
      }
      else
      {
        if (umin > myumin) myumin = umin;
        if (umax < myumax) myumax = umax;
      }
    }

    if (myvmin < vmin || myvmax > vmax)
    {
      if (BS.IsVPeriodic())
      {
        if ((myvmax - myvmin) > (vmax - vmin))
        {
          myvmax = myvmin + (vmax - vmin);
        }
      }
      else
      {
        if (vmin > myvmin) myvmin = vmin;
        if (vmax < myvmax) myvmax = vmax;
      }
    }

    // verif rapide de la validite des bornes calculees. Si fausse,
    // surement un pb de pcurve.

    if (thetype == GeomAbs_BezierSurface &&
        (myumin < -0.5 || myumax > 1.5 || myvmin < -0.5 || myvmax > 1.5))
    {
      MESH_FAILURE(theface);
    }
 
    //define parameters for correct parametrics
    
    deltaX = 1.0;
    deltaY = 1.0;
    Standard_Integer nbVertices = myvemap.Extent();
    const Standard_Real tolclass = Precision::PConfusion(); //0.03*Max(myumax-myumin, myvmax-myvmin);
    
    BRepMesh_Classifier* classifier = 
      new BRepMesh_Classifier(face, tolclass,   internaledges, myvemap, 
                              structure, myumin, myumax, myvmin, myvmax);   
  
    myfacestate = classifier->State();
    if (myfacestate == BRepMesh_SelfIntersectingWire)
    {
      Standard_Integer nbmaill = 0;
      Standard_Real eps = Precision::Confusion();
      while (nbmaill < 5 && myfacestate != BRepMesh_ReMesh)
      {
        nbmaill++;
        
        //clear the structure of links
        const MeshDS_MapOfInteger& aLinks = structure->LinkOfDomain(nbDomains);
        for(MeshDS_MapOfInteger::Iterator anIter(aLinks);anIter.More();anIter.Next())
          structure->ForseRemoveLink(anIter.Key());
        for (i = 1; i <= myvemap.Extent(); i++)
           structure->ForseRemoveNode(myvemap.FindKey(i));
        
        structure->ClearDeleted();
        
        mylistver.Clear();
        myvemap.Clear();
        mylocation2d.Clear();
        internaledges.Clear();

        Standard_Integer j;
        for(j = 1; j <= aShSeq.Length(); j++)
        {
          const TopoDS_Edge& edge = TopoDS::Edge(aShSeq.Value(j));
          if (edges.IsBound(edge))
          {
            edges.UnBind(edge);
            internaledges.UnBind(edge);
          }
        }
        
        
        for( j = 1; j <= aShSeq.Length(); j++)
        {
          const TopoDS_Edge& edge = TopoDS::Edge(aShSeq.Value(j));
          defedge = mapdefle(edge) / 3.;
          defedge = Max(defedge, eps);
          mapdefle.Bind(edge, defedge);
          const Handle(Geom2d_Curve)& C = aCSeq.Value(j);
          Add(edge, face, gFace, C, defedge, aFSeq.Value(j), aLSeq.Value(j));
        }
        
        delete classifier;

        classifier = new BRepMesh_Classifier(face, tolclass, internaledges, myvemap,
                                             structure, myumin, myumax, myvmin, myvmax);

        if (classifier->State() == BRepMesh_NoError)
        {
          myfacestate = BRepMesh_ReMesh;
      }
        nbVertices = myvemap.Extent();
      }
    }
    
    if (myfacestate != BRepMesh_NoError && myfacestate != BRepMesh_ReMesh)
    {
      mynottriangulated.Append(face);
      delete classifier;
      return;
    }
    
    //add internal vertices after self-intersection check
    if(myInternalVerticesMode) {
      for(TopExp_Explorer ex(face,TopAbs_VERTEX ,TopAbs_EDGE); ex.More(); ex.Next()) {
        const TopoDS_Vertex& aVert = TopoDS::Vertex(ex.Current());
        Add(aVert,face,gFace);
      }
      nbVertices = myvemap.Extent();
    }
    

    
    // 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;
    
    Standard_Real du = (myumax-myumin)/20;
    Standard_Real dv = (myvmax-myvmin)/20;
    Standard_Real dfuave=(myumin+myumax)/2, dfucur;
    Standard_Real dfvave=(myvmin+myvmax)/2, dfvcur;
    // U loop
    BS.D0 (myumin, myvmin, P11);
    BS.D0 (myumin, dfvave, P21);
    BS.D0 (myumin, myvmax, P31);
    for (i=0, dfucur=myumin; i <= 20; i++, dfucur+=du)  {
      BS.D0 (dfucur, myvmin, P12);
      BS.D0 (dfucur, dfvave, P22);
      BS.D0 (dfucur, myvmax, P32);
      longu += ( P11.Distance(P12) + P21.Distance(P22) + P31.Distance(P32) );
      P11 = P12;
      P21 = P22;
      P31 = P32;
    }
    // V loop
    BS.D0(myumin, myvmin, P11);
    BS.D0(dfuave, myvmin, P21);
    BS.D0(myumax, myvmin, P31);
    for (i=0, dfvcur=myvmin; i <= 20; i++, dfvcur+=dv) {
      BS.D0 (myumin, dfvcur, P12);
      BS.D0 (dfuave, dfvcur, P22);
      BS.D0 (myumax, dfvcur, P32);
      longv += ( P11.Distance(P12) + P21.Distance(P22) + P31.Distance(P32) );
      P11 = P12;
      P21 = P22;
      P31 = P32;
    }
    longu /= 3.;
    longv /= 3.;
    // akm (bug OCC16) ^^^^^
    
    if (longu <= 1.e-16 || longv <= 1.e-16) {
      //sisi, ca s`est vu!!
#ifdef DEB_MESH_CHRONO
      chMaillEdges.Stop();
      MESH_CHRONO_TSTOP(thetype);
      cout <<"WARNING: la face "<<myNbDomains<<" n'est pas triangulee."<<endl;
#endif
//    delete classifier;
//    MESH_FAILURE(theface);
    }


    if (thetype == GeomAbs_Torus)  {
      gp_Torus Tor = BS.Torus();
      Standard_Real r = Tor.MinorRadius(), R = Tor.MajorRadius();
      Standard_Real Du, Dv;//, pasu, pasv;
      
      Dv = Max(1.0e0 - (defface/r),0.0e0) ;
      Standard_Real oldDv = 2.0 * ACos (Dv);
      oldDv = Min(oldDv, angle);
      Dv  =  0.9*oldDv; //TWOTHIRD * oldDv;
      Dv = oldDv;
      
      Standard_Integer nbV = Max((Standard_Integer)((myvmax-myvmin)/Dv), 2);
      Dv = (myvmax-myvmin)/(nbV+1);
      
      Standard_Real ru = R + r;
      if (ru > 1.e-16) {
        Du = Max(1.0e0 - (defface/ru),0.0e0);
        Du  = (2.0 * ACos (Du));
        Du = Min(Du, angle);
        Standard_Real aa = sqrt(Du*Du + oldDv*oldDv);
        if(aa < gp::Resolution())
          return; 

        Du = Du * Min(oldDv, Du) / aa;
      }
      else Du = Dv;     
      
      Standard_Integer nbU = Max((Standard_Integer)((myumax-myumin)/Du), 2);
      nbU = Max(nbU, (Standard_Integer)(nbV*(myumax-myumin)*R/((myvmax-myvmin)*r)/5.));
      
      Du = (myumax-myumin)/(nbU+1);
      //-- On choisit DeltaX et DeltaY de facon a ce qu on ne saute pas 
      //-- de points sur la grille
      deltaX = Du;
      deltaY = Dv;
    }
    else if (thetype == GeomAbs_Cylinder) {
      Standard_Real aMax = Max(myumax,myvmax);
      deltaX = 0.01;  //1./aMax; //0.01;
      deltaY = 1.0;
    }
    else {
      deltaX = (myumax-myumin)/longu;
      deltaY = (myvmax-myvmin)/longv;
    }    


    myUParam.Clear(); myVParam.Clear();
        
    TColStd_Array1OfInteger tabvert_corr(1, nbVertices);
    
    gp_Pnt2d p2d;
    
    // Check the necessity to fill the map of parameters
    const Standard_Boolean useUVParam = (thetype == GeomAbs_Torus ||
                                         thetype == GeomAbs_BezierSurface ||
                                         thetype == GeomAbs_BSplineSurface);
    
    myUParam.Clear(); myVParam.Clear();
  
    BRepMesh_IDMapOfNodeOfDataStructureOfDelaun aMoveNodes(myvemap.Extent());
    //cout << "UVparam" << endl;
    Standard_Real eps = Precision::PConfusion();
    Standard_Integer nbu1 = 0, nbu2 = 0, nbv1 = 0, nbv2 = 0, nbtot = 0;
    Standard_Boolean IsRectangle = Standard_False;
    TColStd_ListOfReal U1Params, U2Params, V1Params, V2Params;

    for (i = 1; i <= structure->NbNodes(); i++)
    {
      const BRepMesh_Vertex& v = structure->GetNode(i);
      p2d = v.Coord();
      if (useUVParam) {
      //cout << p2d.X() << " " << p2d.Y() << endl;
      Standard_Real u = p2d.X(), v = p2d.Y();
        myUParam.Add(u);
        myVParam.Add(v);

      ++nbtot;
      if(Abs(u - myumin) <= eps) {
        ++nbu1;
        V1Params.Append(v);
      }
      if(Abs(u - myumax) <= eps) {
        ++nbu2;
        V2Params.Append(v);
      }
      if(Abs(v - myvmin) <= eps) {
        ++nbv1;
        U1Params.Append(u);
      }
      if(Abs(v - myvmax) <= eps) {
        ++nbv2;
        U2Params.Append(u);
      }
      }
      gp_XY res;
      res.SetCoord((p2d.X()-minX)/deltaX,(p2d.Y()-minY)/deltaY);
      BRepMesh_Vertex v_new(res,v.Domain(),v.Location3d(),v.Movability());
      const MeshDS_ListOfInteger& alist = structure->GetNodeList(i);
      aMoveNodes.Add(v_new,alist);
      tabvert_corr(i) = i;
    }    
    structure->ReplaceNodes(aMoveNodes);
    
    if(useUVParam) {
      Standard_Integer nb = nbu1 + nbu2 + nbv1 + nbv2 - 4;
      if((nbu1 == nbu2) && (nbv1 == nbv2) && (nbtot == nb)) IsRectangle = Standard_True;
    }


    Standard_Boolean rajout;

     switch (thetype)
     {
     case GeomAbs_Plane:
       rajout = !classifier->NaturalRestriction();
       break;
     case GeomAbs_Sphere:
     case GeomAbs_Torus:
       rajout = Standard_True;
       break;
     default:
       rajout = Standard_False;
     }  

    BRepMesh_Delaun trigu( structure,tabvert_corr, orFace==TopAbs_FORWARD);

    Standard_Boolean isaline;
    isaline = ((myumax-myumin)<1.e-05) || ((myvmax-myvmin)<1.e-05);
    
    Standard_Real aDef = -1;
    if (!isaline && structure->ElemOfDomain(nbDomains).Extent() > 0) {
      TColStd_ListOfInteger badTri, nulTri;
      
      if(!rajout)
      {
        aDef = Control(gFace, defface, mylistver, badTri, nulTri, trigu, Standard_True);
        if( aDef > defface || aDef < 0.)
          rajout = Standard_True;
      }

      if(!rajout) {
      if(useUVParam) {
        if(BS.IsUClosed()) {
          if(myVParam.Extent() > 2) {
            rajout = Standard_True;
          }
        }
        if(BS.IsVClosed()) {
          if(myUParam.Extent() > 2) {
            rajout = Standard_True;
          }
        }
      }
      }

      if(rajout){
      if(IsRectangle) {
        InternalVerticesForRectangle(gFace, mylistver, U1Params, U2Params, V1Params, V2Params);
      }
      else {
        InternalVertices(gFace, mylistver, defface,*classifier);
      }

        if (mylistver.Extent() > 0) {
        BRepMesh_Array1OfVertexOfDelaun verttab(1, mylistver.Extent());
        BRepMesh_ListIteratorOfListOfVertex itVer(mylistver);
        ipn = 1;
        for (; itVer.More(); itVer.Next())
          verttab(ipn++) = itVer.Value();
        trigu.AddVertices(verttab);
      }
        //control internal points
        BRepMesh_ListOfVertex vvlist;
        aDef = Control(gFace, defface, vvlist, badTri, nulTri, trigu, Standard_False);
        mylistver.Append(vvlist);
      }
    }

    //modify structure back
    aMoveNodes.Clear();
    for (i = 1; i <= structure->NbNodes(); i++)
    {
      const BRepMesh_Vertex& v = structure->GetNode(i);
      p2d = v.Coord();
      gp_XY res;
      res.SetCoord(p2d.X()*deltaX+myumin,p2d.Y()*deltaY+myvmin);
      BRepMesh_Vertex v_new(res,v.Domain(),v.Location3d(),v.Movability());
      const MeshDS_ListOfInteger& alist = structure->GetNodeList(i);
      aMoveNodes.Add(v_new,alist);
    }
    structure->ReplaceNodes(aMoveNodes);
    
    if (myinshape && T.IsNull()) AddInShape(face, (aDef < 0.0)? defface : aDef);
  }

  else {
    // recup des triangles pour les mettre dans la structure.
    const TColgp_Array1OfPnt& Nodes = T->Nodes();
    const TColgp_Array1OfPnt2d& UVNodes = T->UVNodes();
    const Poly_Array1OfTriangle& triangles = T->Triangles();
    Standard_Integer nt, ed1, ed2, ed3, tri;
#ifdef DEB
    Standard_Integer nbNodes = 
#endif
      Nodes.Length();
    Standard_Integer nbTriangles = T->NbTriangles();
    BRepMesh_Vertex v1, v2, v3;
    Standard_Integer iv1, iv2, iv3;
    Standard_Integer k, niver, nFree;
    Standard_Integer t[3], n[3];


    Poly_Connect pc(T);

    // count the free edges
    nFree = 0;
    TColStd_MapOfInteger vFree;
    for (i = 1; i <= nbTriangles; i++) {
      pc.Triangles(i,t[0],t[1],t[2]);
      triangles(i).Get(n[0], n[1], n[2]);
      for (j = 0; j < 3; j++) {
      if (t[j] == 0)  {
        k = (j+1) % 3;
        if (!vFree.Contains(n[j])) vFree.Add(n[j]);
        if (!vFree.Contains(n[k])) vFree.Add(n[k]);
        nFree++;
      }
      }
    }

    TColStd_DataMapOfIntegerInteger IDMAP;
    Standard_Integer index;


    // parmi ces noeuds exterieurs, trouver leur indice dans la structure.
    TColStd_MapIteratorOfMapOfInteger it;

    // iteration sur vemap;
    for (i = 1; i <= myvemap.Extent(); i++) {
      index = myvemap.FindKey(i);

      for (it.Initialize(vFree); it.More(); it.Next()) {
      const gp_Pnt2d& P2d = UVNodes(it.Key());
      if (Vertex(index).Coord().IsEqual(P2d.Coord(), Precision::PConfusion())) {
        IDMAP.Bind(it.Key(), index);
        vFree.Remove(it.Key());
        break;
      }
      }
    }
    
    gp_Pnt P1, P2, P3;

    for (nt = 1; nt <= nbTriangles; nt++) {
      triangles(nt).Get(n[0], n[1], n[2]);
      if (!l.IsIdentity()) {
      P1 = Nodes(n[0]).Transformed(l.Transformation());
      P2 = Nodes(n[1]).Transformed(l.Transformation());
      P3 = Nodes(n[2]).Transformed(l.Transformation());
      }
      else {
      P1 = Nodes(n[0]);
      P2 = Nodes(n[1]);
      P3 = Nodes(n[2]);
      }
      const gp_Pnt2d& P12d = UVNodes(n[0]);
      const gp_Pnt2d& P22d = UVNodes(n[1]);
      const gp_Pnt2d& P32d = UVNodes(n[2]);


      if (!IDMAP.IsBound(n[0])) {
      nbLocat++;
      v1.Initialize(P12d.XY(), nbDomains, nbLocat, MeshDS_OnSurface);
      Location3d.Bind(nbLocat, P1);
      iv1 = structure->AddNode(v1);
      IDMAP.Bind(n[0], iv1);
      }
      else {
      iv1 = IDMAP.Find(n[0]);
      }

      if (!IDMAP.IsBound(n[1])) {
      nbLocat++;
      v2.Initialize(P22d.XY(), nbDomains, nbLocat, MeshDS_OnSurface);
      Location3d.Bind(nbLocat, P2);
      iv2 = structure->AddNode(v2);
      IDMAP.Bind(n[1], iv2);
      }
      else {
      iv2 = IDMAP.Find(n[1]);
      }

      if (!IDMAP.IsBound(n[2])) {
      nbLocat++;
      v3.Initialize(P32d.XY(), nbDomains, nbLocat, MeshDS_OnSurface);
      Location3d.Bind(nbLocat, P3);
      iv3 = structure->AddNode(v3);
      IDMAP.Bind(n[2], iv3);
      }
      else {
      iv3 = IDMAP.Find(n[2]);
      }

      if (face.Orientation() == TopAbs_REVERSED) {
      niver=iv2;
      iv2=iv3;
      iv3=niver;
      }

      ed1=structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Free, nbDomains));
      ed2=structure->AddLink(BRepMesh_Edge(iv2,iv3,MeshDS_Free, nbDomains));
      ed3=structure->AddLink(BRepMesh_Edge(iv3,iv1,MeshDS_Free, nbDomains));
      tri = structure->AddElement(BRepMesh_Triangle(Abs(ed1), Abs(ed2), Abs(ed3), 
                                        (ed1>0), (ed2>0), (ed3>0), 
                                        MeshDS_Free, nbDomains));

    }
  }
#ifndef DEB_MESH
  }
  catch(Standard_Failure)
  {
    MESH_FAILURE(theface);
  }
#endif // DEB_MESH
}


//=======================================================================
//function : Add
//purpose  : 
//=======================================================================
void BRepMesh_FastDiscret::Add(const TopoDS_Edge&                  edge, 
                       const TopoDS_Face&                  face, 
                     const Handle(BRepAdaptor_HSurface)& gFace,
                           const Handle(Geom2d_Curve)&         C2d,
                     const Standard_Real                 defedge,
                           const Standard_Real                 first,
                           const Standard_Real                 last)
{
  const TopAbs_Orientation orEdge = edge.Orientation();
  if (orEdge == TopAbs_EXTERNAL) return;

  const Standard_Boolean isEdgeBound = edges.IsBound(edge);

  if (!isEdgeBound)
  {
    if (Update(edge, face, C2d, defedge, first, last))
    {
      return;
    }
  }

  TopoDS_Vertex pBegin, pEnd;
  TopExp::Vertices(edge, pBegin, pEnd);
  if (pBegin.IsNull() || pEnd.IsNull())
  {
    return;
  }

  Standard_Real wFirst, wLast;
  BRep_Tool::Range(edge, face, wFirst, wLast);

  gp_Pnt2d uvFirst, uvLast;
  BRep_Tool::UVPoints(edge, face, uvFirst, uvLast);

  const Standard_Boolean sameUV = uvFirst.IsEqual(uvLast, Precision::PConfusion());

  //Controle vertice tolerances
  gp_Pnt pFirst = gFace->Value(uvFirst.X(), uvFirst.Y());
  gp_Pnt pLast = gFace->Value(uvLast.X(), uvLast.Y());

  Standard_Real mindist = 10. * Max(pFirst.Distance(BRep_Tool::Pnt(pBegin)), 
                            pLast.Distance(BRep_Tool::Pnt(pEnd)));

  if(mindist < BRep_Tool::Tolerance(pBegin) ||
     mindist < BRep_Tool::Tolerance(pEnd) ) mindist = defedge;

  // controle des edges degenerees non codees.

  Standard_Boolean degener = BRep_Tool::Degenerated(edge);

  if (!degener)
  {
    if (pBegin.IsSame(pEnd))
    {
      // calcul de la longueur de l`edge en 3d
      Standard_Real longueur = 0.0;
      Standard_Real du = (wLast-wFirst)/20;
      gp_Pnt P1, P2;
      BRepAdaptor_Curve BC(edge);
      BC.D0(wFirst, P1);
      Standard_Real tolV = BRep_Tool::Tolerance(pBegin);
      Standard_Real tolV2 = 1.2*tolV;
      for (Standard_Integer l = 1; l <= 20; l++) {
      BC.D0(wFirst + l*du, P2);
      longueur += P1.Distance(P2);
      if (longueur > tolV2) break;
      P1 = P2;
      }

      if (longueur < tolV2)
      {
      degener = Standard_True;
      }
    }
  }

  // Correct UV points
  if (sameUV)
  {
    // 1. est-ce vraiment sameUV sans etre denegere
    gp_Pnt2d uvF, uvL;
    C2d->D0(first, uvF);
    C2d->D0(last, uvL);
    if (!uvFirst.IsEqual(uvF, Precision::PConfusion()))
    {
      uvFirst = uvF;
    }
    if (!uvLast.IsEqual(uvL, Precision::PConfusion()))
    {
      uvLast = uvL;
    }
  }

  gp_XY theUV;

  // Process first vertex
  Standard_Integer ipf;
  if (vertices.IsBound(pBegin))
  {
    ipf = vertices.Find(pBegin);
  }
  else
  {
    if (sameUV && vertices.IsBound(pEnd))
    {
      ipf = vertices.Find(pEnd);
    }
    else
    {
      nbLocat++;
      ipf = nbLocat;
      Location3d.Bind(ipf, BRep_Tool::Pnt(pBegin));
    }
    vertices.Bind(pBegin, ipf);
  }
  theUV = FindUV(pBegin, uvFirst, ipf, gFace, mindist);
  BRepMesh_Vertex vf(theUV, nbDomains, ipf, MeshDS_Frontier);
  Standard_Integer ivf = structure->AddNode(vf);

  // Process last vertex
  Standard_Integer ipl;
  if (pEnd.IsSame(pBegin))
  {
    ipl = ipf;
  }
  else
  {
    if (vertices.IsBound(pEnd))
    {
      ipl = vertices.Find(pEnd);
    }
    else
    {
      if (sameUV)
      {
        ipl = ipf;
      }
      else
      {
        nbLocat++;
            ipl = nbLocat;
        Location3d.Bind(ipl, BRep_Tool::Pnt(pEnd));
      }
      vertices.Bind(pEnd,ipl);
    }
  }
  theUV = FindUV(pEnd, uvLast, ipl, gFace, mindist);
  BRepMesh_Vertex vl(theUV, nbDomains, ipl, MeshDS_Frontier);
  Standard_Integer ivl= structure->AddNode(vl);

  Standard_Integer isvf = myvemap.FindIndex(ivf);
  if (isvf == 0) isvf = myvemap.Add(ivf);
  Standard_Integer isvl = myvemap.FindIndex(ivl);
  if (isvl == 0) isvl = myvemap.Add(ivl);
  
  Standard_Real otherdefedge = 0.5*defedge;
  gp_Pnt2d uv;
  BRepMesh_Vertex v2;
  gp_Pnt   P3d;

  Handle(Poly_PolygonOnTriangulation) P1;

  if (!isEdgeBound)
  {
    Handle(Poly_PolygonOnTriangulation) P2;

    if (degener)
    {
      // creation du nouveau:
      TColStd_Array1OfInteger Nodes(1, 2), NodInStruct(1, 2);
      TColStd_Array1OfReal Param(1, 2);
      
      NodInStruct(1) = ipf;
      Nodes(1) = isvf;
      Param(1) = wFirst;
      
      NodInStruct(2) = ipl;
      Nodes(2) = isvl;
      Param(2) = wLast;

      P1 = new Poly_PolygonOnTriangulation(Nodes, Param);
      P2 = new Poly_PolygonOnTriangulation(NodInStruct, Param);
    }
    else
    {
      if (orEdge == TopAbs_INTERNAL) otherdefedge *= 0.5;
      
      BRepAdaptor_Curve cons;
      if (BRep_Tool::SameParameter(edge))
      {
      cons.Initialize(edge);
      }
      else
      {
      cons.Initialize(edge, face);
      }

      TopLoc_Location L;
      Standard_Integer nbpmin = 2;
      if (cons.GetType() == GeomAbs_Circle) nbpmin = 4; //OCC287
      BRepMesh_GeomTool GT(cons, wFirst, wLast, 0.5*angle, otherdefedge, nbpmin);
      // Creation des polygones sur triangulation:
      const Standard_Integer nbnodes = GT.NbPoints();
      TColStd_Array1OfInteger Nodes(1, nbnodes);
      TColStd_Array1OfInteger NodInStruct(1, nbnodes);
      TColStd_Array1OfReal Param(1, nbnodes);
      
      // traitement du 1er point:
      Nodes(1) = isvf;
      NodInStruct(1) = ipf;
      Param(1) = wFirst;
      
      // dernier point:
      Nodes(nbnodes) = isvl;
      NodInStruct(nbnodes) = ipl;
      Param(nbnodes) = wLast;

      Standard_Real puv;
      Standard_Integer i;
      Standard_Integer iv2;
      for (i = 2; i < GT.NbPoints(); i++)
      {
        // Record 3d point
      GT.Value(cons, gFace, i, puv, P3d, uv);
      nbLocat++;
      Location3d.Bind(nbLocat, P3d);
      NodInStruct(i) = nbLocat;
        // Record 2d point
      v2.Initialize(uv.Coord(), nbDomains, nbLocat, MeshDS_OnCurve);
        iv2=structure->AddNode(v2);
        
      Standard_Integer isv = myvemap.FindIndex(iv2);
      if (isv == 0) isv = myvemap.Add(iv2);
      Nodes(i) = isv;
      NodInStruct(i) = nbLocat;
      Param(i) = puv;
      
      if (orEdge == TopAbs_FORWARD)
        structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Frontier, nbDomains));
      else if (orEdge == TopAbs_REVERSED)
        structure->AddLink(BRepMesh_Edge(iv2, ivf, MeshDS_Frontier, nbDomains));
      else if (orEdge == TopAbs_INTERNAL)
        structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Fixed, nbDomains));
      ivf = iv2;
      }

      P1 = new Poly_PolygonOnTriangulation(Nodes, Param);
      P2 = new Poly_PolygonOnTriangulation(NodInStruct, Param);
    }

    P2->Deflection(otherdefedge);
    TColStd_ListOfTransient L;
    L.Append(P2);
    edges.Bind(edge, L);

    if (ivf != ivl) {
      if (orEdge == TopAbs_FORWARD)
      structure->AddLink(BRepMesh_Edge(ivf, ivl, MeshDS_Frontier, nbDomains));
      else if (orEdge == TopAbs_REVERSED)
      structure->AddLink(BRepMesh_Edge(ivl, ivf, MeshDS_Frontier, nbDomains));
      else if (orEdge == TopAbs_INTERNAL)
      structure->AddLink(BRepMesh_Edge(ivf, ivl, MeshDS_Fixed, nbDomains));
    }
      

  }
  // Si on a dİjÒ vu cette Edge et qu`elle n`est pas dİgİnİrİe on reprend
  // les points du polygone calculİs Ò la premiÛre rencontre :
  else
  {
    if (degener)
    {
      // Create 2d polygon for degenerated edge
      TColStd_Array1OfInteger Nodes(1, 2);
      TColStd_Array1OfReal PPar(1, 2);
      
      Nodes(1) = isvf;
      PPar(1) = wFirst;
      
      Nodes(2) = isvl;
      PPar(2) = wLast;
      
      P1 = new Poly_PolygonOnTriangulation(Nodes, PPar);
    }
    else
    {
      // recuperation du polygone:
      const TColStd_ListOfTransient& L = edges.Find(edge);
      const Handle(Poly_PolygonOnTriangulation)& P = 
        *(Handle(Poly_PolygonOnTriangulation)*)&(L.First());
      const TColStd_Array1OfInteger& NOD = P->Nodes();
      Handle(TColStd_HArray1OfReal) Par = P->Parameters();

      // creation du nouveau:
      const Standard_Integer nbnodes = NOD.Length();
      TColStd_Array1OfInteger Nodes(1, nbnodes);
      TColStd_Array1OfReal PPar(1, nbnodes);
      Standard_Integer iv2;
      
      Nodes(1) = isvf;
      PPar(1) = wFirst;
      
      Nodes(nbnodes) = isvl;
      PPar(nbnodes) = wLast;
      
      if (nbnodes > 2)
      {
        Standard_Integer i;
      if (BRep_Tool::SameParameter(edge))
        {
        for (i = 2; i < nbnodes; i++)
          {
            const Standard_Real puv = Par->Value(i);
          C2d->D0(puv, uv);
          v2.Initialize(uv.Coord(), nbDomains, NOD(i), MeshDS_OnCurve);
            iv2 = structure->AddNode(v2);
      
          Standard_Integer isv = myvemap.FindIndex(iv2);
          if (isv == 0) isv = myvemap.Add(iv2);
            Nodes(i) = isv;
          PPar(i) = puv;
            
            if (orEdge==TopAbs_FORWARD)
            structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Frontier, nbDomains));
          else if (orEdge == TopAbs_REVERSED)
            structure->AddLink(BRepMesh_Edge(iv2, ivf, MeshDS_Frontier, nbDomains));
            else if (orEdge == TopAbs_INTERNAL)
            structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Fixed, nbDomains));
        
            ivf = iv2;
        }
      }
        else
        {
        const Standard_Real wFold = Par->Value(Par->Lower());
        const Standard_Real wLold = Par->Value(Par->Upper());

        Standard_Real wKoef = 1.;
        if ((wFold != wFirst || wLold != wLast) && wLold != wFold)
          {
          wKoef = (wLast - wFirst) / (wLold - wFold);
        }
        
        BRepAdaptor_Curve cons(edge, face);
        Extrema_LocateExtPC pcos;
        pcos.Initialize(cons, cons.FirstParameter(), 
                    cons.LastParameter(), Precision::PConfusion());

        Standard_Real wPrev;
        Standard_Real wCur      = wFirst;
        Standard_Real wCurFound = wFirst;
        for (i = 2; i < nbnodes; i++)
          {
          P3d = Location3d(NOD(i));
            // Record 2d point
          wPrev = wCur;
          wCur  = wFirst + wKoef*(Par->Value(i) - wFold);
            wCurFound += (wCur - wPrev);
          pcos.Perform(P3d, wCurFound);
          if (pcos.IsDone()) wCurFound = pcos.Point().Parameter();
          C2d->D0(wCurFound, uv);
          v2.Initialize(uv.Coord(), nbDomains, NOD(i), MeshDS_OnCurve);
            iv2 = structure->AddNode(v2);
            
          Standard_Integer isv = myvemap.FindIndex(iv2);
          if (isv == 0) isv = myvemap.Add(iv2); 
            Nodes(i) = isv;
          PPar(i) = wCurFound;
            
            if (orEdge==TopAbs_FORWARD)
            structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Frontier, nbDomains));
            else if (orEdge == TopAbs_REVERSED)
            structure->AddLink(BRepMesh_Edge(iv2, ivf, MeshDS_Frontier, nbDomains));
            else if (orEdge == TopAbs_INTERNAL)
            structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Fixed, nbDomains));
        
            ivf = iv2;
        }
      }
      }

      P1 = new Poly_PolygonOnTriangulation(Nodes, PPar);
    
    if (ivf != ivl) {
      if (orEdge == TopAbs_FORWARD) 
      structure->AddLink(BRepMesh_Edge(ivf, ivl, MeshDS_Frontier, nbDomains));
      else if (orEdge == TopAbs_REVERSED)
      structure->AddLink(BRepMesh_Edge(ivl, ivf, MeshDS_Frontier, nbDomains));
      else if (orEdge == TopAbs_INTERNAL)
      structure->AddLink(BRepMesh_Edge(ivf, ivl, MeshDS_Fixed, nbDomains));
      }
    }
  }

  P1->Deflection(defedge);
  if (internaledges.IsBound(edge))
  {
    TColStd_ListOfTransient& L = internaledges.ChangeFind(edge);
    if (orEdge == TopAbs_REVERSED)
      L.Append(P1);
    else
      L.Prepend(P1);
  }
  else
  {
    TColStd_ListOfTransient L1;
    L1.Append(P1);
    internaledges.Bind(edge, L1);
  }
}


//=======================================================================
//function : Update(edge)
//purpose  :
//=======================================================================
Standard_Boolean BRepMesh_FastDiscret::Update(const TopoDS_Edge&  edge,
                                const TopoDS_Face&  face,
                                          const Handle(Geom2d_Curve)& C2d,
                                const Standard_Real defedge,
                                          const Standard_Real first,
                                          const Standard_Real last)
{
  TopLoc_Location l;
  Handle(Poly_Triangulation) T, TNull;
  Handle(Poly_PolygonOnTriangulation) Poly, NullPoly;

  Standard_Integer i = 1;
  Standard_Boolean found = Standard_False;
  do
  {
    BRep_Tool::PolygonOnTriangulation(edge,Poly,T,l,i);
    i++;
    if (!found && !T.IsNull() && T->HasUVNodes() && 
      !Poly.IsNull() && Poly->HasParameters())
    {
      if (Poly->Deflection() <= 1.1*defedge)
      {
        // 2d vertex indices
        TopAbs_Orientation orEdge = edge.Orientation();
        Standard_Integer iv1, iv2, ivl;
        Standard_Integer isv1, isv, isvl;
        
        // Get range on 2d curve
      Standard_Real wFirst, wLast;
      BRep_Tool::Range(edge, face, wFirst, wLast);

        // Get end points on 2d curve
      gp_Pnt2d uvFirst, uvLast;
      BRep_Tool::UVPoints(edge, face, uvFirst, uvLast);

        // Get vertices
      TopoDS_Vertex pBegin, pEnd;
      TopExp::Vertices(edge,pBegin,pEnd);

      const Standard_Boolean sameUV =
          uvFirst.IsEqual(uvLast, Precision::PConfusion());

      //Controle vertice tolerances
        BRepAdaptor_Surface  BS(face, Standard_False);
        Handle(BRepAdaptor_HSurface) gFace = new BRepAdaptor_HSurface(BS);


      gp_Pnt pFirst = gFace->Value(uvFirst.X(), uvFirst.Y());
      gp_Pnt pLast  = gFace->Value(uvLast.X(), uvLast.Y());
      
      Standard_Real mindist = 10. * Max(pFirst.Distance(BRep_Tool::Pnt(pBegin)), 
                                pLast.Distance(BRep_Tool::Pnt(pEnd)));

      if (mindist < BRep_Tool::Tolerance(pBegin) ||
          mindist < BRep_Tool::Tolerance(pEnd) ) mindist = defedge;

      if (sameUV)
        {
        // 1. est-ce vraiment sameUV sans etre denegere
        gp_Pnt2d uvF, uvL;
        C2d->D0(first, uvF);
        C2d->D0(last, uvL);
        if (!uvFirst.IsEqual(uvF, Precision::PConfusion())) {
          uvFirst = uvF;
        }
        if (!uvLast.IsEqual(uvL, Precision::PConfusion())) {
          uvLast = uvL; 
        }
      }

      const TColgp_Array1OfPnt&      Nodes   = T->Nodes();
      const TColStd_Array1OfInteger& Indices = Poly->Nodes();
      Handle(TColStd_HArray1OfReal)  Param   = Poly->Parameters();

      const Standard_Integer nbnodes = Indices.Length();
      TColStd_Array1OfInteger NewNodes(1, nbnodes);
        TColStd_Array1OfInteger NewNodInStruct(1, nbnodes);

        gp_Pnt P3d;
      gp_XY theUV;

        // Process first vertex
        Standard_Integer ipf;
      if (vertices.IsBound(pBegin))
        {
        ipf = vertices.Find(pBegin);
      }
        else
        {
        if (sameUV && vertices.IsBound(pEnd))
          {
          ipf = vertices.Find(pEnd);
        }
        else
          {
          P3d = Nodes(Indices(1));
          if (!l.IsIdentity())
              P3d.Transform(l.Transformation());
            nbLocat++;
            Location3d.Bind(nbLocat,P3d);
          ipf = nbLocat;
        }
          vertices.Bind(pBegin,ipf);
      } 
      NewNodInStruct(1) = ipf;
      theUV = FindUV(pBegin, uvFirst, ipf, gFace, mindist);
      BRepMesh_Vertex vf(theUV,nbDomains,ipf,MeshDS_Frontier);
        iv1 = structure->AddNode(vf);
        isv1 = myvemap.FindIndex(iv1);
        if (isv1 == 0) isv1 = myvemap.Add(iv1);
      NewNodes(1) = isv1;

        // Process last vertex
        Standard_Integer ipl;
        if (pEnd.IsSame(pBegin))
        {
          ipl = ipf;
        }
        else
        {
          if (vertices.IsBound(pEnd))
          {
            ipl = vertices.Find(pEnd);
          }
          else
          {
            if (sameUV)
            {
              ipl = ipf;
              ivl = iv1;
              isv1 = isv1;
            }
            else
            {
              nbLocat++;
              Location3d.Bind(nbLocat,Nodes(Indices(nbnodes)).Transformed(l.Transformation()));
              ipl = nbLocat;
            }
            vertices.Bind(pEnd,ipl);
          }
        }
      NewNodInStruct(nbnodes) = ipl;
      theUV = FindUV(pEnd, uvLast, ipl, gFace, mindist);
      BRepMesh_Vertex vl(theUV,nbDomains,ipl,MeshDS_Frontier);
        
        ivl = structure->AddNode(vl);
        isvl = myvemap.FindIndex(ivl);
        if (isvl == 0) isvl = myvemap.Add(ivl);
        
      NewNodes(nbnodes) = isvl;
      
      gp_Pnt2d uv;
      BRepMesh_Vertex v;
      
      if (BRep_Tool::SameParameter(edge))
        {
        for (i = 2; i < Indices.Length(); i++)
          {
            // Record 3d point
          P3d = Nodes(Indices(i));
          if (!l.IsIdentity())
              P3d.Transform(l.Transformation());
          nbLocat++;
          Location3d.Bind(nbLocat, P3d);
          NewNodInStruct(i) = nbLocat;
            // Record 2d point
          uv = C2d->Value(Param->Value(i));
          v.Initialize(uv.Coord(), nbDomains, nbLocat, MeshDS_Frontier);
            iv2 = structure->AddNode(v);
            isv = myvemap.FindIndex(iv2);
            if (isv == 0) isv = myvemap.Add(iv2);
          NewNodes(i) = isv;
            
            //add links
            if (orEdge == TopAbs_FORWARD)
              structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Frontier,nbDomains));
          else if (orEdge == TopAbs_REVERSED)
            structure->AddLink(BRepMesh_Edge(iv2,iv1,MeshDS_Frontier,nbDomains));
          else if (orEdge == TopAbs_INTERNAL)
            structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Fixed,nbDomains));
          iv1 = iv2;  
        }
          
          // last point
          if (iv1 != ivl) {
          if (orEdge == TopAbs_FORWARD)
            structure->AddLink(BRepMesh_Edge(iv1,ivl,MeshDS_Frontier,nbDomains));
          else if (orEdge == TopAbs_REVERSED)
            structure->AddLink(BRepMesh_Edge(ivl,iv1,MeshDS_Frontier,nbDomains));
          else if (orEdge == TopAbs_INTERNAL)
            structure->AddLink(BRepMesh_Edge(iv1,ivl,MeshDS_Fixed,nbDomains));
        }

          
      }
        else
        {
        const Standard_Real wFold = Param->Value(Param->Lower());
        const Standard_Real wLold = Param->Value(Param->Upper());

        Standard_Real wKoef = 1.;
        if ((wFold != wFirst || wLold != wLast) && wLold != wFold)
          {
          wKoef = (wLast - wFirst) / (wLold - wFold);
        }
        
        BRepAdaptor_Curve cons(edge, face);
        Extrema_LocateExtPC pcos;
        pcos.Initialize(cons, cons.FirstParameter(), 
                    cons.LastParameter(), Precision::PConfusion());

        Standard_Real wPrev;
        Standard_Real wCur      = wFirst;
        Standard_Real wCurFound = wFirst;
        for (i = 2; i < Indices.Length(); i++)
          {
            // Record 3d point
          P3d = Nodes(Indices(i));
          if (!l.IsIdentity())
              P3d.Transform(l.Transformation());
          nbLocat++;
          Location3d.Bind(nbLocat, P3d);
          NewNodInStruct(i) = nbLocat;
            // Record 2d point
          wPrev = wCur;
          wCur  = wFirst + wKoef*(Param->Value(i) - wFold);
            wCurFound += (wCur - wPrev);
          pcos.Perform(P3d, wCurFound);
          if (pcos.IsDone()) wCurFound = pcos.Point().Parameter();
          C2d->D0(wCurFound, uv);
          v.Initialize(uv.Coord(), nbDomains, nbLocat, MeshDS_Frontier);
            iv2 = structure->AddNode(v);
            isv = myvemap.FindIndex(iv2);
            if (isv == 0) isv = myvemap.Add(iv2);
          NewNodes(i) = isv;

            
            //add links
            if (orEdge == TopAbs_FORWARD)
              structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Frontier,nbDomains));
          else if (orEdge == TopAbs_REVERSED)
            structure->AddLink(BRepMesh_Edge(iv2,iv1,MeshDS_Frontier,nbDomains));
          else if (orEdge == TopAbs_INTERNAL)
            structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Fixed,nbDomains));
          iv1 = iv2;              
        }
          
          // last point
          if (iv1 != ivl) {
          if (orEdge == TopAbs_FORWARD)
            structure->AddLink(BRepMesh_Edge(iv1,ivl,MeshDS_Frontier,nbDomains));
          else if (orEdge == TopAbs_REVERSED)
            structure->AddLink(BRepMesh_Edge(ivl,iv1,MeshDS_Frontier,nbDomains));
          else if (orEdge == TopAbs_INTERNAL)
            structure->AddLink(BRepMesh_Edge(iv1,ivl,MeshDS_Fixed,nbDomains));
        }
      }

      Handle(Poly_PolygonOnTriangulation) P1 =
          new Poly_PolygonOnTriangulation(NewNodes, Param->Array1());
      P1->Deflection(defedge);
      if (internaledges.IsBound(edge))
        {
        TColStd_ListOfTransient& L = internaledges.ChangeFind(edge);
          if (edge.Orientation() == TopAbs_REVERSED)
            L.Append(P1);
          else
            L.Prepend(P1);
      }
      else
        {
        TColStd_ListOfTransient L1;
        L1.Append(P1);
        internaledges.Bind(edge, L1);
      }

      Handle(Poly_PolygonOnTriangulation) P2 =
          new Poly_PolygonOnTriangulation(NewNodInStruct, Param->Array1());
      P2->Deflection(defedge);
      TColStd_ListOfTransient L;
      L.Append(P2);
      edges.Bind(edge, L);

      found = Standard_True;
      }
      else
      {
      BRep_Builder B;
      B.UpdateEdge(edge,NullPoly,T,l);
      B.UpdateFace(face,TNull);
      }
    }
    else if (!T.IsNull() && !T->HasUVNodes())
    {
      BRep_Builder B;
      B.UpdateEdge(edge,NullPoly,T,l);
      B.UpdateFace(face,TNull);
    }
  }
  while (!Poly.IsNull());

  return found;
}



//=======================================================================
//function : InternalVertices
//purpose  : 
//=======================================================================
void BRepMesh_FastDiscret::InternalVertices ( const Handle(BRepAdaptor_HSurface)& caro,
                                          BRepMesh_ListOfVertex&              InternalV,
                                          const Standard_Real                 defface,
                                          const BRepMesh_Classifier&          classifier)
{
  BRepMesh_Vertex newV;
  gp_Pnt2d p2d;
  gp_Pnt p3d;
  MeshShape_ListOfSurfacePoint pntsOnSurf;
  
  // travail suivant le type de surface

  const BRepAdaptor_Surface& BS = *(BRepAdaptor_Surface*)&(caro->Surface());
  GeomAbs_SurfaceType thetype = caro->GetType();

  if (thetype == GeomAbs_Plane && !classifier.NaturalRestriction())
  {
    // rajout d`un seul point au milieu.
    const Standard_Real U = 0.5*(myumin+myumax);
    const Standard_Real V = 0.5*(myvmin+myvmax);
    if (classifier.Perform(gp_Pnt2d(U, V)) == TopAbs_IN)
    {
      // Record 3d point
      BRepMesh_GeomTool::D0(caro, U, V, p3d);
      nbLocat++;
      Location3d.Bind(nbLocat, p3d);
      // Record 2d point
      p2d.SetCoord((U-myumin)/deltaX, (V-myvmin)/deltaY);
      newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
      InternalV.Append(newV);
    }
  }
  else if (thetype == GeomAbs_Sphere)
  {
    gp_Sphere S = BS.Sphere();
    const Standard_Real R = S.Radius();

    // Calculate parameters for iteration in V direction
    Standard_Real Dv = 1.0 - (defface/R);
    if (Dv < 0.0) Dv = 0.0;
    Standard_Real oldDv = 2.0 * ACos (Dv);
    Dv  =  .7 * oldDv; //.7 ~= sqrt(2.) - Dv is hypotenuse of triangle when oldDv is legs
    const Standard_Real sv = myvmax - myvmin;
    Dv = sv/((Standard_Integer)(sv/Dv) + 1);
    const Standard_Real pasvmax = myvmax-Dv*0.5;

    //Du can be defined from relation: 2*r*Sin(Du/2) = 2*R*Sin(Dv/2), r = R*Cos(v)
    //here approximate relation r*Du = R*Dv is used
    
    Standard_Real Du, pasu, pasv; //, ru;
    const Standard_Real su = myumax-myumin;
    Standard_Boolean Shift = Standard_False;
    for (pasv = myvmin + Dv; pasv < pasvmax; pasv += Dv)
    {
      // Calculate parameters for iteration in U direction
      // 1.-.365*pasv*pasv is simple approximation of Cos(pasv)
      // with condition that it gives ~.1 when pasv = pi/2
      Du = Dv/(1.-.365*pasv*pasv);
      Du = su/((Standard_Integer)(su/Du) + 1);
      Shift = !Shift;
      const Standard_Real d = (Shift)? Du*.5 : 0.;
      const Standard_Real pasumax = myumax-Du*0.5 + d;
      for (pasu = myumin + Du - d; pasu < pasumax; pasu += Du)
      {
      if (classifier.Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
        {
          // Record 3d point
#ifdef DEB_MESH_CHRONO
        D0Internal++;
#endif
        ElSLib::D0(pasu, pasv, S, p3d);
        nbLocat++;
        Location3d.Bind(nbLocat, p3d);
          // Record 2d point
        p2d.SetCoord((pasu-myumin)/deltaX, (pasv-myvmin)/deltaY);
        newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
        InternalV.Append(newV);
      }
      }
    }
  }
  else if (thetype == GeomAbs_Cylinder)
  {
    gp_Cylinder S = BS.Cylinder();
    const Standard_Real R = S.Radius();

    // Calculate parameters for iteration in U direction
    Standard_Real Du = 1.0 - (defface/R);
    if (Du < 0.0) Du = 0.0;
    Du = 2.0 * ACos (Du);
    if (Du > angle) Du = angle;
    const Standard_Real su = myumax - myumin;
    const Standard_Integer nbU = (Standard_Integer)(su/Du);
    Du = su/(nbU+1);

    // Calculate parameters for iteration in V direction
    const Standard_Real sv = myvmax - myvmin;
    Standard_Integer nbV = (Standard_Integer)( nbU*sv/(su*R) );
    nbV = Min(nbV, 100*nbU);
    Standard_Real Dv = sv/(nbV+1);

    Standard_Real pasu, pasv, pasvmax = myvmax-Dv*0.5, pasumax = myumax-Du*0.5;
    for (pasv = myvmin + Dv; pasv < pasvmax; pasv += Dv) {
      for (pasu = myumin + Du; pasu < pasumax; pasu += Du) {
      if (classifier.Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
        {
          // Record 3d point
        ElSLib::D0(pasu, pasv, S, p3d);
          nbLocat++;
          Location3d.Bind(nbLocat, p3d);
          // Record 2d point
        p2d.SetCoord((pasu-myumin)/deltaX, (pasv-myvmin)/deltaY);
        newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
          InternalV.Append(newV);
      }
      }
    }
  }
  else if (thetype == GeomAbs_Cone)
  {
    Standard_Real R, RefR, SAng;
    gp_Cone C = BS.Cone();
    RefR = C.RefRadius();
    SAng = C.SemiAngle();
    R = Max(Abs(RefR+myvmin*Sin(SAng)), Abs(RefR+myvmax*Sin(SAng)));
    Standard_Real Du, Dv, pasu, pasv;
    Du = Max(1.0e0 - (defface/R),0.0e0);
    Du  = (2.0 * ACos (Du));
    Standard_Integer nbU = (Standard_Integer) ( (myumax-myumin)/Du );
    Standard_Integer nbV = (Standard_Integer) ( nbU*(myvmax-myvmin)/((myumax-myumin)*R) );
    Du = (myumax-myumin)/(nbU+1);
    Dv = (myvmax-myvmin)/(nbV+1);

    Standard_Real pasvmax = myvmax-Dv*0.5, pasumax = myumax-Du*0.5;
    for (pasv = myvmin + Dv; pasv < pasvmax; pasv += Dv) {
      for (pasu = myumin + Du; pasu < pasumax; pasu += Du) {
      if (classifier.Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
        {
          // Record 3d point
        ElSLib::D0(pasu, pasv, C, p3d);
          nbLocat++;
          Location3d.Bind(nbLocat, p3d);
          // Record 2d point
          p2d.SetCoord((pasu-myumin)/deltaX, (pasv-myvmin)/deltaY);
        newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
          InternalV.Append(newV);
      }
      }
    }
  }
  else if (thetype == GeomAbs_Torus)
  {
    gp_Torus T = BS.Torus();

    Standard_Boolean insert;
    Standard_Integer i, j, ParamULength, ParamVLength;
    Standard_Real pp, pasu, pasv;
    Standard_Real r = T.MinorRadius(), R = T.MajorRadius();

    TColStd_SequenceOfReal ParamU, ParamV;

    Standard_Real Du, Dv;//, pasu, pasv;
    Dv = Max(1.0e0 - (defface/r),0.0e0) ;
    Standard_Real oldDv = 2.0 * ACos (Dv);
    oldDv = Min(oldDv, angle);
    Dv  =  0.9*oldDv; //TWOTHIRD * oldDv;
    Dv = oldDv;
    
    Standard_Integer nbV = Max((Standard_Integer)((myvmax-myvmin)/Dv), 2);
    Dv = (myvmax-myvmin)/(nbV+1);
    Standard_Real ru = R + r;
    if (ru > 1.e-16)
    {
      Du  = 2.0 * ACos(Max(1.0 - (defface/ru),0.0));
      if (angle < Du) Du = angle;
      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)((myumax-myumin)/Du), 2);
    nbU = Max(nbU , (int)(nbV*(myumax-myumin)*R/((myvmax-myvmin)*r)/5.));
    Du = (myumax-myumin)/(nbU+1);
    
    if (R < r)
    {
      // comme on recupere les points des edges.
      // dans ce cas, les points ne sont pas representatifs.
            
      //-- On choisit DeltaX et DeltaY de facon a ce qu on ne saute pas 
      //-- de points sur la grille
      for (i = 0; i <= nbU; i++) ParamU.Append(myumin + 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, sort array of parameters
      Standard_Real aDU = FUN_CalcAverageDUV(Up,LenU);
      aDU = Max(aDU, Abs(myumax -  myumin) / (Standard_Real) nbU / 2.);
      Standard_Real dUstd = Abs(myumax -  myumin) / (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(myvmax -  myvmin) / (Standard_Real) nbV / 2.);

    Standard_Real dVstd = Abs(myvmax -  myvmin) / (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 umin = myumin+deltaY*0.1;
    Standard_Real vmin = myvmin+deltaX*0.1;
    Standard_Real umax = myumax-deltaY*0.1;
    Standard_Real vmax = myvmax-deltaX*0.1;

    for (i = 1; i <= Lu; i++)
    {
      pasu = ParamU.Value(i);
      if (pasu >= umin && pasu < umax)
      {
      for (j = 1; j <= Lv; j++)
        {
        pasv = ParamV.Value(j);
        if (pasv >= vmin && pasv < vmax)
          {
          if (classifier.Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
            {
              // Record 3d point
            ElSLib::D0(pasu, pasv, T, p3d);
            nbLocat++;
            Location3d.Bind(nbLocat, p3d);
              // Record 2d point
            p2d.SetCoord((pasu-myumin)/deltaX, (pasv-myvmin)/deltaY);
            newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
            InternalV.Append(newV);
          }
        }
      }
      }
    }
  }
  else if (thetype == GeomAbs_BezierSurface || thetype == GeomAbs_BSplineSurface)
  {
    Standard_Integer i, j;
    Standard_Real pp;

    // Sort sequence of U parameters
    TColStd_SequenceOfReal ParamU;
    ParamU.Append(myumin); ParamU.Append(myumax);
    Standard_Integer ParamULength = 2;
    const Standard_Integer LenU = myUParam.Extent();
    const Standard_Real ddu = 0.02*(myumax-myumin);
    
    Standard_Real U1, U2;

    for (j = 1; j <= LenU; j++)
    {
      pp = myUParam(j);

      U1 = ParamU.Value(1);
      for (i = 2; i <= ParamULength; i++)
      {
      U2 = ParamU.Value(i);
      if (pp < U2)
        {
        if ((U2-pp) < ddu)
          {
          ParamU(i) = pp;
        }
        else if ((pp-U1) > ddu)
          {
          ParamU.InsertBefore(i, pp);
          ParamULength++;
        }
        else if (i != 2)
          {
          ParamU(i-1) = pp;
        }
        break;
      }
        U1 = U2;
      }
    }
    
    // Sort sequence of V parameters
    TColStd_SequenceOfReal ParamV;
    ParamV.Append(myvmin); ParamV.Append(myvmax);
    Standard_Integer ParamVLength = 2;
    const Standard_Integer LenV = myVParam.Extent();
    const Standard_Real ddv = 0.02*(myvmax-myvmin);

    Standard_Real V1, V2;
    for (j = 1; j <= LenV; j++)
    {
      pp = myVParam(j);

      V1 = ParamV.Value(1);
      for (i = 2; i <= ParamVLength; i++)
      {
      V2 = ParamV.Value(i);
      if (pp < V2)
        {
        if ((V2-pp) < ddv)
          {
          ParamV(i) = pp;
        }
        else if ((pp-V1) > ddv)
          {
          ParamV.InsertBefore(i, pp);
          ParamVLength++;
        }
        else if (i != 2)
          {
          ParamV(i-1) = pp;
        }
        break;
      }
        V1 = V2;
      }
    }

    // controle des isos U et insertion eventuelle:

    Handle(Geom_Surface) B;
    if (thetype == GeomAbs_BezierSurface) {
      B = BS.Bezier();
    }
    else {
      B = BS.BSpline();
    }

    gp_Pnt P1, P2, PControl;
    Standard_Real u, v, dist;

    // Insert V parameters by deflection criterion
    for (i = 1; i <= ParamULength; i++)
    {
      Handle(Geom_Curve) IsoU = B->UIso(ParamU.Value(i));
      V1 = ParamV.Value(1);
      P1 = IsoU->Value(V1);
      for (j = 2; j <= ParamVLength;)
      {
      V2 = ParamV.Value(j);
      P2 = IsoU->Value(V2);
      v = 0.5*(V1+V2);
      PControl = IsoU->Value(v);
      if (P1.SquareDistance(P2) > 1.e-5)
        {
        gp_Lin L (P1, gp_Dir(gp_Vec(P1, P2)));
        dist = L.Distance(PControl);
      }
      else
        {
        dist = P1.Distance(PControl);
      }
        if (dist > defface)
        {
          // insertion 
          ParamV.InsertBefore(j, v);
          ParamVLength++;
        }
        else
        {
          V1 = V2;
          P1 = P2;
          j++;
        }
      }
    }

    for (i = 2; i < ParamVLength; i++)
    {
      v = ParamV.Value(i);
      Handle(Geom_Curve) IsoV = B->VIso(v);
      U1 = ParamU.Value(1);
      P1 = IsoV->Value(U1);
      for (j = 2; j <= ParamULength;)
      {
      U2 = ParamU.Value(j);
      P2 = IsoV->Value(U2);
      u = 0.5*(U1+U2);
      PControl = IsoV->Value(u);
      if (P1.SquareDistance(P2) > 1.e-5)
        {
        gp_Lin L (P1, gp_Dir(gp_Vec(P1, P2)));
        dist = L.Distance(PControl);
      }
      else
        {
        dist = P1.Distance(PControl);
      }
        if (dist > defface)
        {
          // insertion 
          ParamU.InsertBefore(j, u);
          ParamULength++;
        }
        else
        {
          U1 = U2;
          P1 = P2;
          if (j < ParamULength)
          {
            // Classify intersection point
            if (classifier.Perform(gp_Pnt2d(U1, v)) == TopAbs_IN)
            {
              // Record 3d point
              nbLocat++;
              Location3d.Bind(nbLocat, P1);
              // Record 2d point
              p2d.SetCoord((U1-myumin)/deltaX, (v-myvmin)/deltaY);
              newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
              InternalV.Append(newV);
            }
          }
          j++;
        }
      }
    }
  }
  else
  {
    const Standard_Real theangle = 0.35;

    Standard_Integer i, j, nbpointsU = 10, nbpointsV = 10;
    Adaptor3d_IsoCurve tabu[10], tabv[10];

    TColStd_SequenceOfReal ParamU, ParamV;
    Standard_Real u, v, du, dv;
    Standard_Integer iu, iv;
    Standard_Real f, l;

    du = (myumax-myumin) / (nbpointsU+1); dv = (myvmax-myvmin) / (nbpointsV+1);

    for (iu = 1; iu <= nbpointsU; iu++) {
      u = myumin + iu*du;
      tabu[iu-1].Load(caro);
      tabu[iu-1].Load(GeomAbs_IsoU, u);
    }

    for (iv = 1; iv <= nbpointsV; iv++) {
      v = myvmin + iv*dv;
      tabv[iv-1].Load(caro);
      tabv[iv-1].Load(GeomAbs_IsoV, v);
    }

    Standard_Integer imax = 1, MaxV = 0;
    
    GCPnts_TangentialDeflection* tabGU = new GCPnts_TangentialDeflection[nbpointsU];

    for (i = 0; i <= nbpointsU-1; i++) {
      f = Max(myvmin, tabu[i].FirstParameter());
      l = Min(myvmax, tabu[i].LastParameter());
      GCPnts_TangentialDeflection theDeflection(tabu[i], f, l, theangle, 0.7*defface, 2);
      tabGU[i] = theDeflection;
      if (tabGU[i].NbPoints() > MaxV) {
      MaxV = tabGU[i].NbPoints();
      imax = i;
      }
    }
    
    // recuperation du tableau de parametres V:
    Standard_Integer NV = tabGU[imax].NbPoints();
    for (i = 1; i <= NV; i++) {
      ParamV.Append(tabGU[imax].Parameter(i));
    }
    delete [] tabGU;

    imax = 1;
    Standard_Integer MaxU = 0;

    GCPnts_TangentialDeflection* tabGV = new GCPnts_TangentialDeflection[nbpointsV];

    for (i = 0; i <= nbpointsV-1; i++) {
      f = Max(myumin, tabv[i].FirstParameter());
      l = Min(myumax, tabv[i].LastParameter());
      GCPnts_TangentialDeflection thedeflection2(tabv[i], f, l, theangle, 0.7*defface, 2);
      tabGV[i] = thedeflection2;
      if (tabGV[i].NbPoints() > MaxU) {
      MaxU = tabGV[i].NbPoints();
      imax = i;
      }
    }
    
    // recuperation du tableau de parametres U:
    Standard_Integer NU = tabGV[imax].NbPoints();
    for (i = 1; i <= NU; i++) {
      ParamU.Append(tabGV[imax].Parameter(i));
    }
    delete [] tabGV;
    
    if (ParamU.Length() == 2) {
      ParamU.InsertAfter(1, (myumax+myumin)*0.5);
    }
    if (ParamV.Length() == 2) {
      ParamV.InsertAfter(1, (myvmax+myvmin)*0.5);
    }
    
    TColStd_SequenceOfReal InsertV, InsertU;
    gp_Pnt P1;

    Adaptor3d_IsoCurve IsoV;
    IsoV.Load(caro);

    Standard_Integer Lu = ParamU.Length(), Lv = ParamV.Length();
    
    for (i = 2; i < Lv; i++) {
      v = ParamV.Value(i);
      IsoV.Load(GeomAbs_IsoV, v);
      for (j = 2; j < Lu; j++) {
      u = ParamU.Value(j);
      if (classifier.Perform(gp_Pnt2d(u, v)) == TopAbs_IN)
        {
          // Record 3d point
        P1 = IsoV.Value(u);
        nbLocat++;
          Location3d.Bind(nbLocat, P1);
          // Record 2d point
        p2d.SetCoord((u-myumin)/deltaX, (v-myvmin)/deltaY);
        newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
        InternalV.Append(newV); 
      }
      }
    } 
  }
#ifdef DEB_MESH_CHRONO
  chInternal.Stop();
#endif
}

//=======================================================================
//function : InternalVerticesForRectangle
//purpose  : 
//=======================================================================
void BRepMesh_FastDiscret::
  InternalVerticesForRectangle(const Handle(BRepAdaptor_HSurface)& caro,
                         BRepMesh_ListOfVertex&              InternalV,
                         const TColStd_ListOfReal&           U1Params,
                         const TColStd_ListOfReal&           U2Params,
                         const TColStd_ListOfReal&           V1Params,
                         const TColStd_ListOfReal&           V2Params)
{
  BRepMesh_Vertex newV;
  gp_XY p2d;
  gp_Pnt p3d;


  Standard_Integer i, j, nbu = U1Params.Extent(), 
                         nbv = V1Params.Extent();
  TColStd_Array1OfReal u1prs(1, nbu), u2prs(1, nbu), 
                       v1prs(1, nbv), v2prs(1,nbv);
  TColStd_ListIteratorOfListOfReal Itu1(U1Params), Itu2(U2Params), 
                                   Itv1(V1Params), Itv2(V2Params);
  for(i = 1; Itu1.More(); Itu1.Next(), Itu2.Next(),++i) {
    u1prs(i) = Itu1.Value();
    u2prs(i) = Itu2.Value();
  }
  for(i = 1; Itv1.More(); Itv1.Next(), Itv2.Next(),++i) {
    v1prs(i) = Itv1.Value();
    v2prs(i) = Itv2.Value();
  }
  
  SortTools_ShellSortOfReal aSorter;
  TCollection_CompareOfReal aCompare;
  aSorter.Sort(u1prs, aCompare);
  aSorter.Sort(u2prs, aCompare);
  aSorter.Sort(v1prs, aCompare);
  aSorter.Sort(v2prs, aCompare);

  gp_XY f00(u1prs(1), v1prs(1));
  gp_XY f01(u2prs(1), v1prs(nbv));
  gp_XY f10(u1prs(nbu), v2prs(1));
  gp_XY f11(u2prs(nbu), v2prs(nbv));

  Standard_Real cu = 1./(nbu-1), cv = 1./(nbv-1);

  //Linear Coons-Gordon interpolation for inner points
  Standard_Real x, y, x1, y1;
  for(i = 2, x = cu; i < nbu; ++i, x += cu) {
    x1 = 1.-x;
    gp_XY fx0(u1prs(i), v1prs(1));
    gp_XY fx1(u2prs(i), v1prs(nbv));
 
    for(j = 2, y = cv; j < nbv; ++j, y += cv) {
      y1 = 1.-y;

      gp_XY f0y(u1prs(1), v1prs(j));
      gp_XY f1y(u1prs(nbu), v2prs(j));
            
      gp_XY fxy = y1*fx0 + y*fx1 + x1*f0y + x*f1y -
                  (x1*(y1*f00 + y*f01) + x*(y1*f10 + y*f11));

      // Record 3d point
      BRepMesh_GeomTool::D0(caro, fxy.X(), fxy.Y(), p3d);
      nbLocat++;
      Location3d.Bind(nbLocat, p3d);
      // Record 2d point
      p2d.SetCoord((fxy.X()-myumin)/deltaX, (fxy.Y()-myvmin)/deltaY);
      newV.Initialize(p2d, nbDomains, nbLocat, MeshDS_Free);
      InternalV.Append(newV);
    }
  }
      


}


//=======================================================================
//function : Control
//purpose  : 
//=======================================================================
Standard_Real BRepMesh_FastDiscret::Control
  (const Handle(BRepAdaptor_HSurface)& caro,
   const Standard_Real                 defface,
   BRepMesh_ListOfVertex&              InternalV,
   TColStd_ListOfInteger&              badTriangles,
   TColStd_ListOfInteger&              nulTriangles,
   BRepMesh_Delaun&                    trigu,
   const Standard_Boolean              isfirst)
{
  //IMPORTANT: Constants used in calculations
  const Standard_Real MinimalArea2d = 1.e-9;
  const Standard_Real MinimalSqLength3d = 1.e-12;

  // Define the number of iterations
  Standard_Integer myNbIterations = 11;
  const Standard_Integer nbPasses = (isfirst? 1 : myNbIterations);

  // Initialize stop condition
  Standard_Boolean allDegenerated = Standard_False;
  Standard_Integer nbInserted = 1;

  // Create map of links to skip already processed
  Standard_Integer nbtriangles;

  nbtriangles = structure->ElemOfDomain(nbDomains).Extent();
  if (nbtriangles <= 0) return -1.0;
  MeshShape_MapOfCouple theCouples(3*nbtriangles);

  gp_XY mi2d;
  gp_XYZ vecEd1, vecEd2, vecEd3;
  gp_Pnt pDef;
  Standard_Real dv, defl, maxdef;
  Standard_Integer pass, nf, nl;
  Standard_Integer v1, v2, v3, e1, e2, e3;
  Standard_Boolean o1, o2, o3;
  Standard_Boolean m1, m2, m3;
  BRepMesh_Vertex InsVertex;
  Standard_Boolean caninsert;

  // Perform refinement passes
  for (pass = 1; pass <= nbPasses && nbInserted && !allDegenerated; pass++)
  {
    InternalV.Clear();
    badTriangles.Clear();
    
    // Reset stop condition
    allDegenerated = Standard_True;
    nbInserted = 0;
    maxdef = -1.0;

    // Do not insert nodes in last pass in non-SharedMode
    caninsert = (WithShare || pass < nbPasses);

    // Read mesh size
    nbtriangles = structure->ElemOfDomain(nbDomains).Extent();
    if (nbtriangles <= 0) break;

    // Iterate on current triangles
    BRepMesh_ListOfVertex InsVertices;
    TColStd_MapOfInteger usedEdges;
   // TColStd_MapIteratorOfMapOfInteger triDom;
    MeshDS_MapOfInteger::Iterator triDom;
    const MeshDS_MapOfInteger& TriMap = structure->ElemOfDomain(nbDomains);
    triDom.Initialize(TriMap);
    Standard_Integer aNbPnt = 0;
    for (; triDom.More(); triDom.Next())
    {
      const BRepMesh_Triangle& curTri=Triangle(triDom.Key());
      if (curTri.Movability()==MeshDS_Deleted) continue;
      curTri.Edges(e1, e2, e3, o1, o2, o3);
      
      const BRepMesh_Edge& edg1=Edge(e1);
      const BRepMesh_Edge& edg2=Edge(e2);
      const BRepMesh_Edge& edg3=Edge(e3);
      
      m1 = (edg1.Movability() == MeshDS_Frontier);
      m2 = (edg2.Movability() == MeshDS_Frontier);
      m3 = (edg3.Movability() == MeshDS_Frontier);
      if (o1) {
      v1=edg1.FirstNode();
      v2=edg1.LastNode();
      }
      else {
      v1=edg1.LastNode();
      v2=edg1.FirstNode();
      }
      if (o2)
      v3=edg2.LastNode();
      else
      v3=edg2.FirstNode();

      const BRepMesh_Vertex& vert1=Vertex(v1);
      const BRepMesh_Vertex& vert2=Vertex(v2);
      const BRepMesh_Vertex& vert3=Vertex(v3);

      const gp_XYZ& p1=Location3d(vert1.Location3d()).Coord();
      const gp_XYZ& p2=Location3d(vert2.Location3d()).Coord();
      const gp_XYZ& p3=Location3d(vert3.Location3d()).Coord();

      vecEd1 = p2 - p1;
      vecEd2 = p3 - p2;
      vecEd3 = p1 - p3;

      // Check for degenerated triangle
      if (vecEd1.SquareModulus() < MinimalSqLength3d ||
          vecEd2.SquareModulus() < MinimalSqLength3d ||
          vecEd3.SquareModulus() < MinimalSqLength3d) 
      {
        nulTriangles.Append(triDom.Key());
        continue;
      }

      allDegenerated = Standard_False;

      gp_XY xy1(vert1.Coord().X()*deltaX+myumin,vert1.Coord().Y()*deltaY+myvmin);
      gp_XY xy2(vert2.Coord().X()*deltaX+myumin,vert2.Coord().Y()*deltaY+myvmin);
      gp_XY xy3(vert3.Coord().X()*deltaX+myumin,vert3.Coord().Y()*deltaY+myvmin);

      // Check triangle area in 2d
      if (Abs((xy2-xy1)^(xy3-xy1)) < MinimalArea2d)
      {
            nulTriangles.Append(triDom.Key());
        continue;
      }

      // Check triangle normal
      gp_XYZ normal(vecEd1^vecEd2);
      dv = normal.Modulus();
      if (dv < Precision::Confusion())
      {
            nulTriangles.Append(triDom.Key());
        continue;
      }
      normal /= dv;

      // Check deflection on triangle
      mi2d = (xy1+xy2+xy3)/3.0;
      caro->D0(mi2d.X(), mi2d.Y(), pDef);
      defl = Abs(normal*(pDef.XYZ()-p1));
      if (defl > maxdef) maxdef = defl;
      if (defl > defface)
      {
        if (isfirst) break;
        if (caninsert)
        {
          // Record new vertex
          aNbPnt++;
          nbLocat++;
          Location3d.Bind(nbLocat,pDef);
          mi2d.SetCoord((mi2d.X()-myumin)/deltaX,(mi2d.Y()-myvmin)/deltaY);
          InsVertex.Initialize(mi2d,nbDomains,nbLocat,MeshDS_Free);
          InternalV.Append(InsVertex);
        }
        badTriangles.Append(triDom.Key());
      }
      
      if (!m2) // Not a boundary
      {
        // Check if this link was already processed
        if (v2 < v3) { nf = v2; nl = v3; } else { nf = v3; nl = v2; }
        if (theCouples.Add(MeshShape_Couple(nf,nl)))
        {
          // Check deflection on edge 1
          mi2d = (xy2+xy3)*0.5;
          caro->D0(mi2d.X(), mi2d.Y(), pDef);
          defl = Abs(normal*(pDef.XYZ()-p1));
          if (defl > maxdef) maxdef = defl;
          if (defl > defface)
          {
            if (isfirst) break;
            if (caninsert)
            {
              // Record new vertex
             aNbPnt++;              
              nbLocat++;
              Location3d.Bind(nbLocat,pDef);
              mi2d.SetCoord((mi2d.X()-myumin)/deltaX,(mi2d.Y()-myvmin)/deltaY);
              InsVertex.Initialize(mi2d,nbDomains,nbLocat,MeshDS_Free);
              InternalV.Append(InsVertex);
            }
            badTriangles.Append(triDom.Key());
          }
        }
      }

      if (!m3) // Not a boundary
      {
        // Check if this link was already processed
        if (v1 < v3) { nf = v1; nl = v3; } else { nf = v3; nl = v1; }
        if (theCouples.Add(MeshShape_Couple(nf,nl)))
        {
          // Check deflection on edge 2
          mi2d = (xy3+xy1)*0.5;
          caro->D0(mi2d.X(), mi2d.Y(), pDef);
          defl = Abs(normal*(pDef.XYZ()-p2));
          if (defl > maxdef) maxdef = defl;
          if (defl > defface)
          {
            if (isfirst) break;
            if (caninsert)
            {
              // Record new vertex
              aNbPnt++;
              nbLocat++;
              Location3d.Bind(nbLocat,pDef);
              mi2d.SetCoord((mi2d.X()-myumin)/deltaX,(mi2d.Y()-myvmin)/deltaY);
              InsVertex.Initialize(mi2d,nbDomains,nbLocat,MeshDS_Free);
              InternalV.Append(InsVertex);
            }
            badTriangles.Append(triDom.Key());
          }
        }
      }

      if (!m1) // Not a boundary
      {
        // Check if this link was already processed
        if (v1 < v2) { nf = v1; nl = v2; } else { nf = v2; nl = v1; }
        if (theCouples.Add(MeshShape_Couple(nf,nl)))
        {
          // Check deflection on edge 3
          mi2d = (xy1+xy2)*0.5;
          caro->D0(mi2d.X(), mi2d.Y(), pDef);
          defl = Abs(normal*(pDef.XYZ()-p3));
          if (defl > maxdef) maxdef = defl;
          if (defl > defface)
          {
            if (isfirst) break;
            if (caninsert)
            {
              // Record new vertex
              aNbPnt++;
              nbLocat++;
              Location3d.Bind(nbLocat,pDef);
              mi2d.SetCoord((mi2d.X()-myumin)/deltaX,(mi2d.Y()-myvmin)/deltaY);              
              InsVertex.Initialize(mi2d,nbDomains,nbLocat,MeshDS_Free);
              InternalV.Append(InsVertex);
            }
            badTriangles.Append(triDom.Key());
          }
        }
      }
    }

    if (!isfirst && InternalV.Extent() > 0) 
    {
      BRepMesh_Array1OfVertexOfDelaun verttab(1, InternalV.Extent());
      BRepMesh_ListIteratorOfListOfVertex itVer(InternalV);
      Standard_Integer ipn = 1;
      for (; itVer.More(); itVer.Next())
        verttab(ipn++) = itVer.Value();
      trigu.AddVertices(verttab);
      nbInserted++;
    }
  }

  return maxdef;
}


//=======================================================================
//function : Plan
//purpose  : 
//=======================================================================
static Standard_Boolean Plan(const MeshShape_SurfacePoint& P1,
                       const MeshShape_SurfacePoint& P2,
                       const MeshShape_SurfacePoint& P3,
                       gp_XYZ&        eqPlan,
                       Standard_Real& diPolr)
{
  gp_XYZ v1(P2.Coord()-P1.Coord());
  if (v1.Modulus()>Precision::Confusion()) {
    gp_XYZ v2(P3.Coord()-P2.Coord());
    if (v2.Modulus()>Precision::Confusion()) {
      gp_XYZ v3(P1.Coord()-P3.Coord());
      if (v3.Modulus()>Precision::Confusion()) {
      eqPlan=v1^v2;
      if (eqPlan.Modulus()>0.) {
        eqPlan.Normalize();
        diPolr=P1.Coord()*eqPlan;
        return Standard_True;
      }
      }
    }
  }
  diPolr=0.;
  eqPlan.SetCoord(0., 0., 0.);
  return Standard_False;
}


//=======================================================================
//function : Uindex
//purpose  : 
//=======================================================================

Standard_Integer BRepMesh_FastDiscret::Uindex(const MeshShape_SurfacePoint& P) const
{
//#ifndef DEB 
  return (Standard_Integer) Round((P.UV().X()-myumin)/(myumax-myumin)*LIMITE_TRIANGULATION);
//#else
//  return Round((P.UV().X()-myumin)/(myumax-myumin)*LIMITE_TRIANGULATION);
//#endif
}


//=======================================================================
//function : Vindex
//purpose  : 
//=======================================================================

Standard_Integer BRepMesh_FastDiscret::Vindex(const MeshShape_SurfacePoint& P) const
{
//#ifndef DEB
  return (Standard_Integer) Round((P.UV().Y()-myvmin)/(myvmax-myvmin)*LIMITE_TRIANGULATION);
//#else
//  return Round((P.UV().Y()-myvmin)/(myvmax-myvmin)*LIMITE_TRIANGULATION);
//#endif
}


//=======================================================================
//function : Append
//purpose  : 
//=======================================================================

void BRepMesh_FastDiscret::Append(MeshShape_ListOfSurfacePoint& pntsOnSurf,
                        const MeshShape_SurfacePoint& P)
{
  Standard_Integer IU = Uindex(P), IV = Vindex(P);

  // pas d append dans le cas d une frontiere.
  if ((IU != 0)                    && 
      (IV != 0)                    && 
      (IU != LIMITE_TRIANGULATION) && 
      (IV != LIMITE_TRIANGULATION)) {
    MeshShape_Couple C(IU, IV);
    if (!mymap.Contains(C)) {
      mymap.Add(C);
      pntsOnSurf.Append(P);
    }
  }
}

//=======================================================================
//function : AddInShape
//purpose  : 
//=======================================================================
void BRepMesh_FastDiscret::AddInShape(const TopoDS_Face&  face,
                          const Standard_Real defface)
{
//  gp_Pnt Pt;

  try{
  //TColStd_MapIteratorOfMapOfInteger it;
  MeshDS_MapOfInteger::Iterator it;

  Standard_Integer e1, e2, e3, nTri;
  Standard_Integer v1, v2, v3, iv1, iv2, iv3;
  Standard_Integer i, index;
  Standard_Boolean o1, o2, o3;
  TopAbs_Orientation orFace = face.Orientation();

  //const TColStd_MapOfInteger& TriMap = structure->ElemOfDomain(nbDomains);
  const MeshDS_MapOfInteger& TriMap = structure->ElemOfDomain(nbDomains);
  it.Initialize(TriMap);
    
  nTri = TriMap.Extent();

  if (nTri != 0) {
    
    Poly_Array1OfTriangle Tri(1, nTri);
    
    i = 1;
    
    for (; it.More(); it.Next()) {
      structure->GetElement(it.Key()).Edges(e1, e2, e3, o1, o2, o3);
      
      const BRepMesh_Edge& ve1=structure->GetLink(e1);
      const BRepMesh_Edge& ve2=structure->GetLink(e2);
      const BRepMesh_Edge& ve3=structure->GetLink(e3);
      
      if (o1) {
      v1=ve1.FirstNode();
      }
      else {
      v1=ve1.LastNode();
      }
      if (o2)
      {
        v2=ve2.FirstNode();
      v3=ve2.LastNode();
      }
      else
      {
      v3=ve2.FirstNode();
            v2=ve2.LastNode();
      }
      
      iv1 = myvemap.FindIndex(v1);
      if (iv1 == 0) iv1 = myvemap.Add(v1);
      iv2 = myvemap.FindIndex(v2);
      if (iv2 == 0) iv2 = myvemap.Add(v2);
      iv3 = myvemap.FindIndex(v3);
      if (iv3 == 0) iv3 = myvemap.Add(v3);
      
      if (orFace == TopAbs_REVERSED) Tri(i++).Set(iv1, iv3, iv2);
      else Tri(i++).Set(iv1, iv2, iv3);
    }
    
    Standard_Integer nbVertices = myvemap.Extent();
    Handle(Poly_Triangulation) T = new Poly_Triangulation(nbVertices, nTri, Standard_True);
    Poly_Array1OfTriangle& Trian = T->ChangeTriangles();
    Trian = Tri;
    TColgp_Array1OfPnt&  Nodes = T->ChangeNodes();
    TColgp_Array1OfPnt2d& Nodes2d = T->ChangeUVNodes();
    
    for (i = 1; i <= nbVertices; i++) {
      index = myvemap.FindKey(i);
      Nodes(i) = Pnt(index);
      Nodes2d(i).SetXY(Vertex(index).Coord());
    }
    
    T->Deflection(defface);
    
    // stockage de la triangulation dans la BRep.
    BRep_Builder B;
    TopLoc_Location l = face.Location();
    if (!l.IsIdentity()) {
      gp_Trsf tr = l.Transformation();
      tr.Invert();
      for (i = Nodes.Lower(); i <= Nodes.Upper(); i++) 
      Nodes(i).Transform(tr);
    }
    B.UpdateFace(face, T);

    // mise en place des polygones sur triangulation dans la face:
    MeshShape_DataMapIteratorOfDataMapOfShapeListOfTransient It(internaledges);

    for (; It.More(); It.Next()) {
      if (It.Value().Extent() == 1) {
      const Handle(Poly_PolygonOnTriangulation)& NOD = 
        *((Handle(Poly_PolygonOnTriangulation)*)&(It.Value().First()));
      B.UpdateEdge(TopoDS::Edge(It.Key()), NOD, T, l);
      }
      else {
      const Handle(Poly_PolygonOnTriangulation)& NOD1 = 
        *((Handle(Poly_PolygonOnTriangulation)*)&(It.Value().First()));
      const Handle(Poly_PolygonOnTriangulation)& NOD2 = 
        *((Handle(Poly_PolygonOnTriangulation)*)&(It.Value().Last()));
      B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, NOD2, T, l);
      }
    }
  }
 }
 catch(Standard_Failure)
 {
   MESH_FAILURE(face);
 }
}


//=======================================================================
//function : output
//purpose  : 
//=======================================================================
02742 Standard_Integer BRepMesh_FastDiscret::NbTriangles() const
{
  return structure->NbElements();
}

//=======================================================================
//function : Triangle
//purpose  : 
//=======================================================================

const BRepMesh_Triangle& BRepMesh_FastDiscret::Triangle
02753   (const Standard_Integer Index) const
{
  return structure->GetElement(Index);
}

//=======================================================================
//function : NbEdges
//purpose  : 
//=======================================================================

02763 Standard_Integer BRepMesh_FastDiscret::NbEdges() const
{
  return structure->NbLinks();
}

//=======================================================================
//function : Edge
//purpose  : 
//=======================================================================

02773 const BRepMesh_Edge& BRepMesh_FastDiscret::Edge(const Standard_Integer Index) const
{
  return structure->GetLink(Index);
}

//=======================================================================
//function : NbVertices
//purpose  : 
//=======================================================================

02783 Standard_Integer BRepMesh_FastDiscret::NbVertices() const
{
  return structure->NbNodes();
}

//=======================================================================
//function : Vertex
//purpose  : 
//=======================================================================

const BRepMesh_Vertex& BRepMesh_FastDiscret::Vertex
02794   (const Standard_Integer Index) const
{
  return structure->GetNode(Index);
}

//=======================================================================
//function : Pnt
//purpose  : 
//=======================================================================

02804 const gp_Pnt& BRepMesh_FastDiscret::Pnt(const Standard_Integer Index) const
{
  return Location3d(structure->GetNode(Index).Location3d());
}

//=======================================================================
//function : NbPoint3d
//purpose  : 
//=======================================================================

02814 Standard_Integer BRepMesh_FastDiscret::NbPoint3d() const
{
  return nbLocat;
}

//=======================================================================
//function : Point3d
//purpose  : 
//=======================================================================

02824 const gp_Pnt& BRepMesh_FastDiscret::Point3d(const Standard_Integer Index) const
{
  return Location3d(Index);
}

//=======================================================================
//function : Normal
//purpose  : 
//=======================================================================

02834 void BRepMesh_FastDiscret::Normal(const Standard_Integer Index,
                         gp_Pnt& Pon, gp_Dir& Nor) const
{
  const BRepMesh_Vertex& v=structure->GetNode(Index);
  BRepMesh_GeomTool::Normal(new BRepAdaptor_HSurface(Domains(v.Domain())), 
               v.Coord().X(), v.Coord().Y(), Pon, Nor);
}

//=======================================================================
//function : Result
//purpose  : 
//=======================================================================

Handle(BRepMesh_DataStructureOfDelaun) BRepMesh_FastDiscret::Result() const
{
  return structure;
}


//=======================================================================
//function : NbDomains
//purpose  : 
//=======================================================================

02858 Standard_Integer BRepMesh_FastDiscret::NbDomains() const 
{ return nbDomains;}

//=======================================================================
//function : DomainFace
//purpose  : 
//=======================================================================

02866 const TopoDS_Face& BRepMesh_FastDiscret::DomainFace(const Standard_Integer Index) const 
{ 
  return Domains(Index); 
}

//=======================================================================
//function : VerticesOfDomain
//purpose  : 
//=======================================================================

02876 void BRepMesh_FastDiscret::VerticesOfDomain(const Standard_Integer Index,
                              MeshDS_MapOfInteger&  Indices) const 
{ 
  Indices.Clear();
  
  // recup de la map des edges.
  const MeshDS_MapOfInteger& edmap = structure->LinkOfDomain(Index);

  // iterateur sur les edges.
  //TColStd_MapIteratorOfMapOfInteger iter(edmap);
  MeshDS_MapOfInteger::Iterator iter(edmap);
  
  Standard_Integer ind_edge;
  for (iter.Reset(); iter.More(); iter.Next()) {
    ind_edge = iter.Key();
    const BRepMesh_Edge& Ed = Edge(ind_edge);
    Indices.Add(Ed.FirstNode());
    Indices.Add(Ed.LastNode());
  }
}



//=======================================================================
//function : EdgesOfDomain
//purpose  : 
//=======================================================================

02904 void BRepMesh_FastDiscret::EdgesOfDomain(const Standard_Integer Index,
                             MeshDS_MapOfInteger&  Indices) const 
{ 
  Indices = structure->LinkOfDomain(Index);
}


//=======================================================================
//function : TrianglesOfDomain
//purpose  : 
//=======================================================================

02916 void BRepMesh_FastDiscret::TrianglesOfDomain(const Standard_Integer Index,
                             MeshDS_MapOfInteger&  Indices) const 
{ 
  Indices = structure->ElemOfDomain(Index);
}



//=======================================================================
//function : GetDeflection
//purpose  : 
//=======================================================================

02929 Standard_Real BRepMesh_FastDiscret::GetDeflection() const
{
  return deflection;
}


//=======================================================================
//function : GetAngle
//purpose  : 
//=======================================================================

02940 Standard_Real BRepMesh_FastDiscret::GetAngle() const
{
  return angle;
}


//=======================================================================
//function : FindUV
//purpose  : 
//=======================================================================

gp_XY BRepMesh_FastDiscret::FindUV(const TopoDS_Vertex&   V,
                         const gp_Pnt2d&        XY, 
                         const Standard_Integer ip,
                         const Handle(BRepAdaptor_HSurface)& S,
                         const Standard_Real mindist) 
{
  gp_XY theUV;
  if (mylocation2d.IsBound(ip))
  {
    BRepMesh_ListOfXY& L = mylocation2d.ChangeFind(ip);
    theUV = L.First();
    if (L.Extent() != 1)
    {
      BRepMesh_ListIteratorOfListOfXY it(L);
      it.Next();
      Standard_Real dd, dmin = XY.Distance(gp_Pnt2d(theUV));
      for (; it.More(); it.Next())
      {
      dd = XY.Distance(gp_Pnt2d(it.Value()));
      if (dd < dmin)
        {
        theUV = it.Value();
        dmin = dd;
      }
      }
    }

    const Standard_Real tol = Min(2. * BRep_Tool::Tolerance(V), mindist);

    const Standard_Real Utol2d = .5 * (S->LastUParameter() - S->FirstUParameter());
    const Standard_Real Vtol2d = .5 * (S->LastVParameter() - S->FirstVParameter());

    const gp_Pnt p1 = S->Value(theUV.X(), theUV.Y());
    const gp_Pnt p2 = S->Value(XY.X(), XY.Y());

    if (Abs(theUV.X() - XY.X()) > Utol2d ||
        Abs(theUV.Y() - XY.Y()) > Vtol2d ||
      !p1.IsEqual(p2, tol))
    {
      theUV = XY.Coord();
      L.Append(theUV);
    }
  }
  else
  {
    theUV = XY.Coord();
    BRepMesh_ListOfXY L;
    L.Append(theUV);
    mylocation2d.Bind(ip, L);
  }
  return theUV;
}


BRepMesh_Status BRepMesh_FastDiscret::CurrentFaceStatus() const
{
  return myfacestate;
}

static Standard_Boolean GetVertexParameters(const TopoDS_Vertex& theVert, 
                                         const TopoDS_Face& theFace,
                                         gp_Pnt2d& thePoint)
{
  TopLoc_Location L;
  const Handle(Geom_Surface)& S = BRep_Tool::Surface(theFace,L);
  L = L.Predivided(theVert.Location());
  BRep_ListIteratorOfListOfPointRepresentation itpr
    ((*((Handle(BRep_TVertex)*) &theVert.TShape()))->Points());
  // On regarde dabord si il y des PointRepresentation (cas non Manifold)

  while (itpr.More()) {
    if (itpr.Value()->IsPointOnSurface(S,L)) {
      thePoint.SetCoord(itpr.Value()->Parameter(),
                 itpr.Value()->Parameter2());
      return Standard_True;
    }
    itpr.Next();
  }
  return Standard_False;
}
//=======================================================================
//function : Add
//purpose  : method intended to addition internav vertices in triangulation.
//=======================================================================

void BRepMesh_FastDiscret::Add(const TopoDS_Vertex&                theVert, 
                         const TopoDS_Face&                  theFace, 
                         const Handle(BRepAdaptor_HSurface)& thegFace)
                         
{
  const TopAbs_Orientation anOrient = theVert.Orientation();
  gp_Pnt2d uvXY;
  if( anOrient != TopAbs_INTERNAL || !GetVertexParameters(theVert,theFace,uvXY))
    return;
  Standard_Integer indVert =0;
  if (vertices.IsBound(theVert))
    indVert = vertices.Find(theVert);
  else
  {
    nbLocat++;
    Location3d.Bind(nbLocat, BRep_Tool::Pnt(theVert));
    indVert = nbLocat;
    vertices.Bind(theVert, indVert);
  }
  Standard_Real mindist = BRep_Tool::Tolerance(theVert);
  // gp_Pnt2d uvXY = BRep_Tool::Parameters(theVert,theFace);
  gp_XY anUV = FindUV(theVert, uvXY, indVert, thegFace, mindist);
  BRepMesh_Vertex vf(anUV, nbDomains, indVert, MeshDS_Fixed);
  Standard_Integer ivff = structure->AddNode(vf);
  Standard_Integer isvf = myvemap.FindIndex(ivff);
  if (isvf == 0) isvf = myvemap.Add(ivff);  
}

Generated by  Doxygen 1.6.0   Back to index