Salome HOME
Merge remote-tracking branch 'origin/BR_1064' into BR_HYDRO_IMPS_2016
[modules/hydro.git] / src / HYDROData / HYDROData_DTM.cxx
index e2bbbdeaed916028e1f55157a1d2422810ed6887..df831b619a110a4ba5dd9b61549e7c3be98774fd 100644 (file)
 #include <BRepLib_MakeEdge.hxx>
 #include <BRepLib_MakeWire.hxx>
 #include <BRep_Builder.hxx>
+#include <ShapeAnalysis_Wire.hxx>
+#include <BRepAlgo_NormalProjection.hxx>
+#include <ShapeUpgrade_UnifySameDomain.hxx>
+#include <BRepBuilderAPI_MakePolygon.hxx>
+#include <BOPAlgo_Builder.hxx>
 
-
+#include <TopTools_IndexedDataMapOfShapeShape.hxx>
+#include <TopTools_ListIteratorOfListOfShape.hxx>
+#include <TopTools_SequenceOfShape.hxx>
 
 IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
 
-
 HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ )
   : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ )
 {
@@ -80,7 +86,11 @@ double HYDROData_DTM::CurveUZ::DeltaZ() const
 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
 {
   HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ() );
-  size_t n = size();
+  size_t n = size(), n1 = c.size();
+  if( n!=n1 )
+  {
+    std::cout << "Warning: different number of points in curves: " << n << ", " << n1 << std::endl;
+  }
   res.reserve( n );
   for( int i=0; i<n; i++ )
   {
@@ -126,6 +136,7 @@ HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
 {
   SetReferenceObjects( theProfiles, DataTag_Profiles );
+  Changed( Geom_3d );
 }
 
 double HYDROData_DTM::GetDDZ() const
@@ -136,6 +147,7 @@ double HYDROData_DTM::GetDDZ() const
 void HYDROData_DTM::SetDDZ( double theDDZ )
 {
   SetDouble( DataTag_DDZ, theDDZ );
+  Changed( Geom_3d );
 }
 
 double HYDROData_DTM::GetSpatialStep() const
@@ -146,11 +158,12 @@ double HYDROData_DTM::GetSpatialStep() const
 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
 {
   SetDouble( DataTag_SpatialStep, theSpatialStep );
+  Changed( Geom_3d );
 }
 
 void HYDROData_DTM::PointToWire(const AltitudePoints& pnts, TopoDS_Wire& W )
 {
-  BRepLib_MakeWire WM;
+  /*BRepLib_MakeWire WM;
   if (pnts.empty())
     return;
   for (int i = 0; i < pnts.size() - 1; i++)
@@ -159,8 +172,13 @@ void HYDROData_DTM::PointToWire(const AltitudePoints& pnts, TopoDS_Wire& W )
     gp_Pnt p2(pnts[i+1].X, pnts[i+1].Y, pnts[i+1].Z);    
     WM.Add(BRepLib_MakeEdge(p1, p2).Edge()); 
   }
-  if (WM.IsDone())
-    W = WM.Wire();
+  if (WM.IsDone())*/
+  
+  BRepBuilderAPI_MakePolygon PM;
+  for (int i = 0; i < pnts.size(); i++)
+    PM.Add(gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z));
+  
+  W = PM.Wire();
 }
 
 TopTools_IndexedMapOfOrientedShape HYDROData_DTM::Create3DShape(const AltitudePoints& left,
@@ -207,7 +225,8 @@ void HYDROData_DTM::Update()
   double ddz = GetDDZ();
   double step = GetSpatialStep();
   std::set<int> InvInd;
-  CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, true, true, InvInd, -1 );
+  bool WireIntersections; //__TODO
+  CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, true, true, InvInd, -1, WireIntersections );
   SetAltitudePoints( points );  
 
   SetShape( DataTag_LeftBankShape, OutLeftB);
@@ -235,7 +254,6 @@ void HYDROData_DTM::GetPresentationShapes( TopoDS_Shape& Out3dPres,
   Out3dPres = GetShape( DataTag_3DShape );
   Out2dPres = GetShape( DataTag_2DShape );
 }
-
 void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& InpProfiles,
                                            double ddz,
                                            double step, 
@@ -249,7 +267,8 @@ void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& In
                                            bool Create3dPres,
                                            bool Create2dPres,
                                            std::set<int>& InvInd,
-                                           int thePntsLimit )
+                                           int thePntsLimit,
+                                           bool& WireIntersections)
 {
   int aLower = InpProfiles.Lower(), anUpper = InpProfiles.Upper();
   size_t n = anUpper - aLower + 1;
@@ -276,57 +295,172 @@ void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& In
 
   if( ddz>EPS && step>EPS )
     CreateProfiles(profiles, ddz, step, left, right, points, main_profiles, 
-    Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd );
+    Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections );
 }
 
-void HYDROData_DTM::ProjWireOnPlane(const TopoDS_Wire& inpWire, const Handle_Geom_Plane& RefPlane, TopoDS_Wire& outWire)
+void HYDROData_DTM::ProjWireOnPlane(const TopoDS_Shape& inpWire, const Handle_Geom_Plane& RefPlane, 
+  TopTools_DataMapOfShapeListOfShape* E2PE)
 {
-  //its also possible to use BrepAlgo_NormalProjection here!
-  BRepTools_WireExplorer ex(TopoDS::Wire(inpWire.Oriented(TopAbs_FORWARD)));
-  BRepLib_MakeWire WM;
-  for (;ex.More();ex.Next())
+  BRep_Builder BB;
+
+  //project shape (edges) on planar face
+  TopoDS_Face F;
+  BB.MakeFace(F, RefPlane, Precision::Confusion());
+  BRepAlgo_NormalProjection nproj(F);
+  nproj.Add(inpWire);
+  nproj.SetDefaultParams();
+  nproj.Build();
+  if(!nproj.IsDone())
+    return;
+
+  TopoDS_Shape projRes = nproj.Projection();
+
+  // unite all vertexes/edges from projected result
+  BOPAlgo_Builder anAlgo;
+  TopExp_Explorer exp(projRes, TopAbs_EDGE);
+  for (;exp.More(); exp.Next())
   {
-    const TopoDS_Edge& CE = ex.Current(); 
-    double f, l;
-    Handle(Geom_Curve) C3d = BRep_Tool::Curve(CE, f, l);
-    Handle(Geom_Curve) ProjectedCurve = GeomProjLib::ProjectOnPlane(new Geom_TrimmedCurve(C3d, f, l), RefPlane, RefPlane->Position().Direction(), Standard_True);
-    TopoDS_Edge ProjEdge = BRepLib_MakeEdge(ProjectedCurve);
-    WM.Add(ProjEdge); //auto sharing between edges if vertex is coincident
+    const TopoDS_Edge& E = TopoDS::Edge(exp.Current());
+    if (E.Orientation() != TopAbs_INTERNAL)
+      anAlgo.AddArgument(E);
+  }
+  anAlgo.Perform(); 
+  int stat = anAlgo.ErrorStatus();
+  TopoDS_Shape projResConn = anAlgo.Shape();
+
+  // make wire => vertexes and edges should be untouched after this operation!
+  exp.Init(projResConn, TopAbs_EDGE);
+  //TopTools_ListOfShape llE;
+  //TopoDS_Wire RW;
+  //BRepLib_MakeWire WM;
+
+  //for (;exp.More();exp.Next()) 
+  //  llE.Append(exp.Current());
+  //
+  //WM.Add(llE);
+  //outWire = WM.Wire();
+
+  //outWire.Orientation(inpWire.Orientation()); //take from the original wire
+
+  //history mode: edge to projected edges
+  if (E2PE)
+  {
+    TopExp_Explorer ex(inpWire, TopAbs_EDGE);
+    for (;ex.More();ex.Next())
+    {
+      const TopoDS_Edge& CE = TopoDS::Edge(ex.Current());
+      TopTools_ListOfShape NEL;
+      const TopTools_ListOfShape& LS = nproj.Generated(CE);
+      TopTools_ListIteratorOfListOfShape it(LS);
+      for (;it.More();it.Next())
+      {
+        const TopoDS_Shape& PCE = it.Value();
+        TopTools_ListOfShape PLS = anAlgo.Modified(PCE);
+        if (PLS.IsEmpty())
+          PLS.Append(PCE);
+        TopTools_ListIteratorOfListOfShape itp(PLS);
+        for (;itp.More();itp.Next())
+          NEL.Append(itp.Value());
+      }
+
+      E2PE->Bind(CE, NEL);
+    }
   }
-  outWire = WM.Wire();
-  outWire.Orientation(inpWire.Orientation()); //take from the original wire
 }
-
-
-void HYDROData_DTM::Get2dFaceFrom3dPres(const TopoDS_Compound& cmp, TopoDS_Face& outF )
+#include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
+bool HYDROData_DTM::Get2dFaceFrom3dPres(const TopoDS_Compound& cmp, TopoDS_Face& outF,
+  TopTools_SequenceOfShape* Boundr, std::set<int> ind )
 {
+  //ind : set of indices (starts with 0). index == number of boundary (inlet, outlet, etc..)
+  //in compound cmp. 
+  //if Boundr is not null => this method will return sequence of boundary wires (inlet, outlet...)
+
   Handle_Geom_Plane refpl = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1));
-  BRepLib_MakeWire WM;
+  TopTools_DataMapOfShapeListOfShape E2PE;
+  ProjWireOnPlane(cmp, refpl, &E2PE);
+  TopTools_ListOfShape ELL;
+
   TopoDS_Iterator it(cmp);
-  //TopTools_IndexedMapOfShape IntW;
+  int i = 0;
   for (;it.More(); it.Next())
   {
     const TopoDS_Wire& W = TopoDS::Wire(it.Value());
     if (W.Orientation() != TopAbs_INTERNAL)
     {
-      //use list of edges to protect againts non-manifold cases.
-      //auto sharing between edges will be added automatically
-      TopTools_IndexedMapOfShape ME;
-      TopTools_ListOfShape LE;
-      TopExp::MapShapes(W, TopAbs_EDGE, ME);
-      for (int i = 1; i <= ME.Extent(); i++)
-        LE.Append(ME(i));
-      WM.Add(LE);
+      TopoDS_Wire PW;
+      TopExp_Explorer ex(W, TopAbs_EDGE);   
+      TopTools_ListOfShape LEpW;
+      TopTools_ListOfShape LEpWDummy;
+      for (;ex.More();ex.Next())
+      {
+        const TopoDS_Edge& CE = TopoDS::Edge(ex.Current());
+        TopTools_ListOfShape LS = E2PE.Find(CE);
+        LEpW.Append(LS);
+      }
+      
+      if (ind.count(i) != 0)
+      {
+        BRepLib_MakeWire WM;
+        WM.Add(LEpW);
+        const TopoDS_Wire& WMW = WM.Wire();
+        //assume that wire is a straight line,
+        //take first and last vertex and make simple edge (RE)
+        TopoDS_Vertex VF, VL;
+        TopExp::Vertices(WMW, VF, VL);
+        TopoDS_Edge RE = BRepLib_MakeEdge(VF, VL).Edge();
+        if (RE.IsNull())
+        {
+          LEpWDummy = LEpW; //LEpW will be nullified after appending to ELL
+          ELL.Append(LEpW);
+        }
+        else
+        {
+          LEpWDummy.Append(RE); 
+          ELL.Append(RE);
+        }
+        //TODO: in the new version of OCCT, USD can process separate wires
+        //ShapeUpgrade_UnifySameDomain USD(WMW, 1U, 0U, 1U); //concat bsplines
+        //USD.Build();
+        //const TopoDS_Shape& RSU = USD.Shape();
+        //TopExp_Explorer exp(RSU, TopAbs_EDGE);
+        //TopTools_ListOfShape DummyL;
+        //for (;it.More();it.Next())
+        //  DummyL.Append(it.Value());
+        //if (DummyL.Extent() == 1)
+        //  ELL.Append(DummyL.First()); //if one edge => accept this result 
+        //else
+        //  ELL.Append(LEpW);  //else put 'as is'
+      }
+      else
+      {
+        LEpWDummy = LEpW;
+        ELL.Append(LEpW);
+      }
+
+      if (Boundr)
+      {
+        //make inlet, outlet, left/tight banks [wires]
+        //shouldn't change topology of the edges
+        BRepLib_MakeWire IWM;
+        IWM.Add(LEpWDummy);
+        Boundr->Append(IWM.Wire());
+      }
     }
-    //else
-    //  IntW.Add(W);
+    i++;
   }
 
-  TopoDS_Wire outW;
-  ProjWireOnPlane(WM.Wire(), refpl, outW);
-  BRepBuilderAPI_MakeFace mf(refpl, outW); //check inside is true by def
+  //make primary wire
+  BRepLib_MakeWire WME;
+  WME.Add(ELL);
+
+  const TopoDS_Wire& resW = WME.Wire();
+  BRepBuilderAPI_MakeFace mf(refpl, resW, true); //check inside is true by def
   outF = mf.Face();
 
+  ShapeAnalysis_Wire WA(resW, outF, Precision::Confusion());
+  bool res = WA.CheckSelfIntersection(); //TODO check that this return false if OK
+  return res;
+
   ///!!! the internal wires cant be added with 'internal' ori.
   // it's possible to do with brep builder yet the result will not be correct!
   // more proper way is to use BOP operation here.
@@ -354,7 +488,8 @@ void HYDROData_DTM::CreateProfiles(const std::vector<Handle_HYDROData_Profile>&
                                    TopoDS_Shape& OutOutlet,
                                    bool Create3dPres,
                                    bool Create2dPres,
-                                   std::set<int>& InvInd)
+                                   std::set<int>& InvInd,
+                                   bool& WireIntersections)
 {
   if (theProfiles.empty())
     return;
@@ -383,7 +518,7 @@ void HYDROData_DTM::CreateProfiles(const std::vector<Handle_HYDROData_Profile>&
     if (Create2dPres)
     {
       TopoDS_Face outF;
-      Get2dFaceFrom3dPres(cmp, outF);
+      WireIntersections = Get2dFaceFrom3dPres(cmp, outF); //__TODO
       Out2dPres = outF;
     };
   }
@@ -723,9 +858,11 @@ void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxi
   }
 
   thePoints.reserve( sorted_points.size() );
+  const double EPS = 1E-12;
   std::map<double, AltitudePoint>::const_iterator it = sorted_points.begin(), last = sorted_points.end();
   for( ; it!=last; it++ )
-    thePoints.push_back( it->second );
+    if( thePoints.empty() || thePoints.back().SquareDistance( it->second ) > EPS )
+      thePoints.push_back( it->second );
 }
 
 inline double max( double a, double b )