OpenCascade/occt
1
0
Fork 0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-04-03 17:56:21 +03:00
Code Packages Releases Activity
occt/src/AppDef/AppDef_Variational.cxx
tiv 0423218095 0030895: Coding Rules - specify std namespace explicitly for std::cout and streams 
"endl" manipulator for Message_Messenger is renamed to "Message_EndLine".

The following entities from std namespace are now used
with std:: explicitly specified (from Standard_Stream.hxx):
std::istream,std::ostream,std::ofstream,std::ifstream,std::fstream,
std::filebuf,std::streambuf,std::streampos,std::ios,std::cout,std::cerr,
std::cin,std::endl,std::ends,std::flush,std::setw,std::setprecision,
std::hex,std::dec.
2019-08-16 12:16:38 +03:00

3050 lines
97 KiB
C++
Raw Blame History

// 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.
// Jeannine PANCIATICI le 06/06/96
// Igor FEOKTISTOV 14/12/98 - correction of Approximate() and Init().
// Approximation d une MultiLine de points decrite par le tool MLineTool.
// avec criteres variationnels
#include <AppDef_MultiLine.hxx>
#include <AppDef_SmoothCriterion.hxx>
#include <AppDef_Variational.hxx>
#include <AppParCurves_MultiBSpCurve.hxx>
#include <FEmTool_Assembly.hxx>
#include <FEmTool_Curve.hxx>
#include <gp_VectorWithNullMagnitude.hxx>
#include <math_Matrix.hxx>
#include <PLib_Base.hxx>
#include <Standard_ConstructionError.hxx>
#include <Standard_DimensionError.hxx>
#include <Standard_DomainError.hxx>
#include <Standard_OutOfRange.hxx>
#include <StdFail_NotDone.hxx>
#define No_Standard_RangeError
#define No_Standard_OutOfRange
#define No_Standard_DimensionError
#define No_Standard_ConstructionError
#include <Standard_Stream.hxx>
#include <AppParCurves.hxx>
#include <AppParCurves_Constraint.hxx>
#include <AppParCurves_HArray1OfConstraintCouple.hxx>
#include <AppParCurves_Array1OfMultiPoint.hxx>
#include <AppParCurves_MultiPoint.hxx>
#include <AppParCurves_MultiBSpCurve.hxx>
#include <AppDef_LinearCriteria.hxx>
#include <Convert_CompPolynomialToPoles.hxx>
#include <gp_Pnt.hxx>
#include <gp_Pnt2d.hxx>
#include <gp_Vec.hxx>
#include <gp_Vec2d.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TColgp_Array1OfPnt2d.hxx>
#include <TColgp_Array1OfVec.hxx>
#include <TColgp_Array1OfVec2d.hxx>
#include <TColStd_Array1OfInteger.hxx>
#include <TColStd_HArray1OfInteger.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <TColStd_HArray1OfReal.hxx>
#include <TColStd_HArray2OfReal.hxx>
#include <StdFail_NotDone.hxx>
#include <Standard_SStream.hxx>
#include <Standard_NoSuchObject.hxx>
#include <Precision.hxx>
#include <AppDef_MyLineTool.hxx>
#include <TColStd_HArray2OfInteger.hxx>
#include <TColStd_Array2OfInteger.hxx>
#include <TColStd_Array2OfReal.hxx>
#include <FEmTool_Assembly.hxx>
#include <FEmTool_AssemblyTable.hxx>
#include <FEmTool_Curve.hxx>
#include <math_Matrix.hxx>
#include <math_Vector.hxx>
#include <PLib_Base.hxx>
#include <PLib_JacobiPolynomial.hxx>
#include <PLib_HermitJacobi.hxx>
#include <FEmTool_HAssemblyTable.hxx>
// Add this line:
#include <algorithm>
#if defined(_MSC_VER)
# include <stdio.h>
# include <stdarg.h>
#endif /* _MSC_VER */
//
//=======================================================================
//function : AppDef_Variational
//purpose : Initialization of the fields.
//=======================================================================
//
AppDef_Variational::AppDef_Variational(const AppDef_MultiLine& SSP,
const Standard_Integer FirstPoint,
const Standard_Integer LastPoint,
const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
const Standard_Integer MaxDegree,
const Standard_Integer MaxSegment,
const GeomAbs_Shape Continuity,
const Standard_Boolean WithMinMax,
const Standard_Boolean WithCutting,
const Standard_Real Tolerance,
const Standard_Integer NbIterations):
mySSP(SSP),
myFirstPoint(FirstPoint),
myLastPoint(LastPoint),
myConstraints(TheConstraints),
myMaxDegree(MaxDegree),
myMaxSegment(MaxSegment),
myNbIterations(NbIterations),
myTolerance(Tolerance),
myContinuity(Continuity),
myWithMinMax(WithMinMax),
myWithCutting(WithCutting)
{
// Verifications:
if (myMaxDegree < 1) throw Standard_DomainError();
myMaxDegree = Min (30, myMaxDegree);
//
if (myMaxSegment < 1) throw Standard_DomainError();
//
if (myWithMinMax != 0 && myWithMinMax !=1 ) throw Standard_DomainError();
if (myWithCutting != 0 && myWithCutting !=1 ) throw Standard_DomainError();
//
myIsOverConstr = Standard_False;
myIsCreated = Standard_False;
myIsDone = Standard_False;
switch (myContinuity) {
case GeomAbs_C0:
myNivCont=0;
break ;
case GeomAbs_C1:
myNivCont=1;
break ;
case GeomAbs_C2:
myNivCont=2;
break ;
default:
throw Standard_ConstructionError();
}
//
myNbP2d = AppDef_MyLineTool::NbP2d(SSP);
myNbP3d = AppDef_MyLineTool::NbP3d(SSP);
myDimension = 2 * myNbP2d + 3* myNbP3d ;
//
myPercent[0]=0.4;
myPercent[1]=0.2;
myPercent[2]=0.4;
myKnots= new TColStd_HArray1OfReal(1,2);
myKnots->SetValue(1,0.);
myKnots->SetValue(2,1.);
// Declaration
//
mySmoothCriterion = new AppDef_LinearCriteria(mySSP, myFirstPoint, myLastPoint);
myParameters = new TColStd_HArray1OfReal(myFirstPoint, myLastPoint);
myNbPoints=myLastPoint-myFirstPoint+1;
if (myNbPoints <= 0) throw Standard_ConstructionError();
//
myTabPoints= new TColStd_HArray1OfReal(1,myDimension*myNbPoints);
//
// Table of Points initialization
//
Standard_Integer ipoint,jp2d,jp3d,index;
TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
gp_Pnt2d P2d;
gp_Pnt P3d;
index=1;
for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++)
{
if(myNbP2d !=0 && myNbP3d ==0 ) {
AppDef_MyLineTool::Value(mySSP,ipoint,TabP2d);
for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++)
{ P2d = TabP2d.Value(jp2d);
myTabPoints->SetValue(index++,P2d.X());
myTabPoints->SetValue(index++,P2d.Y());
}
}
if(myNbP3d !=0 && myNbP2d == 0 ) {
AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d);
for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
{ P3d=TabP3d.Value(jp3d);
myTabPoints->SetValue(index++,P3d.X());
myTabPoints->SetValue(index++,P3d.Y());
myTabPoints->SetValue(index++,P3d.Z());
}
}
if(myNbP3d !=0 && myNbP2d != 0 ) {
AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d,TabP2d);
for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
{ P3d=TabP3d.Value(jp3d);
myTabPoints->SetValue(index++,P3d.X());
myTabPoints->SetValue(index++,P3d.Y());
myTabPoints->SetValue(index++,P3d.Z());
}
for ( jp2d = 1 ; jp2d <= myNbP2d ; jp2d++)
{ P2d=TabP2d.Value(jp2d);
myTabPoints->SetValue(index++,P2d.X());
myTabPoints->SetValue(index++,P2d.Y());
}
}
}
Init();
}
//
//=======================================================================
//function : Init
//purpose : Initializes the tables of constraints
// and verifies if the problem is not over-constrained.
// This method is used in the Create and the method SetConstraint.
//=======================================================================
//
void AppDef_Variational::Init()
{
Standard_Integer ipoint,jp2d,jp3d,index,jndex;
Standard_Integer CurMultyPoint;
TColgp_Array1OfVec TabV3d(1, Max(1,myNbP3d));
TColgp_Array1OfVec2d TabV2d(1, Max(1,myNbP2d));
TColgp_Array1OfVec TabV3dcurv(1, Max(1,myNbP3d));
TColgp_Array1OfVec2d TabV2dcurv(1, Max(1,myNbP2d));
gp_Vec Vt3d, Vc3d;
gp_Vec2d Vt2d, Vc2d;
myNbConstraints=myConstraints->Length();
if (myNbConstraints < 0) throw Standard_ConstructionError();
myTypConstraints = new TColStd_HArray1OfInteger(1,Max(1,2*myNbConstraints));
myTabConstraints = new TColStd_HArray1OfReal(1,Max(1,2*myDimension*myNbConstraints));
myTtheta = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints));
myTfthet = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints));
//
// Table of types initialization
Standard_Integer iconstr;
index=1;
jndex=1;
CurMultyPoint = 1;
myNbPassPoints=0;
myNbTangPoints=0;
myNbCurvPoints=0;
AppParCurves_Constraint valcontr;
for ( iconstr = myConstraints->Lower() ; iconstr <= myConstraints->Upper() ; iconstr++)
{
ipoint=(myConstraints->Value(iconstr)).Index();
valcontr=(myConstraints->Value(iconstr)).Constraint();
switch (valcontr) {
case AppParCurves_NoConstraint:
CurMultyPoint -= myNbP3d * 6 + myNbP2d * 2;
break ;
case AppParCurves_PassPoint:
myTypConstraints->SetValue(index++,ipoint);
myTypConstraints->SetValue(index++,0);
myNbPassPoints++;
if(myNbP2d !=0 ) jndex=jndex+4*myNbP2d;
if(myNbP3d !=0 ) jndex=jndex+6*myNbP3d;
break ;
case AppParCurves_TangencyPoint:
myTypConstraints->SetValue(index++,ipoint);
myTypConstraints->SetValue(index++,1);
myNbTangPoints++;
if(myNbP2d !=0 && myNbP3d == 0 )
{
if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV2d) == Standard_False)
throw Standard_ConstructionError();
for (jp2d=1;jp2d<=myNbP2d;jp2d++)
{
Vt2d=TabV2d.Value(jp2d);
Vt2d.Normalize();
myTabConstraints->SetValue(jndex++,Vt2d.X());
myTabConstraints->SetValue(jndex++,Vt2d.Y());
jndex=jndex+2;
InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4);
}
}
if(myNbP3d !=0 && myNbP2d == 0)
{
if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d) == Standard_False)
throw Standard_ConstructionError();
for (jp3d=1;jp3d<=myNbP3d;jp3d++)
{
Vt3d=TabV3d.Value(jp3d);
Vt3d.Normalize();
myTabConstraints->SetValue(jndex++,Vt3d.X());
myTabConstraints->SetValue(jndex++,Vt3d.Y());
myTabConstraints->SetValue(jndex++,Vt3d.Z());
jndex=jndex+3;
InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
}
}
if(myNbP3d !=0 && myNbP2d != 0)
{
if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False)
throw Standard_ConstructionError();
for (jp3d=1;jp3d<=myNbP3d;jp3d++)
{
Vt3d=TabV3d.Value(jp3d);
Vt3d.Normalize();
myTabConstraints->SetValue(jndex++,Vt3d.X());
myTabConstraints->SetValue(jndex++,Vt3d.Y());
myTabConstraints->SetValue(jndex++,Vt3d.Z());
jndex=jndex+3;
InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
}
for (jp2d=1;jp2d<=myNbP2d;jp2d++)
{
Vt2d=TabV2d.Value(jp2d);
Vt2d.Normalize();
myTabConstraints->SetValue(jndex++,Vt2d.X());
myTabConstraints->SetValue(jndex++,Vt2d.Y());
jndex=jndex+2;
InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4);
}
}
break ;
case AppParCurves_CurvaturePoint:
myTypConstraints->SetValue(index++,ipoint);
myTypConstraints->SetValue(index++,2);
myNbCurvPoints++;
if(myNbP2d !=0 && myNbP3d == 0)
{
if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV2d) == Standard_False )
throw Standard_ConstructionError();
if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV2dcurv) == Standard_False)
throw Standard_ConstructionError();
for (jp2d=1;jp2d<=myNbP2d;jp2d++)
{
Vt2d=TabV2d.Value(jp2d);
Vt2d.Normalize();
Vc2d=TabV2dcurv.Value(jp2d);
if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular())
throw Standard_ConstructionError();
myTabConstraints->SetValue(jndex++,Vt2d.X());
myTabConstraints->SetValue(jndex++,Vt2d.Y());
myTabConstraints->SetValue(jndex++,Vc2d.X());
myTabConstraints->SetValue(jndex++,Vc2d.Y());
InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4);
}
}
if(myNbP3d !=0 && myNbP2d == 0 )
{
if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d) == Standard_False )
throw Standard_ConstructionError();
if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV3dcurv) == Standard_False)
throw Standard_ConstructionError();
for (jp3d=1;jp3d<=myNbP3d;jp3d++)
{
Vt3d=TabV3d.Value(jp3d);
Vt3d.Normalize();
Vc3d=TabV3dcurv.Value(jp3d);
if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False)
throw Standard_ConstructionError();
myTabConstraints->SetValue(jndex++,Vt3d.X());
myTabConstraints->SetValue(jndex++,Vt3d.Y());
myTabConstraints->SetValue(jndex++,Vt3d.Z());
myTabConstraints->SetValue(jndex++,Vc3d.X());
myTabConstraints->SetValue(jndex++,Vc3d.Y());
myTabConstraints->SetValue(jndex++,Vc3d.Z());
InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
}
}
if(myNbP3d !=0 && myNbP2d != 0 )
{
if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False )
throw Standard_ConstructionError();
if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV3dcurv,TabV2dcurv) == Standard_False)
throw Standard_ConstructionError();
for (jp3d=1;jp3d<=myNbP3d;jp3d++)
{
Vt3d=TabV3d.Value(jp3d);
Vt3d.Normalize();
Vc3d=TabV3dcurv.Value(jp3d);
if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False)
throw Standard_ConstructionError();
myTabConstraints->SetValue(jndex++,Vt3d.X());
myTabConstraints->SetValue(jndex++,Vt3d.Y());
myTabConstraints->SetValue(jndex++,Vt3d.Z());
myTabConstraints->SetValue(jndex++,Vc3d.X());
myTabConstraints->SetValue(jndex++,Vc3d.Y());
myTabConstraints->SetValue(jndex++,Vc3d.Z());
InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
}
for (jp2d=1;jp2d<=myNbP2d;jp2d++)
{
Vt2d=TabV2d.Value(jp2d);
Vt2d.Normalize();
Vc2d=TabV2dcurv.Value(jp2d);
if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular())
throw Standard_ConstructionError();
myTabConstraints->SetValue(jndex++,Vt2d.X());
myTabConstraints->SetValue(jndex++,Vt2d.Y());
myTabConstraints->SetValue(jndex++,Vc2d.X());
myTabConstraints->SetValue(jndex++,Vc2d.Y());
InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4);
}
}
break ;
default:
throw Standard_ConstructionError();
}
CurMultyPoint += myNbP3d * 6 + myNbP2d * 2;
}
// OverConstraint Detection
Standard_Integer MaxSeg;
if(myWithCutting == Standard_True) MaxSeg = myMaxSegment ;
else MaxSeg = 1;
if (((myMaxDegree-myNivCont)*MaxSeg-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
{
myIsOverConstr =Standard_True;
myIsCreated = Standard_False;
}
else {
InitSmoothCriterion();
myIsCreated = Standard_True;
}
}
//
//=======================================================================
//function : Approximate
//purpose : Makes the approximation with the current fields.
//=======================================================================
//
void AppDef_Variational::Approximate()
{
if (myIsCreated == Standard_False ) throw StdFail_NotDone();
Standard_Real WQuadratic, WQuality;
TColStd_Array1OfReal Ecarts(myFirstPoint, myLastPoint);
mySmoothCriterion->GetWeight(WQuadratic, WQuality);
Handle(FEmTool_Curve) TheCurve;
mySmoothCriterion->GetCurve(TheCurve);
//---------------------------------------------------------------------
TheMotor(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts);
if(myWithMinMax && myTolerance < myMaxError)
Adjusting(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts);
//---------------------------------------------------------------------
Standard_Integer jp2d,jp3d,index,ipole,
NbElem = TheCurve->NbElements();
TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
Standard_Real debfin[2] = {-1., 1};
gp_Pnt2d P2d;
gp_Pnt P3d;
index=0;
{
Handle(TColStd_HArray2OfReal) PolynomialIntervalsPtr =
new TColStd_HArray2OfReal(1,NbElem,1,2) ;
Handle(TColStd_HArray1OfInteger) NbCoeffPtr =
new TColStd_HArray1OfInteger(1, myMaxSegment);
Standard_Integer size = myMaxSegment * (myMaxDegree + 1) * myDimension ;
Handle(TColStd_HArray1OfReal) CoeffPtr = new TColStd_HArray1OfReal(1,size);
CoeffPtr->Init(0.);
Handle(TColStd_HArray1OfReal) IntervallesPtr =
new TColStd_HArray1OfReal(1, NbElem + 1);
IntervallesPtr->ChangeArray1() = TheCurve->Knots();
TheCurve->GetPolynom(CoeffPtr->ChangeArray1());
Standard_Integer ii;
for(ii = 1; ii <= NbElem; ii++)
NbCoeffPtr->SetValue(ii, TheCurve->Degree(ii)+1);
for (ii = PolynomialIntervalsPtr->LowerRow() ;
ii <= PolynomialIntervalsPtr->UpperRow() ;ii++)
{
PolynomialIntervalsPtr->SetValue(ii,1,debfin[0]) ;
PolynomialIntervalsPtr->SetValue(ii,2,debfin[1]) ;
}
/*
printf("\n =========== Parameters for converting\n");
printf("nb_courbes : %d \n", NbElem);
printf("tab_option[4] : %d \n", myNivCont);
printf("myDimension : %d \n", myDimension);
printf("myMaxDegree : %d \n", myMaxDegree);
printf("\n NbcoeffPtr :\n");
for(ii = 1; ii <= NbElem; ii++) printf("NbCoeffPtr(%d) = %d \n", ii, NbCoeffPtr->Value(ii));
size = NbElem*(myMaxDegree + 1) * myDimension;
printf("\n CoeffPtr :\n");
for(ii = 1; ii <= size; ii++) printf("CoeffPtr(%d) = %.8e \n", ii, CoeffPtr->Value(ii));
printf("\n PolinimialIntervalsPtr :\n");
for (ii = PolynomialIntervalsPtr->LowerRow() ;
ii <= PolynomialIntervalsPtr->UpperRow() ;ii++)
{
printf(" %d : [%f, %f] \n", ii, PolynomialIntervalsPtr->Value(ii,1),
PolynomialIntervalsPtr->Value(ii,2)) ;
}
printf("\n IntervallesPtr :\n");
for (ii = IntervallesPtr->Lower();
ii <= IntervallesPtr->Upper() - 1; ii++)
{
printf(" %d : [%f, %f] \n", ii, IntervallesPtr->Value(ii),
IntervallesPtr->Value(ii+1)) ;
}
*/
Convert_CompPolynomialToPoles AConverter(NbElem,
myNivCont,
myDimension,
myMaxDegree,
NbCoeffPtr,
CoeffPtr,
PolynomialIntervalsPtr,
IntervallesPtr) ;
if (AConverter.IsDone())
{
Handle(TColStd_HArray2OfReal) PolesPtr ;
Handle(TColStd_HArray1OfInteger) Mults;
Standard_Integer NbPoles=AConverter.NbPoles();
// Standard_Integer Deg=AConverter.Degree();
AppParCurves_Array1OfMultiPoint TabMU(1,NbPoles);
AConverter.Poles(PolesPtr) ;
AConverter.Knots(myKnots) ;
AConverter.Multiplicities(Mults) ;
for (ipole=PolesPtr->LowerRow();ipole<=PolesPtr->UpperRow();ipole++)
{
index=PolesPtr->LowerCol();
/* if(myNbP2d !=0 )
{
for (jp2d=1;jp2d<=myNbP2d;jp2d++)
{
P2d.SetX(PolesPtr->Value(ipole,index++));
P2d.SetY(PolesPtr->Value(ipole,index++));
TabP2d.SetValue(jp2d,P2d);
}
}*/
if(myNbP3d !=0 )
{
for (jp3d=1;jp3d<=myNbP3d;jp3d++)
{
// std::cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
P3d.SetX(PolesPtr->Value(ipole,index++));
// std::cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
P3d.SetY(PolesPtr->Value(ipole,index++));
// std::cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
P3d.SetZ(PolesPtr->Value(ipole,index++));
TabP3d.SetValue(jp3d,P3d);
}
}
if(myNbP2d !=0 )
{
for (jp2d=1;jp2d<=myNbP2d;jp2d++)
{
P2d.SetX(PolesPtr->Value(ipole,index++));
P2d.SetY(PolesPtr->Value(ipole,index++));
TabP2d.SetValue(jp2d,P2d);
}
}
if(myNbP2d !=0 && myNbP3d !=0)
{
AppParCurves_MultiPoint aMultiPoint(TabP3d,TabP2d);
TabMU.SetValue(ipole,aMultiPoint);
}
else if (myNbP2d !=0)
{
AppParCurves_MultiPoint aMultiPoint(TabP2d);
TabMU.SetValue(ipole,aMultiPoint);
}
else
{
AppParCurves_MultiPoint aMultiPoint(TabP3d);
TabMU.SetValue(ipole,aMultiPoint);
}
}
AppParCurves_MultiBSpCurve aCurve(TabMU,myKnots->Array1(),Mults->Array1());
myMBSpCurve=aCurve;
myIsDone = Standard_True;
}
}
}
//
//=======================================================================
//function : IsCreated
//purpose : returns True if the creation is done
//=======================================================================
//
Standard_Boolean AppDef_Variational::IsCreated() const
{
return myIsCreated;
}
//
//=======================================================================
//function : IsDone
//purpose : returns True if the approximation is ok
//=======================================================================
//
Standard_Boolean AppDef_Variational::IsDone() const
{
return myIsDone;
}
//
//=======================================================================
//function : IsOverConstrained
//purpose : returns True if the problem is overconstrained
// in this case, approximation cannot be done.
//=======================================================================
//
Standard_Boolean AppDef_Variational::IsOverConstrained() const
{
return myIsOverConstr;
}
//
//=======================================================================
//function : Value
//purpose : returns all the BSpline curves approximating the
// MultiLine SSP after minimization of the parameter.
//=======================================================================
//
AppParCurves_MultiBSpCurve AppDef_Variational::Value() const
{
if (myIsDone == Standard_False) throw StdFail_NotDone();
return myMBSpCurve;
}
//
//=======================================================================
//function : MaxError
//purpose : returns the maximum of the distances between
// the points of the multiline and the approximation
// curves.
//=======================================================================
//
Standard_Real AppDef_Variational::MaxError() const
{
if (myIsDone == Standard_False) throw StdFail_NotDone();
return myMaxError;
}
//
//=======================================================================
//function : MaxErrorIndex
//purpose : returns the index of the MultiPoint of ErrorMax
//=======================================================================
//
Standard_Integer AppDef_Variational::MaxErrorIndex() const
{
if (myIsDone == Standard_False) throw StdFail_NotDone();
return myMaxErrorIndex;
}
//
//=======================================================================
//function : QuadraticError
//purpose : returns the quadratic average of the distances between
// the points of the multiline and the approximation
// curves.
//=======================================================================
//
Standard_Real AppDef_Variational::QuadraticError() const
{
if (myIsDone == Standard_False) throw StdFail_NotDone();
return myCriterium[0];
}
//
//=======================================================================
//function : Distance
//purpose : returns the distances between the points of the
// multiline and the approximation curves.
//=======================================================================
//
void AppDef_Variational::Distance(math_Matrix& mat)
{
if (myIsDone == Standard_False) throw StdFail_NotDone();
Standard_Integer ipoint,jp2d,jp3d,index;
TColgp_Array1OfPnt TabP3d(1,Max(1,myNbP3d));
TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
Standard_Integer j0 = mat.LowerCol() - myFirstPoint;
gp_Pnt2d P2d;
gp_Pnt P3d;
gp_Pnt Pt3d;
gp_Pnt2d Pt2d;
for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++)
{
index=1;
if(myNbP3d !=0 ) {
AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d);
for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
{ P3d=TabP3d.Value(jp3d);
myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt3d);
mat(index++, j0 + ipoint)=P3d.Distance(Pt3d);
}
}
if(myNbP2d !=0 ) {
if(myNbP3d == 0 ) AppDef_MyLineTool::Value(mySSP,ipoint,TabP2d);
else AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d,TabP2d);
for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++)
{ P2d = TabP2d.Value(jp2d);
myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt2d);
mat(index++, j0 + ipoint)=P2d.Distance(Pt2d);
}
}
}
}
//
//=======================================================================
//function : AverageError
//purpose : returns the average error between
// the MultiLine and the approximation.
//=======================================================================
//
Standard_Real AppDef_Variational::AverageError() const
{
if (myIsDone == Standard_False) throw StdFail_NotDone();
return myAverageError;
}
//
//=======================================================================
//function : Parameters
//purpose : returns the parameters uses to the approximations
//=======================================================================
//
const Handle(TColStd_HArray1OfReal)& AppDef_Variational::Parameters() const
{
if (myIsDone == Standard_False) throw StdFail_NotDone();
return myParameters;
}
//
//=======================================================================
//function : Knots
//purpose : returns the knots uses to the approximations
//=======================================================================
//
const Handle(TColStd_HArray1OfReal)& AppDef_Variational::Knots() const
{
if (myIsDone == Standard_False) throw StdFail_NotDone();
return myKnots;
}
//
//=======================================================================
//function : Criterium
//purpose : returns the values of the quality criterium.
//=======================================================================
//
void AppDef_Variational::Criterium(Standard_Real& VFirstOrder, Standard_Real& VSecondOrder, Standard_Real& VThirdOrder) const
{
if (myIsDone == Standard_False) throw StdFail_NotDone();
VFirstOrder=myCriterium[1] ;
VSecondOrder=myCriterium[2];
VThirdOrder=myCriterium[3];
}
//
//=======================================================================
//function : CriteriumWeight
//purpose : returns the Weights (as percent) associed to the criterium used in
// the optimization.
//=======================================================================
//
void AppDef_Variational::CriteriumWeight(Standard_Real& Percent1, Standard_Real& Percent2, Standard_Real& Percent3) const
{
Percent1 = myPercent[0];
Percent2 = myPercent[1];
Percent3 = myPercent[2] ;
}
//
//=======================================================================
//function : MaxDegree
//purpose : returns the Maximum Degree used in the approximation
//=======================================================================
//
Standard_Integer AppDef_Variational::MaxDegree() const
{
return myMaxDegree;
}
//
//=======================================================================
//function : MaxSegment
//purpose : returns the Maximum of segment used in the approximation
//=======================================================================
//
Standard_Integer AppDef_Variational::MaxSegment() const
{
return myMaxSegment;
}
//
//=======================================================================
//function : Continuity
//purpose : returns the Continuity used in the approximation
//=======================================================================
//
GeomAbs_Shape AppDef_Variational::Continuity() const
{
return myContinuity;
}
//
//=======================================================================
//function : WithMinMax
//purpose : returns if the approximation search to minimize the
// maximum Error or not.
//=======================================================================
//
Standard_Boolean AppDef_Variational::WithMinMax() const
{
return myWithMinMax;
}
//
//=======================================================================
//function : WithCutting
//purpose : returns if the approximation can insert new Knots or not.
//=======================================================================
//
Standard_Boolean AppDef_Variational::WithCutting() const
{
return myWithCutting;
}
//
//=======================================================================
//function : Tolerance
//purpose : returns the tolerance used in the approximation.
//=======================================================================
//
Standard_Real AppDef_Variational::Tolerance() const
{
return myTolerance;
}
//
//=======================================================================
//function : NbIterations
//purpose : returns the number of iterations used in the approximation.
//=======================================================================
//
Standard_Integer AppDef_Variational::NbIterations() const
{
return myNbIterations;
}
//
//=======================================================================
//function : Dump
//purpose : Prints on the stream o information on the current state
// of the object.
//=======================================================================
//
void AppDef_Variational::Dump(Standard_OStream& o) const
{
o << " \nVariational Smoothing " << std::endl;
o << " Number of multipoints " << myNbPoints << std::endl;
o << " Number of 2d par multipoint " << myNbP2d << std::endl;
o << " Nombre of 3d par multipoint " << myNbP3d << std::endl;
o << " Number of PassagePoint " << myNbPassPoints << std::endl;
o << " Number of TangencyPoints " << myNbTangPoints << std::endl;
o << " Number of CurvaturePoints " << myNbCurvPoints << std::endl;
o << " \nTolerance " << o.setf(std::ios::scientific) << std::setprecision(3) << std::setw(9) << myTolerance;
if ( WithMinMax()) { o << " as Max Error." << std::endl;}
else { o << " as size Error." << std::endl;}
o << "CriteriumWeights : " << myPercent[0] << " , "
<< myPercent[1] << " , " << myPercent[2] << std::endl;
if (myIsDone ) {
o << " MaxError " << std::setprecision(3) << std::setw(9) << myMaxError << std::endl;
o << " Index of MaxError " << myMaxErrorIndex << std::endl;
o << " Average Error " << std::setprecision(3) << std::setw(9) << myAverageError << std::endl;
o << " Quadratic Error " << std::setprecision(3) << std::setw(9) << myCriterium[0] << std::endl;
o << " Tension Criterium " << std::setprecision(3) << std::setw(9) << myCriterium[1] << std::endl;
o << " Flexion Criterium " << std::setprecision(3) << std::setw(9) << myCriterium[2] << std::endl;
o << " Jerk Criterium " << std::setprecision(3) << std::setw(9) << myCriterium[3] << std::endl;
o << " NbSegments " << myKnots->Length()-1 << std::endl;
}
else
{ if (myIsOverConstr) o << "The probleme is overconstraint " << std::endl;
else o << " Erreur dans l''approximation" << std::endl;
}
}
//
//=======================================================================
//function : SetConstraints
//purpose : Define the constraints to approximate
// If this value is incompatible with the others fields
// this method modify nothing and returns false
//=======================================================================
//
Standard_Boolean AppDef_Variational::SetConstraints(const Handle(AppParCurves_HArray1OfConstraintCouple)& aConstraint)
{
myConstraints=aConstraint;
Init();
if (myIsOverConstr ) return Standard_False;
else return Standard_True;
}
//
//=======================================================================
//function : SetParameters
//purpose : Defines the parameters used by the approximations.
//=======================================================================
//
void AppDef_Variational::SetParameters(const Handle(TColStd_HArray1OfReal)& param)
{
myParameters->ChangeArray1() = param->Array1();
}
//
//=======================================================================
//function : SetKnots
//purpose : Defines the knots used by the approximations
// -- If this value is incompatible with the others fields
// -- this method modify nothing and returns false
//=======================================================================
//
Standard_Boolean AppDef_Variational::SetKnots(const Handle(TColStd_HArray1OfReal)& knots)
{
myKnots->ChangeArray1() = knots->Array1();
return Standard_True;
}
//
//=======================================================================
//function : SetMaxDegree
//purpose : Define the Maximum Degree used in the approximation
// If this value is incompatible with the others fields
// this method modify nothing and returns false
//=======================================================================
//
Standard_Boolean AppDef_Variational::SetMaxDegree(const Standard_Integer Degree)
{
if (((Degree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
return Standard_False;
else
{
myMaxDegree=Degree;
InitSmoothCriterion();
return Standard_True;
}
}
//
//=======================================================================
//function : SetMaxSegment
//purpose : Define the maximum number of segments used in the approximation
// If this value is incompatible with the others fields
// this method modify nothing and returns false
//=======================================================================
//
Standard_Boolean AppDef_Variational::SetMaxSegment(const Standard_Integer NbSegment)
{
if ( myWithCutting == Standard_True &&
((myMaxDegree-myNivCont)*NbSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
return Standard_False;
else
{
myMaxSegment=NbSegment;
return Standard_True;
}
}
//
//=======================================================================
//function : SetContinuity
//purpose : Define the Continuity used in the approximation
// If this value is incompatible with the others fields
// this method modify nothing and returns false
//=======================================================================
//
Standard_Boolean AppDef_Variational::SetContinuity(const GeomAbs_Shape C)
{
Standard_Integer NivCont=0;
switch (C) {
case GeomAbs_C0:
NivCont=0;
break ;
case GeomAbs_C1:
NivCont=1;
break ;
case GeomAbs_C2:
NivCont=2;
break ;
default:
throw Standard_ConstructionError();
}
if (((myMaxDegree-NivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
return Standard_False;
else
{
myContinuity= C;
myNivCont=NivCont;
InitSmoothCriterion();
return Standard_True;
}
}
//
//=======================================================================
//function : SetWithMinMax
//purpose : Define if the approximation search to minimize the
// maximum Error or not.
//=======================================================================
//
void AppDef_Variational::SetWithMinMax(const Standard_Boolean MinMax)
{
myWithMinMax=MinMax;
InitSmoothCriterion();
}
//
//=======================================================================
//function : SetWithCutting
//purpose : Define if the approximation can insert new Knots or not.
// If this value is incompatible with the others fields
// this method modify nothing and returns false
//=======================================================================
//
Standard_Boolean AppDef_Variational::SetWithCutting(const Standard_Boolean Cutting)
{
if (Cutting == Standard_False)
{
if (((myMaxDegree-myNivCont)*myKnots->Length()-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
return Standard_False;
else
{
myWithCutting=Cutting;
InitSmoothCriterion();
return Standard_True;
}
}
else
{
if (((myMaxDegree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
return Standard_False;
else
{
myWithCutting=Cutting;
InitSmoothCriterion();
return Standard_True;
}
}
}
//
//=======================================================================
//function : SetCriteriumWeight
//purpose : define the Weights (as percent) associed to the criterium used in
// the optimization.
//=======================================================================
//
void AppDef_Variational::SetCriteriumWeight(const Standard_Real Percent1, const Standard_Real Percent2, const Standard_Real Percent3)
{
if (Percent1 < 0 || Percent2 < 0 || Percent3 < 0 ) throw Standard_DomainError();
Standard_Real Total = Percent1 + Percent2 + Percent3;
myPercent[0] = Percent1/Total;
myPercent[1] = Percent2/Total;
myPercent[2] = Percent3/Total;
InitSmoothCriterion();
}
//
//=======================================================================
//function : SetCriteriumWeight
//purpose : define the Weight (as percent) associed to the
// criterium Order used in the optimization : Others
// weights are updated.
//=======================================================================
//
void AppDef_Variational::SetCriteriumWeight(const Standard_Integer Order, const Standard_Real Percent)
{
if ( Percent < 0 ) throw Standard_DomainError();
if ( Order < 1 || Order > 3 ) throw Standard_ConstructionError();
myPercent[Order-1] = Percent;
Standard_Real Total = myPercent[0] + myPercent[1] + myPercent[2];
myPercent[0] = myPercent[0] / Total;
myPercent[1] = myPercent[1] / Total;
myPercent[2] = myPercent[2] / Total;
InitSmoothCriterion();
}
//
//=======================================================================
//function : SetTolerance
//purpose : define the tolerance used in the approximation.
//=======================================================================
//
void AppDef_Variational::SetTolerance(const Standard_Real Tol)
{
myTolerance=Tol;
InitSmoothCriterion();
}
//
//=======================================================================
//function : SetNbIterations
//purpose : define the number of iterations used in the approximation.
//=======================================================================
//
void AppDef_Variational::SetNbIterations(const Standard_Integer Iter)
{
myNbIterations=Iter;
}
//====================== Private Methodes =============================//
//=======================================================================
//function : TheMotor
//purpose : Smoothing's motor like STRIM routine "MOTLIS"
//=======================================================================
void AppDef_Variational::TheMotor(
Handle(AppDef_SmoothCriterion)& J,
// const Standard_Real WQuadratic,
const Standard_Real ,
const Standard_Real WQuality,
Handle(FEmTool_Curve)& TheCurve,
TColStd_Array1OfReal& Ecarts)
{
// ...
const Standard_Real BigValue = 1.e37, SmallValue = 1.e-6, SmallestValue = 1.e-9;
Handle(TColStd_HArray1OfReal) CurrentTi, NewTi, OldTi;
Handle(TColStd_HArray2OfInteger) Dependence;
Standard_Boolean lestim, lconst, ToOptim, iscut;
Standard_Boolean isnear = Standard_False, again = Standard_True;
Standard_Integer NbEst, ICDANA, NumPnt, Iter;
Standard_Integer MaxNbEst =5;
Standard_Real VOCRI[3] = {BigValue, BigValue, BigValue}, EROLD = BigValue,
VALCRI[3], ERRMAX = BigValue, ERRMOY, ERRQUA;
Standard_Real CBLONG, LNOLD;
Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1;
Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
Handle(FEmTool_Curve) CCurrent, COld, CNew;
Standard_Real EpsLength = SmallValue;
Standard_Real EpsDeg;
Standard_Real e1, e2, e3;
Standard_Real J1min, J2min, J3min;
Standard_Integer iprog;
// (0) Init
J->GetEstimation(e1, e2, e3);
J1min = 1.e-8; J2min = J3min = (e1 + 1.e-8) * 1.e-6;
if(e1 < J1min) e1 = J1min;// Like in
if(e2 < J2min) e2 = J2min;// MOTLIS
if(e3 < J3min) e3 = J3min;
J->SetEstimation(e1, e2, e3);
CCurrent = TheCurve;
CurrentTi = new TColStd_HArray1OfReal(1, myParameters->Length());
CurrentTi->ChangeArray1() = myParameters->Array1();
OldTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
OldTi->ChangeArray1() = CurrentTi->Array1();
COld = CCurrent;
LNOLD = CBLONG = J->EstLength();
Dependence = J->DependenceTable();
J->SetCurve(CCurrent);
FEmTool_Assembly * TheAssembly =
new FEmTool_Assembly (Dependence->Array2(), J->AssemblyTable());
//============ Optimization ============================
// Standard_Integer inagain = 0;
while (again) {
// (1) Loop Optimization / Estimation
lestim = Standard_True;
lconst = Standard_True;
NbEst = 0;
J->SetCurve(CCurrent);
while(lestim) {
// (1.1) Curve's Optimization.
EpsLength = SmallValue * CBLONG / NbrPnt;
CNew = new (FEmTool_Curve) (CCurrent->Dimension(),
CCurrent->NbElements(), CCurrent->Base(), EpsLength);
CNew->Knots() = CCurrent->Knots();
J->SetParameters(CurrentTi);
EpsDeg = Min(WQuality * .1, CBLONG * .001);
Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1());
lconst = Standard_False;
// (1.2) calcul des criteres de qualites et amelioration
// des estimation.
ICDANA = J->QualityValues(J1min, J2min, J3min,
VALCRI[0], VALCRI[1], VALCRI[2]);
if(ICDANA > 0) lconst = Standard_True;
J->ErrorValues(ERRMAX, ERRQUA, ERRMOY);
isnear = ((Sqrt(ERRQUA / NbrPnt) < 2*WQuality) &&
(myNbIterations > 1));
// (1.3) Optimisation des ti par proj orthogonale
// et calcul de l'erreur aux points.
if (isnear) {
NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
Project(CNew, CurrentTi->Array1(),
NewTi->ChangeArray1(),
Ecarts, NumPnt,
ERRMAX, ERRQUA, ERRMOY, 2);
}
else NewTi = CurrentTi;
// (1.4) Progression's test
iprog = 0;
if ((EROLD > WQuality) && (ERRMAX < 0.95*EROLD)) iprog++;
if ((EROLD > WQuality) && (ERRMAX < 0.8*EROLD)) iprog++;
if ((EROLD > WQuality) && (ERRMAX < WQuality)) iprog++;
if ((EROLD > WQuality) && (ERRMAX < 0.99*EROLD)
&& (ERRMAX < 1.1*WQuality)) iprog++;
if ( VALCRI[0] < 0.975 * VOCRI[0]) iprog++;
if ( VALCRI[0] < 0.9 * VOCRI[0]) iprog++;
if ( VALCRI[1] < 0.95 * VOCRI[1]) iprog++;
if ( VALCRI[1] < 0.8 * VOCRI[1]) iprog++;
if ( VALCRI[2] < 0.95 * VOCRI[2]) iprog++;
if ( VALCRI[2] < 0.8 * VOCRI[2]) iprog++;
if ((VOCRI[1]>SmallestValue)&&(VOCRI[2]>SmallestValue)) {
if ((VALCRI[1]/VOCRI[1] + 2*VALCRI[2]/VOCRI[2]) < 2.8) iprog++;
}
if (iprog < 2 && NbEst == 0 ) {
// (1.5) Invalidation of new knots.
VALCRI[0] = VOCRI[0];
VALCRI[1] = VOCRI[1];
VALCRI[2] = VOCRI[2];
ERRMAX = EROLD;
CBLONG = LNOLD;
CCurrent = COld;
CurrentTi = OldTi;
goto L8000; // exit
}
VOCRI[0] = VALCRI[0];
VOCRI[1] = VALCRI[1];
VOCRI[2] = VALCRI[2];
LNOLD = CBLONG;
EROLD = ERRMAX;
CCurrent = CNew;
CurrentTi = NewTi;
// (1.6) Test if the Estimations seems OK, else repeat
NbEst++;
lestim = ( (NbEst<MaxNbEst) && (ICDANA == 2)&& (iprog > 0) );
if (lestim && isnear) {
// (1.7) Optimization of ti by ACR.
std::stable_sort (CurrentTi->begin(), CurrentTi->end());
Standard_Integer Decima = 4;
CCurrent->Length(0., 1., CBLONG);
J->EstLength() = CBLONG;
ACR(CCurrent, CurrentTi->ChangeArray1(), Decima);
lconst = Standard_True;
}
}
// (2) loop of parametric / geometric optimization
Iter = 1;
ToOptim = ((Iter < myNbIterations) && (isnear));
while(ToOptim) {
Iter++;
// (2.1) Save curent results
VOCRI[0] = VALCRI[0];
VOCRI[1] = VALCRI[1];
VOCRI[2] = VALCRI[2];
EROLD = ERRMAX;
LNOLD = CBLONG;
COld = CCurrent;
OldTi->ChangeArray1() = CurrentTi->Array1();
// (2.2) Optimization des ti by ACR.
std::stable_sort (CurrentTi->begin(), CurrentTi->end());
Standard_Integer Decima = 4;
CCurrent->Length(0., 1., CBLONG);
J->EstLength() = CBLONG;
ACR(CCurrent, CurrentTi->ChangeArray1(), Decima);
lconst = Standard_True;
// (2.3) Optimisation des courbes
EpsLength = SmallValue * CBLONG / NbrPnt;
CNew = new (FEmTool_Curve) (CCurrent->Dimension(),
CCurrent->NbElements(), CCurrent->Base(), EpsLength);
CNew->Knots() = CCurrent->Knots();
J->SetParameters(CurrentTi);
EpsDeg = Min(WQuality * .1, CBLONG * .001);
Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1());
CCurrent = CNew;
// (2.4) calcul des criteres de qualites et amelioration
// des estimation.
ICDANA = J->QualityValues(J1min, J2min, J3min, VALCRI[0], VALCRI[1], VALCRI[2]);
if(ICDANA > 0) lconst = Standard_True;
J->GetEstimation(e1, e2, e3);
// (2.5) Optimisation des ti par proj orthogonale
NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
Project(CCurrent, CurrentTi->Array1(),
NewTi->ChangeArray1(),
Ecarts, NumPnt,
ERRMAX, ERRQUA, ERRMOY, 2);
// (2.6) Test de non regression
Standard_Integer iregre = 0;
if (NbrConstraint < NbrPnt) {
if ( (ERRMAX > WQuality) && (ERRMAX > 1.05*EROLD)) iregre++;
if ( (ERRMAX > WQuality) && (ERRMAX > 2*EROLD)) iregre++;
if ( (EROLD > WQuality) && (ERRMAX <= 0.5*EROLD)) iregre--;
}
Standard_Real E1, E2, E3;
J->GetEstimation(E1, E2, E3);
if ( (VALCRI[0] > E1) && (VALCRI[0] > 1.1*VOCRI[0])) iregre++;
if ( (VALCRI[1] > E2) && (VALCRI[1] > 1.1*VOCRI[1])) iregre++;
if ( (VALCRI[2] > E3) && (VALCRI[2] > 1.1*VOCRI[2])) iregre++;
if (iregre >= 2) {
// if (iregre >= 1) {
// (2.7) on restaure l'iteration precedente
VALCRI[0] = VOCRI[0];
VALCRI[1] = VOCRI[1];
VALCRI[2] = VOCRI[2];
ERRMAX = EROLD;
CBLONG = LNOLD;
CCurrent = COld;
CurrentTi->ChangeArray1() = OldTi->Array1();
ToOptim = Standard_False;
}
else { // Iteration is Ok.
CCurrent = CNew;
CurrentTi = NewTi;
}
if (Iter >= myNbIterations) ToOptim = Standard_False;
}
// (3) Decoupe eventuelle
if ( (CCurrent->NbElements() < myMaxSegment) && myWithCutting ) {
// (3.1) Sauvgarde de l'etat precedent
VOCRI[0] = VALCRI[0];
VOCRI[1] = VALCRI[1];
VOCRI[2] = VALCRI[2];
EROLD = ERRMAX;
COld = CCurrent;
OldTi->ChangeArray1() = CurrentTi->Array1();
// (3.2) On arrange les ti : Trie + recadrage sur (0,1)
// ---> On trie, afin d'assurer l'ordre par la suite.
std::stable_sort (CurrentTi->begin(), CurrentTi->end());
if ((CurrentTi->Value(1)!= 0.) ||
(CurrentTi->Value(NbrPnt)!= 1.)) {
Standard_Real t, DelatT =
1.0 /(CurrentTi->Value(NbrPnt)-CurrentTi->Value(1));
for (Standard_Integer ii=2; ii<NbrPnt; ii++) {
t = (CurrentTi->Value(ii)-CurrentTi->Value(1))*DelatT;
CurrentTi->SetValue(ii, t);
}
CurrentTi->SetValue(1, 0.);
CurrentTi->SetValue(NbrPnt, 1.);
}
// (3.3) Insert new Knots
SplitCurve(CCurrent, CurrentTi->Array1(), EpsLength, CNew, iscut);
if (!iscut) again = Standard_False;
else {
CCurrent = CNew;
// New Knots => New Assembly.
J->SetCurve(CNew);
delete TheAssembly;
TheAssembly = new FEmTool_Assembly (Dependence->Array2(),
J->AssemblyTable());
}
}
else { again = Standard_False;}
}
// ================ Great loop end ===================
L8000:
// (4) Compute the best Error.
NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
Project(CCurrent, CurrentTi->Array1(),
NewTi->ChangeArray1(),
Ecarts, NumPnt,
ERRMAX, ERRQUA, ERRMOY, 10);
// (5) field's update
TheCurve = CCurrent;
J->EstLength() = CBLONG;
myParameters->ChangeArray1() = NewTi->Array1();
myCriterium[0] = ERRQUA;
myCriterium[1] = Sqrt(VALCRI[0]);
myCriterium[2] = Sqrt(VALCRI[1]);
myCriterium[3] = Sqrt(VALCRI[2]);
myMaxError = ERRMAX;
myMaxErrorIndex = NumPnt;
if(NbrPnt > NbrConstraint)
myAverageError = ERRMOY / (NbrPnt - NbrConstraint);
else
myAverageError = ERRMOY / NbrConstraint;
delete TheAssembly;
}
//=======================================================================
//function : Optimization
//purpose : (like FORTRAN subroutine MINIMI)
//=======================================================================
void AppDef_Variational::Optimization(Handle(AppDef_SmoothCriterion)& J,
FEmTool_Assembly& A,
const Standard_Boolean ToAssemble,
const Standard_Real EpsDeg,
Handle(FEmTool_Curve)& Curve,
const TColStd_Array1OfReal& Parameters) const
{
Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
NbElm = Curve->NbElements(),
NbDim = Curve->Dimension();
Handle(FEmTool_HAssemblyTable) AssTable;
math_Matrix H(0, MxDeg, 0, MxDeg);
math_Vector G(0, MxDeg), Sol(1, A.NbGlobVar());
Standard_Integer el, dim;
A.GetAssemblyTable(AssTable);
Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
Standard_Real CBLONG = J->EstLength();
// Updating Assembly
if (ToAssemble)
A.NullifyMatrix();
A.NullifyVector();
for (el = 1; el <= NbElm; el++) {
if (ToAssemble) {
J->Hessian(el, 1, 1, H);
for(dim = 1; dim <= NbDim; dim++)
A.AddMatrix(el, dim, dim, H);
}
for(dim = 1; dim <= NbDim; dim++) {
J->Gradient(el, dim, G);
A.AddVector(el, dim, G);
}
}
// Solution of system
if (ToAssemble) {
if(NbConstr != 0) { //Treatment of constraints
AssemblingConstraints(Curve, Parameters, CBLONG, A);
}
A.Solve();
}
A.Solution(Sol);
// Updating J
J->SetCurve(Curve);
J->InputVector(Sol, AssTable);
// Updating Curve and reduction of degree
Standard_Integer Newdeg;
Standard_Real MaxError;
if(NbConstr == 0) {
for(el = 1; el <= NbElm; el++) {
Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError);
}
}
else {
TColStd_Array1OfReal& TabInt = Curve->Knots();
Standard_Integer Icnt = 1, p0 = Parameters.Lower() - myFirstPoint, point;
for(el = 1; el <= NbElm; el++) {
while((Icnt < NbConstr) &&
(Parameters(p0 + myTypConstraints->Value(2 * Icnt - 1)) <= TabInt(el))) Icnt++;
point = p0 + myTypConstraints->Value(2 * Icnt - 1);
if(Parameters(point) <= TabInt(el) || Parameters(point) >= TabInt(el + 1))
Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError);
else
if(Curve->Degree(el) < MxDeg) Curve->SetDegree(el, MxDeg);
}
}
}
void AppDef_Variational::Project(const Handle(FEmTool_Curve)& C,
const TColStd_Array1OfReal& Ti,
TColStd_Array1OfReal& ProjTi,
TColStd_Array1OfReal& Distance,
Standard_Integer& NumPoints,
Standard_Real& MaxErr,
Standard_Real& QuaErr,
Standard_Real& AveErr,
const Standard_Integer NbIterations) const
{
// Initialisation
const Standard_Real Seuil = 1.e-9, Eps = 1.e-12;
MaxErr = QuaErr = AveErr = 0.;
Standard_Integer Ipnt, NItCv, Iter, i, i0 = -myDimension, d0 = Distance.Lower() - 1;
Standard_Real TNew, Dist, T0, Dist0, F1, F2, Aux, DF, Ecart;
Standard_Boolean EnCour;
TColStd_Array1OfReal ValOfC(1, myDimension), FirstDerOfC(1, myDimension),
SecndDerOfC(1, myDimension);
for(Ipnt = 1; Ipnt <= ProjTi.Length(); Ipnt++) {
i0 += myDimension;
TNew = Ti(Ipnt);
EnCour = Standard_True;
NItCv = 0;
Iter = 0;
C->D0(TNew, ValOfC);
Dist = 0;
for(i = 1; i <= myDimension; i++) {
Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
Dist += Aux * Aux;
}
Dist = Sqrt(Dist);
// ------- Newton's method for solving (C'(t),C(t) - P) = 0
while( EnCour ) {
Iter++;
T0 = TNew;
Dist0 = Dist;
C->D2(TNew, SecndDerOfC);
C->D1(TNew, FirstDerOfC);
F1 = F2 = 0.;
for(i = 1; i <= myDimension; i++) {
Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
DF = FirstDerOfC(i);
F1 += Aux*DF; // (C'(t),C(t) - P)
F2 += DF*DF + Aux * SecndDerOfC(i); // ((C'(t),C(t) - P))'
}
if(Abs(F2) < Eps)
EnCour = Standard_False;
else {
// Formula of Newton x(k+1) = x(k) - F(x(k))/F'(x(k))
TNew -= F1 / F2;
if(TNew < 0.) TNew = 0.;
if(TNew > 1.) TNew = 1.;
// Analysis of result
C->D0(TNew, ValOfC);
Dist = 0;
for(i = 1; i <= myDimension; i++) {
Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
Dist += Aux * Aux;
}
Dist = Sqrt(Dist);
Ecart = Dist0 - Dist;
if(Ecart <= -Seuil) {
// Pas d'amelioration on s'arrete
EnCour = Standard_False;
TNew = T0;
Dist = Dist0;
}
else if(Ecart <= Seuil)
// Convergence
NItCv++;
else
NItCv = 0;
if((NItCv >= 2) || (Iter >= NbIterations)) EnCour = Standard_False;
}
}
ProjTi(Ipnt) = TNew;
Distance(d0 + Ipnt) = Dist;
if(Dist > MaxErr) {
MaxErr = Dist;
NumPoints = Ipnt;
}
QuaErr += Dist * Dist;
AveErr += Dist;
}
NumPoints = NumPoints + myFirstPoint - 1;// Setting NumPoints to interval [myFirstPoint, myLastPoint]
}
void AppDef_Variational::ACR(Handle(FEmTool_Curve)& Curve,
TColStd_Array1OfReal& Ti,
const Standard_Integer Decima) const
{
const Standard_Real Eps = 1.e-8;
TColStd_Array1OfReal& Knots = Curve->Knots();
Standard_Integer NbrPnt = Ti.Length(), TiFirst = Ti.Lower(), TiLast = Ti.Upper(),
KFirst = Knots.Lower(), KLast = Knots.Upper();
Standard_Real CbLong, DeltaT, VTest, UNew, UOld, DU, TPara, TOld, DTInv, Ratio;
Standard_Integer ipnt, ii, IElm, IOld, POld, PCnt, ICnt=0;
Standard_Integer NbCntr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
// (1) Calcul de la longueur de courbe
Curve->Length(Ti(TiFirst), Ti(TiLast), CbLong);
// (2) Mise de l'acr dans Ti
if(NbrPnt >= 2) {
// (2.0) Initialisation
DeltaT = (Ti(TiLast) - Ti(TiFirst)) / Decima;
VTest = Ti(TiFirst) + DeltaT;
if(NbCntr > 0) {
PCnt = myTypConstraints->Value(1) - myFirstPoint + TiFirst;
ICnt = 1;
}
else
PCnt = TiLast + 1;
UOld = 0.;
TOld = Ti(TiFirst);
POld = TiFirst;
IElm = KFirst;
IOld = IElm;
Ti(TiFirst) = 0.;
for(ipnt = TiFirst + 1; ipnt <= TiLast; ipnt++) {
while((ICnt <= NbCntr) && (PCnt < ipnt)) {
ICnt++;
PCnt = myTypConstraints->Value(2*ICnt-1) - myFirstPoint + TiFirst;
}
TPara = Ti(ipnt);
if(TPara >= VTest || PCnt == ipnt) {
if ( Ti(TiLast) - TPara <= 1.e-2*DeltaT) {
ipnt = TiLast;
TPara = Ti(ipnt);
}
// (2.2), (2.3) Cacul la longueur de courbe
Curve->Length(Ti(TiFirst), TPara, UNew);
UNew /= CbLong;
while(Knots(IElm + 1) < TPara && IElm < KLast - 1) IElm++;
// (2.4) Mise a jours des parametres de decoupe
DTInv = 1. / (TPara - TOld);
DU = UNew - UOld;
for(ii = IOld+1; ii <= IElm; ii++) {
Ratio = (Knots(ii) - TOld) * DTInv;
Knots(ii) = UOld + Ratio * DU;
}
// (2.5) Mise a jours des parametres de points.
//Very strange loop, because it never works (POld+1 > ipnt-1)
for(ii = POld+1; ii <= ipnt-1; ii++) {
Ratio = ( Ti(ii) - TOld ) * DTInv;
Ti(ii) = UOld + Ratio * DU;
}
Ti(ipnt) = UNew;
UOld = UNew;
IOld = IElm;
TOld = TPara;
POld = ipnt;
}
// --> Nouveau seuil parametrique pour le decimage
if(TPara >= VTest) {
// ii = RealToInt((TPara - VTest + Eps) / DeltaT);
// VTest += (ii + 1) * DeltaT;
VTest += Ceiling((TPara - VTest + Eps) / DeltaT) * DeltaT;
if(VTest > 1. - Eps) VTest = 1.;
}
}
}
// --- On ajuste les valeurs extremes
Ti(TiFirst) = 0.;
Ti(TiLast) = 1.;
ii = TiLast - 1;
while ( Ti(ii) > Knots(KLast) ) {
Ti(ii) = 1.;
--ii;
}
Knots(KFirst) = 0.;
Knots(KLast) = 1.;
}
//----------------------------------------------------------//
// Standard_Integer NearIndex //
// Purpose: searching nearest index of TabPar corresponding //
// given value T (is similar to fortran routine MSRRE2) //
//----------------------------------------------------------//
static Standard_Integer NearIndex(const Standard_Real T,
const TColStd_Array1OfReal& TabPar,
const Standard_Real Eps, Standard_Integer& Flag)
{
Standard_Integer Loi = TabPar.Lower(), Upi = TabPar.Upper();
Flag = 0;
if(T < TabPar(Loi)) { Flag = -1; return Loi;}
if(T > TabPar(Upi)) { Flag = 1; return Upi;}
Standard_Integer Ibeg = Loi, Ifin = Upi, Imidl;
while(Ibeg + 1 != Ifin) {
Imidl = (Ibeg + Ifin) / 2;
if((T >= TabPar(Ibeg)) && (T <= TabPar(Imidl)))
Ifin = Imidl;
else
Ibeg = Imidl;
}
if(Abs(T - TabPar(Ifin)) < Eps) return Ifin;
return Ibeg;
}
//----------------------------------------------------------//
// void GettingKnots //
// Purpose: calculating values of new Knots for elements //
// with degree that is equal Deg //
//----------------------------------------------------------//
static void GettingKnots(const TColStd_Array1OfReal& TabPar,
const Handle(FEmTool_Curve)& InCurve,
const Standard_Integer Deg,
Standard_Integer& NbElm,
TColStd_Array1OfReal& NewKnots)
{
const Standard_Real Eps = 1.e-12;
TColStd_Array1OfReal& OldKnots = InCurve->Knots();
Standard_Integer NbMaxOld = InCurve->NbElements();
Standard_Integer NbMax = NewKnots.Upper(), Ipt, Ipt1, Ipt2;
Standard_Integer el = 0, i1 = OldKnots.Lower(), i0 = i1 - 1, Flag;
Standard_Real TPar;
while((NbElm < NbMax) && (el < NbMaxOld)) {
el++; i0++; i1++; // i0, i1 are indexes of left and right knots of element el
if(InCurve->Degree(el) == Deg) {
NbElm++;
Ipt1 = NearIndex(OldKnots(i0), TabPar, Eps, Flag);
if(Flag != 0) Ipt1 = TabPar.Lower();
Ipt2 = NearIndex(OldKnots(i1), TabPar, Eps, Flag);
if(Flag != 0) Ipt2 = TabPar.Upper();
if(Ipt2 - Ipt1 >= 1) {
Ipt = (Ipt1 + Ipt2) / 2;
if(2 * Ipt == Ipt1 + Ipt2)
TPar = 2. * TabPar(Ipt);
else
TPar = TabPar(Ipt) + TabPar(Ipt + 1);
NewKnots(NbElm) = (OldKnots(i0) + OldKnots(i1) + TPar) / 4.;
}
else
NewKnots(NbElm) = (OldKnots(i0) + OldKnots(i1)) / 2.;
}
}
}
void AppDef_Variational::SplitCurve(const Handle(FEmTool_Curve)& InCurve,
const TColStd_Array1OfReal& Ti,
const Standard_Real CurveTol,
Handle(FEmTool_Curve)& OutCurve,
Standard_Boolean& iscut) const
{
Standard_Integer NbElmOld = InCurve->NbElements();
if(NbElmOld >= myMaxSegment) {iscut = Standard_False; return;}
#ifdef OCCT_DEBUG
Standard_Integer MaxDegree =
#endif
InCurve->Base()->WorkDegree();
Standard_Integer NbElm = NbElmOld;
TColStd_Array1OfReal NewKnots(NbElm + 1, myMaxSegment);
#ifndef OCCT_DEBUG
GettingKnots(Ti, InCurve, InCurve->Base()->WorkDegree(), NbElm, NewKnots);
GettingKnots(Ti, InCurve, InCurve->Base()->WorkDegree() - 1, NbElm, NewKnots);
#else
GettingKnots(Ti, InCurve, MaxDegree, NbElm, NewKnots);
GettingKnots(Ti, InCurve, MaxDegree - 1, NbElm, NewKnots);
#endif
if(NbElm > NbElmOld) {
iscut = Standard_True;
OutCurve = new FEmTool_Curve(InCurve->Dimension(), NbElm, InCurve->Base(), CurveTol);
TColStd_Array1OfReal& OutKnots = OutCurve->Knots();
TColStd_Array1OfReal& InKnots = InCurve->Knots();
Standard_Integer i, i0 = OutKnots.Lower();
for(i = InKnots.Lower(); i <= InKnots.Upper(); i++) OutKnots(i) = InKnots(i);
for(i = NbElmOld + 1; i <= NbElm; i++) OutKnots(i + i0) = NewKnots(i);
std::sort (OutKnots.begin(), OutKnots.end());
}
else
iscut = Standard_False;
}
//=======================================================================
//function : InitSmoothCriterion
//purpose : Initializes the SmoothCriterion
//=======================================================================
void AppDef_Variational::InitSmoothCriterion()
{
const Standard_Real Eps2 = 1.e-6, Eps3 = 1.e-9;
// const Standard_Real J1 = .01, J2 = .001, J3 = .001;
Standard_Real Length;
InitParameters(Length);
mySmoothCriterion->SetParameters(myParameters);
Standard_Real E1, E2, E3;
InitCriterionEstimations(Length, E1, E2, E3);
/*
J1 = 1.e-8; J2 = J3 = (E1 + 1.e-8) * 1.e-6;
if(E1 < J1) E1 = J1;
if(E2 < J2) E2 = J2;
if(E3 < J3) E3 = J3;
*/
mySmoothCriterion->EstLength() = Length;
mySmoothCriterion->SetEstimation(E1, E2, E3);
Standard_Real WQuadratic, WQuality;
if(!myWithMinMax && myTolerance != 0.)
WQuality = myTolerance;
else if (myTolerance == 0.)
WQuality = 1.;
else
WQuality = Max(myTolerance, Eps2 * Length);
Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
WQuadratic = Sqrt((Standard_Real)(myNbPoints - NbConstr)) * WQuality;
if(WQuadratic > Eps3) WQuadratic = 1./ WQuadratic;
if(WQuadratic == 0.) WQuadratic = Max(Sqrt(E1), 1.);
mySmoothCriterion->SetWeight(WQuadratic, WQuality,
myPercent[0], myPercent[1], myPercent[2]);
Handle(PLib_Base) TheBase = new PLib_HermitJacobi(myMaxDegree, myContinuity);
Handle(FEmTool_Curve) TheCurve;
Standard_Integer NbElem;
Standard_Real CurvTol = Eps2 * Length / myNbPoints;
// Decoupe de l'intervalle en fonction des contraintes
if(myWithCutting == Standard_True && NbConstr != 0) {
InitCutting(TheBase, CurvTol, TheCurve);
}
else {
NbElem = 1;
TheCurve = new FEmTool_Curve(myDimension, NbElem, TheBase, CurvTol);
TheCurve->Knots().SetValue(TheCurve->Knots().Lower(),
myParameters->Value(myFirstPoint));
TheCurve->Knots().SetValue(TheCurve->Knots().Upper(),
myParameters->Value(myLastPoint));
}
mySmoothCriterion->SetCurve(TheCurve);
return;
}
//
//=======================================================================
//function : InitParameters
//purpose : Calculation of initial estimation of the multicurve's length
// and parameters for multipoints.
//=======================================================================
//
void AppDef_Variational::InitParameters(Standard_Real& Length)
{
const Standard_Real Eps1 = Precision::Confusion() * .01;
Standard_Real aux, dist;
Standard_Integer i, i0, i1 = 0, ipoint;
Length = 0.;
myParameters->SetValue(myFirstPoint, Length);
for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint; ipoint++) {
i0 = i1; i1 += myDimension;
dist = 0;
for(i = 1; i <= myDimension; i++) {
aux = myTabPoints->Value(i1 + i) - myTabPoints->Value(i0 + i);
dist += aux * aux;
}
Length += Sqrt(dist);
myParameters->SetValue(ipoint, Length);
}
if(Length <= Eps1)
throw Standard_ConstructionError("AppDef_Variational::InitParameters");
for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint - 1; ipoint++)
myParameters->SetValue(ipoint, myParameters->Value(ipoint) / Length);
myParameters->SetValue(myLastPoint, 1.);
// Avec peu de point il y a sans doute sous estimation ...
if(myNbPoints < 10) Length *= (1. + 0.1 / (myNbPoints - 1));
}
//
//=======================================================================
//function : InitCriterionEstimations
//purpose : Calculation of initial estimation of the linear criteria
//
//=======================================================================
//
void AppDef_Variational::InitCriterionEstimations(const Standard_Real Length,
Standard_Real& E1,
Standard_Real& E2,
Standard_Real& E3) const
{
E1 = Length * Length;
const Standard_Real Eps1 = Precision::Confusion() * .01;
math_Vector VTang1(1, myDimension), VTang2(1, myDimension), VTang3(1, myDimension),
VScnd1(1, myDimension), VScnd2(1, myDimension), VScnd3(1, myDimension);
// ========== Treatment of first point =================
Standard_Integer ipnt = myFirstPoint;
EstTangent(ipnt, VTang1);
ipnt++;
EstTangent(ipnt, VTang2);
ipnt++;
EstTangent(ipnt, VTang3);
ipnt = myFirstPoint;
EstSecnd(ipnt, VTang1, VTang2, Length, VScnd1);
ipnt++;
EstSecnd(ipnt, VTang1, VTang3, Length, VScnd2);
// Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 Begin
// Standard_Real Delta = .5 * (myParameters->Value(ipnt) - myParameters->Value(--ipnt));
Standard_Integer anInd = ipnt;
Standard_Real Delta = .5 * (myParameters->Value(anInd) - myParameters->Value(--ipnt));
// Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 End
if(Delta <= Eps1) Delta = 1.;
E2 = VScnd1.Norm2() * Delta;
E3 = (Delta > Eps1) ? VScnd2.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
// ========== Treatment of internal points =================
Standard_Integer CurrPoint = 2;
for(ipnt = myFirstPoint + 1; ipnt < myLastPoint; ipnt++) {
Delta = .5 * (myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1));
if(CurrPoint == 1) {
if(ipnt + 1 != myLastPoint) {
EstTangent(ipnt + 2, VTang3);
EstSecnd(ipnt + 1, VTang1, VTang3, Length, VScnd2);
}
else
EstSecnd(ipnt + 1, VTang1, VTang2, Length, VScnd2);
E2 += VScnd1.Norm2() * Delta;
E3 += (Delta > Eps1) ? VScnd2.Subtracted(VScnd3).Norm2() / (4. * Delta) : 0.;
}
else if(CurrPoint == 2) {
if(ipnt + 1 != myLastPoint) {
EstTangent(ipnt + 2, VTang1);
EstSecnd(ipnt + 1, VTang2, VTang1, Length, VScnd3);
}
else
EstSecnd(ipnt + 1, VTang2, VTang3, Length, VScnd3);
E2 += VScnd2.Norm2() * Delta;
E3 += (Delta > Eps1) ? VScnd3.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
}
else {
if(ipnt + 1 != myLastPoint) {
EstTangent(ipnt + 2, VTang2);
EstSecnd(ipnt + 1, VTang3, VTang2, Length, VScnd1);
}
else
EstSecnd(ipnt + 1, VTang3, VTang1, Length, VScnd1);
E2 += VScnd3.Norm2() * Delta;
E3 += (Delta > Eps1) ? VScnd1.Subtracted(VScnd2).Norm2() / (4. * Delta) : 0.;
}
CurrPoint++; if(CurrPoint == 4) CurrPoint = 1;
}
// ========== Treatment of last point =================
Delta = .5 * (myParameters->Value(myLastPoint) - myParameters->Value(myLastPoint - 1));
if(Delta <= Eps1) Delta = 1.;
Standard_Real aux;
if(CurrPoint == 1) {
E2 += VScnd1.Norm2() * Delta;
aux = VScnd1.Subtracted(VScnd3).Norm2();
E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
}
else if(CurrPoint == 2) {
E2 += VScnd2.Norm2() * Delta;
aux = VScnd2.Subtracted(VScnd1).Norm2();
E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
}
else {
E2 += VScnd3.Norm2() * Delta;
aux = VScnd3.Subtracted(VScnd2).Norm2();
E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
}
aux = Length * Length;
E2 *= aux; E3 *= aux;
}
//
//=======================================================================
//function : EstTangent
//purpose : Calculation of estimation of the Tangent
// (see fortran routine MMLIPRI)
//=======================================================================
//
void AppDef_Variational::EstTangent(const Standard_Integer ipnt,
math_Vector& VTang) const
{
Standard_Integer i ;
const Standard_Real Eps1 = Precision::Confusion() * .01;
const Standard_Real EpsNorm = 1.e-9;
Standard_Real Wpnt = 1.;
if(ipnt == myFirstPoint) {
// Estimation at first point
if(myNbPoints < 3)
Wpnt = 0.;
else {
Standard_Integer adr1 = 1, adr2 = adr1 + myDimension,
adr3 = adr2 + myDimension;
math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
// Parabolic interpolation
// if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
// first derivative for t=0 is A1 = ((d2-1)*P1 + P2 - d2*P3)/(d*(1-d))
// d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*d
Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
Standard_Real V2 = 0.;
if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
if(V2 > Eps1) {
Standard_Real d = V1 / (V1 + V2), d1;
d1 = 1. / (d * (1 - d)); d *= d;
VTang = ((d - 1.) * Pnt1 + Pnt2 - d * Pnt3) * d1;
}
else {
// Simple 2-point estimation
VTang = Pnt2 - Pnt1;
}
}
}
else if (ipnt == myLastPoint) {
// Estimation at last point
if(myNbPoints < 3)
Wpnt = 0.;
else {
Standard_Integer adr1 = (myLastPoint - 3) * myDimension + 1,
adr2 = adr1 + myDimension,
adr3 = adr2 + myDimension;
math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
// Parabolic interpolation
// if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
// first derivative for t=1 is 2*A2 + A1 = ((d2+1)*P1 - P2 - d2*P3)/(d*(1-d))
// d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*(d-2)
Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
Standard_Real V2 = 0.;
if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
if(V2 > Eps1) {
Standard_Real d = V1 / (V1 + V2), d1;
d1 = 1. / (d * (1 - d)); d *= d - 2;
VTang = ((d + 1.) * Pnt1 - Pnt2 - d * Pnt3) * d1;
}
else {
// Simple 2-point estimation
VTang = Pnt3 - Pnt2;
}
}
}
else {
Standard_Integer adr1 = (ipnt - myFirstPoint - 1) * myDimension + 1,
adr2 = adr1 + 2 * myDimension;
math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
VTang = Pnt2 - Pnt1;
}
Standard_Real Vnorm = VTang.Norm();
if(Vnorm <= EpsNorm)
VTang.Init(0.);
else
VTang /= Vnorm;
// Estimation with constraints
Standard_Real Wcnt = 0.;
Standard_Integer IdCnt = 1;
// Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
math_Vector VCnt(1, myDimension, 0.);
if(NbConstr > 0) {
while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
IdCnt <= NbConstr) IdCnt++;
if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
(myTypConstraints->Value(2 * IdCnt) >= 1)) {
Wcnt = 1.;
Standard_Integer i0 = 2 * myDimension * (IdCnt - 1), k = 0;
for( i = 1; i <= myNbP3d; i++) {
for(Standard_Integer j = 1; j <= 3; j++)
VCnt(++k) = myTabConstraints->Value(++i0);
i0 += 3;
}
for(i = 1; i <= myNbP2d; i++) {
for(Standard_Integer j = 1; j <= 2; j++)
VCnt(++k) = myTabConstraints->Value(++i0);
i0 += 2;
}
}
}
// Averaging of estimation
Standard_Real Denom = Wpnt + Wcnt;
if(Denom == 0.) Denom = 1.;
else Denom = 1. / Denom;
VTang = (Wpnt * VTang + Wcnt * VCnt) * Denom;
Vnorm = VTang.Norm();
if(Vnorm <= EpsNorm)
VTang.Init(0.);
else
VTang /= Vnorm;
}
//
//=======================================================================
//function : EstSecnd
//purpose : Calculation of estimation of the second derivative
// (see fortran routine MLIMSCN)
//=======================================================================
//
void AppDef_Variational::EstSecnd(const Standard_Integer ipnt,
const math_Vector& VTang1,
const math_Vector& VTang2,
const Standard_Real Length,
math_Vector& VScnd) const
{
Standard_Integer i ;
const Standard_Real Eps = 1.e-9;
Standard_Real Wpnt = 1.;
Standard_Real aux;
if(ipnt == myFirstPoint)
aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt);
else if(ipnt == myLastPoint)
aux = myParameters->Value(ipnt) - myParameters->Value(ipnt - 1);
else
aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1);
if(aux <= Eps)
aux = 1.;
else
aux = 1. / aux;
VScnd = (VTang2 - VTang1) * aux;
// Estimation with constraints
Standard_Real Wcnt = 0.;
Standard_Integer IdCnt = 1;
// Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
math_Vector VCnt(1, myDimension, 0.);
if(NbConstr > 0) {
while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
IdCnt <= NbConstr) IdCnt++;
if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
(myTypConstraints->Value(2 * IdCnt) >= 2)) {
Wcnt = 1.;
Standard_Integer i0 = 2 * myDimension * (IdCnt - 1) + 3, k = 0;
for( i = 1; i <= myNbP3d; i++) {
for(Standard_Integer j = 1; j <= 3; j++)
VCnt(++k) = myTabConstraints->Value(++i0);
i0 += 3;
}
i0--;
for(i = 1; i <= myNbP2d; i++) {
for(Standard_Integer j = 1; j <= 2; j++)
VCnt(++k) = myTabConstraints->Value(++i0);
i0 += 2;
}
}
}
// Averaging of estimation
Standard_Real Denom = Wpnt + Wcnt;
if(Denom == 0.) Denom = 1.;
else Denom = 1. / Denom;
VScnd = (Wpnt * VScnd + (Wcnt * Length) * VCnt) * Denom;
}
//
//=======================================================================
//function : InitCutting
//purpose : Realisation of curve's cutting a priory accordingly to
// constraints (see fortran routine MLICUT)
//=======================================================================
//
void AppDef_Variational::InitCutting(const Handle(PLib_Base)& aBase,
const Standard_Real CurvTol,
Handle(FEmTool_Curve)& aCurve) const
{
// Definition of number of elements
Standard_Integer ORCMx = -1, NCont = 0, i, kk, NbElem;
Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
for(i = 1; i <= NbConstr; i++) {
kk = Abs(myTypConstraints->Value(2 * i)) + 1;
ORCMx = Max(ORCMx, kk);
NCont += kk;
}
if(ORCMx > myMaxDegree - myNivCont)
throw Standard_ConstructionError("AppDef_Variational::InitCutting");
Standard_Integer NLibre = Max(myMaxDegree - myNivCont - (myMaxDegree + 1) / 4,
myNivCont + 1);
NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
while((NbElem > myMaxSegment) && (NLibre < myMaxDegree - myNivCont)) {
NLibre++;
NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
}
if(NbElem > myMaxSegment)
throw Standard_ConstructionError("AppDef_Variational::InitCutting");
aCurve = new FEmTool_Curve(myDimension, NbElem, aBase, CurvTol);
Standard_Integer NCnt = (NCont - 1) / NbElem + 1;
Standard_Integer NPlus = NbElem - (NCnt * NbElem - NCont);
TColStd_Array1OfReal& Knot = aCurve->Knots();
Standard_Integer IDeb = 0, IFin = NbConstr + 1,
NDeb = 0, NFin = 0,
IndEl = Knot.Lower(), IUpper = Knot.Upper(),
NbEl = 0;
Knot(IndEl) = myParameters->Value(myFirstPoint);
Knot(IUpper) = myParameters->Value(myLastPoint);
while(NbElem - NbEl > 1) {
IndEl++; NbEl++;
if(NPlus == 0) NCnt--;
while(NDeb < NCnt && IDeb < IFin) {
IDeb++;
NDeb += Abs(myTypConstraints->Value(2 * IDeb)) + 1;
}
if(NDeb == NCnt) {
NDeb = 0;
if(NPlus == 1 &&
myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) > Knot(IndEl - 1))
Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
else
Knot(IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) +
myParameters->Value(myTypConstraints->Value(2 * IDeb + 1))) / 2;
}
else {
NDeb -= NCnt;
Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
}
NPlus--;
if(NPlus == 0) NCnt--;
if(NbElem - NbEl == 1) break;
NbEl++;
while(NFin < NCnt && IDeb < IFin) {
IFin--;
NFin += Abs(myTypConstraints->Value(2 * IFin)) + 1;
}
if(NFin == NCnt) {
NFin = 0;
Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
}
else {
NFin -= NCnt;
if(myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) < Knot(IUpper - IndEl + 1))
Knot(IUpper + 1 - IndEl) = myParameters->Value(myTypConstraints->Value(2 * IFin - 1));
else
Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
}
}
}
//=======================================================================
//function : Adjusting
//purpose : Smoothing's adjusting like STRIM routine "MAJLIS"
//=======================================================================
void AppDef_Variational::Adjusting(
Handle(AppDef_SmoothCriterion)& J,
Standard_Real& WQuadratic,
Standard_Real& WQuality,
Handle(FEmTool_Curve)& TheCurve,
TColStd_Array1OfReal& Ecarts)
{
// std::cout << "=========== Adjusting =============" << std::endl;
/* Initialized data */
const Standard_Integer mxiter = 2;
const Standard_Real eps1 = 1e-6;
Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1;
Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
Standard_Real CurvTol = eps1 * J->EstLength() / NbrPnt;
/* Local variables */
Standard_Integer iter, ipnt;
Standard_Real ecart, emold, erold, tpara;
Standard_Real vocri[4], j1cibl, vtest, vseuil;
Standard_Integer i, numint, flag;
TColStd_Array1OfReal tbpoid(myFirstPoint, myLastPoint);
Standard_Boolean loptim, lrejet;
Handle(AppDef_SmoothCriterion) JNew;
Handle(FEmTool_Curve) CNew;
Standard_Real E1, E2, E3;
/* (0.b) Initialisations */
loptim = Standard_True;
iter = 0;
tbpoid.Init(1.);
/* ============ boucle sur le moteur de lissage ============== */
vtest = WQuality * .9;
j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint));
while(loptim) {
++iter;
/* (1) Sauvegarde de l'etat precedents */
vocri[0] = myCriterium[0];
vocri[1] = myCriterium[1];
vocri[2] = myCriterium[2];
vocri[3] = myCriterium[3];
erold = myMaxError;
emold = myAverageError;
/* (2) Augmentation du poids des moindre carre */
if (j1cibl > vtest) {
WQuadratic = j1cibl / vtest * WQuadratic;
}
/* (3) Augmentation du poid associe aux points a problemes */
vseuil = WQuality * .88;
for (ipnt = myFirstPoint; ipnt <= myLastPoint; ++ipnt) {
if (Ecarts(ipnt) > vtest) {
ecart = (Ecarts(ipnt) - vseuil) / WQuality;
tbpoid(ipnt) = (ecart * 3 + 1.) * tbpoid(ipnt);
}
}
/* (4) Decoupe force */
if (TheCurve->NbElements() < myMaxSegment && myWithCutting) {
numint = NearIndex(myParameters->Value(myMaxErrorIndex), TheCurve->Knots(), 0, flag);
tpara = (TheCurve->Knots()(numint) + TheCurve->Knots()(numint + 1) +
myParameters->Value(myMaxErrorIndex) * 2) / 4;
CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements() + 1,
TheCurve->Base(), CurvTol);
for(i = 1; i <= numint; i++) CNew->Knots()(i) = TheCurve->Knots()(i);
for(i = numint + 1; i <= TheCurve->Knots().Length(); i++)
CNew->Knots()(i + 1) = TheCurve->Knots()(i);
CNew->Knots()(numint + 1) = tpara;
} else {
CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements(), TheCurve->Base(), CurvTol);
CNew->Knots() = TheCurve->Knots();
}
JNew = new AppDef_LinearCriteria(mySSP, myFirstPoint, myLastPoint);
JNew->EstLength() = J->EstLength();
J->GetEstimation(E1, E2, E3);
JNew->SetEstimation(E1, E2, E3);
JNew->SetParameters(myParameters);
JNew->SetWeight(WQuadratic, WQuality, myPercent[0], myPercent[1], myPercent[2]);
JNew->SetWeight(tbpoid);
JNew->SetCurve(CNew);
/* (5) Relissage */
TheMotor(JNew, WQuadratic, WQuality, CNew, Ecarts);
/* (6) Tests de rejet */
j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint));
vseuil = Sqrt(vocri[1]) + (erold - myMaxError) * 4;
lrejet = ((myMaxError > WQuality) && (myMaxError > erold * 1.01)) || (Sqrt(myCriterium[1]) > vseuil * 1.05);
if (lrejet) {
myCriterium[0] = vocri[0];
myCriterium[1] = vocri[1];
myCriterium[2] = vocri[2];
myCriterium[3] = vocri[3];
myMaxError = erold;
myAverageError = emold;
loptim = Standard_False;
}
else {
J = JNew;
TheCurve = CNew;
J->SetCurve(TheCurve);
}
/* (7) Test de convergence */
if (((iter >= mxiter) && (myMaxSegment == CNew->NbElements())) || myMaxError < WQuality) {
loptim = Standard_False;
}
}
}
static Standard_Boolean NotParallel(gp_Vec& T, gp_Vec& V)
{
V = T;
V.SetX(V.X() + 1.);
if (V.CrossMagnitude(T) > 1.e-12)
return Standard_True;
V.SetY(V.Y() + 1.);
if (V.CrossMagnitude(T) > 1.e-12)
return Standard_True;
V.SetZ(V.Z() + 1.);
if (V.CrossMagnitude(T) > 1.e-12)
return Standard_True;
return Standard_False;
}
void AppDef_Variational::AssemblingConstraints(const Handle(FEmTool_Curve)& Curve,
const TColStd_Array1OfReal& Parameters,
const Standard_Real CBLONG,
FEmTool_Assembly& A ) const
{
Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
NbElm = Curve->NbElements(),
NbDim = Curve->Dimension();
TColStd_Array1OfReal G0(0, MxDeg), G1(0, MxDeg), G2(0, MxDeg);
math_Vector V0((Standard_Real*)&G0(0), 0, MxDeg),
V1((Standard_Real*)&G1(0), 0, MxDeg),
V2((Standard_Real*)&G2(0), 0, MxDeg);
Standard_Integer IndexOfConstraint, Ng3d, Ng2d, NBeg2d, NPass, NgPC1,
NTang3d, NTang2d,
Point, TypOfConstr,
p0 = Parameters.Lower() - myFirstPoint,
curel = 1, el, i, ipnt, ityp, j, k, pnt, curdim,
jt, Ntheta = 6 * myNbP3d + 2 * myNbP2d;
Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
// Ng3d = 3 * NbConstr + 2 * myNbTangPoints + 5 * myNbCurvPoints;
// Ng2d = 2 * NbConstr + 1 * myNbTangPoints + 3 * myNbCurvPoints;
Ng3d = 3 * NbConstr + 3 * myNbTangPoints + 5 * myNbCurvPoints;
Ng2d = 2 * NbConstr + 2 * myNbTangPoints + 3 * myNbCurvPoints;
NBeg2d = Ng3d * myNbP3d;
// NgPC1 = NbConstr + myNbCurvPoints;
NgPC1 = NbConstr + myNbTangPoints + myNbCurvPoints;
NPass = 0;
NTang3d = 3 * NgPC1;
NTang2d = 2 * NgPC1;
TColStd_Array1OfReal& Intervals = Curve->Knots();
Standard_Real t, R1, R2;
Handle(PLib_Base) myBase = Curve->Base();
Handle(PLib_HermitJacobi) myHermitJacobi = Handle(PLib_HermitJacobi)::DownCast (myBase);
Standard_Integer Order = myHermitJacobi->NivConstr() + 1;
Standard_Real UFirst, ULast, coeff, c0, mfact, mfact1;
A.NullifyConstraint();
ipnt = -1;
ityp = 0;
for(i = 1; i <= NbConstr; i++) {
ipnt += 2; ityp += 2;
Point = myTypConstraints->Value(ipnt);
TypOfConstr = myTypConstraints->Value(ityp);
t = Parameters(p0 + Point);
for(el = curel; el <= NbElm; ) {
if( t <= Intervals(++el) ) {
curel = el - 1;
break;
}
}
UFirst = Intervals(curel); ULast = Intervals(curel + 1);
coeff = (ULast - UFirst)/2.; c0 = (ULast + UFirst)/2.;
t = (t - c0) / coeff;
if(TypOfConstr == 0) {
myBase->D0(t, G0);
for(k = 1; k < Order; k++) {
mfact = Pow(coeff, k);
G0(k) *= mfact;
G0(k + Order) *= mfact;
}
}
else if(TypOfConstr == 1) {
myBase->D1(t, G0, G1);
for(k = 1; k < Order; k++) {
mfact = Pow(coeff, k);
G0(k) *= mfact;
G0(k + Order) *= mfact;
G1(k) *= mfact;
G1(k + Order) *= mfact;
}
mfact = 1./coeff;
for(k = 0; k <= MxDeg; k++) {
G1(k) *= mfact;
}
}
else {
myBase->D2(t, G0, G1, G2);
for(k = 1; k < Order; k++) {
mfact = Pow(coeff, k);
G0(k) *= mfact;
G0(k + Order) *= mfact;
G1(k) *= mfact;
G1(k + Order) *= mfact;
G2(k) *= mfact;
G2(k + Order) *= mfact;
}
mfact = 1. / coeff;
mfact1 = mfact / coeff;
for(k = 0; k <= MxDeg; k++) {
G1(k) *= mfact;
G2(k) *= mfact1;
}
}
NPass++;
j = NbDim * (Point - myFirstPoint);
Standard_Integer n0 = NPass;
curdim = 0;
for(pnt = 1; pnt <= myNbP3d; pnt++) {
IndexOfConstraint = n0;
for(k = 1; k <= 3; k++) {
curdim++;
A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
IndexOfConstraint += NgPC1;
}
j +=3;
n0 += Ng3d;
}
n0 = NPass + NBeg2d;
for(pnt = 1; pnt <= myNbP2d; pnt++) {
IndexOfConstraint = n0;
for(k = 1; k <= 2; k++) {
curdim++;
A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
IndexOfConstraint += NgPC1;
}
j +=2;
n0 += Ng2d;
}
/* if(TypOfConstr == 1) {
IndexOfConstraint = NTang3d + 1;
jt = Ntheta * (i - 1);
for(pnt = 1; pnt <= myNbP3d; pnt++) {
for(k = 1; k <= 3; k++) {
A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.);
A.AddConstraint(IndexOfConstraint + 1, curel, k, myTtheta->Value(jt + 3 + k) * V1, 0.);
}
IndexOfConstraint += Ng3d;
jt += 6;
}
IndexOfConstraint = NBeg2d + NTang2d + 1;
for(pnt = 1; pnt <= myNbP2d; pnt++) {
for(k = 1; k <= 2; k++) {
A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.);
}
IndexOfConstraint += Ng2d;
jt += 2;
}
NTang3d += 2;
NTang2d += 1;
} */
if(TypOfConstr == 1) {
NPass++;
n0 = NPass;
j = 2 * NbDim * (i - 1);
curdim = 0;
for(pnt = 1; pnt <= myNbP3d; pnt++) {
IndexOfConstraint = n0;
for(k = 1; k <= 3; k++) {
curdim++;
A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
IndexOfConstraint += NgPC1;
}
n0 += Ng3d;
j += 6;
}
n0 = NPass + NBeg2d;
for(pnt = 1; pnt <= myNbP2d; pnt++) {
IndexOfConstraint = n0;
for(k = 1; k <= 2; k++) {
curdim++;
A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
IndexOfConstraint += NgPC1;
}
n0 += Ng2d;
j += 4;
}
}
if(TypOfConstr == 2) {
NPass++;
n0 = NPass;
j = 2 * NbDim * (i - 1);
curdim = 0;
for(pnt = 1; pnt <= myNbP3d; pnt++) {
IndexOfConstraint = n0;
for(k = 1; k <= 3; k++) {
curdim++;
A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
IndexOfConstraint += NgPC1;
}
n0 += Ng3d;
j += 6;
}
n0 = NPass + NBeg2d;
for(pnt = 1; pnt <= myNbP2d; pnt++) {
IndexOfConstraint = n0;
for(k = 1; k <= 2; k++) {
curdim++;
A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
IndexOfConstraint += NgPC1;
}
n0 += Ng2d;
j += 4;
}
j = 2 * NbDim * (i - 1) + 3;
jt = Ntheta * (i - 1);
IndexOfConstraint = NTang3d + 1;
curdim = 0;
for(pnt = 1; pnt <= myNbP3d; pnt++) {
R1 = 0.; R2 = 0.;
for(k = 1; k <= 3; k++) {
R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
R2 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + 3 + k);
}
R1 *= CBLONG * CBLONG;
R2 *= CBLONG * CBLONG;
for(k = 1; k <= 3; k++) {
curdim++;
if(k > 1) R1 = R2 = 0.;
A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
A.AddConstraint(IndexOfConstraint + 1, curel, curdim, myTfthet->Value(jt + 3 + k) * V2, R2);
}
IndexOfConstraint += Ng3d;
j += 6;
jt += 6;
}
j--;
IndexOfConstraint = NBeg2d + NTang2d + 1;
for(pnt = 1; pnt <= myNbP2d; pnt++) {
R1 = 0.;
for(k = 1; k <= 2; k++) {
R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
}
R1 *= CBLONG * CBLONG;
for(k = 1; k <= 2; k++) {
curdim++;
if(k > 1) R1 = 0.;
A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
}
IndexOfConstraint += Ng2d;
j += 4;
jt += 2;
}
NTang3d += 2;
NTang2d += 1;
}
}
}
Standard_Boolean AppDef_Variational::InitTthetaF(const Standard_Integer ndimen,
const AppParCurves_Constraint typcon,
const Standard_Integer begin,
const Standard_Integer jndex)
{
if ((ndimen < 2)||(ndimen >3))
return Standard_False;
gp_Vec T, V;
gp_Vec theta1, theta2;
gp_Vec F;
Standard_Real XX, XY, YY, XZ, YZ, ZZ;
if ((typcon == AppParCurves_TangencyPoint)||(typcon == AppParCurves_CurvaturePoint))
{
T.SetX(myTabConstraints->Value(jndex));
T.SetY(myTabConstraints->Value(jndex + 1));
if (ndimen == 3)
T.SetZ(myTabConstraints->Value(jndex + 2));
else
T.SetZ(0.);
if (ndimen == 2)
{
V.SetX(0.);
V.SetY(0.);
V.SetZ(1.);
}
if (ndimen == 3)
if (!NotParallel(T, V))
return Standard_False;
theta1 = V ^ T;
theta1.Normalize();
myTtheta->SetValue(begin, theta1.X());
myTtheta->SetValue(begin + 1, theta1.Y());
if (ndimen == 3)
{
theta2 = T ^ theta1;
theta2.Normalize();
myTtheta->SetValue(begin + 2, theta1.Z());
myTtheta->SetValue(begin + 3, theta2.X());
myTtheta->SetValue(begin + 4, theta2.Y());
myTtheta->SetValue(begin + 5, theta2.Z());
}
// Calculation of myTfthet
if (typcon == AppParCurves_CurvaturePoint)
{
XX = Pow(T.X(), 2);
XY = T.X() * T.Y();
YY = Pow(T.Y(), 2);
if (ndimen == 2)
{
F.SetX(YY * theta1.X() - XY * theta1.Y());
F.SetY(XX * theta1.Y() - XY * theta1.X());
myTfthet->SetValue(begin, F.X());
myTfthet->SetValue(begin + 1, F.Y());
}
if (ndimen == 3)
{
XZ = T.X() * T.Z();
YZ = T.Y() * T.Z();
ZZ = Pow(T.Z(), 2);
F.SetX((ZZ + YY) * theta1.X() - XY * theta1.Y() - XZ * theta1.Z());
F.SetY((XX + ZZ) * theta1.Y() - XY * theta1.X() - YZ * theta1.Z());
F.SetZ((XX + YY) * theta1.Z() - XZ * theta1.X() - YZ * theta1.Y());
myTfthet->SetValue(begin, F.X());
myTfthet->SetValue(begin + 1, F.Y());
myTfthet->SetValue(begin + 2, F.Z());
F.SetX((ZZ + YY) * theta2.X() - XY * theta2.Y() - XZ * theta2.Z());
F.SetY((XX + ZZ) * theta2.Y() - XY * theta2.X() - YZ * theta2.Z());
F.SetZ((XX + YY) * theta2.Z() - XZ * theta2.X() - YZ * theta2.Y());
myTfthet->SetValue(begin + 3, F.X());
myTfthet->SetValue(begin + 4, F.Y());
myTfthet->SetValue(begin + 5, F.Z());
}
}
}
return Standard_True;
}
View Git Blame Copy Permalink