241c92c45696230636c68416becdbc5d447efe71
[occt.git] / src / IVtkOCC / IVtkOCC_ShapeMesher.cxx
1 // Created on: 2011-10-14 
2 // Created by: Roman KOZLOV
3 // Copyright (c) 2011-2014 OPEN CASCADE SAS 
4 // 
5 // This file is part of Open CASCADE Technology software library.
6 //
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
12 //
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
15
16 #include <Adaptor3d_IsoCurve.hxx>
17 #include <Bnd_Box.hxx>
18 #include <BRep_Tool.hxx>
19 #include <BRepMesh_IncrementalMesh.hxx>
20 #include <BRepBndLib.hxx>
21 #include <BRepMesh_DiscretFactory.hxx>
22 #include <BRepMesh_DiscretRoot.hxx>
23 #include <BRepTools.hxx>
24 #include <Hatch_Hatcher.hxx>
25 #include <GCPnts_QuasiUniformDeflection.hxx>
26 #include <GCPnts_TangentialDeflection.hxx>
27 #include <Geom_BezierSurface.hxx>
28 #include <Geom_BSplineSurface.hxx>
29 #include <Geom2dAdaptor_Curve.hxx>
30 #include <GeomAdaptor_Curve.hxx>
31 #include <gp_Dir2d.hxx>
32 #include <gp_Pnt2d.hxx>
33 #include <IVtkOCC_ShapeMesher.hxx>
34 #include <NCollection_Array1.hxx>
35 #include <Poly_Polygon3D.hxx>
36 #include <Poly_PolygonOnTriangulation.hxx>
37 #include <Poly_Triangulation.hxx>
38 #include <Precision.hxx>
39 #include <Prs3d.hxx>
40 #include <Prs3d_Drawer.hxx>
41 #include <Quantity_Length.hxx>
42 #include <Standard_ErrorHandler.hxx>
43 #include <TColgp_SequenceOfPnt2d.hxx>
44 #include <TColStd_Array1OfReal.hxx>
45 #include <TopExp.hxx>
46 #include <TopExp_Explorer.hxx>
47
48 IMPLEMENT_STANDARD_RTTIEXT(IVtkOCC_ShapeMesher,IVtk_IShapeMesher)
49
50 // Handle implementation
51
52
53 //================================================================
54 // Function : internalBuild
55 // Purpose  : 
56 //================================================================
57 void IVtkOCC_ShapeMesher::internalBuild()
58 {
59   // TODO: do we need any protection here so as not to triangualte
60   // the shape twice??? This can be done e.g. by checking if
61   // triangulation exists for TopoDS_Shape..
62   meshShape();
63
64   // Free vertices and free edges should always be shown.
65   // Shared edges are needed in WF representation only.
66   // TODO: how to filter free edges at visualization level????
67   addFreeVertices();
68   addEdges();
69
70   // Build wireframe points and cells (lines for isolines)
71   addWireFrameFaces();
72
73   // Build shaded representation (based on Poly_Triangulation)
74   addShadedFaces();
75 }
76
77 //================================================================
78 // Function : GetShapeObj
79 // Purpose  : 
80 //================================================================
81 const IVtkOCC_Shape::Handle IVtkOCC_ShapeMesher::GetShapeObj() const
82 {
83   return (IVtkOCC_Shape::Handle::DownCast(myShapeObj));
84 }
85
86 //================================================================
87 // Function : GetDeflection
88 // Purpose  : Returns absolute deflection used by this algorithm.
89 //================================================================
90 Standard_Real IVtkOCC_ShapeMesher::GetDeflection() const
91 {
92   if (myDeflection < Precision::Confusion()) // if not yet initialized
93   {
94     Handle(Prs3d_Drawer) aDefDrawer = new Prs3d_Drawer();
95     aDefDrawer->SetTypeOfDeflection (Aspect_TOD_RELATIVE);
96     aDefDrawer->SetDeviationCoefficient (GetDeviationCoeff());
97     myDeflection = Prs3d::GetDeflection (GetShapeObj()->GetShape(), aDefDrawer);
98   }
99
100   return myDeflection;
101 }
102
103 //================================================================
104 // Function : meshShape
105 // Purpose  : 
106 //================================================================
107 void IVtkOCC_ShapeMesher::meshShape()
108 {
109   const TopoDS_Shape& anOcctShape = GetShapeObj()->GetShape();
110   if (anOcctShape.IsNull())
111   {
112     return;
113   }
114
115   //Clean triangulation before compute incremental mesh
116   BRepTools::Clean (anOcctShape);
117
118   //Compute triangulation
119   Standard_Real aDeflection = GetDeflection();
120   if (aDeflection < Precision::Confusion())
121   {
122     return;
123   }
124
125   try
126   {
127     OCC_CATCH_SIGNALS
128
129     Handle(BRepMesh_DiscretRoot) anAlgo;
130     anAlgo = BRepMesh_DiscretFactory::Get().Discret (anOcctShape,
131                                                      aDeflection,
132                                                      GetDeviationAngle());
133     if (!anAlgo.IsNull())
134     {
135       anAlgo->Perform();
136     }
137   }
138   catch (Standard_Failure)
139   { }
140 }
141
142 //================================================================
143 // Function : addFreeVertices
144 // Purpose  : 
145 //================================================================
146 void IVtkOCC_ShapeMesher::addFreeVertices()
147 {
148   TopTools_IndexedDataMapOfShapeListOfShape aVertexMap;
149   TopExp::MapShapesAndAncestors (GetShapeObj()->GetShape(),
150                                  TopAbs_VERTEX,
151                                  TopAbs_EDGE,
152                                  aVertexMap);
153
154   Standard_Integer aVertNum = aVertexMap.Extent();
155   IVtk_MeshType aType;
156   for (Standard_Integer anIt = 1; anIt <= aVertNum; anIt++)
157   {
158     if (aVertexMap.FindFromIndex(anIt).IsEmpty())
159     {
160       aType = MT_FreeVertex;
161     }
162     else
163     {
164       aType = MT_SharedVertex;
165     }
166     const TopoDS_Vertex& aVertex = TopoDS::Vertex (aVertexMap.FindKey (anIt));
167     addVertex (aVertex, GetShapeObj()->GetSubShapeId (aVertex), aType);
168   }
169 }
170
171 //================================================================
172 // Function : addEdges
173 // Purpose  : 
174 //================================================================
175 void IVtkOCC_ShapeMesher::addEdges()
176 {
177   TopTools_IndexedDataMapOfShapeListOfShape anEdgesMap;
178   TopExp::MapShapesAndAncestors (GetShapeObj()->GetShape(), 
179                                  TopAbs_EDGE, 
180                                  TopAbs_FACE,
181                                  anEdgesMap);
182
183   int aNbFaces;
184   IVtk_MeshType aType;
185   myEdgesTypes.Clear();
186
187   TopExp_Explorer anEdgeIter (GetShapeObj()->GetShape(), TopAbs_EDGE);
188   for (; anEdgeIter.More(); anEdgeIter.Next())
189   {
190     const TopoDS_Edge& anOcctEdge = TopoDS::Edge (anEdgeIter.Current());
191     aNbFaces = anEdgesMap.FindFromKey (anOcctEdge).Extent();
192     if (aNbFaces == 0)
193     {
194       aType = MT_FreeEdge;
195     }
196     else if (aNbFaces == 1)
197     {
198       aType = MT_BoundaryEdge;
199     }
200     else
201     {
202       aType = MT_SharedEdge;
203     }
204     addEdge (anOcctEdge, GetShapeObj()->GetSubShapeId (anOcctEdge), aType);
205     myEdgesTypes.Bind (anOcctEdge, aType);
206   }
207 }
208
209 //================================================================
210 // Function : addWireFrameFaces
211 // Purpose  : 
212 //================================================================
213 void IVtkOCC_ShapeMesher::addWireFrameFaces()
214 {
215   // Check the deflection value once for all faces
216   if (GetDeflection() < Precision::Confusion())
217   {
218     return;
219   }
220
221   TopExp_Explorer aFaceIter (GetShapeObj()->GetShape(), TopAbs_FACE);
222   for (; aFaceIter.More(); aFaceIter.Next())
223   {
224     const TopoDS_Face& anOcctFace = TopoDS::Face (aFaceIter.Current());
225     try
226     {
227       OCC_CATCH_SIGNALS
228       addWFFace (anOcctFace, 
229                  GetShapeObj()->GetSubShapeId (anOcctFace));
230     }
231     catch (Standard_Failure)
232     { }
233   }
234 }
235
236 //================================================================
237 // Function : addShadedFaces
238 // Purpose  : 
239 //================================================================
240 void IVtkOCC_ShapeMesher::addShadedFaces()
241 {
242   TopExp_Explorer aFaceIter (GetShapeObj()->GetShape(), TopAbs_FACE);
243   for (; aFaceIter.More(); aFaceIter.Next())
244   {
245     const TopoDS_Face& anOcctFace = TopoDS::Face (aFaceIter.Current());
246     addShadedFace (anOcctFace,
247                    GetShapeObj()->GetSubShapeId (anOcctFace));
248   }
249 }
250
251 //================================================================
252 // Function : addVertex
253 // Purpose  : 
254 //================================================================
255 void IVtkOCC_ShapeMesher::addVertex (const TopoDS_Vertex& theVertex,
256                                      const IVtk_IdType    theShapeId,
257                                      const IVtk_MeshType  theMeshType)
258 {
259   if (theVertex.IsNull())
260   {
261     return;
262   }
263
264   gp_Pnt aPnt3d = BRep_Tool::Pnt (theVertex);
265
266   IVtk_PointId anId = 
267     myShapeData->InsertCoordinate (aPnt3d.X(), aPnt3d.Y(), aPnt3d.Z());
268   myShapeData->InsertVertex (theShapeId, anId, theMeshType);
269
270 }
271
272 //================================================================
273 // Function : processPolyline
274 // Purpose  : 
275 //================================================================
276 void IVtkOCC_ShapeMesher::processPolyline (Standard_Integer          theNbNodes,
277                                       const TColgp_Array1OfPnt&      thePoints,
278                                       const TColStd_Array1OfInteger& thePointIds,
279                                       const IVtk_IdType              theOcctId,
280                                       bool                           theNoTransform,
281                                       gp_Trsf                        theTransformation,
282                                       const IVtk_MeshType            theMeshType)
283 {
284   if (theNbNodes < 2)
285   {
286     return;
287   }
288
289   IVtk_PointIdList aPolyPointIds;
290
291   IVtk_PointId anId;
292   for (Standard_Integer aJ = 0; aJ < theNbNodes; aJ++)
293   {
294     Standard_Integer aPntId = thePointIds (aJ + 1);
295     gp_Pnt point = thePoints (aPntId);
296
297     if (!theNoTransform)
298     {
299       // Apply the transformation to points
300       point.Transform (theTransformation);
301     }
302
303     anId = myShapeData->InsertCoordinate (point.X(), point.Y(), point.Z());
304     aPolyPointIds.Append (anId);
305   }
306
307   myShapeData->InsertLine (theOcctId, &aPolyPointIds, theMeshType);
308 }
309
310 //================================================================
311 // Function : addEdge
312 // Purpose  : 
313 //================================================================
314 void IVtkOCC_ShapeMesher::addEdge (const TopoDS_Edge&  theEdge,
315                                    const IVtk_IdType   theShapeId,
316                                    const IVtk_MeshType theMeshType)
317 {
318   if (theEdge.IsNull() || BRep_Tool::Degenerated (theEdge))
319   {
320     return;
321   }
322
323   // Two discrete representations of an OCCT edge are possible:
324   // 1. Polygon on trinagulation - holds Ids of points
325   // contained in Poly_Triangulation object
326   Handle(Poly_PolygonOnTriangulation) aPolyOnTriangulation;
327   Handle(Poly_Triangulation) aTriangulation;
328   TopLoc_Location aLocation;
329   BRep_Tool::PolygonOnTriangulation (theEdge,
330                                      aPolyOnTriangulation,
331                                      aTriangulation,
332                                      aLocation,
333                                      1);
334
335   // 2. 3D polygon - holds 3D points
336   Handle(Poly_Polygon3D) aPoly3d;
337   if (aPolyOnTriangulation.IsNull())
338   {
339     aPoly3d = BRep_Tool::Polygon3D (theEdge, aLocation);
340   }
341
342   if (aPoly3d.IsNull() && aPolyOnTriangulation.IsNull())
343   {
344     return;
345   }
346
347   // Handle a non-identity transofmation applied to the edge
348   gp_Trsf anEdgeTransf;
349   bool noTransform = true;
350   if (!aLocation.IsIdentity())
351   {
352     noTransform = false;
353     anEdgeTransf = aLocation.Transformation();
354   }
355
356   if (!aPoly3d.IsNull())
357   {
358     Standard_Integer aNbNodes = aPoly3d->NbNodes();
359     const TColgp_Array1OfPnt& aPoints = aPoly3d->Nodes();
360     TColStd_Array1OfInteger aPointIds (1, aNbNodes);
361
362     for (Standard_Integer anI = 1; anI <= aNbNodes; anI++)
363     {
364       aPointIds.SetValue (anI, anI);
365     }
366
367     processPolyline (aNbNodes,
368                      aPoints,
369                      aPointIds,
370                      theShapeId,
371                      noTransform,
372                      anEdgeTransf,
373                      theMeshType);
374   }
375   else
376   {
377     Standard_Integer aNbNodes = aPolyOnTriangulation->NbNodes();
378     const TColStd_Array1OfInteger& aPointIds = aPolyOnTriangulation->Nodes();
379     const TColgp_Array1OfPnt& aPoints = aTriangulation->Nodes();
380
381     processPolyline (aNbNodes,
382                      aPoints,
383                      aPointIds,
384                      theShapeId,
385                      noTransform,
386                      anEdgeTransf,
387                      theMeshType);
388   }
389 }
390
391
392 //================================================================
393 // Function : FindLimits
394 // Purpose  : Static internal function, finds parametrical limits of the curve.
395 //! @param [in] theCurve 3D curve adaptor used to retrieve the curve geometry
396 //! @param [in] theLimit maximum allowed absolute parameter value
397 //! @param [out] theFirst minimum parameter value for the curve
398 //! @param [out] theLast maximum parameter value for the curve
399 //================================================================
400 static void FindLimits (const Adaptor3d_Curve& theCurve,
401                         const Standard_Real&   theLimit,
402                         Standard_Real&         theFirst,
403                         Standard_Real&         theLast)
404 {
405   theFirst = Max(theCurve.FirstParameter(), theFirst);
406   theLast  = Min(theCurve.LastParameter(), theLast);
407   Standard_Boolean isFirstInf = Precision::IsNegativeInfinite (theFirst);
408   Standard_Boolean isLastInf  = Precision::IsPositiveInfinite (theLast);
409
410   if (isFirstInf || isLastInf)
411   {
412     gp_Pnt aP1, aP2;
413     Standard_Real aDelta = 1;
414     if (isFirstInf && isLastInf)
415     {
416       do
417       {
418         aDelta *= 2;
419         theFirst = - aDelta;
420         theLast  =   aDelta;
421         theCurve.D0 (theFirst, aP1);
422         theCurve.D0 (theLast, aP2);
423       } while (aP1.Distance (aP2) < theLimit);
424     }
425     else if (isFirstInf)
426     {
427       theCurve.D0 (theLast, aP2);
428       do {
429         aDelta *= 2;
430         theFirst = theLast - aDelta;
431         theCurve.D0 (theFirst, aP1);
432       } while (aP1.Distance(aP2) < theLimit);
433     }
434     else if (isLastInf)
435     {
436       theCurve.D0 (theFirst, aP1);
437       do
438       {
439         aDelta *= 2;
440         theLast = theFirst + aDelta;
441         theCurve.D0 (theLast, aP2);
442       } while (aP1.Distance (aP2) < theLimit);
443     }
444   }
445 }
446
447 //================================================================
448 // Function : FindLimits
449 // Purpose  :Static helper function, builds a discrete representation
450 //! (sequence of points) for the given curve.
451 //!
452 //! @param [in] theCurve 3D curve adaptor used to retrieve the curve geometry
453 //! @param [in] theDeflection absolute deflection value
454 //! @param [in] theAngle deviation angle value
455 //! @param [in] theU1 minimal curve parameter value
456 //! @param [in] theU2 maximal curve parameter value
457 //! @param [out] thePoints the container for generated polyline
458 //================================================================
459 static void DrawCurve (Adaptor3d_Curve&      theCurve,
460                        const Quantity_Length theDeflection,
461                        const Standard_Real   theAngle,
462                        const Standard_Real   theU1,
463                        const Standard_Real   theU2,
464                        IVtk_Polyline&        thePoints)
465 {
466   switch (theCurve.GetType())
467   {
468   case GeomAbs_Line:
469     {
470       gp_Pnt aPnt = theCurve.Value(theU1);
471       thePoints.Append (aPnt);
472
473       aPnt = theCurve.Value(0.5 * (theU1 + theU2));
474       thePoints.Append (aPnt);
475
476       aPnt = theCurve.Value (theU2);
477       thePoints.Append(aPnt);
478    }
479     break;
480   default:
481     {
482       Standard_Integer aNbInter = theCurve.NbIntervals(GeomAbs_C1);
483       Standard_Integer anI, aJ;
484       TColStd_Array1OfReal aParams(1, aNbInter+1);
485       theCurve.Intervals(aParams, GeomAbs_C1);
486       Standard_Real anU1, anU2;
487       Standard_Integer NumberOfPoints;
488
489       for (aJ = 1; aJ <= aNbInter; aJ++)
490       {
491         anU1 = aParams (aJ); anU2 = aParams (aJ + 1);
492         if (anU2 > anU1 && anU1 < anU2)
493         {
494           anU1 = Max(anU1, anU1);
495           anU2 = Min(anU2, anU2);
496
497           GCPnts_TangentialDeflection anAlgo (theCurve, anU1, anU2, theAngle, theDeflection);
498           NumberOfPoints = anAlgo.NbPoints();
499
500           if (NumberOfPoints > 0)
501           {
502             for (anI = 1; anI < NumberOfPoints; anI++)
503             {
504               thePoints.Append(anAlgo.Value (anI));
505             }
506             if (aJ == aNbInter)
507             {
508               thePoints.Append (anAlgo.Value (NumberOfPoints));
509             }
510           }
511         }
512       }
513     }
514   }
515 }
516
517 //================================================================
518 // Function : buildIsoLines
519 // Purpose  : 
520 //================================================================
521 void IVtkOCC_ShapeMesher::buildIsoLines (const Handle(BRepAdaptor_HSurface)& theFace,
522                                          const Standard_Boolean          theIsDrawUIso,
523                                          const Standard_Boolean          theIsDrawVIso,
524                                          const Standard_Integer          theNBUiso,
525                                          const Standard_Integer          theNBViso,
526                                          IVtk_PolylineList&              thePolylines)
527 {
528   Standard_Real anUF, anUL, aVF, aVL;
529   anUF = theFace->FirstUParameter();
530   anUL = theFace->LastUParameter();
531   aVF = theFace->FirstVParameter();
532   aVL = theFace->LastVParameter();
533
534   // Restrict maximal parameter value
535   // in OCCT it's 5e+5 by default
536   const Standard_Real aLimit = 5e+5;
537
538   // compute bounds of the restriction
539   Standard_Real anUMin, anUMax, aVMin, aVMax;
540   Standard_Integer anI;
541
542   anUMin = Max (anUF, -aLimit);
543   anUMax = Min (anUL, aLimit);
544   aVMin = Max (aVF, -aLimit);
545   aVMax = Min (aVL, aLimit);
546
547   // update min max for the hatcher.
548   gp_Pnt2d aP1,aP2;
549   gp_Pnt aDummyPnt;
550
551   Standard_Real aDdefle = Max (anUMax - anUMin, aVMax - aVMin) * GetDeviationCoeff();
552   TColgp_SequenceOfPnt2d aTabPoints;
553
554   anUMin = aVMin = 1.e100;
555   anUMax = aVMax = -1.e100;
556
557   // Process the edges
558   TopExp_Explorer aToolRst;
559   TopoDS_Face aTopoFace (((BRepAdaptor_Surface*)&(theFace->Surface()))->Face());
560   for (aToolRst.Init (aTopoFace, TopAbs_EDGE); aToolRst.More(); aToolRst.Next())
561   {
562     TopAbs_Orientation anOrient = aToolRst.Current().Orientation();
563     // Skip INTERNAL and EXTERNAL edges
564     if (anOrient == TopAbs_FORWARD || anOrient == TopAbs_REVERSED)
565     {
566       Standard_Real anU1, anU2;
567       const Handle(Geom2d_Curve)& aCurve = 
568         BRep_Tool::CurveOnSurface (TopoDS::Edge (aToolRst.Current()),
569                                    aTopoFace,
570                                    anU1, anU2);
571       if (aCurve.IsNull())
572       {
573         continue;
574       }
575
576       Geom2dAdaptor_Curve aRCurve;
577       aRCurve.Load (aCurve, anU1, anU2);
578       if (aRCurve.GetType() != GeomAbs_Line)
579       {
580         GCPnts_QuasiUniformDeflection aUDP(aRCurve, aDdefle);
581         if (aUDP.IsDone())
582         {
583           Standard_Integer NumberOfPoints = aUDP.NbPoints();
584           if ( NumberOfPoints >= 2 )
585           {
586             aDummyPnt = aUDP.Value (1);
587             aP2.SetCoord (aDummyPnt.X(), aDummyPnt.Y());
588             anUMin = Min (aP2.X(), anUMin);
589             anUMax = Max (aP2.X(), anUMax);
590             aVMin = Min (aP2.Y(), aVMin);
591             aVMax = Max (aP2.Y(), aVMax);
592             for (anI = 2; anI <= NumberOfPoints; anI++)
593             {
594               aP1 = aP2;
595               aDummyPnt = aUDP.Value (anI);
596               aP2.SetCoord (aDummyPnt.X(), aDummyPnt.Y());
597               anUMin = Min(aP2.X(), anUMin);
598               anUMax = Max(aP2.X(), anUMax);
599               aVMin = Min(aP2.Y(), aVMin);
600               aVMax = Max(aP2.Y(), aVMax);
601
602               if(anOrient == TopAbs_FORWARD )
603               {
604                 //isobuild.Trim(P1,P2);
605                 aTabPoints.Append (aP1);
606                 aTabPoints.Append (aP2);
607               }
608               else
609               {
610                 //isobuild.Trim(P2,P1);
611                 aTabPoints.Append (aP2);
612                 aTabPoints.Append (aP1);
613               }
614             }
615           }
616         }
617         else
618         {
619           cout << "Cannot evaluate curve on surface"<<endl;
620         }
621       }
622       else
623       {
624         anU1 = aRCurve.FirstParameter();
625         anU2 = aRCurve.LastParameter();
626         // MSV 17.08.06 OCC13144: U2 occured less than U1, to overcome it
627         // ensure that distance U2-U1 is not greater than aLimit*2,
628         // if greater then choose an origin and use aLimit to define
629         // U1 and U2 anew
630         Standard_Real aOrigin = 0.;
631         if (!Precision::IsNegativeInfinite(anU1) || !Precision::IsPositiveInfinite (anU2))
632         {
633           if (Precision::IsNegativeInfinite (anU1))
634           {
635             aOrigin = anU2 - aLimit;
636           }
637           else if (Precision::IsPositiveInfinite (anU2))
638           {
639             aOrigin = anU1 + aLimit;
640           }
641           else
642           {
643             aOrigin = (anU1 + anU2) * 0.5;
644           }
645         }
646
647         anU1 = Max (aOrigin - aLimit, anU1);
648         anU2 = Min (aOrigin + aLimit, anU2);
649         aP1 = aRCurve.Value (anU1);
650         aP2 = aRCurve.Value (anU2);
651         anUMin = Min(aP1.X(), anUMin);
652         anUMax = Max(aP1.X(), anUMax);
653         aVMin = Min(aP1.Y(), aVMin);
654         aVMax = Max(aP1.Y(), aVMax);
655         anUMin = Min(aP2.X(), anUMin);
656         anUMax = Max(aP2.X(), anUMax);
657         aVMin = Min(aP2.Y(), aVMin);
658         aVMax = Max(aP2.Y(), aVMax);
659         if(anOrient == TopAbs_FORWARD )
660         {
661          // isobuild.Trim(P1,P2);
662           aTabPoints.Append (aP1);
663           aTabPoints.Append (aP2);
664         }
665         else
666         {
667           //isobuild.Trim(P2,P1);
668           aTabPoints.Append (aP2);
669           aTabPoints.Append (aP1);
670         }
671       }
672     }
673   }
674
675   // load the isos
676   const Standard_Real anIntersectionTolerance = 1.e-5;
677   Hatch_Hatcher anIsoBuild (anIntersectionTolerance, Standard_True );
678
679   Standard_Boolean isUClosed = theFace->IsUClosed();
680   Standard_Boolean isVClosed = theFace->IsVClosed();
681
682   if (!isUClosed)
683   {
684     anUMin = anUMin + (anUMax - anUMin) / 1000.0;
685     anUMax = anUMax - (anUMax - anUMin) /1000.0;
686   }
687
688   if (!isVClosed)
689   {
690     aVMin = aVMin + (aVMax - aVMin) /1000.0;
691     aVMax = aVMax - (aVMax - aVMin) /1000.0;
692   }
693
694   if (theIsDrawUIso)
695   {
696     if (theNBUiso > 0)
697     {
698       isUClosed = Standard_False;
699       Standard_Real aDu= isUClosed ? (anUMax - anUMin) / theNBUiso : (anUMax - anUMin) / (1 + theNBUiso);
700       for (anI = 1; anI <= theNBUiso; anI++)
701       {
702         anIsoBuild.AddXLine (anUMin + aDu*anI);
703       }
704     }
705   }
706   if (theIsDrawVIso)
707   {
708     if (theNBViso > 0)
709     {
710       isVClosed = Standard_False;
711       Standard_Real aDv= isVClosed ? (aVMax - aVMin) / theNBViso : (aVMax - aVMin) / (1 + theNBViso);
712       for (anI = 1; anI <= theNBViso; anI++)
713       {
714         anIsoBuild.AddYLine (aVMin + aDv*anI);
715       }
716     }
717   }
718
719   Standard_Integer aLength = aTabPoints.Length();
720   for (anI = 1; anI <= aLength; anI += 2)
721   {
722     anIsoBuild.Trim (aTabPoints (anI),aTabPoints (anI + 1));
723   }
724
725   // Create the polylines for isos
726   Adaptor3d_IsoCurve anIso;
727   anIso.Load(theFace);
728   Handle(Geom_Curve) aBCurve;
729   const BRepAdaptor_Surface& aBSurf = *(BRepAdaptor_Surface*)&(theFace->Surface());
730   GeomAbs_SurfaceType aType = theFace->GetType();
731
732   Standard_Integer aNumberOfLines = anIsoBuild.NbLines();
733   Handle(Geom_Surface) aGeomSurf;
734   if (aType == GeomAbs_BezierSurface)
735   {
736     aGeomSurf = aBSurf.Bezier();
737   }
738   else if (aType == GeomAbs_BSplineSurface)
739   {
740     aGeomSurf = aBSurf.BSpline();
741   }
742
743   Standard_Real aDeflection = GetDeflection();
744   Standard_Real anAngle     = GetDeviationAngle();
745   for (anI = 1; anI <= aNumberOfLines; anI++)
746   {
747     Standard_Integer aNumberOfIntervals = anIsoBuild.NbIntervals(anI);
748     Standard_Real aCoord = anIsoBuild.Coordinate(anI);
749     for (Standard_Integer aJ = 1; aJ <= aNumberOfIntervals; aJ++)
750     {
751       Standard_Real aB1 = anIsoBuild.Start (anI, aJ);
752       Standard_Real aB2 = anIsoBuild.End(anI, aJ);
753
754       if (!aGeomSurf.IsNull())
755       {
756         if (anIsoBuild.IsXLine (anI))
757         {
758           aBCurve = aGeomSurf->UIso (aCoord);
759         }
760         else
761         {
762           aBCurve = aGeomSurf->VIso (aCoord);
763         }
764
765         GeomAdaptor_Curve aGeomCurve (aBCurve);
766         FindLimits (aGeomCurve, aLimit, aB1, aB2);
767         if (aB2 - aB1 > Precision::Confusion())
768         {
769           IVtk_Polyline aPoints;
770           DrawCurve (aGeomCurve, aDeflection, anAngle, aB1, aB2, aPoints);
771           thePolylines.Append (aPoints);
772         }
773       }
774       else
775       {
776         if (anIsoBuild.IsXLine (anI))
777         {
778           anIso.Load (GeomAbs_IsoU, aCoord, aB1, aB2);
779         }
780         else
781         {
782           anIso.Load (GeomAbs_IsoV, aCoord, aB1, aB2);
783         }
784         FindLimits (anIso, aLimit, aB1, aB2);
785         if (aB2 - aB1>Precision::Confusion())
786         {
787           IVtk_Polyline aPoints;
788           DrawCurve (anIso, aDeflection, anAngle, aB1, aB2, aPoints);
789           thePolylines.Append (aPoints);
790         }
791       }
792     }
793   }
794 }
795
796 //================================================================
797 // Function : addWFFace
798 // Purpose  : 
799 //================================================================
800 void IVtkOCC_ShapeMesher::addWFFace (const TopoDS_Face& theFace,
801                                      const IVtk_IdType  theShapeId)
802 {
803   if (theFace.IsNull())
804   {
805     return;
806   }
807
808   TopoDS_Face aFaceToMesh = theFace;
809   aFaceToMesh.Orientation (TopAbs_FORWARD);
810
811   // The code that builds wireframe representation for a TopoDS_Face
812   // has been adapted from some OCCT 6.5.1 methods:
813   // - Prs3d_WFShape::Add()
814   // - StdPrs_WFDeflectionRestrictedFace::Add()
815   // - StdPrs_DeflectionCurve::Add()
816
817   // Add face's edges here but with the face ID
818   TopExp_Explorer anEdgeIter (aFaceToMesh, TopAbs_EDGE );
819   for (; anEdgeIter.More(); anEdgeIter.Next())
820   {
821     const TopoDS_Edge& anOcctEdge = TopoDS::Edge (anEdgeIter.Current());
822     addEdge (anOcctEdge, theShapeId, myEdgesTypes (anOcctEdge));
823   }
824
825   TopLoc_Location aLoc;
826   const Handle(Geom_Surface)& aGeomSurf = BRep_Tool::Surface (aFaceToMesh, aLoc);
827   if (aGeomSurf.IsNull())
828   {
829     return;
830   }
831
832   BRepAdaptor_Surface aSurf;
833   aSurf.Initialize (aFaceToMesh);
834   Handle(BRepAdaptor_HSurface) aSurfAdaptor = new BRepAdaptor_HSurface (aSurf);
835
836   IVtk_PolylineList aPolylines;
837   gp_Trsf aDummyTrsf;
838
839   // Building U isolines
840   // Introducing a local scope here to simplify variable naming
841   {
842     buildIsoLines (aSurfAdaptor,
843                    myNbIsos[0] != 0,
844                    Standard_False,
845                    myNbIsos[0],
846                    0,
847                    aPolylines);
848
849     IVtk_PolylineList::Iterator anIt (aPolylines);
850     for (; anIt.More(); anIt.Next())
851     {
852       const IVtk_Polyline& aPntSeq = anIt.Value();
853       Standard_Integer aNbNodes = aPntSeq.Length();
854       TColgp_Array1OfPnt aPoints (1, aNbNodes);
855       for (Standard_Integer aJ = 1; aJ <= aNbNodes; aJ++)
856       {
857         aPoints.SetValue (aJ, aPntSeq.Value(aJ));
858       }
859
860       TColStd_Array1OfInteger aPointIds (1, aNbNodes);
861       for (Standard_Integer anI = 1; anI <= aNbNodes; anI++)
862       {
863         aPointIds.SetValue (anI, anI);
864       }
865
866       processPolyline (aNbNodes,
867                        aPoints,
868                        aPointIds,
869                        theShapeId,
870                        true,
871                        aDummyTrsf, 
872                        MT_IsoLine);
873     }
874   }
875
876   // Building V isolines
877   {
878     aPolylines.Clear();
879     buildIsoLines (aSurfAdaptor,
880                    Standard_False,
881                    myNbIsos[1] != 0,
882                    0,
883                    myNbIsos[1],
884                    aPolylines);
885
886     IVtk_PolylineList::Iterator anIt (aPolylines);
887     for (; anIt.More(); anIt.Next())
888     {
889       const IVtk_Polyline& aPntSeq = anIt.Value();
890       Standard_Integer aNbNodes = aPntSeq.Length();
891       TColgp_Array1OfPnt aPoints (1, aNbNodes);
892       for (int aJ = 1; aJ <= aNbNodes; aJ++)
893       {
894         aPoints.SetValue (aJ, aPntSeq.Value (aJ));
895       }
896
897       TColStd_Array1OfInteger aPointIds (1, aNbNodes);
898       for (Standard_Integer anI = 1; anI <= aNbNodes; anI++)
899       {
900         aPointIds.SetValue (anI, anI);
901       }
902
903       processPolyline (aNbNodes,
904                        aPoints,
905                        aPointIds,
906                        theShapeId,
907                        true,
908                        aDummyTrsf, 
909                        MT_IsoLine);
910     }
911   }
912 }
913
914 //================================================================
915 // Function : addShadedFace
916 // Purpose  : 
917 //================================================================
918 void IVtkOCC_ShapeMesher::addShadedFace (const TopoDS_Face& theFace,
919                                          const IVtk_IdType  theShapeId)
920 {
921   if (theFace.IsNull())
922   {
923     return;
924   }
925
926   // Build triangulation of the face.
927   TopLoc_Location aLoc;
928   Handle(Poly_Triangulation) anOcctTriangulation = BRep_Tool::Triangulation (theFace, aLoc);
929   if (anOcctTriangulation.IsNull())
930   {
931     return;
932   }
933
934   gp_Trsf aPntTransform;
935   Standard_Boolean noTransform = Standard_True;
936   if (!aLoc.IsIdentity())
937   {
938     noTransform = Standard_False;
939     aPntTransform = aLoc.Transformation();
940   }
941
942   // Get triangulation points.
943   const TColgp_Array1OfPnt& aPoints = anOcctTriangulation->Nodes();
944   Standard_Integer aNbPoints = anOcctTriangulation->NbNodes();
945
946   // Keep inserted points id's of triangulation in an array.
947   NCollection_Array1<IVtk_PointId> aPointIds (1, aNbPoints);
948   IVtk_PointId anId;
949
950   Standard_Integer anI;
951   for (anI = 1; anI <= aNbPoints; anI++)
952   {
953     gp_Pnt aPoint = aPoints (anI);
954
955     if (!noTransform)
956     {
957       aPoint.Transform (aPntTransform);
958     }
959
960     // Add a point into output shape data and keep its id in the array.
961     anId = myShapeData->InsertCoordinate (aPoint.X(), aPoint.Y(), aPoint.Z());
962     aPointIds.SetValue (anI, anId);
963   }
964
965   // Create triangles on the created triangulation points.
966   const Poly_Array1OfTriangle& aTriangles = anOcctTriangulation->Triangles();
967   Standard_Integer aNbTriangles = anOcctTriangulation->NbTriangles();
968   Standard_Integer aN1, aN2, aN3;
969   for (anI = 1; anI <= aNbTriangles; anI++)
970   {
971     aTriangles(anI).Get (aN1, aN2, aN3); // get indexes of triangle's points
972     // Insert new triangle on these points into output shape data.
973     myShapeData->InsertTriangle (
974       theShapeId, aPointIds(aN1), aPointIds(aN2), aPointIds(aN3), MT_ShadedFace);
975   }
976 }