Loading files
3d097a032This file is pinned to the selected release. Download the original to open it in your design tool.
// Created on: 1997-05-05
// Created by: Jerome LEMONIER
// Copyright (c) 1996-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.
#include <GeomPlate_BuildPlateSurface.hxx>
#include <Adaptor2d_Curve2d.hxx>
#include <Adaptor3d_Curve.hxx>
#include <Adaptor3d_CurveOnSurface.hxx>
#include <Approx_CurveOnSurface.hxx>
#include <Extrema_POnSurf.hxx>
#include <GCPnts_AbscissaPoint.hxx>
#include <Geom2d_BezierCurve.hxx>
#include <Geom2d_BSplineCurve.hxx>
#include <Geom2d_Curve.hxx>
#include <Geom2dAdaptor_Curve.hxx>
#include <Geom2dInt_GInter.hxx>
#include <Geom_BSplineSurface.hxx>
#include <Geom_Plane.hxx>
#include <Geom_RectangularTrimmedSurface.hxx>
#include <Geom_Surface.hxx>
#include <GeomAdaptor.hxx>
#include <GeomAdaptor_Surface.hxx>
#include <GeomLProp_SLProps.hxx>
#include <GeomPlate_BuildAveragePlane.hxx>
#include <GeomPlate_CurveConstraint.hxx>
#include <NCollection_Sequence.hxx>
#include <NCollection_Array1.hxx>
#include <NCollection_HArray1.hxx>
#include <GeomPlate_MakeApprox.hxx>
#include <GeomPlate_PointConstraint.hxx>
#include <GeomPlate_Aij.hxx>
#include <GeomPlate_Surface.hxx>
#include <gp_Pnt.hxx>
#include <gp_Pnt2d.hxx>
#include <gp_Vec.hxx>
#include <gp_Vec2d.hxx>
#include <IntRes2d_IntersectionPoint.hxx>
#include <Law_Interpol.hxx>
#include <LocalAnalysis_SurfaceContinuity.hxx>
#include <Plate_D2.hxx>
#include <Plate_FreeGtoCConstraint.hxx>
#include <Plate_GtoCConstraint.hxx>
#include <Plate_PinpointConstraint.hxx>
#include <Plate_Plate.hxx>
#include <Precision.hxx>
#include <ProjLib_HCompProjectedCurve.hxx>
#include <Standard_ConstructionError.hxx>
#include <Standard_Integer.hxx>
#include <Message_ProgressScope.hxx>
#include <cstdio>
#ifdef OCCT_DEBUG
#include <OSD_Chronometer.hxx>
static int Affich = 0;
#endif
//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
// =========================================================
// C O N S T R U C T O R S
// =========================================================
//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
//---------------------------------------------------------
// Constructor compatible with the old version
//---------------------------------------------------------
GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface(
const occ::handle<NCollection_HArray1<int>>& NPoints,
const occ::handle<NCollection_HArray1<occ::handle<Adaptor3d_Curve>>>& TabCurve,
const occ::handle<NCollection_HArray1<int>>& Tang,
const int Degree,
const int NbIter,
const double Tol2d,
const double Tol3d,
const double TolAng,
const double,
const bool Anisotropie)
: myAnisotropie(Anisotropie),
myDegree(Degree),
myNbIter(NbIter),
myTol2d(Tol2d),
myTol3d(Tol3d),
myTolAng(TolAng),
myNbBounds(0)
{
int NTCurve = TabCurve->Length(); // Number of linear constraints
myNbPtsOnCur = 0; // Different calculation of the number of points depending on the length
myLinCont = new NCollection_HSequence<occ::handle<GeomPlate_CurveConstraint>>;
myPntCont = new NCollection_HSequence<occ::handle<GeomPlate_PointConstraint>>;
if (myNbIter < 1)
{
throw Standard_ConstructionError("GeomPlate : Number of iteration must be >= 1");
}
if (NTCurve == 0)
{
throw Standard_ConstructionError("GeomPlate : the bounds Array is null");
}
if (Tang->Length() == 0)
{
throw Standard_ConstructionError("GeomPlate : the constraints Array is null");
}
int nbp = 0;
int i;
for (i = 1; i <= NTCurve; i++)
{
nbp += NPoints->Value(i);
}
if (nbp == 0)
{
throw Standard_ConstructionError(
"GeomPlate : the resolution is impossible if the number of constraints points is 0");
}
if (myDegree < 2)
{
throw Standard_ConstructionError("GeomPlate ; the degree resolution must be upper of 2");
}
// Filling fields passing from the old constructor to the new one
for (i = 1; i <= NTCurve; i++)
{
occ::handle<GeomPlate_CurveConstraint> Cont =
new GeomPlate_CurveConstraint(TabCurve->Value(i), Tang->Value(i), NPoints->Value(i));
myLinCont->Append(Cont);
}
mySurfInitIsGive = false;
myIsLinear = true;
myFree = false;
}
//------------------------------------------------------------------
// Constructor with initial surface and degree
//------------------------------------------------------------------
GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface(const occ::handle<Geom_Surface>& Surf,
const int Degree,
const int NbPtsOnCur,
const int NbIter,
const double Tol2d,
const double Tol3d,
const double TolAng,
const double /*TolCurv*/,
const bool Anisotropie)
: mySurfInit(Surf),
myAnisotropie(Anisotropie),
myDegree(Degree),
myNbPtsOnCur(NbPtsOnCur),
myNbIter(NbIter),
myTol2d(Tol2d),
myTol3d(Tol3d),
myTolAng(TolAng),
myNbBounds(0)
{
if (myNbIter < 1)
{
throw Standard_ConstructionError("GeomPlate : Number of iteration must be >= 1");
}
if (myDegree < 2)
{
throw Standard_ConstructionError("GeomPlate : the degree must be above 2");
}
myLinCont = new NCollection_HSequence<occ::handle<GeomPlate_CurveConstraint>>;
myPntCont = new NCollection_HSequence<occ::handle<GeomPlate_PointConstraint>>;
mySurfInitIsGive = true;
myIsLinear = true;
myFree = false;
}
//---------------------------------------------------------
// Constructor with degree
//---------------------------------------------------------
GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface(const int Degree,
const int NbPtsOnCur,
const int NbIter,
const double Tol2d,
const double Tol3d,
const double TolAng,
const double /*TolCurv*/,
const bool Anisotropie)
: myAnisotropie(Anisotropie),
myDegree(Degree),
myNbPtsOnCur(NbPtsOnCur),
myNbIter(NbIter),
myTol2d(Tol2d),
myTol3d(Tol3d),
myTolAng(TolAng),
myNbBounds(0)
{
if (myNbIter < 1)
{
throw Standard_ConstructionError("GeomPlate : Number of iteration must be >= 1");
}
if (myDegree < 2)
{
throw Standard_ConstructionError("GeomPlate : the degree resolution must be upper of 2");
}
myLinCont = new NCollection_HSequence<occ::handle<GeomPlate_CurveConstraint>>;
myPntCont = new NCollection_HSequence<occ::handle<GeomPlate_PointConstraint>>;
mySurfInitIsGive = false;
myIsLinear = true;
myFree = false;
}
//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
// =========================================================
// P U B L I C M E T H O D S
// =========================================================
//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
//=================================================================================================
// Reorder the table of transformations
// After the call of CourbeJointive the order of curves is modified
// Ex : initial order of curves ==> A B C D E F
// In TabInit we note ==> 1 2 3 4 5 6
// after CourbeJointive ==> A E C B D F
// TabInit ==> 1 5 3 2 4 6
// after TrierTab the Table contains ==> 1 4 3 5 2 6
// It is also possible to access the 2nd curve by taking TabInit[2]
// i.e. the 4th from the table of classified curves
//-------------------------------------------------------------------------
static void TrierTab(occ::handle<NCollection_HArray1<int>>& Tab)
{
// Parse the table of transformations to find the initial order
int Nb = Tab->Length();
NCollection_Array1<int> TabTri(1, Nb);
for (int i = 1; i <= Nb; i++)
{
TabTri.SetValue(Tab->Value(i), i);
}
Tab->ChangeArray1() = TabTri;
}
//---------------------------------------------------------
// Function : ProjectCurve
//---------------------------------------------------------
occ::handle<Geom2d_Curve> GeomPlate_BuildPlateSurface::ProjectCurve(
const occ::handle<Adaptor3d_Curve>& Curv)
{
// Project a curve on a plane
occ::handle<Geom2d_Curve> Curve2d;
occ::handle<GeomAdaptor_Surface> hsur = new GeomAdaptor_Surface(mySurfInit);
gp_Pnt2d P2d;
occ::handle<ProjLib_HCompProjectedCurve> HProjector =
new ProjLib_HCompProjectedCurve(hsur, Curv, myTol3d / 10, myTol3d / 10);
double UdebCheck, UfinCheck, ProjUdeb, ProjUfin;
UdebCheck = Curv->FirstParameter();
UfinCheck = Curv->LastParameter();
HProjector->Bounds(1, ProjUdeb, ProjUfin);
if (HProjector->NbCurves() != 1 || std::abs(UdebCheck - ProjUdeb) > Precision::PConfusion()
|| std::abs(UfinCheck - ProjUfin) > Precision::PConfusion())
{
if (HProjector->IsSinglePnt(1, P2d))
{
// solution in a point
NCollection_Array1<gp_Pnt2d> poles(1, 2);
poles.Init(P2d);
Curve2d = new (Geom2d_BezierCurve)(poles);
}
else
{
Curve2d.Nullify(); // No continuous solution
#ifdef OCCT_DEBUG
std::cout << "BuildPlateSurace :: No continuous projection" << std::endl;
#endif
}
}
else
{
GeomAbs_Shape Continuity = GeomAbs_C1;
int MaxDegree = 10, MaxSeg;
double Udeb, Ufin;
HProjector->Bounds(1, Udeb, Ufin);
MaxSeg = 20 + HProjector->NbIntervals(GeomAbs_C3);
Approx_CurveOnSurface appr(HProjector, hsur, Udeb, Ufin, myTol3d);
appr.Perform(MaxSeg, MaxDegree, Continuity, false, true);
Curve2d = appr.Curve2d();
}
return Curve2d;
}
//---------------------------------------------------------
// Function : ProjectedCurve
//---------------------------------------------------------
occ::handle<Adaptor2d_Curve2d> GeomPlate_BuildPlateSurface::ProjectedCurve(
occ::handle<Adaptor3d_Curve>& Curv)
{
// Projection of a curve on the initial surface
occ::handle<GeomAdaptor_Surface> hsur = new GeomAdaptor_Surface(mySurfInit);
occ::handle<ProjLib_HCompProjectedCurve> HProjector =
new ProjLib_HCompProjectedCurve(hsur, Curv, myTolU / 10, myTolV / 10);
if (HProjector->NbCurves() != 1)
{
HProjector.Nullify(); // No continuous solution
#ifdef OCCT_DEBUG
std::cout << "BuildPlateSurace :: No continuous projection" << std::endl;
#endif
}
else
{
double First1, Last1, First2, Last2;
First1 = Curv->FirstParameter();
Last1 = Curv->LastParameter();
HProjector->Bounds(1, First2, Last2);
if (std::abs(First1 - First2) <= std::max(myTolU, myTolV)
&& std::abs(Last1 - Last2) <= std::max(myTolU, myTolV))
{
HProjector = occ::down_cast<ProjLib_HCompProjectedCurve>(
HProjector->Trim(First2, Last2, Precision::PConfusion()));
}
else
{
HProjector.Nullify(); // No continuous solution
#ifdef OCCT_DEBUG
std::cout << "BuildPlateSurace :: No complete projection" << std::endl;
#endif
}
}
return HProjector;
}
//---------------------------------------------------------
// Function : ProjectPoint
//---------------------------------------------------------
// Projects a point on the initial surface
//---------------------------------------------------------
gp_Pnt2d GeomPlate_BuildPlateSurface::ProjectPoint(const gp_Pnt& p3d)
{
Extrema_POnSurf P;
myProj.Perform(p3d);
int nearest = 1;
if (myProj.NbExt() > 1)
{
double dist2mini = myProj.SquareDistance(1);
for (int i = 2; i <= myProj.NbExt(); i++)
{
if (myProj.SquareDistance(i) < dist2mini)
{
dist2mini = myProj.SquareDistance(i);
nearest = i;
}
}
}
P = myProj.Point(nearest);
double u, v;
P.Parameter(u, v);
gp_Pnt2d p2d;
p2d.SetCoord(u, v);
return p2d;
}
//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
// =========================================================
// P U B L I C M E T H O D S
// =========================================================
//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
//---------------------------------------------------------
// Function : Init
//---------------------------------------------------------
// Initializes linear and point constraints
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::Init()
{
myLinCont->Clear();
myPntCont->Clear();
myPntCont = new NCollection_HSequence<occ::handle<GeomPlate_PointConstraint>>;
myLinCont = new NCollection_HSequence<occ::handle<GeomPlate_CurveConstraint>>;
}
//---------------------------------------------------------
// Function : LoadInitSurface
//---------------------------------------------------------
// Loads the initial surface
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::LoadInitSurface(const occ::handle<Geom_Surface>& Surf)
{
mySurfInit = Surf;
mySurfInitIsGive = true;
}
//---------------------------------------------------------
// Function : Add
//---------------------------------------------------------
//---------------------------------------------------------
// Adds a linear constraint
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::Add(const occ::handle<GeomPlate_CurveConstraint>& Cont)
{
myLinCont->Append(Cont);
}
void GeomPlate_BuildPlateSurface::SetNbBounds(const int NbBounds)
{
myNbBounds = NbBounds;
}
//---------------------------------------------------------
// Function : Add
//---------------------------------------------------------
//---------------------------------------------------------
// Adds a point constraint
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::Add(const occ::handle<GeomPlate_PointConstraint>& Cont)
{
myPntCont->Append(Cont);
}
//---------------------------------------------------------
// Function : Perform
// Calculates the surface filled with loaded constraints
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::Perform(const Message_ProgressRange& theProgress)
{
myGeomPlateSurface.Nullify();
#ifdef OCCT_DEBUG
// Timing
OSD_Chronometer Chrono;
Chrono.Reset();
Chrono.Start();
#endif
if (myNbBounds == 0)
{
myNbBounds = myLinCont->Length();
}
myPlate.Init();
//=====================================================================
// Declaration of variables.
//=====================================================================
int NTLinCont = myLinCont->Length(), NTPntCont = myPntCont->Length(), NbBoucle = 0;
bool Fini = true;
if ((NTLinCont + NTPntCont) == 0)
{
#ifdef OCCT_DEBUG
std::cout << "WARNING : GeomPlate : The number of constraints is null." << std::endl;
#endif
return;
}
//======================================================================
// Initial Surface
//======================================================================
Message_ProgressScope aPS(theProgress, "Calculating the surface filled", 100, true);
if (!mySurfInitIsGive)
{
ComputeSurfInit(aPS.Next(10));
if (aPS.UserBreak())
{
return;
}
}
else
{
if (NTLinCont >= 2)
{ // Table of transformations to preserve the initial order, see TrierTab
myInitOrder = new NCollection_HArray1<int>(1, NTLinCont);
for (int l = 1; l <= NTLinCont; l++)
{
myInitOrder->SetValue(l, l);
}
if (!CourbeJointive(myTol3d))
{ // throw Standard_Failure("Curves are not joined");
#ifdef OCCT_DEBUG
std::cout << "WARNING : Courbes non jointives a " << myTol3d << " pres" << std::endl;
#endif
}
TrierTab(myInitOrder); // Reorder the table of transformations
}
else if (NTLinCont > 0) // Patch
{
mySense = new NCollection_HArray1<int>(1, NTLinCont, 0);
myInitOrder = new NCollection_HArray1<int>(1, NTLinCont, 1);
}
}
if (mySurfInit.IsNull())
{
return;
}
double u1, v1, u2, v2;
mySurfInit->Bounds(u1, v1, u2, v2);
GeomAdaptor_Surface aSurfInit(mySurfInit);
myTolU = aSurfInit.UResolution(myTol3d);
myTolV = aSurfInit.VResolution(myTol3d);
myProj.Initialize(aSurfInit, u1, v1, u2, v2, myTolU, myTolV);
//======================================================================
// Projection of curves
//======================================================================
bool Ok = true;
for (int i = 1; i <= NTLinCont; i++)
{
if (myLinCont->Value(i)->Curve2dOnSurf().IsNull())
{
occ::handle<Geom2d_Curve> Curve2d = ProjectCurve(myLinCont->Value(i)->Curve3d());
if (Curve2d.IsNull())
{
Ok = false;
break;
}
myLinCont->ChangeValue(i)->SetCurve2dOnSurf(Curve2d);
}
}
if (!Ok)
{
GeomPlate_MakeApprox App(myGeomPlateSurface, myTol3d, 1, 3, 15 * myTol3d, -1, GeomAbs_C0, 1.3);
mySurfInit = App.Surface();
mySurfInit->Bounds(u1, v1, u2, v2);
GeomAdaptor_Surface Surf(mySurfInit);
myTolU = Surf.UResolution(myTol3d);
myTolV = Surf.VResolution(myTol3d);
myProj.Initialize(Surf, u1, v1, u2, v2, myTolU, myTolV);
Ok = true;
for (int i = 1; i <= NTLinCont; i++)
{
occ::handle<Geom2d_Curve> Curve2d = ProjectCurve(myLinCont->Value(i)->Curve3d());
if (Curve2d.IsNull())
{
Ok = false;
break;
}
myLinCont->ChangeValue(i)->SetCurve2dOnSurf(Curve2d);
}
if (!Ok)
{
mySurfInit = myPlanarSurfInit;
mySurfInit->Bounds(u1, v1, u2, v2);
GeomAdaptor_Surface SurfNew(mySurfInit);
myTolU = SurfNew.UResolution(myTol3d);
myTolV = SurfNew.VResolution(myTol3d);
myProj.Initialize(SurfNew, u1, v1, u2, v2, myTolU, myTolV);
for (int i = 1; i <= NTLinCont; i++)
{
myLinCont->ChangeValue(i)->SetCurve2dOnSurf(ProjectCurve(myLinCont->Value(i)->Curve3d()));
}
}
else
{ // Project the points
for (int i = 1; i <= NTPntCont; i++)
{
gp_Pnt P;
myPntCont->Value(i)->D0(P);
myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P));
}
}
}
//======================================================================
// Projection of points
//======================================================================
for (int i = 1; i <= NTPntCont; i++)
{
if (!myPntCont->Value(i)->HasPnt2dOnSurf())
{
gp_Pnt P;
myPntCont->Value(i)->D0(P);
myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P));
}
}
//======================================================================
// Number of points by curve
//======================================================================
if ((NTLinCont != 0) && (myNbPtsOnCur != 0))
{
CalculNbPtsInit();
}
//======================================================================
// Management of incompatibilites between curves
//======================================================================
occ::handle<NCollection_HArray1<NCollection_Sequence<double>>> PntInter;
occ::handle<NCollection_HArray1<NCollection_Sequence<double>>> PntG1G1;
if (NTLinCont != 0)
{
PntInter = new NCollection_HArray1<NCollection_Sequence<double>>(1, NTLinCont);
PntG1G1 = new NCollection_HArray1<NCollection_Sequence<double>>(1, NTLinCont);
Intersect(PntInter, PntG1G1);
}
//======================================================================
// Loop to obtain a better surface
//======================================================================
myFree = !myIsLinear;
do
{
#ifdef OCCT_DEBUG
if (Affich && NbBoucle)
{
std::cout << "Resultats boucle" << NbBoucle << std::endl;
std::cout << "DistMax=" << myG0Error << std::endl;
if (myG1Error != 0)
std::cout << "AngleMax=" << myG1Error << std::endl;
if (myG2Error != 0)
std::cout << "CourbMax=" << myG2Error << std::endl;
}
#endif
NbBoucle++;
if (NTLinCont != 0)
{ //====================================================================
// Calculate the total number of points and the maximum of points by curve
//====================================================================
int NPointMax = 0;
for (int i = 1; i <= NTLinCont; i++)
{
if ((myLinCont->Value(i)->NbPoints()) > NPointMax)
{
NPointMax = (int)(myLinCont->Value(i)->NbPoints());
}
}
//====================================================================
// Discretization of curves
//====================================================================
Discretise(PntInter, PntG1G1);
//====================================================================
// Preparation of constraint points for plate
//====================================================================
LoadCurve(NbBoucle);
if (myPntCont->Length() != 0)
{
LoadPoint(NbBoucle);
}
//====================================================================
// Construction of the surface
//====================================================================
myPlate.SolveTI(myDegree, ComputeAnisotropie(), aPS.Next(90));
if (aPS.UserBreak())
{
return;
}
if (!myPlate.IsDone())
{
#ifdef OCCT_DEBUG
std::cout << "WARNING : GeomPlate : calculation of Plate failed" << std::endl;
#endif
return;
}
myGeomPlateSurface = new GeomPlate_Surface(mySurfInit, myPlate);
double Umin, Umax, Vmin, Vmax;
myPlate.UVBox(Umin, Umax, Vmin, Vmax);
myGeomPlateSurface->SetBounds(Umin, Umax, Vmin, Vmax);
Fini = VerifSurface(NbBoucle);
if ((NbBoucle >= myNbIter) && (!Fini))
{
#ifdef OCCT_DEBUG
std::cout << "Warning: objective was not reached" << std::endl;
#endif
Fini = true;
}
if ((NTPntCont != 0) && (Fini))
{
double di, an, cu;
VerifPoints(di, an, cu);
}
}
else
{
LoadPoint(NbBoucle);
//====================================================================
// Construction of the surface
//====================================================================
myPlate.SolveTI(myDegree, ComputeAnisotropie(), aPS.Next(90));
if (aPS.UserBreak())
{
return;
}
if (!myPlate.IsDone())
{
#ifdef OCCT_DEBUG
std::cout << "WARNING : GeomPlate : calculation of Plate failed" << std::endl;
#endif
return;
}
myGeomPlateSurface = new GeomPlate_Surface(mySurfInit, myPlate);
double Umin, Umax, Vmin, Vmax;
myPlate.UVBox(Umin, Umax, Vmin, Vmax);
myGeomPlateSurface->SetBounds(Umin, Umax, Vmin, Vmax);
Fini = true;
double di, an, cu;
VerifPoints(di, an, cu);
}
} while (!Fini); // End loop for better surface
#ifdef OCCT_DEBUG
if (NTLinCont != 0)
{
std::cout << "======== Global results ===========" << std::endl;
std::cout << "DistMax=" << myG0Error << std::endl;
if (myG1Error != 0)
std::cout << "AngleMax=" << myG1Error << std::endl;
if (myG2Error != 0)
std::cout << "CourbMax=" << myG2Error << std::endl;
}
Chrono.Stop();
double Tps;
Chrono.Show(Tps);
std::cout << "*** END OF GEOMPLATE ***" << std::endl;
std::cout << "Time of calculation : " << Tps << std::endl;
std::cout << "Number of loops : " << NbBoucle << std::endl;
#endif
}
//---------------------------------------------------------
// Function : EcartContraintesMIL
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::EcartContraintesMil(
const int c,
occ::handle<NCollection_HArray1<double>>& d,
occ::handle<NCollection_HArray1<double>>& an,
occ::handle<NCollection_HArray1<double>>& courb)
{
int NbPt = myParCont->Value(c).Length();
double U;
if (NbPt < 3)
{
NbPt = 4;
}
else
{
NbPt = myParCont->Value(c).Length();
}
gp_Vec v1i, v1f, v2i, v2f, v3i, v3f;
gp_Pnt Pi, Pf;
gp_Pnt2d P2d;
int i;
occ::handle<GeomPlate_CurveConstraint> LinCont = myLinCont->Value(c);
switch (LinCont->Order())
{
case 0:
for (i = 1; i < NbPt; i++)
{
U = (myParCont->Value(c).Value(i) + myParCont->Value(c).Value(i + 1)) / 2;
LinCont->D0(U, Pi);
if (!LinCont->ProjectedCurve().IsNull())
{
P2d = LinCont->ProjectedCurve()->Value(U);
}
else
{
if (!LinCont->Curve2dOnSurf().IsNull())
{
P2d = LinCont->Curve2dOnSurf()->Value(U);
}
else
{
P2d = ProjectPoint(Pi);
}
}
myGeomPlateSurface->D0(P2d.Coord(1), P2d.Coord(2), Pf);
an->Init(0);
courb->Init(0);
d->ChangeValue(i) = Pf.Distance(Pi);
}
break;
case 1:
for (i = 1; i < NbPt; i++)
{
U = (myParCont->Value(c).Value(i) + myParCont->Value(c).Value(i + 1)) / 2;
LinCont->D1(U, Pi, v1i, v2i);
if (!LinCont->ProjectedCurve().IsNull())
{
P2d = LinCont->ProjectedCurve()->Value(U);
}
else
{
if (!LinCont->Curve2dOnSurf().IsNull())
{
P2d = LinCont->Curve2dOnSurf()->Value(U);
}
else
{
P2d = ProjectPoint(Pi);
}
}
myGeomPlateSurface->D1(P2d.Coord(1), P2d.Coord(2), Pf, v1f, v2f);
d->ChangeValue(i) = Pf.Distance(Pi);
v3i = v1i ^ v2i;
v3f = v1f ^ v2f;
double angle = v3f.Angle(v3i);
if (angle > (M_PI / 2))
{
an->ChangeValue(i) = M_PI - angle;
}
else
{
an->ChangeValue(i) = angle;
}
courb->Init(0);
}
break;
case 2:
occ::handle<Geom_Surface> Splate(myGeomPlateSurface);
LocalAnalysis_SurfaceContinuity CG2;
for (i = 1; i < NbPt; i++)
{
U = (myParCont->Value(c).Value(i) + myParCont->Value(c).Value(i + 1)) / 2;
LinCont->D0(U, Pi);
if (!LinCont->ProjectedCurve().IsNull())
{
P2d = LinCont->ProjectedCurve()->Value(U);
}
else
{
if (!LinCont->Curve2dOnSurf().IsNull())
{
P2d = LinCont->Curve2dOnSurf()->Value(U);
}
else
{
P2d = ProjectPoint(Pi);
}
}
GeomLProp_SLProps Prop(Splate, P2d.Coord(1), P2d.Coord(2), 2, 0.001);
CG2.ComputeAnalysis(Prop, myLinCont->Value(c)->LPropSurf(U), GeomAbs_G2);
d->ChangeValue(i) = CG2.C0Value();
an->ChangeValue(i) = CG2.G1Angle();
courb->ChangeValue(i) = CG2.G2CurvatureGap();
}
break;
}
}
//---------------------------------------------------------
// Function : Disc2dContour
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::Disc2dContour(const int /*nbp*/,
NCollection_Sequence<gp_XY>& Seq2d)
{
#ifdef OCCT_DEBUG
if (Seq2d.Length() != 4)
std::cout << "Number of points should be equal to 4 for Disc2dContour" << std::endl;
#endif
// initialization
Seq2d.Clear();
// sampling in "cosine" + 3 points on each interval
int NTCurve = myLinCont->Length();
int NTPntCont = myPntCont->Length();
gp_Pnt2d P2d;
gp_XY UV;
gp_Pnt PP;
double u1, v1, u2, v2;
int i;
mySurfInit->Bounds(u1, v1, u2, v2);
GeomAdaptor_Surface Surf(mySurfInit);
myProj.Initialize(Surf, u1, v1, u2, v2, myTolU, myTolV);
for (i = 1; i <= NTPntCont; i++)
{
if (myPntCont->Value(i)->Order() != -1)
{
P2d = myPntCont->Value(i)->Pnt2dOnSurf();
UV.SetX(P2d.Coord(1));
UV.SetY(P2d.Coord(2));
Seq2d.Append(UV);
}
}
for (i = 1; i <= NTCurve; i++)
{
occ::handle<GeomPlate_CurveConstraint> LinCont = myLinCont->Value(i);
if (LinCont->Order() != -1)
{
int NbPt = myParCont->Value(i).Length();
// first point of constraint (j=0)
if (!LinCont->ProjectedCurve().IsNull())
{
P2d = LinCont->ProjectedCurve()->Value(myParCont->Value(i).Value(1));
}
else
{
if (!LinCont->Curve2dOnSurf().IsNull())
{
P2d = LinCont->Curve2dOnSurf()->Value(myParCont->Value(i).Value(1));
}
else
{
LinCont->D0(myParCont->Value(i).Value(1), PP);
P2d = ProjectPoint(PP);
}
}
UV.SetX(P2d.Coord(1));
UV.SetY(P2d.Coord(2));
Seq2d.Append(UV);
for (int j = 2; j < NbPt; j++)
{
double Uj = myParCont->Value(i).Value(j), Ujp1 = myParCont->Value(i).Value(j + 1);
if (!LinCont->ProjectedCurve().IsNull())
{
P2d = LinCont->ProjectedCurve()->Value((Ujp1 + 3 * Uj) / 4);
}
else
{
if (!LinCont->Curve2dOnSurf().IsNull())
{
P2d = LinCont->Curve2dOnSurf()->Value((Ujp1 + 3 * Uj) / 4);
}
else
{
LinCont->D0((Ujp1 + 3 * Uj) / 4, PP);
P2d = ProjectPoint(PP);
}
}
UV.SetX(P2d.Coord(1));
UV.SetY(P2d.Coord(2));
Seq2d.Append(UV);
// point 1/2 previous
if (!LinCont->ProjectedCurve().IsNull())
{
P2d = LinCont->ProjectedCurve()->Value((Ujp1 + Uj) / 2);
}
else
{
if (!LinCont->Curve2dOnSurf().IsNull())
{
P2d = LinCont->Curve2dOnSurf()->Value((Ujp1 + Uj) / 2);
}
else
{
LinCont->D0((Ujp1 + Uj) / 2, PP);
P2d = ProjectPoint(PP);
}
}
UV.SetX(P2d.Coord(1));
UV.SetY(P2d.Coord(2));
Seq2d.Append(UV);
// point 3/4 previous
if (!LinCont->ProjectedCurve().IsNull())
{
P2d = LinCont->ProjectedCurve()->Value((3 * Ujp1 + Uj) / 4);
}
else
{
if (!LinCont->Curve2dOnSurf().IsNull())
{
P2d = LinCont->Curve2dOnSurf()->Value((3 * Ujp1 + Uj) / 4);
}
else
{
LinCont->D0((3 * Ujp1 + Uj) / 4, PP);
P2d = ProjectPoint(PP);
}
}
UV.SetX(P2d.Coord(1));
UV.SetY(P2d.Coord(2));
Seq2d.Append(UV);
// current constraint point
if (!LinCont->ProjectedCurve().IsNull())
{
P2d = LinCont->ProjectedCurve()->Value(Ujp1);
}
else
{
if (!LinCont->Curve2dOnSurf().IsNull())
{
P2d = LinCont->Curve2dOnSurf()->Value(Ujp1);
}
else
{
LinCont->D0(Ujp1, PP);
P2d = ProjectPoint(PP);
}
}
UV.SetX(P2d.Coord(1));
UV.SetY(P2d.Coord(2));
Seq2d.Append(UV);
}
}
}
}
//---------------------------------------------------------
// Function : Disc3dContour
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::Disc3dContour(const int /*nbp*/,
const int iordre,
NCollection_Sequence<gp_XYZ>& Seq3d)
{
#ifdef OCCT_DEBUG
if (Seq3d.Length() != 4)
std::cout << "nbp should be equal to 4 for Disc3dContour" << std::endl;
if (iordre != 0 && iordre != 1)
std::cout << "incorrect order for Disc3dContour" << std::endl;
#endif
// initialization
Seq3d.Clear();
// sampling in "cosine" + 3 points on each interval
double u1, v1, u2, v2;
mySurfInit->Bounds(u1, v1, u2, v2);
GeomAdaptor_Surface Surf(mySurfInit);
myProj.Initialize(Surf, u1, v1, u2, v2, Surf.UResolution(myTol3d), Surf.VResolution(myTol3d));
int NTCurve = myLinCont->Length();
int NTPntCont = myPntCont->Length();
// gp_Pnt2d P2d;
gp_Pnt P3d;
gp_Vec v1h, v2h, v3h;
gp_XYZ Pos;
int i;
for (i = 1; i <= NTPntCont; i++)
{
if (myPntCont->Value(i)->Order() != -1)
{
if (iordre == 0)
{
myPntCont->Value(i)->D0(P3d);
Pos.SetX(P3d.X());
Pos.SetY(P3d.Y());
Pos.SetZ(P3d.Z());
Seq3d.Append(Pos);
}
else
{
myPntCont->Value(i)->D1(P3d, v1h, v2h);
v3h = v1h ^ v2h;
Pos.SetX(v3h.X());
Pos.SetY(v3h.Y());
Pos.SetZ(v3h.Z());
Seq3d.Append(Pos);
}
}
}
for (i = 1; i <= NTCurve; i++)
{
if (myLinCont->Value(i)->Order() != -1)
{
int NbPt = myParCont->Value(i).Length();
// first constraint point (j=0)
// int NbPt=myParCont->Length();
if (iordre == 0)
{
myLinCont->Value(i)->D0(myParCont->Value(i).Value(1), P3d);
Pos.SetX(P3d.X());
Pos.SetY(P3d.Y());
Pos.SetZ(P3d.Z());
Seq3d.Append(Pos);
}
else
{
myLinCont->Value(i)->D1(myParCont->Value(i).Value(1), P3d, v1h, v2h);
v3h = v1h ^ v2h;
Pos.SetX(v3h.X());
Pos.SetY(v3h.Y());
Pos.SetZ(v3h.Z());
Seq3d.Append(Pos);
}
for (int j = 2; j < NbPt; j++)
{
double Uj = myParCont->Value(i).Value(j), Ujp1 = myParCont->Value(i).Value(j + 1);
if (iordre == 0)
{
// point 1/4 previous
myLinCont->Value(i)->D0((Ujp1 + 3 * Uj) / 4, P3d);
Pos.SetX(P3d.X());
Pos.SetY(P3d.Y());
Pos.SetZ(P3d.Z());
Seq3d.Append(Pos);
// point 1/2 previous
myLinCont->Value(i)->D0((Ujp1 + Uj) / 2, P3d);
Pos.SetX(P3d.X());
Pos.SetY(P3d.Y());
Pos.SetZ(P3d.Z());
Seq3d.Append(Pos);
// point 3/4 previous
myLinCont->Value(i)->D0((3 * Ujp1 + Uj) / 4, P3d);
Pos.SetX(P3d.X());
Pos.SetY(P3d.Y());
Pos.SetZ(P3d.Z());
Seq3d.Append(Pos);
// current constraint point
myLinCont->Value(i)->D0(Ujp1, P3d);
Pos.SetX(P3d.X());
Pos.SetY(P3d.Y());
Pos.SetZ(P3d.Z());
Seq3d.Append(Pos);
}
else
{
// point 1/4 previous
myLinCont->Value(i)->D1((Ujp1 + 3 * Uj) / 4, P3d, v1h, v2h);
v3h = v1h ^ v2h;
Pos.SetX(v3h.X());
Pos.SetY(v3h.Y());
Pos.SetZ(v3h.Z());
Seq3d.Append(Pos);
// point 1/2 previous
myLinCont->Value(i)->D1((Ujp1 + Uj) / 2, P3d, v1h, v2h);
v3h = v1h ^ v2h;
Pos.SetX(v3h.X());
Pos.SetY(v3h.Y());
Pos.SetZ(v3h.Z());
Seq3d.Append(Pos);
// point 3/4 previous
myLinCont->Value(i)->D1((3 * Ujp1 + Uj) / 4, P3d, v1h, v2h);
v3h = v1h ^ v2h;
Pos.SetX(v3h.X());
Pos.SetY(v3h.Y());
Pos.SetZ(v3h.Z());
Seq3d.Append(Pos);
// current constraint point
myLinCont->Value(i)->D1(Ujp1, P3d, v1h, v2h);
v3h = v1h ^ v2h;
Pos.SetX(v3h.X());
Pos.SetY(v3h.Y());
Pos.SetZ(v3h.Z());
Seq3d.Append(Pos);
}
}
}
}
}
//---------------------------------------------------------
// Function : IsDone
//---------------------------------------------------------
bool GeomPlate_BuildPlateSurface::IsDone() const
{
return myPlate.IsDone();
}
//---------------------------------------------------------
// Function : Surface
//---------------------------------------------------------
occ::handle<GeomPlate_Surface> GeomPlate_BuildPlateSurface::Surface() const
{
return myGeomPlateSurface;
}
//---------------------------------------------------------
// Function : SurfInit
//---------------------------------------------------------
occ::handle<Geom_Surface> GeomPlate_BuildPlateSurface::SurfInit() const
{
return mySurfInit;
}
//---------------------------------------------------------
// Function : Sense
//---------------------------------------------------------
occ::handle<NCollection_HArray1<int>> GeomPlate_BuildPlateSurface::Sense() const
{
int NTCurve = myLinCont->Length();
occ::handle<NCollection_HArray1<int>> Sens = new NCollection_HArray1<int>(1, NTCurve);
for (int i = 1; i <= NTCurve; i++)
{
Sens->SetValue(i, mySense->Value(myInitOrder->Value(i)));
}
return Sens;
}
//---------------------------------------------------------
// Function : Curve2d
//---------------------------------------------------------
occ::handle<NCollection_HArray1<occ::handle<Geom2d_Curve>>> GeomPlate_BuildPlateSurface::Curves2d()
const
{
int NTCurve = myLinCont->Length();
occ::handle<NCollection_HArray1<occ::handle<Geom2d_Curve>>> C2dfin =
new NCollection_HArray1<occ::handle<Geom2d_Curve>>(1, NTCurve);
for (int i = 1; i <= NTCurve; i++)
{
C2dfin->SetValue(i, myLinCont->Value(myInitOrder->Value(i))->Curve2dOnSurf());
}
return C2dfin;
}
//---------------------------------------------------------
// Function : Order
//---------------------------------------------------------
occ::handle<NCollection_HArray1<int>> GeomPlate_BuildPlateSurface::Order() const
{
occ::handle<NCollection_HArray1<int>> result =
new NCollection_HArray1<int>(1, myLinCont->Length());
for (int i = 1; i <= myLinCont->Length(); i++)
{
result->SetValue(myInitOrder->Value(i), i);
}
return result;
}
//---------------------------------------------------------
// Function : G0Error
//---------------------------------------------------------
double GeomPlate_BuildPlateSurface::G0Error() const
{
return myG0Error;
}
//---------------------------------------------------------
// Function : G1Error
//---------------------------------------------------------
double GeomPlate_BuildPlateSurface::G1Error() const
{
return myG1Error;
}
//---------------------------------------------------------
// Function : G2Error
//---------------------------------------------------------
double GeomPlate_BuildPlateSurface::G2Error() const
{
return myG2Error;
}
//=================================================================================================
double GeomPlate_BuildPlateSurface::G0Error(const int Index)
{
occ::handle<NCollection_HArray1<double>> tdistance =
new NCollection_HArray1<double>(1, myNbPtsOnCur);
occ::handle<NCollection_HArray1<double>> tangle =
new NCollection_HArray1<double>(1, myNbPtsOnCur);
occ::handle<NCollection_HArray1<double>> tcurvature =
new NCollection_HArray1<double>(1, myNbPtsOnCur);
EcartContraintesMil(Index, tdistance, tangle, tcurvature);
double MaxDistance = 0.;
for (int i = 1; i <= myNbPtsOnCur; i++)
{
if (tdistance->Value(i) > MaxDistance)
{
MaxDistance = tdistance->Value(i);
}
}
return MaxDistance;
}
//=================================================================================================
double GeomPlate_BuildPlateSurface::G1Error(const int Index)
{
occ::handle<NCollection_HArray1<double>> tdistance =
new NCollection_HArray1<double>(1, myNbPtsOnCur);
occ::handle<NCollection_HArray1<double>> tangle =
new NCollection_HArray1<double>(1, myNbPtsOnCur);
occ::handle<NCollection_HArray1<double>> tcurvature =
new NCollection_HArray1<double>(1, myNbPtsOnCur);
EcartContraintesMil(Index, tdistance, tangle, tcurvature);
double MaxAngle = 0.;
for (int i = 1; i <= myNbPtsOnCur; i++)
{
if (tangle->Value(i) > MaxAngle)
{
MaxAngle = tangle->Value(i);
}
}
return MaxAngle;
}
//=================================================================================================
double GeomPlate_BuildPlateSurface::G2Error(const int Index)
{
occ::handle<NCollection_HArray1<double>> tdistance =
new NCollection_HArray1<double>(1, myNbPtsOnCur);
occ::handle<NCollection_HArray1<double>> tangle =
new NCollection_HArray1<double>(1, myNbPtsOnCur);
occ::handle<NCollection_HArray1<double>> tcurvature =
new NCollection_HArray1<double>(1, myNbPtsOnCur);
EcartContraintesMil(Index, tdistance, tangle, tcurvature);
double MaxCurvature = 0.;
for (int i = 1; i <= myNbPtsOnCur; i++)
{
if (tcurvature->Value(i) > MaxCurvature)
{
MaxCurvature = tcurvature->Value(i);
}
}
return MaxCurvature;
}
occ::handle<GeomPlate_CurveConstraint> GeomPlate_BuildPlateSurface::CurveConstraint(
const int order) const
{
return myLinCont->Value(order);
}
occ::handle<GeomPlate_PointConstraint> GeomPlate_BuildPlateSurface::PointConstraint(
const int order) const
{
return myPntCont->Value(order);
}
//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
// =========================================================
// P R I V A T E M E T H O D S
// =========================================================
//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
//=======================================================================
// Function : CourbeJointive
// Purpose : Create a chain of curves to calculate the
// initial surface with the method of max flow.
// Return true if it is a closed contour.
//=======================================================================
bool GeomPlate_BuildPlateSurface::CourbeJointive(const double tolerance)
{
int nbf = myLinCont->Length();
double Ufinal1, Uinit1, Ufinal2, Uinit2;
mySense = new NCollection_HArray1<int>(1, nbf, 0);
bool result = true;
int j = 1, i;
gp_Pnt P1, P2;
while (j <= (myNbBounds - 1))
{
int a = 0;
i = j + 1;
if (i > myNbBounds)
{
result = false;
a = 2;
}
while (a < 1)
{
if (i > myNbBounds)
{
result = false;
a = 2;
}
else
{
Uinit1 = myLinCont->Value(j)->FirstParameter();
Ufinal1 = myLinCont->Value(j)->LastParameter();
Uinit2 = myLinCont->Value(i)->FirstParameter();
Ufinal2 = myLinCont->Value(i)->LastParameter();
if (mySense->Value(j) == 1)
{
Ufinal1 = Uinit1;
}
myLinCont->Value(j)->D0(Ufinal1, P1);
myLinCont->Value(i)->D0(Uinit2, P2);
if (P1.Distance(P2) < tolerance)
{
if (i != j + 1)
{
occ::handle<GeomPlate_CurveConstraint> tampon = myLinCont->Value(j + 1);
myLinCont->SetValue(j + 1, myLinCont->Value(i));
myLinCont->SetValue(i, tampon);
int Tmp = myInitOrder->Value(j + 1);
// See function TrierTab for the functioning of myInitOrder
myInitOrder->SetValue(j + 1, myInitOrder->Value(i));
myInitOrder->SetValue(i, Tmp);
};
a = 2;
mySense->SetValue(j + 1, 0);
}
else
{
myLinCont->Value(i)->D0(Ufinal2, P2);
if (P1.Distance(P2) < tolerance)
{
if (i != j + 1)
{
occ::handle<GeomPlate_CurveConstraint> tampon = myLinCont->Value(j + 1);
myLinCont->SetValue(j + 1, myLinCont->Value(i));
myLinCont->SetValue(i, tampon);
int Tmp = myInitOrder->Value(j + 1);
// See function TrierTab for the functioning of myInitOrder
myInitOrder->SetValue(j + 1, myInitOrder->Value(i));
myInitOrder->SetValue(i, Tmp);
};
a = 2;
mySense->SetValue(j + 1, 1);
}
}
}
i++;
}
j++;
}
Uinit1 = myLinCont->Value(myNbBounds)->FirstParameter();
Ufinal1 = myLinCont->Value(myNbBounds)->LastParameter();
Uinit2 = myLinCont->Value(1)->FirstParameter();
Ufinal2 = myLinCont->Value(1)->LastParameter();
myLinCont->Value(myNbBounds)->D0(Ufinal1, P1);
myLinCont->Value(1)->D0(Uinit2, P2);
if ((mySense->Value(myNbBounds) == 0) && (P1.Distance(P2) < tolerance))
{
return ((result));
}
myLinCont->Value(myNbBounds)->D0(Uinit1, P1);
if ((mySense->Value(myNbBounds) == 1) && (P1.Distance(P2) < tolerance))
{
return ((result));
}
else
{
return false;
}
}
//=================================================================================================
// Calculate the initial surface either by the method of max flow or by
// the method of the plane of inertia if the contour is not closed or if
// there are point constraints.
//-------------------------------------------------------------------------
void GeomPlate_BuildPlateSurface::ComputeSurfInit(const Message_ProgressRange& theProgress)
{
int nopt = 2, popt = 2, Np = 1;
bool isHalfSpace = true;
double LinTol = 0.001, AngTol = 0.001; // AngTol = 0.0001; //LinTol = 0.0001
// Option to calculate the initial plane
int NTLinCont = myLinCont->Length(), NTPntCont = myPntCont->Length();
// Table of transformation to preserve the initial order see TrierTab
if (NTLinCont != 0)
{
myInitOrder = new NCollection_HArray1<int>(1, NTLinCont);
for (int i = 1; i <= NTLinCont; i++)
{
myInitOrder->SetValue(i, i);
}
}
bool CourbeJoint = (NTLinCont != 0) && CourbeJointive(myTol3d);
if (CourbeJoint && IsOrderG1())
{
nopt = 3;
// Table contains the cloud of points for calculation of the plane
int NbPoint = 20, Discr = NbPoint / 4, pnum = 0;
occ::handle<NCollection_HArray1<gp_Pnt>> Pts =
new NCollection_HArray1<gp_Pnt>(1, (NbPoint + 1) * NTLinCont + NTPntCont);
NCollection_Sequence<gp_Vec> Vecs, NewVecs;
NCollection_Sequence<GeomPlate_Aij> Aset;
double Uinit, Ufinal, Uif;
gp_Vec LastVec;
int i;
for (i = 1; i <= NTLinCont; i++)
{
int Order = myLinCont->Value(i)->Order();
NewVecs.Clear();
Uinit = myLinCont->Value(i)->FirstParameter();
Ufinal = myLinCont->Value(i)->LastParameter();
Uif = Ufinal - Uinit;
if (mySense->Value(i) == 1)
{
Uinit = Ufinal;
Uif = -Uif;
}
gp_Vec Vec1, Vec2, Normal;
bool ToReverse = false;
if (i > 1 && Order >= GeomAbs_G1)
{
gp_Pnt P;
myLinCont->Value(i)->D1(Uinit, P, Vec1, Vec2);
Normal = Vec1 ^ Vec2;
if (LastVec.IsOpposite(Normal, AngTol))
{
ToReverse = true;
}
}
for (int j = 0; j <= NbPoint; j++)
{ // Number of points per curve = 20
// Linear distribution
double Inter = j * Uif / (NbPoint);
if (Order < GeomAbs_G1 || j % Discr != 0)
{
myLinCont->Value(i)->D0(Uinit + Inter, Pts->ChangeValue(++pnum));
}
else
{
myLinCont->Value(i)->D1(Uinit + Inter, Pts->ChangeValue(++pnum), Vec1, Vec2);
Normal = Vec1 ^ Vec2;
Normal.Normalize();
if (ToReverse)
{
Normal.Reverse();
}
bool isNew = true;
int k;
for (k = 1; k <= Vecs.Length(); k++)
{
if (Vecs(k).IsEqual(Normal, LinTol, AngTol))
{
isNew = false;
break;
}
}
if (isNew)
{
for (k = 1; k <= NewVecs.Length(); k++)
{
if (NewVecs(k).IsEqual(Normal, LinTol, AngTol))
{
isNew = false;
break;
}
}
}
if (isNew)
{
NewVecs.Append(Normal);
}
}
}
if (Order >= GeomAbs_G1)
{
isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace(NewVecs, Vecs, Aset, LinTol, AngTol);
if (!isHalfSpace)
{
break;
}
LastVec = Normal;
}
} // for (i = 1; i <= NTLinCont; i++)
if (isHalfSpace)
{
for (i = 1; i <= NTPntCont; i++)
{
int Order = myPntCont->Value(i)->Order();
NewVecs.Clear();
gp_Vec Vec1, Vec2, Normal;
if (Order < GeomAbs_G1)
{
myPntCont->Value(i)->D0(Pts->ChangeValue(++pnum));
}
else
{
myPntCont->Value(i)->D1(Pts->ChangeValue(++pnum), Vec1, Vec2);
Normal = Vec1 ^ Vec2;
Normal.Normalize();
bool isNew = true;
for (int k = 1; k <= Vecs.Length(); k++)
{
if (Vecs(k).IsEqual(Normal, LinTol, AngTol))
{
isNew = false;
break;
}
}
if (isNew)
{
NewVecs.Append(Normal);
isHalfSpace =
GeomPlate_BuildAveragePlane::HalfSpace(NewVecs, Vecs, Aset, LinTol, AngTol);
if (!isHalfSpace)
{
NewVecs(1).Reverse();
isHalfSpace =
GeomPlate_BuildAveragePlane::HalfSpace(NewVecs, Vecs, Aset, LinTol, AngTol);
}
if (!isHalfSpace)
{
break;
}
}
}
} // for (i = 1; i <= NTPntCont; i++)
if (isHalfSpace)
{
bool NullExist = true;
while (NullExist)
{
NullExist = false;
for (i = 1; i <= Vecs.Length(); i++)
{
if (Vecs(i).SquareMagnitude() == 0.)
{
NullExist = true;
Vecs.Remove(i);
break;
}
}
}
GeomPlate_BuildAveragePlane BAP(Vecs, Pts);
double u1, u2, v1, v2;
BAP.MinMaxBox(u1, u2, v1, v2);
// The space is greater for projections
double du = u2 - u1;
double dv = v2 - v1;
u1 -= du;
u2 += du;
v1 -= dv;
v2 += dv;
mySurfInit = new Geom_RectangularTrimmedSurface(BAP.Plane(), u1, u2, v1, v2);
}
} // if (isHalfSpace)
if (!isHalfSpace)
{
#ifdef OCCT_DEBUG
std::cout << std::endl << "Normals are not in half space" << std::endl << std::endl;
#endif
myIsLinear = false;
nopt = 2;
}
} // if (NTLinCont != 0 && (CourbeJoint = CourbeJointive( myTol3d )) && IsOrderG1())
if (NTLinCont != 0)
{
TrierTab(myInitOrder); // Reorder the table of transformations
}
if (nopt != 3)
{
if (NTPntCont != 0)
{
nopt = 1; // Calculate by the method of plane of inertia
}
else if (!CourbeJoint || NTLinCont != myNbBounds)
{ // throw Standard_Failure("Curves are not joined");
#ifdef OCCT_DEBUG
std::cout << "WARNING : Curves are non-adjacent with tolerance " << myTol3d << std::endl;
#endif
nopt = 1;
}
double LenT = 0;
int Npt = 0;
int NTPoint = 20 * NTLinCont;
int i;
for (i = 1; i <= NTLinCont; i++)
{
LenT += myLinCont->Value(i)->Length();
}
for (i = 1; i <= NTLinCont; i++)
{
int NbPoint = (int)(NTPoint * (myLinCont->Value(i)->Length()) / LenT);
if (NbPoint < 10)
{
NbPoint = 10;
}
(void)Npt; // unused but set for debug
Npt += NbPoint;
}
// Table containing a cloud of points for calculation of the plane
occ::handle<NCollection_HArray1<gp_Pnt>> Pts =
new NCollection_HArray1<gp_Pnt>(1, 20 * NTLinCont + NTPntCont);
int NbPoint = 20;
double Uinit, Ufinal, Uif;
for (i = 1; i <= NTLinCont; i++)
{
Uinit = myLinCont->Value(i)->FirstParameter();
Ufinal = myLinCont->Value(i)->LastParameter();
Uif = Ufinal - Uinit;
if (mySense->Value(i) == 1)
{
Uinit = Ufinal;
Uif = -Uif;
}
for (int j = 0; j < NbPoint; j++)
{ // Number of points per curve = 20
// Linear distribution
double Inter = j * Uif / (NbPoint);
gp_Pnt P;
myLinCont->Value(i)->D0(Uinit + Inter, P);
Pts->SetValue(Np++, P);
}
}
for (i = 1; i <= NTPntCont; i++)
{
gp_Pnt P;
myPntCont->Value(i)->D0(P);
Pts->SetValue(Np++, P);
}
if (!CourbeJoint)
{
myNbBounds = 0;
}
GeomPlate_BuildAveragePlane BAP(Pts, NbPoint * myNbBounds, myTol3d / 1000, popt, nopt);
if (!BAP.IsPlane())
{
#ifdef OCCT_DEBUG
std::cout << "WARNING : GeomPlate : the initial surface is not a plane." << std::endl;
#endif
return;
}
double u1, u2, v1, v2;
BAP.MinMaxBox(u1, u2, v1, v2);
// The space is greater for projections
double du = u2 - u1;
double dv = v2 - v1;
u1 -= du;
u2 += du;
v1 -= dv;
v2 += dv;
mySurfInit = new Geom_RectangularTrimmedSurface(BAP.Plane(), u1, u2, v1, v2);
} // if (nopt != 3)
// Comparing metrics of curves and projected curves
if (NTLinCont != 0 && myIsLinear)
{
occ::handle<Geom_Surface> InitPlane =
(occ::down_cast<Geom_RectangularTrimmedSurface>(mySurfInit))->BasisSurface();
// clang-format off
double Ratio = 0., R1 = 2., R2 = 0.6; //R1 = 3, R2 = 0.5;//R1 = 1.4, R2 = 0.8; //R1 = 5., R2 = 0.2;
// clang-format on
occ::handle<GeomAdaptor_Surface> hsur = new GeomAdaptor_Surface(InitPlane);
int NbPoint = 20;
// gp_Pnt P;
// gp_Vec DerC, DerCproj, DU, DV;
// gp_Pnt2d P2d;
// gp_Vec2d DProj;
for (int i = 1; i <= NTLinCont && myIsLinear; i++)
{
double FirstPar = myLinCont->Value(i)->FirstParameter();
double LastPar = myLinCont->Value(i)->LastParameter();
double Uif = (LastPar - FirstPar) / (NbPoint);
occ::handle<Adaptor3d_Curve> Curve = myLinCont->Value(i)->Curve3d();
occ::handle<ProjLib_HCompProjectedCurve> ProjCurve =
new ProjLib_HCompProjectedCurve(hsur, Curve, myTol3d, myTol3d);
Adaptor3d_CurveOnSurface AProj(ProjCurve, hsur);
gp_Pnt P;
gp_Vec DerC, DerCproj;
for (int j = 1; j < NbPoint && myIsLinear; j++)
{
double Inter = FirstPar + j * Uif;
Curve->D1(Inter, P, DerC);
AProj.D1(Inter, P, DerCproj);
double A1 = DerC.Magnitude();
double A2 = DerCproj.Magnitude();
if (A2 <= 1.e-20)
{
Ratio = 1.e20;
}
else
{
Ratio = A1 / A2;
}
if (Ratio > R1 || Ratio < R2)
{
myIsLinear = false;
break;
}
}
}
#ifdef OCCT_DEBUG
if (!myIsLinear)
std::cout << "Metrics are too different :" << Ratio << std::endl;
#endif
// myIsLinear = true; // !!
} // comparing metrics of curves and projected curves
if (!myIsLinear)
{
myPlanarSurfInit = mySurfInit;
double u1, v1, u2, v2;
mySurfInit->Bounds(u1, v1, u2, v2);
GeomAdaptor_Surface Surf(mySurfInit);
myTolU = Surf.UResolution(myTol3d);
myTolV = Surf.VResolution(myTol3d);
myProj.Initialize(Surf, u1, v1, u2, v2, myTolU, myTolV);
//======================================================================
// Projection of curves
//======================================================================
int i;
for (i = 1; i <= NTLinCont; i++)
{
if (myLinCont->Value(i)->Curve2dOnSurf().IsNull())
{
myLinCont->ChangeValue(i)->SetCurve2dOnSurf(ProjectCurve(myLinCont->Value(i)->Curve3d()));
}
}
//======================================================================
// Projection of points
//======================================================================
for (i = 1; i <= NTPntCont; i++)
{
gp_Pnt P;
myPntCont->Value(i)->D0(P);
if (!myPntCont->Value(i)->HasPnt2dOnSurf())
{
myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P));
}
}
//======================================================================
// Number of points by curve
//======================================================================
if ((NTLinCont != 0) && (myNbPtsOnCur != 0))
{
CalculNbPtsInit();
}
//======================================================================
// Management of incompatibilities between curves
//======================================================================
occ::handle<NCollection_HArray1<NCollection_Sequence<double>>> PntInter;
occ::handle<NCollection_HArray1<NCollection_Sequence<double>>> PntG1G1;
if (NTLinCont != 0)
{
PntInter = new NCollection_HArray1<NCollection_Sequence<double>>(1, NTLinCont);
PntG1G1 = new NCollection_HArray1<NCollection_Sequence<double>>(1, NTLinCont);
Intersect(PntInter, PntG1G1);
}
//====================================================================
// Discretization of curves
//====================================================================
Discretise(PntInter, PntG1G1);
//====================================================================
// Preparation of points of constraint for plate
//====================================================================
LoadCurve(0, 0);
if (myPntCont->Length() != 0)
{
LoadPoint(0, 0);
}
//====================================================================
// Construction of the surface
//====================================================================
Message_ProgressScope aPS(theProgress, "ComputeSurfInit", 1);
myPlate.SolveTI(2, ComputeAnisotropie(), aPS.Next());
if (theProgress.UserBreak())
{
return;
}
if (!myPlate.IsDone())
{
#ifdef OCCT_DEBUG
std::cout << "WARNING : GeomPlate : calculation of Plate failed" << std::endl;
#endif
return;
}
myGeomPlateSurface = new GeomPlate_Surface(mySurfInit, myPlate);
GeomPlate_MakeApprox App(myGeomPlateSurface, myTol3d, 1, 3, 15 * myTol3d, -1, GeomAbs_C0);
mySurfInit = App.Surface();
mySurfInitIsGive = true;
myPlate.Init(); // Reset
for (i = 1; i <= NTLinCont; i++)
{
occ::handle<Geom2d_Curve> NullCurve;
NullCurve.Nullify();
myLinCont->ChangeValue(i)->SetCurve2dOnSurf(NullCurve);
}
}
}
//---------------------------------------------------------
// Function : Intersect
//---------------------------------------------------------
// Find intersections between 2d curves
// If the intersection is compatible (in cases G1-G1)
// remove the point on one of two curves
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::Intersect(
occ::handle<NCollection_HArray1<NCollection_Sequence<double>>>& PntInter,
occ::handle<NCollection_HArray1<NCollection_Sequence<double>>>& PntG1G1)
{
int NTLinCont = myLinCont->Length();
Geom2dInt_GInter Intersection;
Geom2dAdaptor_Curve Ci, Cj;
IntRes2d_IntersectionPoint int2d;
gp_Pnt P1, P2;
gp_Pnt2d P2d;
gp_Vec2d V2d;
// if (!mySurfInitIsGive)
for (int i = 1; i <= NTLinCont; i++)
{
// double NbPnt_i=myLinCont->Value(i)->NbPoints();
// Find the intersection with each curve including the curve itself
Ci.Load(myLinCont->Value(i)->Curve2dOnSurf());
for (int j = i; j <= NTLinCont; j++)
{
Cj.Load(myLinCont->Value(j)->Curve2dOnSurf());
if (i == j)
{
Intersection.Perform(Ci, myTol2d * 10, myTol2d * 10);
}
else
{
Intersection.Perform(Ci, Cj, myTol2d * 10, myTol2d * 10);
}
if (!Intersection.IsEmpty())
{ // there is one intersection
int nbpt = Intersection.NbPoints();
// number of points of intersection
for (int k = 1; k <= nbpt; k++)
{
int2d = Intersection.Point(k);
myLinCont->Value(i)->D0(int2d.ParamOnFirst(), P1);
myLinCont->Value(j)->D0(int2d.ParamOnSecond(), P2);
#ifdef OCCT_DEBUG
if (Affich > 1)
{
std::cout << " Intersection " << k << " entre " << i << " &" << j << std::endl;
std::cout << " Distance = " << P1.Distance(P2) << std::endl;
}
#endif
if (P1.Distance(P2) < myTol3d)
{ // 2D intersection corresponds to close 3D points.
// Note the interval, in which the point needs to be removed
// to avoid duplications, which cause
// error in plate. The point on curve i is removed;
// the point on curve j is preserved;
// the length of interval is a length 2d
// corresponding in 3d to myTol3d
double tolint = Ci.Resolution(myTol3d);
Ci.D1(int2d.ParamOnFirst(), P2d, V2d);
double aux = V2d.Magnitude();
if (aux > 1.e-7)
{
aux = myTol3d / aux;
if (aux > 100 * tolint)
{
tolint *= 100;
}
else
{
tolint = aux;
}
}
else
{
tolint *= 100;
}
PntInter->ChangeValue(i).Append(int2d.ParamOnFirst() - tolint);
PntInter->ChangeValue(i).Append(int2d.ParamOnFirst() + tolint);
// If G1-G1
if ((myLinCont->Value(i)->Order() == 1) && (myLinCont->Value(j)->Order() == 1))
{
gp_Vec v11, v12, v13, v14, v15, v16, v21, v22, v23, v24, v25, v26;
myLinCont->Value(i)->D2(int2d.ParamOnFirst(), P1, v11, v12, v13, v14, v15);
myLinCont->Value(j)->D2(int2d.ParamOnSecond(), P2, v21, v22, v23, v24, v25);
v16 = v11 ^ v12;
v26 = v21 ^ v22;
double ant = v16.Angle(v26);
if (ant > (M_PI / 2))
{
ant = M_PI - ant;
}
if ((std::abs(v16 * v15 - v16 * v25) > (myTol3d / 1000))
|| (std::abs(ant) > myTol3d / 1000))
// Non-compatible ==> remove zone in constraint G1
// corresponding to 3D tolerance of 0.01
{
double coin;
double Tol = 100 * myTol3d;
double A1;
gp_Pnt2d P1temp, P2temp;
gp_Vec2d V1, V2;
myLinCont->Value(i)->Curve2dOnSurf()->D1(int2d.ParamOnFirst(), P1temp, V1);
myLinCont->Value(j)->Curve2dOnSurf()->D1(int2d.ParamOnSecond(), P2temp, V2);
A1 = V1.Angle(V2);
if (A1 > (M_PI / 2))
{
A1 = M_PI - A1;
}
if (std::abs(std::abs(A1) - M_PI) < myTolAng)
{
Tol = 100000 * myTol3d;
}
#ifdef OCCT_DEBUG
if (Affich)
std::cout << "Angle between curves " << i << "," << j << " "
<< std::abs(std::abs(A1) - M_PI) << std::endl;
#endif
coin = Ci.Resolution(Tol);
double Par1 = int2d.ParamOnFirst() - coin, Par2 = int2d.ParamOnFirst() + coin;
// Storage of the interval for curve i
PntG1G1->ChangeValue(i).Append(Par1);
PntG1G1->ChangeValue(i).Append(Par2);
coin = Cj.Resolution(Tol);
Par1 = int2d.ParamOnSecond() - coin;
Par2 = int2d.ParamOnSecond() + coin;
// Storage of the interval for curve j
PntG1G1->ChangeValue(j).Append(Par1);
PntG1G1->ChangeValue(j).Append(Par2);
}
}
// If G0-G1
if ((myLinCont->Value(i)->Order() == 0 && myLinCont->Value(j)->Order() == 1)
|| (myLinCont->Value(i)->Order() == 1 && myLinCont->Value(j)->Order() == 0))
{
gp_Vec vec, vecU, vecV, N;
if (myLinCont->Value(i)->Order() == 0)
{
occ::handle<Adaptor3d_Curve> theCurve = myLinCont->Value(i)->Curve3d();
theCurve->D1(int2d.ParamOnFirst(), P1, vec);
myLinCont->Value(j)->D1(int2d.ParamOnSecond(), P2, vecU, vecV);
}
else
{
occ::handle<Adaptor3d_Curve> theCurve = myLinCont->Value(j)->Curve3d();
theCurve->D1(int2d.ParamOnSecond(), P2, vec);
myLinCont->Value(i)->D1(int2d.ParamOnFirst(), P1, vecU, vecV);
}
N = vecU ^ vecV;
double Angle = vec.Angle(N);
Angle = std::abs(M_PI / 2 - Angle);
if (Angle > myTolAng / 10.) //????????? //if (std::abs( scal ) > myTol3d/100)
{
// Non-compatible ==> one removes zone in constraint G0 and G1
// corresponding to 3D tolerance of 0.01
double coin;
double Tol = 100 * myTol3d;
double A1;
gp_Pnt2d P1temp, P2temp;
gp_Vec2d V1, V2;
myLinCont->Value(i)->Curve2dOnSurf()->D1(int2d.ParamOnFirst(), P1temp, V1);
myLinCont->Value(j)->Curve2dOnSurf()->D1(int2d.ParamOnSecond(), P2temp, V2);
A1 = V1.Angle(V2);
if (A1 > M_PI / 2)
{
A1 = M_PI - A1;
}
if (std::abs(std::abs(A1) - M_PI) < myTolAng)
{
Tol = 100000 * myTol3d;
}
#ifdef OCCT_DEBUG
if (Affich)
std::cout << "Angle entre Courbe " << i << "," << j << " "
<< std::abs(std::abs(A1) - M_PI) << std::endl;
#endif
if (myLinCont->Value(i)->Order() == 1)
{
coin = Ci.Resolution(Tol);
coin *= Angle / myTolAng * 10.;
#ifdef OCCT_DEBUG
std::cout << std::endl << "coin = " << coin << std::endl;
#endif
double Par1 = int2d.ParamOnFirst() - coin;
double Par2 = int2d.ParamOnFirst() + coin;
// Storage of the interval for curve i
PntG1G1->ChangeValue(i).Append(Par1);
PntG1G1->ChangeValue(i).Append(Par2);
}
else
{
coin = Cj.Resolution(Tol);
coin *= Angle / myTolAng * 10.;
#ifdef OCCT_DEBUG
std::cout << std::endl << "coin = " << coin << std::endl;
#endif
double Par1 = int2d.ParamOnSecond() - coin;
double Par2 = int2d.ParamOnSecond() + coin;
// Storage of the interval for curve j
PntG1G1->ChangeValue(j).Append(Par1);
PntG1G1->ChangeValue(j).Append(Par2);
}
}
}
} // if (P1.Distance( P2 ) < myTol3d)
else
{
// 2D intersection corresponds to extended 3D points.
// Note the interval where it is necessary to remove
// the points to avoid duplications causing
// error in plate. The point on curve i is removed,
// the point on curve j is preserved.
// The length of interval is 2D length
// corresponding to the distance of points in 3D to myTol3d
double tolint, Dist;
Dist = P1.Distance(P2);
tolint = Ci.Resolution(Dist);
PntInter->ChangeValue(i).Append(int2d.ParamOnFirst() - tolint);
PntInter->ChangeValue(i).Append(int2d.ParamOnFirst() + tolint);
if (j != i)
{
tolint = Cj.Resolution(Dist);
PntInter->ChangeValue(j).Append(int2d.ParamOnSecond() - tolint);
PntInter->ChangeValue(j).Append(int2d.ParamOnSecond() + tolint);
}
#ifdef OCCT_DEBUG
std::cout << "Attention: Two points 3d have the same projection dist = " << Dist
<< std::endl;
#endif
}
}
}
}
}
}
//---------------------------------------------------------
// Function : Discretize
//---------------------------------------------------------
// Discretize curves according to parameters
// the table of sequences Parcont contains all
// parameter of points on curves
// Field myPlateCont contains parameter of points on a plate;
// it excludes duplicate points and incompatible zones.
// The first part corresponds to verification of compatibility
// and to removal of duplicate points.
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::Discretise(
const occ::handle<NCollection_HArray1<NCollection_Sequence<double>>>& PntInter,
const occ::handle<NCollection_HArray1<NCollection_Sequence<double>>>& PntG1G1)
{
int NTLinCont = myLinCont->Length();
bool ACR;
occ::handle<Geom2d_Curve> C2d;
Geom2dAdaptor_Curve AC2d;
// occ::handle<Adaptor_HCurve2d> HC2d;
occ::handle<Law_Interpol> acrlaw = new (Law_Interpol)();
myPlateCont = new NCollection_HArray1<NCollection_Sequence<double>>(1, NTLinCont);
myParCont = new NCollection_HArray1<NCollection_Sequence<double>>(1, NTLinCont);
//===========================================================================
// Construction of the table containing parameters of constraint points
//===========================================================================
double Uinit, Ufinal, Length2d = 0, Inter;
double CurLength;
int NbPnt_i, NbPtInter, NbPtG1G1;
occ::handle<GeomPlate_CurveConstraint> LinCont;
for (int i = 1; i <= NTLinCont; i++)
{
LinCont = myLinCont->Value(i);
Uinit = LinCont->FirstParameter();
Ufinal = LinCont->LastParameter();
// HC2d=LinCont->ProjectedCurve();
// if(HC2d.IsNull())
// ACR = (!HC2d.IsNull() || !C2d.IsNull());
C2d = LinCont->Curve2dOnSurf();
ACR = (!C2d.IsNull());
if (ACR)
{
// Construct a law close to curvilinear abscissa
if (!C2d.IsNull())
{
AC2d.Load(C2d);
}
// AC2d.Load(LinCont->Curve2dOnSurf());
int ii, Nbint = 20;
double U;
NCollection_Array1<gp_Pnt2d> tabP2d(1, Nbint + 1);
tabP2d(1).SetY(Uinit);
tabP2d(1).SetX(0.);
tabP2d(Nbint + 1).SetY(Ufinal);
/* if (!HC2d.IsNull())
Length2d = GCPnts_AbscissaPoint::Length(HC2d->Curve2d(), Uinit, Ufinal);
else*/
Length2d = GCPnts_AbscissaPoint::Length(AC2d, Uinit, Ufinal);
tabP2d(Nbint + 1).SetX(Length2d);
for (ii = 2; ii <= Nbint; ii++)
{
U = Uinit + (Ufinal - Uinit) * ((1 - std::cos((ii - 1) * M_PI / (Nbint))) / 2);
tabP2d(ii).SetY(U);
/* if (!HC2d.IsNull()) {
double L = GCPnts_AbscissaPoint::Length(HC2d->Curve2d(), Uinit, U);
tabP2d(ii).SetX(L);
}
else*/
tabP2d(ii).SetX(GCPnts_AbscissaPoint::Length(AC2d, Uinit, U));
}
acrlaw->Set(tabP2d);
}
NbPnt_i = (int)(LinCont->NbPoints());
NbPtInter = PntInter->Value(i).Length();
NbPtG1G1 = PntG1G1->Value(i).Length();
#ifdef OCCT_DEBUG
if (Affich > 1)
{
std::cout << "Courbe : " << i << std::endl;
std::cout << " NbPnt, NbPtInter, NbPtG1G1 :" << NbPnt_i << ", " << NbPtInter << ", "
<< NbPtG1G1 << std::endl;
}
#endif
for (int j = 1; j <= NbPnt_i; j++)
{
// Distribution of points in cosine following ACR 2D
// To avoid points of accumulation in 2D
// Inter=Uinit+(Uif)*((-cos(M_PI*((j-1)/(NbPnt_i-1)))+1)/2);
if (j == NbPnt_i)
{
Inter = Ufinal; // to avoid bug on Sun
}
else if (ACR)
{
CurLength = Length2d * (1 - std::cos((j - 1) * M_PI / (NbPnt_i - 1))) / 2;
Inter = acrlaw->Value(CurLength);
}
else
{
Inter = Uinit + (Ufinal - Uinit) * ((1 - std::cos((j - 1) * M_PI / (NbPnt_i - 1))) / 2);
}
myParCont->ChangeValue(i).Append(Inter); // add a point
if (NbPtInter != 0)
{
for (int l = 1; l <= NbPtInter; l += 2)
{
// check if the point Inter is in the interval
// PntInter[i] PntInter[i+1]
// in which case it is not necessary to store it (problem with duplicates)
if ((Inter > PntInter->Value(i).Value(l)) && (Inter < PntInter->Value(i).Value(l + 1)))
{
l = NbPtInter + 2;
// leave the loop without storing the point
}
else
{
if (l + 1 >= NbPtInter)
{
// one has parsed the entire table : the point
// does not belong to a common point interval
if (NbPtG1G1 != 0)
{
// if there exists an incompatible interval
for (int k = 1; k <= NbPtG1G1; k += 2)
{
if ((Inter > PntG1G1->Value(i).Value(k))
&& (Inter < PntG1G1->Value(i).Value(k + 1)))
{
k = NbPtG1G1 + 2; // to leave the loop
// Add points of constraint G0
gp_Pnt P3d, PP, Pdif;
gp_Pnt2d P2d;
AC2d.D0(Inter, P2d);
LinCont->D0(Inter, P3d);
mySurfInit->D0(P2d.Coord(1), P2d.Coord(2), PP);
Pdif.SetCoord(-PP.Coord(1) + P3d.Coord(1),
-PP.Coord(2) + P3d.Coord(2),
-PP.Coord(3) + P3d.Coord(3));
Plate_PinpointConstraint PC(P2d.XY(), Pdif.XYZ(), 0, 0);
myPlate.Load(PC);
}
else // the point does not belong to interval G1
{
if (k + 1 >= NbPtG1G1)
{
myPlateCont->ChangeValue(i).Append(Inter);
// add the point
}
}
}
}
else
{
myPlateCont->ChangeValue(i).Append(Inter);
// add the point
}
}
}
}
}
else
{
if (NbPtG1G1 != 0) // there exist an incompatible interval
{
for (int k = 1; k <= NbPtG1G1; k += 2)
{
if ((Inter > PntG1G1->Value(i).Value(k)) && (Inter < PntG1G1->Value(i).Value(k + 1)))
{
k = NbPtG1G1 + 2; // to leave the loop
// Add points of constraint G0
gp_Pnt P3d, PP, Pdif;
gp_Pnt2d P2d;
AC2d.D0(Inter, P2d);
LinCont->D0(Inter, P3d);
mySurfInit->D0(P2d.Coord(1), P2d.Coord(2), PP);
Pdif.SetCoord(-PP.Coord(1) + P3d.Coord(1),
-PP.Coord(2) + P3d.Coord(2),
-PP.Coord(3) + P3d.Coord(3));
Plate_PinpointConstraint PC(P2d.XY(), Pdif.XYZ(), 0, 0);
myPlate.Load(PC);
}
else // the point does not belong to interval G1
{
if (k + 1 >= NbPtG1G1)
{
myPlateCont->ChangeValue(i).Append(Inter);
// add the point
}
}
}
}
else
{
if (((!mySurfInitIsGive)
&& (Geom2dAdaptor_Curve(LinCont->Curve2dOnSurf()).GetType() != GeomAbs_Circle))
|| ((j > 1) && (j < NbPnt_i)))
{ // exclude extremities
myPlateCont->ChangeValue(i).Append(Inter); // add the point
}
}
}
}
}
}
//---------------------------------------------------------
// Function : CalculNbPtsInit
//---------------------------------------------------------
// Calculate the number of points by curve depending on the
// length for the first iteration
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::CalculNbPtsInit()
{
double LenT = 0;
int NTLinCont = myLinCont->Length();
int NTPoint = (int)(myNbPtsOnCur * NTLinCont);
int i;
for (i = 1; i <= NTLinCont; i++)
{
LenT += myLinCont->Value(i)->Length();
}
for (i = 1; i <= NTLinCont; i++)
{
int Cont = myLinCont->Value(i)->Order();
switch (Cont)
{
case 0: // Case G0 *1.2
myLinCont->ChangeValue(i)->SetNbPoints(
int(1.2 * NTPoint * (myLinCont->Value(i)->Length()) / LenT));
break;
case 1: // Case G1 *1
myLinCont->ChangeValue(i)->SetNbPoints(
int(NTPoint * (myLinCont->Value(i)->Length()) / LenT));
break;
case 2: // Case G2 *0.7
myLinCont->ChangeValue(i)->SetNbPoints(
int(0.7 * NTPoint * (myLinCont->Value(i)->Length()) / LenT));
break;
}
if (myLinCont->Value(i)->NbPoints() < 3)
{
myLinCont->ChangeValue(i)->SetNbPoints(3);
}
}
}
//---------------------------------------------------------
// Function : LoadCurve
//---------------------------------------------------------
// Starting from table myParCont load all the points noted in plate
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::LoadCurve(const int NbBoucle, const int OrderMax)
{
gp_Pnt P3d, Pdif, PP;
gp_Pnt2d P2d;
int NTLinCont = myLinCont->Length(), i, j;
int Tang, Nt;
for (i = 1; i <= NTLinCont; i++)
{
occ::handle<GeomPlate_CurveConstraint> CC = myLinCont->Value(i);
if (CC->Order() != -1)
{
Tang = std::min(CC->Order(), OrderMax);
Nt = myPlateCont->Value(i).Length();
if (Tang != -1)
{
for (j = 1; j <= Nt; j++)
{ // Loading of points G0 on boundaries
CC->D0(myPlateCont->Value(i).Value(j), P3d);
if (!CC->ProjectedCurve().IsNull())
{
P2d = CC->ProjectedCurve()->Value(myPlateCont->Value(i).Value(j));
}
else
{
if (!CC->Curve2dOnSurf().IsNull())
{
P2d = CC->Curve2dOnSurf()->Value(myPlateCont->Value(i).Value(j));
}
else
{
P2d = ProjectPoint(P3d);
}
}
mySurfInit->D0(P2d.Coord(1), P2d.Coord(2), PP);
Pdif.SetCoord(-PP.Coord(1) + P3d.Coord(1),
-PP.Coord(2) + P3d.Coord(2),
-PP.Coord(3) + P3d.Coord(3));
Plate_PinpointConstraint PC(P2d.XY(), Pdif.XYZ(), 0, 0);
myPlate.Load(PC);
// Loading of points G1
if (Tang == 1)
{ // ==1
gp_Vec V1, V2, V3, V4;
CC->D1(myPlateCont->Value(i).Value(j), PP, V1, V2);
mySurfInit->D1(P2d.Coord(1), P2d.Coord(2), PP, V3, V4);
Plate_D1 D1final(V1.XYZ(), V2.XYZ());
Plate_D1 D1init(V3.XYZ(), V4.XYZ());
if (!myFree)
{
Plate_GtoCConstraint GCC(P2d.XY(), D1init, D1final);
myPlate.Load(GCC);
}
else if (NbBoucle == 1)
{
Plate_FreeGtoCConstraint FreeGCC(P2d.XY(), D1init, D1final);
myPlate.Load(FreeGCC);
}
else
{
gp_Vec DU, DV, Normal, DerPlateU, DerPlateV;
Normal = V1 ^ V2;
// Normal.Normalize();
double norm = Normal.Magnitude();
if (norm > 1.e-12)
{
Normal /= norm;
}
DerPlateU = myPrevPlate.EvaluateDerivative(P2d.XY(), 1, 0);
DerPlateV = myPrevPlate.EvaluateDerivative(P2d.XY(), 0, 1);
DU.SetLinearForm(-(V3 + DerPlateU).Dot(Normal), Normal, DerPlateU);
DV.SetLinearForm(-(V4 + DerPlateV).Dot(Normal), Normal, DerPlateV);
Plate_PinpointConstraint PinU(P2d.XY(), DU.XYZ(), 1, 0);
Plate_PinpointConstraint PinV(P2d.XY(), DV.XYZ(), 0, 1);
myPlate.Load(PinU);
myPlate.Load(PinV);
}
}
// Loading of points G2
if (Tang == 2) // ==2
{
gp_Vec V1, V2, V3, V4, V5, V6, V7, V8, V9, V10;
CC->D2(myPlateCont->Value(i).Value(j), PP, V1, V2, V5, V6, V7);
mySurfInit->D2(P2d.Coord(1), P2d.Coord(2), PP, V3, V4, V8, V9, V10);
Plate_D1 D1final(V1.XYZ(), V2.XYZ());
Plate_D1 D1init(V3.XYZ(), V4.XYZ());
Plate_D2 D2final(V5.XYZ(), V6.XYZ(), V7.XYZ());
Plate_D2 D2init(V8.XYZ(), V9.XYZ(), V10.XYZ());
// if (! myFree)
// {
Plate_GtoCConstraint GCC(P2d.XY(), D1init, D1final, D2init, D2final);
myPlate.Load(GCC);
// }
// else // Good but too expansive
// {
// Plate_FreeGtoCConstraint FreeGCC( P2d.XY(),
// D1init, D1final, D2init, D2final );
// myPlate.Load( FreeGCC );
// }
}
}
}
}
}
}
//---------------------------------------------------------
// Function : LoadPoint
//---------------------------------------------------------
// void GeomPlate_BuildPlateSurface::LoadPoint(const int NbBoucle,
void GeomPlate_BuildPlateSurface::LoadPoint(const int, const int OrderMax)
{
gp_Pnt P3d, Pdif, PP;
gp_Pnt2d P2d;
int NTPntCont = myPntCont->Length();
int Tang, i;
// gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10;
// Loading of points of point constraints
for (i = 1; i <= NTPntCont; i++)
{
myPntCont->Value(i)->D0(P3d);
P2d = myPntCont->Value(i)->Pnt2dOnSurf();
mySurfInit->D0(P2d.Coord(1), P2d.Coord(2), PP);
Pdif.SetCoord(-PP.Coord(1) + P3d.Coord(1),
-PP.Coord(2) + P3d.Coord(2),
-PP.Coord(3) + P3d.Coord(3));
Plate_PinpointConstraint PC(P2d.XY(), Pdif.XYZ(), 0, 0);
myPlate.Load(PC);
Tang = std::min(myPntCont->Value(i)->Order(), OrderMax);
if (Tang == 1)
{ // ==1
gp_Vec V1, V2, V3, V4;
myPntCont->Value(i)->D1(PP, V1, V2);
mySurfInit->D1(P2d.Coord(1), P2d.Coord(2), PP, V3, V4);
Plate_D1 D1final(V1.XYZ(), V2.XYZ());
Plate_D1 D1init(V3.XYZ(), V4.XYZ());
if (!myFree)
{
Plate_GtoCConstraint GCC(P2d.XY(), D1init, D1final);
myPlate.Load(GCC);
}
else
{
Plate_FreeGtoCConstraint FreeGCC(P2d.XY(), D1init, D1final);
myPlate.Load(FreeGCC);
}
}
// Loading of points G2 GeomPlate_PlateG0Criterion
if (Tang == 2) // ==2
{
gp_Vec V1, V2, V3, V4, V5, V6, V7, V8, V9, V10;
myPntCont->Value(i)->D2(PP, V1, V2, V5, V6, V7);
// gp_Vec Tv2 = V1^V2;
mySurfInit->D2(P2d.Coord(1), P2d.Coord(2), PP, V3, V4, V8, V9, V10);
Plate_D1 D1final(V1.XYZ(), V2.XYZ());
Plate_D1 D1init(V3.XYZ(), V4.XYZ());
Plate_D2 D2final(V5.XYZ(), V6.XYZ(), V7.XYZ());
Plate_D2 D2init(V8.XYZ(), V9.XYZ(), V10.XYZ());
// if (! myFree)
// {
Plate_GtoCConstraint GCC(P2d.XY(), D1init, D1final, D2init, D2final);
myPlate.Load(GCC);
// }
// else // Good but too expansive
// {
// Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final, D2init//, D2final );
// myPlate.Load( FreeGCC );
// }
}
}
}
//---------------------------------------------------------
// Function : VerifSurface
//---------------------------------------------------------
bool GeomPlate_BuildPlateSurface::VerifSurface(const int NbBoucle)
{
//======================================================================
// Calculate errors
//======================================================================
int NTLinCont = myLinCont->Length();
bool Result = true;
// variable for error calculation
myG0Error = 0, myG1Error = 0, myG2Error = 0;
for (int i = 1; i <= NTLinCont; i++)
{
occ::handle<GeomPlate_CurveConstraint> LinCont;
LinCont = myLinCont->Value(i);
if (LinCont->Order() != -1)
{
int NbPts_i = myParCont->Value(i).Length();
if (NbPts_i < 3)
{
NbPts_i = 4;
}
occ::handle<NCollection_HArray1<double>> tdist =
new NCollection_HArray1<double>(1, NbPts_i - 1);
occ::handle<NCollection_HArray1<double>> tang =
new NCollection_HArray1<double>(1, NbPts_i - 1);
occ::handle<NCollection_HArray1<double>> tcourb =
new NCollection_HArray1<double>(1, NbPts_i - 1);
EcartContraintesMil(i, tdist, tang, tcourb);
double diffDistMax = 0, diffAngMax = 0;
// double SdiffDist=0, SdiffAng=0;
int NdiffDist = 0, NdiffAng = 0;
for (int j = 1; j < NbPts_i; j++)
{
if (tdist->Value(j) > myG0Error)
{
myG0Error = tdist->Value(j);
}
if (tang->Value(j) > myG1Error)
{
myG1Error = tang->Value(j);
}
if (tcourb->Value(j) > myG2Error)
{
myG2Error = tcourb->Value(j);
}
double U;
if (myParCont->Value(i).Length() > 3)
{
U = (myParCont->Value(i).Value(j) + myParCont->Value(i).Value(j + 1)) / 2;
}
else
{
U = LinCont->FirstParameter()
+ (LinCont->LastParameter() - LinCont->FirstParameter()) * (j - 1) / (NbPts_i - 2);
}
double diffDist = tdist->Value(j) - LinCont->G0Criterion(U), diffAng;
if (LinCont->Order() > 0)
{
diffAng = tang->Value(j) - LinCont->G1Criterion(U);
}
else
{
diffAng = 0;
}
// find the maximum variation of error and calculate the average
if (diffDist > 0)
{
diffDist = diffDist / LinCont->G0Criterion(U);
if (diffDist > diffDistMax)
{
diffDistMax = diffDist;
}
// SdiffDist+=diffDist;
NdiffDist++;
}
else if ((diffAng > 0) && (LinCont->Order() == 1))
{
diffAng = diffAng / myLinCont->Value(i)->G1Criterion(U);
if (diffAng > diffAngMax)
{
diffAngMax = diffAng;
}
// SdiffAng+=diffAng;
NdiffAng++;
}
}
if (NdiffDist > 0)
{ // at least one point is not acceptable in G0
double Coef;
if (LinCont->Order() == 0)
{
Coef = 0.6 * std::log(diffDistMax + 7.4);
// 7.4 corresponds to the calculation of min. coefficient = 1.2 is e^1.2/0.6
}
else
{
Coef = std::log(diffDistMax + 3.3);
}
// 3.3 corresponds to calculation of min. coefficient = 1.2 donc e^1.2
if (Coef > 3)
{
Coef = 3;
}
// experimentally after the coefficient becomes bad for L cases
if ((NbBoucle > 1) && (diffDistMax > 2))
{
Coef = 1.6;
}
if (LinCont->NbPoints() >= std::floor(LinCont->NbPoints() * Coef))
{
Coef = 2; // to provide increase of the number of points
}
LinCont->SetNbPoints(int(LinCont->NbPoints() * Coef));
Result = false;
}
else if (NdiffAng > 0) // at least 1 point is not acceptable in G1
{
double Coef = 1.5;
if ((LinCont->NbPoints() + 1) >= std::floor(LinCont->NbPoints() * Coef))
{
Coef = 2;
}
LinCont->SetNbPoints(int(LinCont->NbPoints() * Coef));
Result = false;
}
}
}
if (!Result)
{
if (myFree && NbBoucle == 1)
{
myPrevPlate = myPlate;
}
myPlate.Init();
}
return Result;
}
//---------------------------------------------------------
// Function : VerifPoint
//---------------------------------------------------------
void GeomPlate_BuildPlateSurface::VerifPoints(double& Dist, double& Ang, double& Curv) const
{
int NTPntCont = myPntCont->Length();
gp_Pnt Pi, Pf;
gp_Pnt2d P2d;
gp_Vec v1i, v1f, v2i, v2f, v3i, v3f;
Ang = 0;
Dist = 0, Curv = 0;
occ::handle<GeomPlate_PointConstraint> PntCont;
for (int i = 1; i <= NTPntCont; i++)
{
PntCont = myPntCont->Value(i);
switch (PntCont->Order())
{
case 0:
P2d = PntCont->Pnt2dOnSurf();
PntCont->D0(Pi);
myGeomPlateSurface->D0(P2d.Coord(1), P2d.Coord(2), Pf);
Dist = Pf.Distance(Pi);
break;
case 1:
PntCont->D1(Pi, v1i, v2i);
P2d = PntCont->Pnt2dOnSurf();
myGeomPlateSurface->D1(P2d.Coord(1), P2d.Coord(2), Pf, v1f, v2f);
Dist = Pf.Distance(Pi);
v3i = v1i ^ v2i;
v3f = v1f ^ v2f;
Ang = v3f.Angle(v3i);
if (Ang > (M_PI / 2))
{
Ang = M_PI - Ang;
}
break;
case 2:
occ::handle<Geom_Surface> Splate(myGeomPlateSurface);
LocalAnalysis_SurfaceContinuity CG2;
P2d = PntCont->Pnt2dOnSurf();
GeomLProp_SLProps Prop(Splate, P2d.Coord(1), P2d.Coord(2), 2, 0.001);
CG2.ComputeAnalysis(Prop, PntCont->LPropSurf(), GeomAbs_G2);
Dist = CG2.C0Value();
Ang = CG2.G1Angle();
Curv = CG2.G2CurvatureGap();
break;
}
}
}
double GeomPlate_BuildPlateSurface::ComputeAnisotropie() const
{
if (myAnisotropie)
{
// Temporary
return 1.0;
}
else
{
return 1.0;
}
}
bool GeomPlate_BuildPlateSurface::IsOrderG1() const
{
bool result = true;
for (int i = 1; i <= myLinCont->Length(); i++)
{
if (myLinCont->Value(i)->Order() < 1)
{
result = false;
break;
}
}
return result;
}