1 // Created by: Ekaterina SMIRNOVA
2 // Copyright (c) 2008-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
15 #include <BRepMesh_FastDiscretFace.hxx>
17 #include <BRepMesh_PairOfPolygon.hxx>
18 #include <BRepMesh_ShapeTool.hxx>
19 #include <Poly_PolygonOnTriangulation.hxx>
20 #include <Poly_Triangulation.hxx>
22 #include <BRepAdaptor_Surface.hxx>
23 #include <BRepAdaptor_HSurface.hxx>
24 #include <BRepAdaptor_Curve.hxx>
25 #include <Adaptor3d_IsoCurve.hxx>
26 #include <Adaptor3d_HCurve.hxx>
28 #include <BRep_ListIteratorOfListOfPointRepresentation.hxx>
29 #include <BRep_PointRepresentation.hxx>
30 #include <BRep_TVertex.hxx>
31 #include <BRep_Tool.hxx>
33 #include <GeomLib.hxx>
34 #include <Geom_Surface.hxx>
35 #include <Geom_BSplineSurface.hxx>
36 #include <Geom_BezierSurface.hxx>
37 #include <GCPnts_TangentialDeflection.hxx>
38 #include <GCPnts_AbscissaPoint.hxx>
40 #include <Standard_ErrorHandler.hxx>
41 #include <Standard_Failure.hxx>
42 #include <TColStd_Array1OfReal.hxx>
43 #include <TColStd_ListOfInteger.hxx>
44 #include <TColStd_SequenceOfReal.hxx>
45 #include <TColStd_Array1OfInteger.hxx>
46 #include <TColStd_HArray1OfReal.hxx>
47 #include <TColgp_Array1OfPnt2d.hxx>
48 #include <TopTools_DataMapOfShapeReal.hxx>
50 #include <TopExp_Explorer.hxx>
52 #include <TopoDS_Vertex.hxx>
55 #include <NCollection_Map.hxx>
56 #include <Bnd_Box2d.hxx>
60 //#define DEBUG_MESH "mesh.tcl"
66 IMPLEMENT_STANDARD_RTTIEXT(BRepMesh_FastDiscretFace,Standard_Transient)
68 static Standard_Real FUN_CalcAverageDUV(TColStd_Array1OfReal& P, const Standard_Integer PLen)
70 Standard_Integer i, j, n = 0;
71 Standard_Real p, result = 0.;
73 for(i = 1; i <= PLen; i++)
76 for(j = i + 1; j <= PLen; j++)
88 p = Abs(P(i) - P(i-1));
96 return (n? (result / (Standard_Real) n) : -1.);
101 Standard_Real deflectionOfSegment (
102 const gp_Pnt& theFirstPoint,
103 const gp_Pnt& theLastPoint,
104 const gp_Pnt& theMidPoint)
106 // 23.03.2010 skl for OCC21645 - change precision for comparison
107 if (theFirstPoint.SquareDistance (theLastPoint) > Precision::SquareConfusion ())
109 gp_Lin aLin (theFirstPoint, gp_Dir (gp_Vec (theFirstPoint, theLastPoint)));
110 return aLin.Distance (theMidPoint);
113 return theFirstPoint.Distance (theMidPoint);
116 Standard_Boolean IsCompexSurface (const GeomAbs_SurfaceType theType)
119 theType != GeomAbs_Sphere &&
120 theType != GeomAbs_Cylinder &&
121 theType != GeomAbs_Cone &&
122 theType != GeomAbs_Torus);
125 //! Auxiliary class used to extract geometrical parameters of fixed TopoDS_Vertex.
130 DEFINE_STANDARD_ALLOC
132 FixedVExplorer(const TopoDS_Vertex& theVertex)
133 : myVertex(theVertex)
137 const TopoDS_Vertex& Vertex() const
142 Standard_Boolean IsSameUV() const
144 return Standard_False;
147 TopoDS_Vertex SameVertex() const
149 return TopoDS_Vertex();
154 return BRep_Tool::Pnt(myVertex);
159 void operator =(const FixedVExplorer& /*theOther*/)
164 const TopoDS_Vertex& myVertex;
167 void ComputeErrFactors (const Standard_Real theDeflection,
168 const Handle(Adaptor3d_HSurface)& theFace,
169 Standard_Real& theErrFactorU,
170 Standard_Real& theErrFactorV)
172 theErrFactorU = theDeflection * 10.;
173 theErrFactorV = theDeflection * 10.;
175 switch (theFace->GetType ())
177 case GeomAbs_Cylinder:
183 case GeomAbs_SurfaceOfExtrusion:
184 case GeomAbs_SurfaceOfRevolution:
186 Handle (Adaptor3d_HCurve) aCurve = theFace->BasisCurve ();
187 if (aCurve->GetType () == GeomAbs_BSplineCurve && aCurve->Degree () > 2)
189 theErrFactorV /= (aCurve->Degree () * aCurve->NbKnots ());
193 case GeomAbs_BezierSurface:
195 if (theFace->UDegree () > 2)
197 theErrFactorU /= (theFace->UDegree ());
199 if (theFace->VDegree () > 2)
201 theErrFactorV /= (theFace->VDegree ());
205 case GeomAbs_BSplineSurface:
207 if (theFace->UDegree () > 2)
209 theErrFactorU /= (theFace->UDegree () * theFace->NbUKnots ());
211 if (theFace->VDegree () > 2)
213 theErrFactorV /= (theFace->VDegree () * theFace->NbVKnots ());
220 theErrFactorU = theErrFactorV = 1.;
224 void AdjustCellsCounts (const Handle(Adaptor3d_HSurface)& theFace,
225 const Standard_Integer theNbVertices,
226 Standard_Integer& theCellsCountU,
227 Standard_Integer& theCellsCountV)
229 const GeomAbs_SurfaceType aType = theFace->GetType ();
230 if (aType == GeomAbs_OtherSurface)
232 // fallback to the default behavior
233 theCellsCountU = theCellsCountV = -1;
237 Standard_Real aSqNbVert = theNbVertices;
238 if (aType == GeomAbs_Plane)
240 theCellsCountU = theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
242 else if (aType == GeomAbs_Cylinder || aType == GeomAbs_Cone)
244 theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
246 else if (aType == GeomAbs_SurfaceOfExtrusion || aType == GeomAbs_SurfaceOfRevolution)
248 Handle (Adaptor3d_HCurve) aCurve = theFace->BasisCurve ();
249 if (aCurve->GetType () == GeomAbs_Line ||
250 (aCurve->GetType () == GeomAbs_BSplineCurve && aCurve->Degree () < 2))
252 // planar, cylindrical, conical cases
253 if (aType == GeomAbs_SurfaceOfExtrusion)
254 theCellsCountU = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
256 theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
258 if (aType == GeomAbs_SurfaceOfExtrusion)
260 // V is always a line
261 theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
264 else if (aType == GeomAbs_BezierSurface || aType == GeomAbs_BSplineSurface)
266 if (theFace->UDegree () < 2)
268 theCellsCountU = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
270 if (theFace->VDegree () < 2)
272 theCellsCountV = (Standard_Integer)Ceiling (Pow (2, Log10 (aSqNbVert)));
276 theCellsCountU = Max (theCellsCountU, 2);
277 theCellsCountV = Max (theCellsCountV, 2);
282 //=======================================================================
283 //function : BRepMesh_FastDiscretFace
285 //=======================================================================
286 BRepMesh_FastDiscretFace::BRepMesh_FastDiscretFace(
287 const Standard_Real theAngle,
288 const Standard_Real theMinSize,
289 const Standard_Boolean isInternalVerticesMode,
290 const Standard_Boolean isControlSurfaceDeflection)
292 myInternalVerticesMode(isInternalVerticesMode),
293 myMinSize(theMinSize),
294 myIsControlSurfaceDeflection(isControlSurfaceDeflection)
298 //=======================================================================
301 //=======================================================================
302 void BRepMesh_FastDiscretFace::Perform(const Handle(BRepMesh_FaceAttribute)& theAttribute)
305 commitSurfaceTriangulation();
308 //=======================================================================
309 //function : initDataStructure
311 //=======================================================================
312 void BRepMesh_FastDiscretFace::initDataStructure()
314 const Standard_Real aTolU = myAttribute->ToleranceU();
315 const Standard_Real aTolV = myAttribute->ToleranceV();
316 const Standard_Real uCellSize = 14.0 * aTolU;
317 const Standard_Real vCellSize = 14.0 * aTolV;
319 const Standard_Real deltaX = myAttribute->GetDeltaX();
320 const Standard_Real deltaY = myAttribute->GetDeltaY();
322 Handle(NCollection_IncAllocator) aAllocator =
323 new NCollection_IncAllocator(BRepMesh::MEMORY_BLOCK_SIZE_HUGE);
324 myStructure = new BRepMesh_DataStructureOfDelaun(aAllocator);
325 myStructure->Data()->SetCellSize ( uCellSize / deltaX, vCellSize / deltaY);
326 myStructure->Data()->SetTolerance( aTolU / deltaX, aTolV / deltaY);
328 myAttribute->ChangeStructure() = myStructure;
329 myAttribute->ChangeSurfacePoints() = new BRepMesh::DMapOfIntegerPnt(1, aAllocator);
330 myAttribute->ChangeSurfaceVertices()= new BRepMesh::DMapOfVertexInteger(1, aAllocator);
332 // Check the necessity to fill the map of parameters
333 const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
334 GeomAbs_SurfaceType thetype = gFace->GetType();
335 const Standard_Boolean isBSpline = (thetype == GeomAbs_BezierSurface ||
336 thetype == GeomAbs_BSplineSurface);
337 const Standard_Boolean useUVParam = (thetype == GeomAbs_Torus || IsCompexSurface (thetype));
339 Handle(NCollection_BaseAllocator) aBaseAlloc = aAllocator;
340 myUParam.Clear (aBaseAlloc);
341 myVParam.Clear (aBaseAlloc);
343 // essai de determination de la longueur vraie:
344 // akm (bug OCC16) : We must calculate these measures in non-singular
345 // parts of face. Let`s try to compute average value of three
346 // (umin, (umin+umax)/2, umax), and respectively for v.
348 //Standard_Real longu = 0.0, longv = 0.0; //, last , first;
349 //gp_Pnt P11, P12, P21, P22, P31, P32;
350 BRepMesh::HVectorOfVertex& aBoundaryNodes = myAttribute->ChangeMeshNodes();
351 BRepMesh::VectorOfVertex::Iterator aNodesIt(*aBoundaryNodes);
352 for (; aNodesIt.More(); aNodesIt.Next())
354 BRepMesh_Vertex& aNode = aNodesIt.ChangeValue();
355 gp_XY aPnt2d = aNode.Coord();
359 myUParam.Add(aPnt2d.X());
360 myVParam.Add(aPnt2d.Y());
363 aNode.ChangeCoord() = myAttribute->Scale(aPnt2d, Standard_True);
364 myStructure->AddNode(aNode, Standard_True);
366 aBoundaryNodes.Nullify();
370 const Standard_Real aRange[2][2] = {
371 {myAttribute->GetUMin(), myAttribute->GetUMax()},
372 {myAttribute->GetVMin(), myAttribute->GetVMax()}
375 const GeomAbs_Shape aContinuity = GeomAbs_CN;
376 for (Standard_Integer i = 0; i < 2; ++i)
378 const Standard_Boolean isU = (i == 0);
379 const Standard_Integer aIntervalsNb = isU ?
380 gFace->NbUIntervals(aContinuity) :
381 gFace->NbVIntervals(aContinuity);
383 BRepMesh::IMapOfReal& aParams = isU ? myUParam : myVParam;
384 if (aIntervalsNb < aParams.Size())
387 TColStd_Array1OfReal aIntervals(1, aIntervalsNb + 1);
389 gFace->UIntervals(aIntervals, aContinuity);
391 gFace->VIntervals(aIntervals, aContinuity);
393 for (Standard_Integer j = 1; j <= aIntervals.Upper(); ++j)
395 const Standard_Real aParam = aIntervals(j);
396 if (aParam > aRange[i][0] && aParam < aRange[i][1])
402 ////////////////////////////////////////////////////////////
403 //add internal vertices after self-intersection check
404 if ( myInternalVerticesMode )
406 TopExp_Explorer anExplorer(myAttribute->Face(), TopAbs_VERTEX, TopAbs_EDGE);
407 for ( ; anExplorer.More(); anExplorer.Next() )
408 add(TopoDS::Vertex(anExplorer.Current()));
411 const BRepMesh::HDMapOfShapePairOfPolygon& aEdges = myAttribute->ChangeInternalEdges();
412 TopExp_Explorer aWireIt(myAttribute->Face(), TopAbs_WIRE);
413 for (; aWireIt.More(); aWireIt.Next())
415 TopExp_Explorer aEdgeIt(aWireIt.Current(), TopAbs_EDGE);
416 for (; aEdgeIt.More(); aEdgeIt.Next())
418 const TopoDS_Edge& aEdge = TopoDS::Edge(aEdgeIt.Current());
419 BRepMesh_PairOfPolygon aPair;
420 if (!aEdges->Find(aEdge, aPair))
423 TopAbs_Orientation aOri = aEdge.Orientation();
424 const Handle(Poly_PolygonOnTriangulation)& aPolygon =
425 aOri == TopAbs_REVERSED ? aPair.Last() : aPair.First();
427 const TColStd_Array1OfInteger& aIndices = aPolygon->Nodes();
428 const Standard_Integer aNodesNb = aPolygon->NbNodes();
430 Standard_Integer aPrevId = aIndices(1);
431 for (Standard_Integer i = 2; i <= aNodesNb; ++i)
433 const Standard_Integer aCurId = aIndices(i);
434 addLinkToMesh(aPrevId, aCurId, aOri);
441 //=======================================================================
442 //function : addLinkToMesh
444 //=======================================================================
445 void BRepMesh_FastDiscretFace::addLinkToMesh(
446 const Standard_Integer theFirstNodeId,
447 const Standard_Integer theLastNodeId,
448 const TopAbs_Orientation theOrientation)
450 if (theOrientation == TopAbs_FORWARD)
451 myStructure->AddLink(BRepMesh_Edge(theFirstNodeId, theLastNodeId, BRepMesh_Frontier));
452 else if (theOrientation == TopAbs_REVERSED)
453 myStructure->AddLink(BRepMesh_Edge(theLastNodeId, theFirstNodeId, BRepMesh_Frontier));
454 else if (theOrientation == TopAbs_INTERNAL)
455 myStructure->AddLink(BRepMesh_Edge(theFirstNodeId, theLastNodeId, BRepMesh_Fixed));
458 //=======================================================================
461 //=======================================================================
462 void BRepMesh_FastDiscretFace::add(const Handle(BRepMesh_FaceAttribute)& theAttribute)
464 if (!theAttribute->IsValid() || theAttribute->ChangeMeshNodes()->IsEmpty())
467 myAttribute = theAttribute;
470 BRepMesh::HIMapOfInteger& aVertexEdgeMap = myAttribute->ChangeVertexEdgeMap();
471 Standard_Integer nbVertices = aVertexEdgeMap->Extent();
472 BRepMesh::Array1OfInteger tabvert_corr(1, nbVertices);
473 for ( Standard_Integer i = 1; i <= nbVertices; ++i )
476 Handle (Adaptor3d_HSurface) aSurface (myAttribute->Surface ());
477 GeomAbs_SurfaceType aType = aSurface->GetType ();
479 while (aType == GeomAbs_OffsetSurface)
481 aSurface = aSurface->BasisSurface ();
482 aType = aSurface->GetType ();
485 // For better meshing performance we try to estimate the acceleration circles grid structure sizes:
486 // For each parametric direction (U, V) we estimate firstly an approximate distance between the future points -
487 // this estimation takes into account the required face deflection and the complexity of the face.
488 // Particularly, the complexity of the faces based on BSpline curves and surfaces requires much more points.
489 // At the same time, for planar faces and linear parts of the arbitrary surfaces usually no intermediate points
491 // The general idea for each parametric direction:
492 // cells_count = 2 ^ log10 ( estimated_points_count )
493 // For linear parametric direction we fall back to the initial vertex count:
494 // cells_count = 2 ^ log10 ( initial_vertex_count )
496 Standard_Real anErrFactorU, anErrFactorV;
497 ComputeErrFactors(myAttribute->GetDefFace (), aSurface, anErrFactorU, anErrFactorV);
499 Standard_Integer aCellsCountU, aCellsCountV;
500 if (aType == GeomAbs_Torus)
502 aCellsCountU = (Standard_Integer)Ceiling(Pow(2, Log10(
503 (myAttribute->GetUMax() - myAttribute->GetUMin()) / myAttribute->GetDeltaX())));
504 aCellsCountV = (Standard_Integer)Ceiling(Pow(2, Log10(
505 (myAttribute->GetVMax() - myAttribute->GetVMin()) / myAttribute->GetDeltaY())));
507 else if (aType == GeomAbs_Cylinder)
509 aCellsCountU = (Standard_Integer)Ceiling(Pow(2, Log10(
510 (myAttribute->GetUMax() - myAttribute->GetUMin()) / myAttribute->GetDeltaX() /
511 (myAttribute->GetVMax() - myAttribute->GetVMin()))));
512 aCellsCountV = (Standard_Integer)Ceiling(Pow(2, Log10(
513 (myAttribute->GetVMax() - myAttribute->GetVMin()) / anErrFactorV)));
517 aCellsCountU = (Standard_Integer)Ceiling(Pow(2, Log10(
518 (myAttribute->GetUMax() - myAttribute->GetUMin()) / myAttribute->GetDeltaX() / anErrFactorU)));
519 aCellsCountV = (Standard_Integer)Ceiling(Pow(2, Log10(
520 (myAttribute->GetVMax() - myAttribute->GetVMin()) / myAttribute->GetDeltaY() / anErrFactorV)));
523 AdjustCellsCounts(aSurface, nbVertices, aCellsCountU, aCellsCountV);
525 BRepMesh_Delaun trigu(myStructure, tabvert_corr, aCellsCountU, aCellsCountV);
527 //removed all free edges from triangulation
528 const Standard_Integer nbLinks = myStructure->NbLinks();
529 for( Standard_Integer i = 1; i <= nbLinks; i++ )
531 if( myStructure->ElementsConnectedTo(i).Extent() < 1 )
533 BRepMesh_Edge& anEdge = (BRepMesh_Edge&)trigu.GetEdge(i);
534 if ( anEdge.Movability() == BRepMesh_Deleted )
537 anEdge.SetMovability(BRepMesh_Free);
538 myStructure->RemoveLink(i);
542 Standard_Boolean rajout =
543 (aType == GeomAbs_Sphere || aType == GeomAbs_Torus);
545 // Check the necessity to fill the map of parameters
546 const Standard_Boolean useUVParam = (aType == GeomAbs_Torus ||
547 aType == GeomAbs_BezierSurface ||
548 aType == GeomAbs_BSplineSurface);
550 const Standard_Real umax = myAttribute->GetUMax();
551 const Standard_Real umin = myAttribute->GetUMin();
552 const Standard_Real vmax = myAttribute->GetVMax();
553 const Standard_Real vmin = myAttribute->GetVMin();
555 Standard_Boolean isaline =
556 ((umax - umin) < Precision::PConfusion() ||
557 (vmax - vmin) < Precision::PConfusion());
559 Standard_Real aDef = -1;
560 if ( !isaline && myStructure->ElementsOfDomain().Extent() > 0 )
564 // compute maximal deflection
565 aDef = control(trigu, Standard_True);
566 rajout = (aDef > myAttribute->GetDefFace() || aDef < 0.);
568 if (aType != GeomAbs_Plane)
570 if (!rajout && useUVParam)
572 rajout = (myVParam.Extent() > 2 &&
573 (aSurface->IsUClosed() || aSurface->IsVClosed()));
578 insertInternalVertices(trigu);
580 //control internal points
581 if (myIsControlSurfaceDeflection)
582 aDef = control(trigu, Standard_False);
587 //modify myStructure back
588 BRepMesh::HVectorOfVertex& aMeshNodes = myStructure->Data()->ChangeVertices();
589 for ( Standard_Integer i = 1; i <= myStructure->NbNodes(); ++i )
591 BRepMesh_Vertex& aNode = aMeshNodes->ChangeValue(i - 1);
592 aNode.ChangeCoord() = myAttribute->Scale(aNode.Coord(), Standard_False);
594 const BRepMesh::ListOfInteger& alist = myStructure->LinksConnectedTo(i);
595 // Register internal nodes used in triangulation
596 if (!alist.IsEmpty() && aVertexEdgeMap->FindIndex(i) == 0)
597 aVertexEdgeMap->Add(i);
601 myAttribute->SetDefFace(aDef);
604 //=======================================================================
605 //function : addVerticesToMesh
607 //=======================================================================
608 Standard_Boolean BRepMesh_FastDiscretFace::addVerticesToMesh(
609 const BRepMesh::ListOfVertex& theVertices,
610 BRepMesh_Delaun& theMeshBuilder)
612 if (theVertices.IsEmpty())
613 return Standard_False;
615 BRepMesh::Array1OfVertexOfDelaun aArrayOfNewVertices(1, theVertices.Extent());
616 BRepMesh::ListOfVertex::Iterator aVertexIt(theVertices);
617 for (Standard_Integer aVertexId = 0; aVertexIt.More(); aVertexIt.Next())
618 aArrayOfNewVertices(++aVertexId) = aVertexIt.Value();
620 theMeshBuilder.AddVertices(aArrayOfNewVertices);
621 return Standard_True;
624 //=======================================================================
625 //function : filterParameters
627 //=======================================================================
628 static void filterParameters(const BRepMesh::IMapOfReal& theParams,
629 const Standard_Real theMinDist,
630 const Standard_Real theFilterDist,
631 BRepMesh::SequenceOfReal& theResult)
633 // Sort sequence of parameters
634 const Standard_Integer anInitLen = theParams.Extent();
636 TColStd_Array1OfReal aParamArray(1, anInitLen);
638 for (j = 1; j <= anInitLen; j++)
639 aParamArray(j) = theParams(j);
641 std::sort (aParamArray.begin(), aParamArray.end());
643 // mandatory pre-filtering using the first (minimal) filter value
644 Standard_Integer aParamLength = 1;
645 for (j = 2; j <= anInitLen; j++)
647 if ((aParamArray(j)-aParamArray(aParamLength)) > theMinDist)
649 if (++aParamLength < j)
650 aParamArray(aParamLength) = aParamArray(j);
654 //perform filtering on series
655 Standard_Real aLastAdded, aLastCandidate;
656 Standard_Boolean isCandidateDefined = Standard_False;
657 aLastAdded = aParamArray(1);
658 aLastCandidate = aLastAdded;
659 theResult.Append(aLastAdded);
661 for(j=2; j < aParamLength; j++)
663 Standard_Real aVal = aParamArray(j);
664 if(aVal-aLastAdded > theFilterDist)
667 if(isCandidateDefined) {
668 aLastAdded = aLastCandidate;
669 isCandidateDefined = Standard_False;
676 theResult.Append(aLastAdded);
680 aLastCandidate = aVal;
681 isCandidateDefined = Standard_True;
683 theResult.Append(aParamArray(aParamLength));
686 //=======================================================================
687 //function : insertInternalVertices
689 //=======================================================================
690 void BRepMesh_FastDiscretFace::insertInternalVertices(BRepMesh_Delaun& theMeshBuilder)
692 Handle(NCollection_IncAllocator) anAlloc = new NCollection_IncAllocator;
693 BRepMesh::ListOfVertex aNewVertices(anAlloc);
694 const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
695 switch (gFace->GetType())
698 insertInternalVerticesSphere(aNewVertices);
701 case GeomAbs_Cylinder:
702 insertInternalVerticesCylinder(aNewVertices);
706 insertInternalVerticesCone(aNewVertices);
710 insertInternalVerticesTorus(aNewVertices);
714 insertInternalVerticesOther(aNewVertices);
718 addVerticesToMesh(aNewVertices, theMeshBuilder);
721 //=======================================================================
722 //function : insertInternalVerticesSphere
724 //=======================================================================
725 void BRepMesh_FastDiscretFace::insertInternalVerticesSphere(
726 BRepMesh::ListOfVertex& theNewVertices)
728 Standard_Real aRange[] = {
729 myAttribute->GetVMin(), myAttribute->GetVMax(),
730 myAttribute->GetUMin(), myAttribute->GetUMax()
733 gp_Sphere aSphere = myAttribute->Surface()->Sphere();
735 // Calculate parameters for iteration in V direction
736 Standard_Real aStep = 0.7 * GCPnts_TangentialDeflection::ArcAngularStep(
737 aSphere.Radius(), myAttribute->GetDefFace(), myAngle, myMinSize);
739 Standard_Real aDd[2] = {aStep, aStep};
740 Standard_Real aPasMax[2] = {0., 0.};
741 for (Standard_Integer i = 0; i < 2; ++i)
743 const Standard_Real aMax = aRange[2 * i + 1];
744 const Standard_Real aDiff = aMax - aRange[2 * i + 0];
745 aDd[i] = aDiff / ((Standard_Integer)(aDiff / aDd[i]) + 1);
746 aPasMax[i] = aMax - Precision::PConfusion();
749 const Standard_Real aHalfDu = aDd[1] * 0.5;
750 Standard_Boolean Shift = Standard_False;
751 Standard_Real aPasV = aRange[0] + aDd[0];
752 for (; aPasV < aPasMax[0]; aPasV += aDd[0])
755 const Standard_Real d = (Shift) ? aHalfDu : 0.;
756 Standard_Real aPasU = aRange[2] + d;
757 for (; aPasU < aPasMax[1]; aPasU += aDd[1])
759 tryToInsertAnalyticVertex(gp_Pnt2d(aPasU, aPasV), aSphere, theNewVertices);
764 //=======================================================================
765 //function : insertInternalVerticesCylinder
767 //=======================================================================
768 void BRepMesh_FastDiscretFace::insertInternalVerticesCylinder(
769 BRepMesh::ListOfVertex& theNewVertices)
771 const Standard_Real umax = myAttribute->GetUMax();
772 const Standard_Real umin = myAttribute->GetUMin();
773 const Standard_Real vmax = myAttribute->GetVMax();
774 const Standard_Real vmin = myAttribute->GetVMin();
776 gp_Cylinder aCylinder = myAttribute->Surface()->Cylinder();
777 const Standard_Real aRadius = aCylinder.Radius();
779 Standard_Integer nbU = 0;
780 Standard_Integer nbV = 0;
781 const Standard_Real su = umax - umin;
782 const Standard_Real sv = vmax - vmin;
783 const Standard_Real aArcLen = su * aRadius;
784 if (aArcLen > myAttribute->GetDefFace ())
786 // Calculate parameters for iteration in U direction
787 const Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep (
788 aRadius, myAttribute->GetDefFace (), myAngle, myMinSize);
789 nbU = (Standard_Integer)(su / Du);
791 // Calculate parameters for iteration in V direction
792 const Standard_Real aDv = nbU*sv / aArcLen;
793 // Protection against overflow during casting to int in case
794 // of long cylinder with small radius.
795 nbV = aDv > static_cast<Standard_Real> (IntegerLast ()) ?
796 0 : (Standard_Integer)(aDv);
797 nbV = Min (nbV, 100 * nbU);
800 const Standard_Real Du = su / (nbU + 1);
801 const Standard_Real Dv = sv / (nbV + 1);
803 Standard_Real pasu, pasv, pasvmax = vmax - Dv*0.5, pasumax = umax - Du*0.5;
804 for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv)
806 for (pasu = umin + Du; pasu < pasumax; pasu += Du)
808 tryToInsertAnalyticVertex(gp_Pnt2d(pasu, pasv), aCylinder, theNewVertices);
813 //=======================================================================
814 //function : insertInternalVerticesCone
816 //=======================================================================
817 void BRepMesh_FastDiscretFace::insertInternalVerticesCone(
818 BRepMesh::ListOfVertex& theNewVertices)
820 const Standard_Real umax = myAttribute->GetUMax();
821 const Standard_Real umin = myAttribute->GetUMin();
822 const Standard_Real vmax = myAttribute->GetVMax();
823 const Standard_Real vmin = myAttribute->GetVMin();
825 gp_Cone aCone = myAttribute->Surface()->Cone();
826 Standard_Real RefR = aCone.RefRadius();
827 Standard_Real SAng = aCone.SemiAngle();
828 Standard_Real aRadius = Max(Abs(RefR + vmin*Sin(SAng)), Abs(RefR + vmax*Sin(SAng)));
830 Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep(
831 aRadius, myAttribute->GetDefFace(), myAngle, myMinSize);
833 Standard_Real Dv, pasu, pasv;
834 Standard_Integer nbU = (Standard_Integer)((umax - umin) / Du);
835 Standard_Integer nbV = (Standard_Integer)(nbU*(vmax - vmin) / ((umax - umin)*aRadius));
836 Du = (umax - umin) / (nbU + 1);
837 Dv = (vmax - vmin) / (nbV + 1);
839 Standard_Real pasvmax = vmax - Dv*0.5, pasumax = umax - Du*0.5;
840 for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv)
842 for (pasu = umin + Du; pasu < pasumax; pasu += Du)
844 tryToInsertAnalyticVertex(gp_Pnt2d(pasu, pasv), aCone, theNewVertices);
849 //=======================================================================
850 //function : insertInternalVerticesTorus
852 //=======================================================================
853 void BRepMesh_FastDiscretFace::insertInternalVerticesTorus(
854 BRepMesh::ListOfVertex& theNewVertices)
856 const Standard_Real umax = myAttribute->GetUMax();
857 const Standard_Real umin = myAttribute->GetUMin();
858 const Standard_Real vmax = myAttribute->GetVMax();
859 const Standard_Real vmin = myAttribute->GetVMin();
860 const Standard_Real deltaX = myAttribute->GetDeltaX();
861 const Standard_Real deltaY = myAttribute->GetDeltaY();
862 const Standard_Real aDefFace = myAttribute->GetDefFace();
864 gp_Torus T = myAttribute->Surface()->Torus();
866 Standard_Boolean insert;
867 Standard_Integer i, j, ParamULength, ParamVLength;
868 Standard_Real pp, pasu, pasv;
869 Standard_Real r = T.MinorRadius(), R = T.MajorRadius();
871 BRepMesh::SequenceOfReal ParamU, ParamV;
873 Standard_Real oldDv = GCPnts_TangentialDeflection::ArcAngularStep(
874 r, aDefFace, myAngle, myMinSize);
876 Standard_Real Dv = 0.9*oldDv; //TWOTHIRD * oldDv;
880 Standard_Integer nbV = Max((Standard_Integer)((vmax - vmin) / Dv), 2);
881 Dv = (vmax - vmin) / (nbV + 1);
882 Standard_Real ru = R + r;
885 Du = GCPnts_TangentialDeflection::ArcAngularStep(
886 ru, aDefFace, myAngle, myMinSize);
888 Standard_Real aa = sqrt(Du*Du + oldDv*oldDv);
889 if (aa < gp::Resolution())
891 Du *= Min(oldDv, Du) / aa;
895 Standard_Integer nbU = Max((Standard_Integer)((umax - umin) / Du), 2);
896 nbU = Max(nbU, (int)(nbV*(umax - umin)*R / ((vmax - vmin)*r) / 5.));
897 Du = (umax - umin) / (nbU + 1);
901 // As the points of edges are returned.
902 // in this case, the points are not representative.
904 //-- Choose DeltaX and DeltaY so that to avoid skipping points on the grid
905 for (i = 0; i <= nbU; i++) ParamU.Append(umin + i* Du);
910 // Number of mapped U parameters
911 const Standard_Integer LenU = myUParam.Extent();
912 // Fill array of U parameters
913 TColStd_Array1OfReal Up(1, LenU);
914 for (j = 1; j <= LenU; j++) Up(j) = myUParam(j);
916 // Calculate DU, leave array of parameters
917 Standard_Real aDU = FUN_CalcAverageDUV(Up, LenU);
918 aDU = Max(aDU, Abs(umax - umin) / (Standard_Real)nbU / 2.);
919 Standard_Real dUstd = Abs(umax - umin) / (Standard_Real)LenU;
920 if (aDU > dUstd) dUstd = aDU;
922 for (j = 1; j <= LenU; j++)
925 insert = Standard_True;
926 ParamULength = ParamU.Length();
927 for (i = 1; i <= ParamULength && insert; i++)
929 insert = (Abs(ParamU.Value(i) - pp) > (0.5*dUstd));
931 if (insert) ParamU.Append(pp);
936 // Number of mapped V parameters
937 const Standard_Integer LenV = myVParam.Extent();
938 // Fill array of V parameters
939 TColStd_Array1OfReal Vp(1, LenV);
940 for (j = 1; j <= LenV; j++) Vp(j) = myVParam(j);
941 // Calculate DV, sort array of parameters
942 Standard_Real aDV = FUN_CalcAverageDUV(Vp, LenV);
943 aDV = Max(aDV, Abs(vmax - vmin) / (Standard_Real)nbV / 2.);
945 Standard_Real dVstd = Abs(vmax - vmin) / (Standard_Real)LenV;
946 if (aDV > dVstd) dVstd = aDV;
948 for (j = 1; j <= LenV; j++)
952 insert = Standard_True;
953 ParamVLength = ParamV.Length();
954 for (i = 1; i <= ParamVLength && insert; i++)
956 insert = (Abs(ParamV.Value(i) - pp) > (dVstd*2. / 3.));
958 if (insert) ParamV.Append(pp);
961 Standard_Integer Lu = ParamU.Length(), Lv = ParamV.Length();
962 Standard_Real uminnew = umin + deltaY*0.1;
963 Standard_Real vminnew = vmin + deltaX*0.1;
964 Standard_Real umaxnew = umax - deltaY*0.1;
965 Standard_Real vmaxnew = vmax - deltaX*0.1;
967 for (i = 1; i <= Lu; i++)
969 pasu = ParamU.Value(i);
970 if (pasu >= uminnew && pasu < umaxnew)
972 for (j = 1; j <= Lv; j++)
974 pasv = ParamV.Value(j);
975 if (pasv >= vminnew && pasv < vmaxnew)
977 tryToInsertAnalyticVertex(gp_Pnt2d(pasu, pasv), T, theNewVertices);
984 //=======================================================================
985 //function : insertInternalVerticesOther
987 //=======================================================================
988 void BRepMesh_FastDiscretFace::insertInternalVerticesOther(
989 BRepMesh::ListOfVertex& theNewVertices)
991 const Standard_Real aRange[2][2] = {
992 { myAttribute->GetUMax(), myAttribute->GetUMin() },
993 { myAttribute->GetVMax(), myAttribute->GetVMin() }
996 const Standard_Real aDelta[2] = {
997 myAttribute->GetDeltaX(),
998 myAttribute->GetDeltaY()
1001 const Standard_Real aDefFace = myAttribute->GetDefFace();
1002 const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
1004 Handle(NCollection_IncAllocator) anAlloc = new NCollection_IncAllocator;
1005 BRepMesh::SequenceOfReal aParams[2] = { BRepMesh::SequenceOfReal(anAlloc),
1006 BRepMesh::SequenceOfReal(anAlloc) };
1007 for (Standard_Integer i = 0; i < 2; ++i)
1009 Standard_Boolean isU = (i == 0);
1010 Standard_Real aRes = isU ?
1011 gFace->UResolution(aDefFace) :
1012 gFace->VResolution(aDefFace);
1014 // Sort and filter sequence of parameters
1015 Standard_Real aMinDiff = Precision::PConfusion();
1017 aMinDiff /= aDelta[i];
1019 aMinDiff = Max(myMinSize, aMinDiff);
1021 Standard_Real aRangeDiff = aRange[i][0] - aRange[i][1];
1022 Standard_Real aDiffMaxLim = 0.1 * aRangeDiff;
1023 Standard_Real aDiffMinLim = Max(0.005 * aRangeDiff, 2. * aRes);
1024 Standard_Real aDiff = Max(myMinSize, Min(aDiffMaxLim, aDiffMinLim));
1025 filterParameters(isU ? myUParam : myVParam, aMinDiff, aDiff, aParams[i]);
1028 // check intermediate isolines
1029 Handle (Geom_Surface) aSurface = gFace->ChangeSurface ().Surface ().Surface ();
1031 BRepMesh::MapOfReal aParamsToRemove[2] = { BRepMesh::MapOfReal(1, anAlloc),
1032 BRepMesh::MapOfReal(1, anAlloc) };
1033 BRepMesh::MapOfReal aParamsForbiddenToRemove[2] = { BRepMesh::MapOfReal(1, anAlloc),
1034 BRepMesh::MapOfReal(1, anAlloc) };
1036 // insert additional points where it is needed to conform criteria.
1037 // precision for normals calculation
1038 const Standard_Real aNormPrec = Precision::Confusion();
1039 for (Standard_Integer k = 0; k < 2; ++k)
1041 const Standard_Integer aOtherIndex = (k + 1) % 2;
1042 BRepMesh::SequenceOfReal& aParams1 = aParams[k];
1043 BRepMesh::SequenceOfReal& aParams2 = aParams[aOtherIndex];
1044 const Standard_Boolean isU = (k == 0);
1045 for (Standard_Integer i = 2; i < aParams1.Length(); ++i)
1047 const Standard_Real aParam1 = aParams1(i);
1048 GeomAdaptor_Curve aIso(isU ?
1049 aSurface->UIso(aParam1) : aSurface->VIso(aParam1));
1051 Standard_Real aPrevParam2 = aParams2(1);
1054 aIso.D1(aPrevParam2, aPrevPnt2, aPrevVec2);
1055 for (Standard_Integer j = 2; j <= aParams2.Length();)
1057 Standard_Real aParam2 = aParams2(j);
1060 aIso.D1(aParam2, aPnt2, aVec2);
1062 Standard_Real aMidParam = 0.5 * (aPrevParam2 + aParam2);
1063 gp_Pnt aMidPnt = aIso.Value(aMidParam);
1065 Standard_Real aDist = deflectionOfSegment(aPrevPnt2, aPnt2, aMidPnt);
1066 if (aDist > aDefFace && aDist > myMinSize)
1069 aParams2.InsertBefore(j, aMidParam);
1072 //put regular grig for normals
1073 gp_Pnt2d aStPnt1, aStPnt2;
1076 aStPnt1 = gp_Pnt2d(aParam1, aPrevParam2);
1077 aStPnt2 = gp_Pnt2d(aParam1, aMidParam);
1081 aStPnt1 = gp_Pnt2d(aPrevParam2, aParam1);
1082 aStPnt2 = gp_Pnt2d(aMidParam, aParam1);
1085 gp_Dir N1(0, 0, 1), N2(0, 0, 1);
1086 Standard_Integer aSt1 = GeomLib::NormEstim(aSurface, aStPnt1, aNormPrec, N1);
1087 Standard_Integer aSt2 = GeomLib::NormEstim(aSurface, aStPnt2, aNormPrec, N2);
1089 const Standard_Real aAngle = N2.Angle(N1);
1090 if (aSt1 < 1 && aSt2 < 1 && aAngle > myAngle)
1092 const Standard_Real aLen = GCPnts_AbscissaPoint::Length(
1093 aIso, aPrevParam2, aMidParam, aDefFace);
1095 if (aLen > myMinSize)
1098 aParams2.InsertBefore(j, aMidParam);
1102 aPrevParam2 = aParam2;
1110 #ifdef DEBUG_InsertInternal
1111 // output numbers of parameters along U and V
1112 cout << "aParams: " << aParams[0].Length() << " " << aParams[1].Length() << endl;
1114 // try to reduce number of points evaluating of isolines sampling
1115 for (Standard_Integer k = 0; k < 2; ++k)
1117 const Standard_Integer aOtherIndex = (k + 1) % 2;
1118 BRepMesh::SequenceOfReal& aParams1 = aParams[k];
1119 BRepMesh::SequenceOfReal& aParams2 = aParams[aOtherIndex];
1120 const Standard_Boolean isU = (k == 0);
1121 BRepMesh::MapOfReal& aToRemove2 = aParamsToRemove[aOtherIndex];
1122 BRepMesh::MapOfReal& aForbiddenToRemove1 = aParamsForbiddenToRemove[k];
1123 BRepMesh::MapOfReal& aForbiddenToRemove2 = aParamsForbiddenToRemove[aOtherIndex];
1124 for (Standard_Integer i = 2; i < aParams1.Length(); ++i)
1126 const Standard_Real aParam1 = aParams1(i);
1127 GeomAdaptor_Curve aIso(isU ?
1128 aSurface->UIso (aParam1) : aSurface->VIso (aParam1));
1129 #ifdef DEBUG_InsertInternal
1130 // write polyline containing initial parameters to the file
1132 ofstream ff(DEBUG_InsertInternal, std::ios_base::app);
1133 ff << "polyline " << (k == 0 ? "u" : "v") << i << " ";
1134 for (Standard_Integer j = 1; j <= aParams2.Length(); j++)
1137 aIso.D0(aParams2(j), aP);
1138 ff << aP.X() << " " << aP.Y() << " " << aP.Z() << " ";
1144 Standard_Real aPrevParam2 = aParams2(1);
1147 aIso.D1 (aPrevParam2, aPrevPnt2, aPrevVec2);
1148 for (Standard_Integer j = 2; j <= aParams2.Length();)
1150 Standard_Real aParam2 = aParams2(j);
1153 aIso.D1 (aParam2, aPnt2, aVec2);
1155 // Here we should leave at least 3 parameters as far as
1156 // we must have at least one parameter related to surface
1157 // internals in order to prevent movement of triangle body
1158 // outside the surface in case of highly curved ones, e.g.
1160 if (aParams2.Length() > 3 && j < aParams2.Length())
1162 // Remove too dense points
1163 const Standard_Real aNextParam = aParams2(j + 1);
1166 aIso.D1(aNextParam, aNextPnt, aNextVec);
1168 // Lets check current parameter.
1169 // If it fits deflection, we can remove it.
1170 Standard_Real aDist = deflectionOfSegment(aPrevPnt2, aNextPnt, aPnt2);
1171 if (aDist < aDefFace)
1173 // Lets check parameters for angular deflection.
1174 if (aPrevVec2.SquareMagnitude() > gp::Resolution() &&
1175 aNextVec.SquareMagnitude() > gp::Resolution() &&
1176 aPrevVec2.Angle(aNextVec) < myAngle)
1178 // For current Iso line we can remove this parameter.
1179 #ifdef DEBUG_InsertInternal
1180 // write point of removed parameter
1182 ofstream ff(DEBUG_InsertInternal, std::ios_base::app);
1183 ff << "point " << (k == 0 ? "u" : "v") << i << "r_" << j << " "
1184 << aPnt2.X() << " " << aPnt2.Y() << " " << aPnt2.Z() << endl;
1187 aToRemove2.Add(aParam2);
1188 aPrevParam2 = aNextParam;
1189 aPrevPnt2 = aNextPnt;
1190 aPrevVec2 = aNextVec;
1195 // We have found a place on the surface refusing
1196 // removement of this parameter.
1197 #ifdef DEBUG_InsertInternal
1198 // write point of forbidden to remove parameter
1200 ofstream ff(DEBUG_InsertInternal, std::ios_base::app);
1201 ff << "vertex " << (k == 0 ? "u" : "v") << i << "f_" << j << " "
1202 << aPnt2.X() << " " << aPnt2.Y() << " " << aPnt2.Z() << endl;
1205 aForbiddenToRemove1.Add(aParam1);
1206 aForbiddenToRemove2.Add(aParam2);
1210 aPrevParam2 = aParam2;
1217 // remove filtered out parameters
1218 for (Standard_Integer k = 0; k < 2; ++k)
1220 BRepMesh::SequenceOfReal& aParamsk = aParams[k];
1221 for (Standard_Integer i = 1; i <= aParamsk.Length();)
1223 const Standard_Real aParam = aParamsk.Value(i);
1224 if (aParamsToRemove[k].Contains(aParam) &&
1225 !aParamsForbiddenToRemove[k].Contains(aParam))
1233 #ifdef DEBUG_InsertInternal
1234 // write polylines containing remaining parameters
1236 ofstream ff(DEBUG_InsertInternal, std::ios_base::app);
1237 for (Standard_Integer k = 0; k < 2; ++k)
1239 for (Standard_Integer i = 1; i <= aParams[k].Length(); i++)
1241 ff << "polyline " << (k == 0 ? "uo" : "vo") << i << " ";
1242 for (Standard_Integer j = 1; j <= aParams[1 - k].Length(); j++)
1246 gFace->D0(aParams[k](i), aParams[1 - k](j), aP);
1248 gFace->D0(aParams[1 - k](j), aParams[k](i), aP);
1249 ff << aP.X() << " " << aP.Y() << " " << aP.Z() << " ";
1257 // insert nodes of the regular grid
1258 const BRepMesh::HClassifier& aClassifier = myAttribute->ChangeClassifier();
1259 for (Standard_Integer i = 1; i <= aParams[0].Length(); ++i)
1261 const Standard_Real aParam1 = aParams[0].Value (i);
1262 for (Standard_Integer j = 1; j <= aParams[1].Length(); ++j)
1264 const Standard_Real aParam2 = aParams[1].Value (j);
1265 gp_Pnt2d aPnt2d(aParam1, aParam2);
1267 // Classify intersection point
1268 if (aClassifier->Perform(aPnt2d) == TopAbs_IN)
1271 gFace->D0(aPnt2d.X(), aPnt2d.Y(), aPnt);
1272 insertVertex(aPnt, aPnt2d.Coord(), theNewVertices);
1278 //=======================================================================
1279 //function : checkDeflectionAndInsert
1281 //=======================================================================
1282 Standard_Boolean BRepMesh_FastDiscretFace::checkDeflectionAndInsert(
1283 const gp_Pnt& thePnt3d,
1285 const Standard_Boolean isDeflectionCheckOnly,
1286 const Standard_Real theTriangleDeflection,
1287 const Standard_Real theFaceDeflection,
1288 const BRepMesh_CircleTool& theCircleTool,
1289 BRepMesh::ListOfVertex& theVertices,
1290 Standard_Real& theMaxTriangleDeflection,
1291 const Handle(NCollection_IncAllocator)& theTempAlloc)
1293 if (theTriangleDeflection > theMaxTriangleDeflection)
1294 theMaxTriangleDeflection = theTriangleDeflection;
1296 if (theTriangleDeflection < theFaceDeflection)
1297 return Standard_True;
1299 if (myMinSize > Precision::Confusion())
1301 // Iterator in the list of indexes of circles containing the node
1302 BRepMesh::ListOfInteger& aCirclesList =
1303 const_cast<BRepMesh_CircleTool&>(theCircleTool).Select(
1304 myAttribute->Scale(theUV, Standard_True));
1306 BRepMesh::MapOfInteger aUsedNodes(10, theTempAlloc);
1307 BRepMesh::ListOfInteger::Iterator aCircleIt(aCirclesList);
1308 for (; aCircleIt.More(); aCircleIt.Next())
1310 const BRepMesh_Triangle& aTriangle =
1311 myStructure->GetElement(aCircleIt.Value());
1313 Standard_Integer aNodes[3];
1314 myStructure->ElementNodes(aTriangle, aNodes);
1316 for (Standard_Integer i = 0; i < 3; ++i)
1318 const Standard_Integer aNodeId = aNodes[i];
1319 if (aUsedNodes.Contains(aNodeId))
1322 aUsedNodes.Add(aNodeId);
1323 const BRepMesh_Vertex& aNode = myStructure->GetNode(aNodeId);
1324 const gp_Pnt& aPoint = myAttribute->GetPoint(aNode);
1326 if (thePnt3d.SquareDistance(aPoint) < myMinSize * myMinSize)
1327 return Standard_True;
1332 if (isDeflectionCheckOnly)
1333 return Standard_False;
1335 insertVertex(thePnt3d, theUV, theVertices);
1336 return Standard_True;
1339 //=======================================================================
1340 //function : control
1342 //=======================================================================
1343 Standard_Real BRepMesh_FastDiscretFace::control(
1344 BRepMesh_Delaun& theTrigu,
1345 const Standard_Boolean theIsFirst)
1347 Standard_Integer aTrianglesNb = myStructure->ElementsOfDomain().Extent();
1348 if (aTrianglesNb < 1)
1351 //IMPORTANT: Constants used in calculations
1352 const Standard_Real MinimalArea2d = 1.e-9;
1353 const Standard_Real MinimalSqLength3d = 1.e-12;
1354 const Standard_Real aSqDefFace = myAttribute->GetDefFace() * myAttribute->GetDefFace();
1356 const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
1358 Handle(Geom_Surface) aBSpline;
1359 const GeomAbs_SurfaceType aSurfType = gFace->GetType ();
1360 if (IsCompexSurface (aSurfType) && aSurfType != GeomAbs_SurfaceOfExtrusion)
1361 aBSpline = gFace->ChangeSurface ().Surface().Surface();
1363 Handle(NCollection_IncAllocator) anAlloc =
1364 new NCollection_IncAllocator(BRepMesh::MEMORY_BLOCK_SIZE_HUGE);
1365 NCollection_DataMap<Standard_Integer, gp_Dir> aNorMap(1, anAlloc);
1366 BRepMesh::MapOfIntegerInteger aStatMap(1, anAlloc);
1367 NCollection_Map<BRepMesh_Edge> aCouples(3 * aTrianglesNb, anAlloc);
1368 const BRepMesh_CircleTool& aCircles = theTrigu.Circles();
1370 // Perform refinement passes
1371 // Define the number of iterations
1372 Standard_Integer aIterationsNb = 11;
1373 const Standard_Integer aPassesNb = (theIsFirst ? 1 : aIterationsNb);
1374 // Initialize stop condition
1375 Standard_Real aMaxSqDef = -1.;
1376 Standard_Integer aPass = 1, aInsertedNb = 1;
1377 Standard_Boolean isAllDegenerated = Standard_False;
1378 Handle(NCollection_IncAllocator) aTempAlloc =
1379 new NCollection_IncAllocator(BRepMesh::MEMORY_BLOCK_SIZE_HUGE);
1380 for (; aPass <= aPassesNb && aInsertedNb && !isAllDegenerated; ++aPass)
1382 aTempAlloc->Reset(Standard_False);
1383 BRepMesh::ListOfVertex aNewVertices(aTempAlloc);
1385 // Reset stop condition
1388 isAllDegenerated = Standard_True;
1390 aTrianglesNb = myStructure->ElementsOfDomain().Extent();
1391 if (aTrianglesNb < 1)
1394 // Iterate on current triangles
1395 const BRepMesh::MapOfInteger& aTriangles = myStructure->ElementsOfDomain();
1396 BRepMesh::MapOfInteger::Iterator aTriangleIt(aTriangles);
1397 for (; aTriangleIt.More(); aTriangleIt.Next())
1399 const Standard_Integer aTriangleId = aTriangleIt.Key();
1400 const BRepMesh_Triangle& aCurrentTriangle = myStructure->GetElement(aTriangleId);
1402 if (aCurrentTriangle.Movability() == BRepMesh_Deleted)
1405 Standard_Integer v[3];
1406 myStructure->ElementNodes(aCurrentTriangle, v);
1408 Standard_Integer e[3];
1409 Standard_Boolean o[3];
1410 aCurrentTriangle.Edges(e, o);
1414 Standard_Boolean m[3];
1415 for (Standard_Integer i = 0; i < 3; ++i)
1417 m[i] = (myStructure->GetLink(e[i]).Movability() == BRepMesh_Frontier);
1419 const BRepMesh_Vertex& aVertex = myStructure->GetNode(v[i]);
1420 xy[i] = myAttribute->Scale(aVertex.Coord(), Standard_False);
1421 p [i] = myAttribute->GetPoint(aVertex).Coord();
1425 Standard_Boolean isDegeneratedTri = Standard_False;
1426 for (Standard_Integer i = 0; i < 3 && !isDegeneratedTri; ++i)
1428 aLinkVec[i] = p[(i + 1) % 3] - p[i];
1429 isDegeneratedTri = (aLinkVec[i].SquareModulus() < MinimalSqLength3d);
1432 if (isDegeneratedTri)
1435 isAllDegenerated = Standard_False;
1437 // Check triangle area in 2d
1438 if (Abs((xy[1]-xy[0])^(xy[2]-xy[1])) < MinimalArea2d)
1441 // Check triangle normal
1443 Standard_Real aSqDef = -1.;
1444 Standard_Boolean isSkipped = Standard_False;
1445 gp_XYZ normal(aLinkVec[0] ^ aLinkVec[1]);
1446 if (normal.SquareModulus () < gp::Resolution())
1451 // Check deflection at triangle centroid
1452 gp_XY aCenter2d = (xy[0] + xy[1] + xy[2]) / 3.0;
1453 gFace->D0(aCenter2d.X(), aCenter2d.Y(), pDef);
1454 aSqDef = Abs(normal * (pDef.XYZ() - p[0]));
1457 isSkipped = !checkDeflectionAndInsert(pDef, aCenter2d, theIsFirst,
1458 aSqDef, aSqDefFace, aCircles, aNewVertices, aMaxSqDef, aTempAlloc);
1463 // Check deflection at triangle links
1464 for (Standard_Integer i = 0; i < 3 && !isSkipped; ++i)
1466 if (m[i]) // is a boundary
1469 // Check if this link was already processed
1470 if (aCouples.Add(myStructure->GetLink(e[i])))
1472 // Check deflection on edge 1
1473 Standard_Integer j = (i + 1) % 3;
1474 gp_XY mi2d = (xy[i] + xy[j]) * 0.5;
1475 gFace->D0(mi2d.X(), mi2d.Y(), pDef);
1476 gp_Lin aLin(p[i], gp_Vec(p[i], p[j]));
1477 aSqDef = aLin.SquareDistance(pDef);
1479 isSkipped = !checkDeflectionAndInsert(pDef, mi2d, theIsFirst,
1480 aSqDef, aSqDefFace, aCircles, aNewVertices, aMaxSqDef, aTempAlloc);
1487 //check normal on bsplines
1488 if (theIsFirst && !aBSpline.IsNull())
1490 gp_Dir N[3] = { gp::DZ(), gp::DZ(), gp::DZ() };
1491 Standard_Integer aSt[3];
1493 for (Standard_Integer i = 0; i < 3; ++i)
1495 if (aNorMap.IsBound(v[i]))
1497 aSt[i] = aStatMap.Find(v[i]);
1498 N[i] = aNorMap.Find(v[i]);
1502 aSt[i] = GeomLib::NormEstim(aBSpline, gp_Pnt2d(xy[i]), Precision::Confusion(), N[i]);
1503 aStatMap.Bind(v[i], aSt[i]);
1504 aNorMap.Bind(v[i], N[i]);
1508 Standard_Real aAngle[3];
1509 for (Standard_Integer i = 0; i < 3; ++i)
1510 aAngle[i] = N[(i + 1) % 3].Angle(N[i]);
1512 if (aSt[0] < 1 && aSt[1] < 1 && aSt[2] < 1)
1514 if (aAngle[0] > myAngle || aAngle[1] > myAngle || aAngle[2] > myAngle)
1527 // append to the file triangles in the form of polyline commands;
1528 // so the file can be sourced in draw to analyze triangles on each pass of the algorithm.
1530 ofstream ftt(DEBUG_MESH, std::ios_base::app);
1531 Standard_Integer nbElem = myStructure->NbElements();
1532 for (Standard_Integer i = 1; i <= nbElem; i++)
1534 const BRepMesh_Triangle& aTri = myStructure->GetElement(i);
1535 if (aTri.Movability() == BRepMesh_Deleted)
1537 Standard_Integer n[3];
1538 myStructure->ElementNodes(aTri, n);
1539 const BRepMesh_Vertex& aV1 = myStructure->GetNode(n[0]);
1540 const BRepMesh_Vertex& aV2 = myStructure->GetNode(n[1]);
1541 const BRepMesh_Vertex& aV3 = myStructure->GetNode(n[2]);
1542 const gp_Pnt& aP1 = myAttribute->GetPoint(aV1);
1543 const gp_Pnt& aP2 = myAttribute->GetPoint(aV2);
1544 const gp_Pnt& aP3 = myAttribute->GetPoint(aV3);
1545 ftt << "polyline t" << aPass << "_" << i << " "
1546 << aP1.X() << " " << aP1.Y() << " " << aP1.Z() << " "
1547 << aP2.X() << " " << aP2.Y() << " " << aP2.Z() << " "
1548 << aP3.X() << " " << aP3.Y() << " " << aP3.Z() << " "
1549 << aP1.X() << " " << aP1.Y() << " " << aP1.Z() << endl;
1551 // write points to insert on the current pass
1552 BRepMesh::ListOfVertex::Iterator it(aNewVertices);
1553 for (Standard_Integer i=1; it.More(); it.Next(), i++)
1555 const BRepMesh_Vertex& aVertex = it.Value();
1556 const gp_Pnt& aP = myAttribute->GetPoint(aVertex);
1557 ftt << "vertex vt" << aPass << "_" << i << " "
1558 << aP.X() << " " << aP.Y() << " " << aP.Z() << endl;
1562 if (addVerticesToMesh(aNewVertices, theTrigu))
1566 return (aMaxSqDef < 0) ? aMaxSqDef : Sqrt(aMaxSqDef);
1569 //=======================================================================
1572 //=======================================================================
1573 void BRepMesh_FastDiscretFace::add(const TopoDS_Vertex& theVertex)
1575 if (theVertex.Orientation() != TopAbs_INTERNAL)
1582 gp_Pnt2d aPnt2d = BRep_Tool::Parameters(theVertex, myAttribute->Face());
1583 // check UV values for internal vertices
1584 if (myAttribute->ChangeClassifier()->Perform(aPnt2d) != TopAbs_IN)
1587 NCollection_Handle<FixedVExplorer> aFixedVExplorer = new FixedVExplorer(theVertex);
1588 Standard_Integer aIndex = myAttribute->GetVertexIndex(aFixedVExplorer);
1589 gp_XY anUV = BRepMesh_ShapeTool::FindUV(aIndex, aPnt2d,
1590 BRep_Tool::Tolerance(theVertex), myAttribute);
1592 Standard_Integer aTmpId1, aTmpId2;
1593 anUV = myAttribute->Scale(anUV, Standard_True);
1594 myAttribute->AddNode(aIndex, anUV, BRepMesh_Fixed, aTmpId1, aTmpId2);
1596 catch (Standard_Failure)
1601 //=======================================================================
1602 //function : insertVertex
1604 //=======================================================================
1605 void BRepMesh_FastDiscretFace::insertVertex(
1606 const gp_Pnt& thePnt3d,
1608 BRepMesh::ListOfVertex& theVertices)
1610 Standard_Integer aNbLocat = myAttribute->LastPointId();
1611 myAttribute->ChangeSurfacePoints()->Bind(++aNbLocat, thePnt3d);
1613 gp_XY aPnt2d = myAttribute->Scale(theUV, Standard_True);
1614 BRepMesh_Vertex aVertex(aPnt2d, aNbLocat, BRepMesh_Free);
1615 theVertices.Append(aVertex);
1618 //=======================================================================
1619 //function : commitSurfaceTriangulation
1621 //=======================================================================
1622 void BRepMesh_FastDiscretFace::commitSurfaceTriangulation()
1624 if (myAttribute.IsNull() || !myAttribute->IsValid())
1627 const TopoDS_Face& aFace = myAttribute->Face();
1628 BRepMesh_ShapeTool::NullifyFace(aFace);
1630 Handle(BRepMesh_DataStructureOfDelaun)& aStructure = myAttribute->ChangeStructure();
1631 const BRepMesh::MapOfInteger& aTriangles = aStructure->ElementsOfDomain();
1633 if (aTriangles.IsEmpty())
1635 myAttribute->SetStatus(BRepMesh_Failure);
1639 BRepMesh::HIMapOfInteger& aVetrexEdgeMap = myAttribute->ChangeVertexEdgeMap();
1642 Standard_Integer aVerticesNb = aVetrexEdgeMap->Extent();
1643 Standard_Integer aTrianglesNb = aTriangles.Extent();
1644 Handle(Poly_Triangulation) aNewTriangulation =
1645 new Poly_Triangulation(aVerticesNb, aTrianglesNb, Standard_True);
1647 Poly_Array1OfTriangle& aPolyTrianges = aNewTriangulation->ChangeTriangles();
1649 Standard_Integer aTriangeId = 1;
1650 BRepMesh::MapOfInteger::Iterator aTriIt(aTriangles);
1651 for (; aTriIt.More(); aTriIt.Next())
1653 const BRepMesh_Triangle& aCurElem = aStructure->GetElement(aTriIt.Key());
1655 Standard_Integer aNode[3];
1656 aStructure->ElementNodes(aCurElem, aNode);
1658 Standard_Integer aNodeId[3];
1659 for (Standard_Integer i = 0; i < 3; ++i)
1660 aNodeId[i] = aVetrexEdgeMap->FindIndex(aNode[i]);
1662 aPolyTrianges(aTriangeId++).Set(aNodeId[0], aNodeId[1], aNodeId[2]);
1666 TColgp_Array1OfPnt& aNodes = aNewTriangulation->ChangeNodes();
1667 TColgp_Array1OfPnt2d& aNodes2d = aNewTriangulation->ChangeUVNodes();
1669 for (Standard_Integer i = 1; i <= aVerticesNb; ++i)
1671 Standard_Integer aVertexId = aVetrexEdgeMap->FindKey(i);
1672 const BRepMesh_Vertex& aVertex = aStructure->GetNode(aVertexId);
1673 const gp_Pnt& aPoint = myAttribute->GetPoint(aVertex);
1676 aNodes2d(i) = aVertex.Coord();
1679 aNewTriangulation->Deflection(myAttribute->GetDefFace());
1680 BRepMesh_ShapeTool::AddInFace(aFace, aNewTriangulation);
1682 // Delete unused data
1685 myAttribute->ChangeStructure().Nullify();
1686 myAttribute->ChangeSurfacePoints().Nullify();
1687 myAttribute->ChangeSurfaceVertices().Nullify();
1689 myAttribute->ChangeClassifier().Nullify();
1690 myAttribute->ChangeLocation2D().Nullify();
1691 myAttribute->ChangeVertexEdgeMap().Nullify();