#include <TopTools_IndexedMapOfShape.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <TopExp.hxx>
+#include <TopTools_IndexedMapOfOrientedShape.hxx>
+#include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
#include <BRepLib_MakeEdge.hxx>
#include <BRepLib_MakeWire.hxx>
#include <ShapeUpgrade_UnifySameDomain.hxx>
#include <BRepBuilderAPI_MakePolygon.hxx>
#include <BOPAlgo_Builder.hxx>
-
+#include <BRepAdaptor_Curve.hxx>
+#include <GeomProjLib.hxx>
+#include <gp_Pln.hxx>
#include <TopTools_IndexedDataMapOfShapeShape.hxx>
#include <TopTools_ListIteratorOfListOfShape.hxx>
+#include <TopTools_SequenceOfShape.hxx>
+#include <assert.h>
+#include <float.h>
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 )
+HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ, double theMaxZ )
+ : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ ), myMaxZ (theMaxZ)
{
}
return myDeltaZ;
}
+double HYDROData_DTM::CurveUZ::MaxZ() const
+{
+ return myMaxZ;
+}
+
HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
{
- HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ() );
+ HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ(), MaxZ() + c.MaxZ() );
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++ )
+ int q = n < n1 ? n : n1;
+ res.reserve( q );
+ for( int i=0; i<q; i++ )
{
PointUZ p;
p.U = operator[]( i ).U + c[i].U;
HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
{
- HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d, DeltaZ()*d );
+ HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d, DeltaZ()*d, MaxZ()*d );
size_t n = size();
res.reserve( n );
for( int i=0; i<n; i++ )
Changed( Geom_3d );
}
-void HYDROData_DTM::PointToWire(const AltitudePoints& pnts, TopoDS_Wire& W )
+void HYDROData_DTM::PointsToWire(const AltitudePoints& pnts, TopoDS_Wire& W )
{
- /*BRepLib_MakeWire WM;
- if (pnts.empty())
- return;
- for (int i = 0; i < pnts.size() - 1; i++)
- {
- gp_Pnt p1(pnts[i].X, pnts[i].Y, pnts[i].Z);
- 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())*/
-
+
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();
}
+void HYDROData_DTM::PointsToEdge(const AltitudePoints& pnts, TopoDS_Edge& E )
+{
+ Handle(TColgp_HArray1OfPnt) gpPoints = new TColgp_HArray1OfPnt(1, (int)pnts.size());
+
+ for (int i = 0; i < pnts.size(); i++)
+ gpPoints->SetValue(i+1, gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z));
+
+ GeomAPI_Interpolate anInterpolator(gpPoints, Standard_False,1.0e-6);
+ anInterpolator.Perform() ;
+ if (anInterpolator.IsDone())
+ {
+ Handle(Geom_BSplineCurve) C = anInterpolator.Curve();
+ E = BRepBuilderAPI_MakeEdge(C).Edge();
+ }
+}
+
TopTools_IndexedMapOfOrientedShape HYDROData_DTM::Create3DShape(const AltitudePoints& left,
const AltitudePoints& right,
const std::vector<AltitudePoints>& main_profiles)
{
TopTools_IndexedMapOfOrientedShape ll;
- TopoDS_Wire LWire, RWire;
- PointToWire(left, LWire);
- PointToWire(right, RWire);
- if (!LWire.IsNull())
- ll.Add(LWire.Oriented(TopAbs_FORWARD));
+ //TopoDS_Wire LWire, RWire;
+ //PointsToWire(left, LWire);
+ //PointsToWire(right, RWire);
+ TopoDS_Edge LEdge, REdge;
+ PointsToEdge(left, LEdge);
+ PointsToEdge(right, REdge);
+ if (!LEdge.IsNull())
+ ll.Add(LEdge.Oriented(TopAbs_FORWARD));
for (int k = 0; k < main_profiles.size(); k++)
{
TopoDS_Wire W;
- PointToWire(main_profiles[k], W);
+ PointsToWire(main_profiles[k], W);
TopAbs_Orientation Ori = TopAbs_INTERNAL;
if (k == 0 || k == main_profiles.size() - 1)
Ori = TopAbs_FORWARD;
ll.Add(W.Oriented(Ori));
}
- if (!RWire.IsNull())
- ll.Add(RWire.Oriented(TopAbs_FORWARD));
+ if (!REdge.IsNull())
+ ll.Add(REdge.Oriented(TopAbs_FORWARD));
//yes, add subshapes in this order (left + profiles + right)
//otherwise the projected wire will be non-manifold
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);
SetShape( DataTag_RightBankShape, OutRightB);
SetShape( DataTag_InletShape, OutInlet);
Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections );
}
-void HYDROData_DTM::ProjWireOnPlane(const TopoDS_Shape& inpWire, const Handle(Geom_Plane)& RefPlane,
- TopTools_DataMapOfShapeListOfShape* E2PE)
+bool HYDROData_DTM::GetPlanarFaceFromBanks( const TopoDS_Edge& LB, const TopoDS_Edge& RB, TopoDS_Face& outF,
+ TopTools_SequenceOfShape* Boundr)
{
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;
+ Handle_Geom_Plane refpl = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1));
- 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& 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);
- }
- }
-}
-#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));
- TopTools_DataMapOfShapeListOfShape E2PE;
- ProjWireOnPlane(cmp, refpl, &E2PE);
- TopTools_ListOfShape ELL;
-
- TopoDS_Iterator it(cmp);
- int i = 0;
- for (;it.More(); it.Next())
- {
- const TopoDS_Wire& W = TopoDS::Wire(it.Value());
- if (W.Orientation() != TopAbs_INTERNAL)
- {
- 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);
- }
+ TopoDS_Vertex VFI, VLI, VFO, VLO;
+ TopoDS_Edge prLB;
+ TopoDS_Edge prRB;
- 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());
- }
- }
- i++;
- }
- //make primary wire
- BRepLib_MakeWire WME;
- WME.Add(ELL);
+ BRepAdaptor_Curve LBAD(LB);
+ Handle_Geom_Curve LBPC = GeomProjLib::ProjectOnPlane(LBAD.Curve().Curve(), refpl, gp_Dir(0, 0, -1), 1 );
+ prLB = BRepLib_MakeEdge(LBPC).Edge();
- const TopoDS_Wire& resW = WME.Wire();
- BRepBuilderAPI_MakeFace mf(refpl, resW, true); //check inside is true by def
- outF = mf.Face();
+ BRepAdaptor_Curve RBAD(RB);
+ Handle_Geom_Curve RBPC = GeomProjLib::ProjectOnPlane(RBAD.Curve().Curve(), refpl, gp_Dir(0, 0, -1), 1 );
+ prRB = BRepLib_MakeEdge(RBPC).Edge();
- ShapeAnalysis_Wire WA(resW, outF, Precision::Confusion());
- bool res = WA.CheckSelfIntersection(); //TODO check that this return false if OK
- return res;
+ TopExp::Vertices(prLB, VFI, VFO, 1);
+ TopExp::Vertices(prRB, VLI, VLO, 1);
+ TopoDS_Edge prIL = BRepLib_MakeEdge(VFI, VLI).Edge();
+ TopoDS_Edge prOL = BRepLib_MakeEdge(VFO, VLO).Edge();
+ TopoDS_Wire prW = BRepLib_MakeWire(prLB, prIL, prOL, prRB).Wire();
+ outF = BRepBuilderAPI_MakeFace(refpl->Pln(), prW, 1).Face();
- ///!!! 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.
- /*for (int i = 1; i <= IntW.Extent(); i++)
+ if (Boundr)
{
- TopoDS_Wire outIW;
- const TopoDS_Wire& W = TopoDS::Wire(IntW(i));
- ProjWireOnPlane(W, refpl, outIW);
- BB.Add(outF, outIW);
- }*/
+ Boundr->Append(prLB);
+ Boundr->Append(prIL);
+ Boundr->Append(prOL);
+ Boundr->Append(prRB);
+ }
+
+ ShapeAnalysis_Wire WA(prW, outF, Precision::Confusion());
+ bool res = WA.CheckSelfIntersection();
+ return !res;
}
void HYDROData_DTM::CreateProfiles(const std::vector<Handle(HYDROData_Profile)>& theProfiles,
bool Create3dPres,
bool Create2dPres,
std::set<int>& InvInd,
- bool& WireIntersections)
+ bool& ProjStat)
{
if (theProfiles.empty())
return;
if (Create2dPres)
{
TopoDS_Face outF;
- WireIntersections = Get2dFaceFrom3dPres(cmp, outF); //__TODO
+ ProjStat = GetPlanarFaceFromBanks(TopoDS::Edge(OutLeftB), TopoDS::Edge(OutRightB), outF, NULL);
Out2dPres = outF;
};
}
HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
int lo = points.Lower();
int up = points.Upper();
- theZMin = std::numeric_limits<double>::max();
+ theZMin = DBL_MAX;
theZMax = -theZMin;
for( int i=lo; i<=up; i++ )
{
TColgp_Array1OfVec2d tangents( 1, (int)n );
Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
- for( size_t i = 1; i <= n; i++ )
+ // Stage 1. Orient all profiles to be co-oriented with the first profile
+ theProfiles[0]->Update();
+ for( size_t i = 1; i < n; i++ )
{
- Handle(HYDROData_Profile) aProfile = theProfiles[i-1];
- Handle(HYDROData_Profile) aPrevProfile = i==1 ? theProfiles[i-1] : theProfiles[i-2];
- Handle(HYDROData_Profile) aNextProfile = i==n ? theProfiles[i-1] : theProfiles[i];
+ Handle(HYDROData_Profile) aProfile = theProfiles[i];
+ Handle(HYDROData_Profile) aPrevProfile = theProfiles[i-1];
- if( !IsCooriented( aProfile, aNextProfile ) )
+ if( !IsCooriented( aProfile, aPrevProfile ) )
{
gp_XY lp, rp;
aProfile->GetLeftPoint( lp, true );
aProfile->SetRightPoint( lp, true );
}
aProfile->Update();
+ }
+
+ // Stage 2. Calculate normals so that each normal "points" to the next profile
+ for( size_t i = 0; i < n; i++ )
+ {
+ Handle(HYDROData_Profile) aProfile = theProfiles[i];
+ Handle(HYDROData_Profile) aNextProfile = i==n-1 ? theProfiles[i-1] : theProfiles[i+1];
gp_Pnt aLowest;
gp_Vec2d aNormal;
double zmin, zmax;
- gp_XYZ curP = aProfile->GetBottomPoint();
+ gp_XYZ curP = aProfile->GetBottomPoint(true);
gp_XY curP2d = gp_XY(curP.X(), curP.Y());
- gp_XYZ nextP;
- if( i==n )
- nextP = aPrevProfile->GetBottomPoint(true);
- else
- nextP = aNextProfile->GetBottomPoint(true);
-
+ gp_XYZ nextP = aNextProfile->GetBottomPoint(true);
gp_XY nextP2d = gp_XY(nextP.X(), nextP.Y());
gp_Vec2d aPrTangent;
aNormal.SetCoord( -aPrTangent.Y(), aPrTangent.X() );
gp_Vec2d aDirToNextProfile(nextP2d.X() - curP2d.X(), nextP2d.Y() - curP2d.Y() );
- if( i==n )
+ if( i==n-1 )
aDirToNextProfile.Reverse();
+
if (aNormal.Dot(aDirToNextProfile) < 0)
aNormal.Reverse();
aNormal.Normalize();
- points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
- tangents.SetValue( (int)i, aNormal );
- flags->SetValue( (int)i, Standard_True );
+ points->SetValue( (int)(i+1), gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
+ tangents.SetValue( (int)(i+1), aNormal );
+ flags->SetValue( (int)(i+1), Standard_True );
}
Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
{
- double umin = std::numeric_limits<double>::max(),
+ double umin = DBL_MAX,
umax = -umin;
size_t n = intersections.size();
}
void HYDROData_DTM::ProfileDiscretization( const Handle(HYDROData_Profile)& theProfile,
- double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
+ double theXCurv, double theMinZ, double theMaxZ, double theTopZ, double theDDZ,
CurveUZ& theMidPointCurve,
CurveUZ& theWidthCurve,
int& intersection_nb,
double theTolerance)
{
- double aDblMax = std::numeric_limits<double>::max(),
+ double aDblMax = DBL_MAX,
aUMin = aDblMax,
aUMax = -aUMin,
aVMax = 1000000;
curves.push_back( aT2 );
int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
- theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
+ theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ);
theMidPointCurve.reserve( psize );
- theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
+ theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ );
theWidthCurve.reserve( psize );
n = curves.size();
size_t n = theMidCurve.size();
std::map<double, AltitudePoint> sorted_points;
- for( size_t i=0; i<n; i++ )
+ bool isOnTop = false;
+ for( size_t i=0; i<n; i++ ) // build the two banks of the interpolated profile, from bottom to top
{
double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
- double z = theMidCurve[i].Z + theMidCurve.DeltaZ();
-
- AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
+ bool arrivedOnTop = false;
+ double z = 0;
+ if (theMidCurve[i].Z <= theMidCurve.MaxZ())
+ z = theMidCurve[i].Z + theMidCurve.DeltaZ();
+ else
+ {
+ z = theMidCurve.MaxZ() + theMidCurve.DeltaZ(); // limit z to linear interpolation between maxima on extremity profiles
+ arrivedOnTop = true; // do not keep points after this one
+ }
+ if (!isOnTop)
+ {
+ AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
- sorted_points[param1] = p3d_1;
- sorted_points[param2] = p3d_2;
+ sorted_points[param1] = p3d_1;
+ sorted_points[param2] = p3d_2;
+ }
+ //if (arrivedOnTop)
+ // isOnTop =true; // do not keep points after this one (commented: leads to strange limits of 2D shape)
}
thePoints.reserve( sorted_points.size() );
thePoints.push_back( it->second );
}
-inline double max( double a, double b )
-{
- if( a>b )
- return a;
- else
- return b;
-}
-
-inline double min( double a, double b )
-{
- if( a<b )
- return a;
- else
- return b;
-}
+//inline double max( double a, double b )
+//{
+// if( a>b )
+// return a;
+// else
+// return b;
+//}
+//
+//inline double min( double a, double b )
+//{
+// if( a<b )
+// return a;
+// else
+// return b;
+//}
#include <BRepLib_MakeWire.hxx>
GetProperties( theProfileB, lowestB, dirB, zminB, zmaxB );
- double hmax = max( zmaxA-zminA, zmaxB-zminB );
+ double hmax = zmaxA-zminA > zmaxB-zminB ? zmaxA-zminA : zmaxB-zminB;
//double dz = zminB - zminA;
//double zmin = min( zminA, zminB );
//double zmax = max( zmaxA, zmaxB );
- CurveUZ midA(0, gp_Vec2d(), 0), midB(0, gp_Vec2d(), 0);
- CurveUZ widA(0, gp_Vec2d(), 0), widB(0, gp_Vec2d(), 0);
+ CurveUZ midA(0, gp_Vec2d(), 0, 0), midB(0, gp_Vec2d(), 0, 0);
+ CurveUZ widA(0, gp_Vec2d(), 0, 0), widB(0, gp_Vec2d(), 0, 0);
- ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, theDDZ, midA, widA, inter_nb_1 );
- ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, theDDZ, midB, widB, inter_nb_2 );
+ ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, zmaxA-zminA, theDDZ, midA, widA, inter_nb_1 );
+ ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, zmaxB-zminB, theDDZ, midB, widB, inter_nb_2 );
std::vector<CurveUZ> mid, wid;
Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
