OpenCascade/occt
1
0
Fork 0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-08-29 14:00:49 +03:00
Code Activity

0024708: Convertation of the generic classes to the non-generic. Part 2

Browse Source 
Generic classes from "AppParCurves" package:
    "AppDef_SmoothCriterion", "AppDef_LinearCriteria" and "AppDef_Variational" moved to the corresponding non-generic classes "AppDef_SmoothCriterion", "AppDef_LinearCriteria" and "AppDef_Variational" to "AppDef" package. Also several "*.cxx" files of "AppDef_Variational" class merged to one ".cxx".
Generic class from "IntImp" package:
    "IntImp_ZerCOnSSParFunc" moved to the corresponding non-generic class "IntPatch_CSFunction" to "IntPatch" package.
Next unused generic classes were removed:

- IntCurveSurface_SurfaceTool
- Intf_InterferencePolygon3d
And some other minor changes.
This commit is contained in:
dln
2014-03-05 18:22:43 +04:00
committed by abv
parent 84c71f29e4
commit f62de37212
37 changed files with 3482 additions and 4806 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
src/AppParCurves/AppParCurves.cdl
View File

@@ -21,8 +21,7 @@ package AppParCurves
-- described by the tool MLineTool.
-- The result of the approximation will be a MultiCurve.
uses math, FEmTool, gp, TColgp, StdFail, TColStd, TCollection, Standard, MMgt, GeomAbs,
PLib
uses math, FEmTool, gp, TColgp, StdFail, TColStd, TCollection, Standard, MMgt, GeomAbs, PLib
is
@@ -82,18 +81,6 @@ is
-- sum||C(ui)-Qi||2 with a gradient method.
-- The Result is a Bspline set.
deferred class SmoothCriterion;
generic class LinearCriteria;
class MyCriterion;
---Level: Internal
generic class Variational;
---Purpose: computes the approximation of a Multiline by
-- Variational optimization.
--- instantiate classes:
--
@@ -105,7 +92,6 @@ is
(ConstraintCouple,Array1OfConstraintCouple);
class Array1OfMultiPoint instantiates Array1 from TCollection
(MultiPoint);

141
src/AppParCurves/AppParCurves_LinearCriteria.cdl
View File

@@ -1,141 +0,0 @@
-- Created on: 1997-09-11
-- Created by: Philippe MANGIN
-- Copyright (c) 1997-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.
generic class LinearCriteria from AppParCurves
(MultiLine as any;
ToolLine as any) -- as ToolLine(MultiLine)
inherits SmoothCriterion from AppParCurves
---Purpose: defined an Linear Criteria to used in variational
-- Smoothing of points.
uses
Vector from math,
Matrix from math,
Curve from FEmTool,
HAssemblyTable from FEmTool,
ElementaryCriterion from FEmTool,
HArray2OfInteger from TColStd,
HArray1OfReal from TColStd,
Array1OfReal from TColStd
raises
NotImplemented,
DomainError
is
Create(SSP: MultiLine;
FirstPoint, LastPoint: Integer) returns LinearCriteria;
SetParameters(me : mutable; Parameters : HArray1OfReal);
SetCurve(me : mutable; C :Curve from FEmTool)
is static;
GetCurve(me; C : out Curve from FEmTool)
is static;
SetEstimation(me : mutable; E1, E2, E3 : Real)
is static;
EstLength(me : mutable)
---C++: return &
returns Real is static;
GetEstimation(me; E1, E2, E3 : out Real)
is static;
AssemblyTable(me)
returns HAssemblyTable from FEmTool
is static;
DependenceTable(me)
returns HArray2OfInteger from TColStd
is static;
QualityValues (me : mutable; J1min, J2min, J3min : Real;
J1, J2, J3 : out Real)
returns Integer is static;
ErrorValues(me : mutable;
MaxError, QuadraticError, AverageError : out Real)
is static;
Hessian(me : mutable ;
Element : Integer;
Dimension1 : Integer;
Dimension2 : Integer;
H : out Matrix from math)
raises DomainError -- If DependenceTable(Dimension1,Dimension2) is False
is static;
Gradient(me : mutable;
Element : Integer;
Dimension : Integer;
G : out Vector from math)
is static;
InputVector(me : mutable; X : Vector from math;
AssTable : HAssemblyTable from FEmTool)
---Purpose: Convert the assembly Vector in an Curve;
--
raises DomainError;
SetWeight(me: mutable;
QuadraticWeight, QualityWeight : Real;
percentJ1, percentJ2, percentJ3 : Real)
is static;
GetWeight(me; QuadraticWeight, QualityWeight : out Real)
is static;
SetWeight(me: mutable;
Weight : Array1OfReal)
is static;
BuildCache(me: mutable; E : Integer) is private;
fields
mySSP : MultiLine;
myParameters : HArray1OfReal;
myCache : HArray1OfReal;
myCriteria : ElementaryCriterion from FEmTool[3];
myEstimation: Real[3];
myQuadraticWeight, myQualityWeight : Real;
myPercent : Real[3];
myPntWeight : Array1OfReal;
myCurve : Curve from FEmTool;
myLength : Real;
myE : Integer;
IF, IL : Integer;
end LinearCriteria;

757
src/AppParCurves/AppParCurves_LinearCriteria.gxx
View File

@@ -1,757 +0,0 @@
// Created on: 1998-11-30
// Created by: Igor FEOKTISTOV
// Copyright (c) 1998-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 <PLib_Base.hxx>
#include <PLib_JacobiPolynomial.hxx>
#include <PLib_HermitJacobi.hxx>
#include <GeomAbs_Shape.hxx>
#include <TColStd_HArray2OfReal.hxx>
#include <FEmTool_LinearTension.hxx>
#include <FEmTool_LinearFlexion.hxx>
#include <FEmTool_LinearJerk.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TColgp_Array1OfPnt2d.hxx>
#include <gp_Pnt2d.hxx>
#include <gp_Pnt.hxx>
#include <math_Matrix.hxx>
#include <math_Gauss.hxx>
static Standard_Integer order(const Handle(PLib_Base)& B)
{
return (*( Handle(PLib_HermitJacobi)*)&B)->NivConstr();
}
//=======================================================================
//function :
//purpose :
//=======================================================================
AppParCurves_LinearCriteria::AppParCurves_LinearCriteria(const MultiLine& SSP,
const Standard_Integer FirstPoint,
const Standard_Integer LastPoint):
mySSP(SSP),
myPntWeight(FirstPoint, LastPoint),
myE(0)
{
myPntWeight.Init(1.);
}
//=======================================================================
//function :
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::SetParameters(const Handle(TColStd_HArray1OfReal)& Parameters)
{
myParameters = Parameters;
myE = 0; // Cache become invalid.
}
//=======================================================================
//function : SetCurve
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::SetCurve(const Handle(FEmTool_Curve)& C)
{
if(myCurve.IsNull()) {
myCurve = C;
Standard_Integer MxDeg = myCurve->Base()->WorkDegree(),
NbDim = myCurve->Dimension(),
Order = order(myCurve->Base());
GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
switch (Order) {
case 0 : ConstraintOrder = GeomAbs_C0;
break;
case 1 : ConstraintOrder = GeomAbs_C1;
break;
case 2 : ConstraintOrder = GeomAbs_C2;
}
myCriteria[0] = new FEmTool_LinearTension(MxDeg, ConstraintOrder);
myCriteria[1] = new FEmTool_LinearFlexion(MxDeg, ConstraintOrder);
myCriteria[2] = new FEmTool_LinearJerk (MxDeg, ConstraintOrder);
Handle(TColStd_HArray2OfReal) Coeff = new TColStd_HArray2OfReal(0, 0, 1, NbDim);
myCriteria[0]->Set(Coeff);
myCriteria[1]->Set(Coeff);
myCriteria[2]->Set(Coeff);
}
else if (myCurve != C) {
Standard_Integer OldMxDeg = myCurve->Base()->WorkDegree(),
OldNbDim = myCurve->Dimension(),
OldOrder = order(myCurve->Base());
myCurve = C;
Standard_Integer MxDeg = myCurve->Base()->WorkDegree(),
NbDim = myCurve->Dimension(),
Order = order(myCurve->Base());
if(MxDeg != OldMxDeg || Order != OldOrder) {
GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
switch (Order) {
case 0 : ConstraintOrder = GeomAbs_C0;
break;
case 1 : ConstraintOrder = GeomAbs_C1;
break;
case 2 : ConstraintOrder = GeomAbs_C2;
}
myCriteria[0] = new FEmTool_LinearTension(MxDeg, ConstraintOrder);
myCriteria[1] = new FEmTool_LinearFlexion(MxDeg, ConstraintOrder);
myCriteria[2] = new FEmTool_LinearJerk (MxDeg, ConstraintOrder);
Handle(TColStd_HArray2OfReal) Coeff = new TColStd_HArray2OfReal(0, 0, 1, NbDim);
myCriteria[0]->Set(Coeff);
myCriteria[1]->Set(Coeff);
myCriteria[2]->Set(Coeff);
}
else if(NbDim != OldNbDim) {
Handle(TColStd_HArray2OfReal) Coeff = new TColStd_HArray2OfReal(0, 0, 1, NbDim);
myCriteria[0]->Set(Coeff);
myCriteria[1]->Set(Coeff);
myCriteria[2]->Set(Coeff);
}
}
}
//=======================================================================
//function : GetCurve
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::GetCurve(Handle(FEmTool_Curve)& C) const
{
C = myCurve;
}
//=======================================================================
//function : SetEstimation
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::SetEstimation(const Standard_Real E1,
const Standard_Real E2,
const Standard_Real E3)
{
myEstimation[0] = E1;
myEstimation[1] = E2;
myEstimation[2] = E3;
}
Standard_Real& AppParCurves_LinearCriteria::EstLength()
{
return myLength;
}
//=======================================================================
//function : GetEstimation
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::GetEstimation(Standard_Real& E1,
Standard_Real& E2,
Standard_Real& E3) const
{
E1 = myEstimation[0];
E2 = myEstimation[1];
E3 = myEstimation[2];
}
//=======================================================================
//function : AssemblyTable
//purpose :
//=======================================================================
Handle(FEmTool_HAssemblyTable) AppParCurves_LinearCriteria::AssemblyTable() const
{
if(myCurve.IsNull()) Standard_DomainError::Raise("AppParCurves_LinearCriteria::AssemblyTable");
Standard_Integer NbDim = myCurve->Dimension(),
NbElm = myCurve->NbElements(),
nc1 = order(myCurve->Base()) + 1;
Standard_Integer MxDeg = myCurve->Base()->WorkDegree() ;
Handle(FEmTool_HAssemblyTable) AssTable = new FEmTool_HAssemblyTable(1, NbDim, 1, NbElm);
Handle(TColStd_HArray1OfInteger) GlobIndex, Aux;
Standard_Integer i, el = 1, dim = 1, NbGlobVar = 0, gi0;
// For dim = 1
// For first element (el = 1)
GlobIndex = new TColStd_HArray1OfInteger(0, MxDeg);
for(i = 0; i < nc1; i++) {
NbGlobVar++;
GlobIndex->SetValue(i, NbGlobVar);
}
gi0 = MxDeg - 2 * nc1 + 1;
for(i = nc1; i < 2*nc1; i++) {
NbGlobVar++;
GlobIndex->SetValue(i, NbGlobVar + gi0);
}
for(i = 2*nc1; i <= MxDeg; i++) {
NbGlobVar++;
GlobIndex->SetValue(i, NbGlobVar - nc1);
}
gi0 = NbGlobVar - nc1 + 1;
AssTable->SetValue(dim, el, GlobIndex);
// For rest elements
for(el = 2; el <= NbElm; el++) {
GlobIndex = new TColStd_HArray1OfInteger(0, MxDeg);
for(i = 0; i < nc1; i++) GlobIndex->SetValue(i, gi0 + i);
gi0 = MxDeg - 2 * nc1 + 1;
for(i = nc1; i < 2*nc1; i++) {
NbGlobVar++;
GlobIndex->SetValue(i, NbGlobVar + gi0);
}
for(i = 2*nc1; i <= MxDeg; i++) {
NbGlobVar++;
GlobIndex->SetValue(i, NbGlobVar - nc1);
}
gi0 = NbGlobVar - nc1 + 1;
AssTable->SetValue(dim, el, GlobIndex);
}
// For other dimensions
gi0 = NbGlobVar;
for(dim = 2; dim <= NbDim; dim++) {
for(el = 1; el <= NbElm; el++) {
Aux = AssTable->Value(1, el);
GlobIndex = new TColStd_HArray1OfInteger(0, MxDeg);
for(i = 0; i <= MxDeg; i++) GlobIndex->SetValue(i, Aux->Value(i) + NbGlobVar);
AssTable->SetValue(dim, el, GlobIndex);
}
NbGlobVar += gi0;
}
return AssTable;
}
//=======================================================================
//function :
//purpose :
//=======================================================================
Handle(TColStd_HArray2OfInteger) AppParCurves_LinearCriteria::DependenceTable() const
{
if(myCurve.IsNull()) Standard_DomainError::Raise("AppParCurves_LinearCriteria::DependenceTable");
Standard_Integer Dim = myCurve->Dimension();
Handle(TColStd_HArray2OfInteger) DepTab =
new TColStd_HArray2OfInteger(1, Dim, 1, Dim, 0);
Standard_Integer i;
for(i=1; i <= Dim; i++) DepTab->SetValue(i,i,1);
return DepTab;
}
//=======================================================================
//function : QualityValues
//purpose :
//=======================================================================
Standard_Integer AppParCurves_LinearCriteria::QualityValues(const Standard_Real J1min,
const Standard_Real J2min,
const Standard_Real J3min,
Standard_Real& J1,
Standard_Real& J2,
Standard_Real& J3)
{
if(myCurve.IsNull()) Standard_DomainError::Raise("AppParCurves_LinearCriteria::QualityValues");
Standard_Integer NbDim = myCurve->Dimension(),
NbElm = myCurve->NbElements();
TColStd_Array1OfReal& Knots = myCurve->Knots();
Handle(TColStd_HArray2OfReal) Coeff;
Standard_Integer el, deg = 0, curdeg, i;
Standard_Real UFirst, ULast;
J1 = J2 = J3 = 0.;
for(el = 1; el <= NbElm; el++) {
curdeg = myCurve->Degree(el);
if(deg != curdeg) {
deg = curdeg;
Coeff = new TColStd_HArray2OfReal(0, deg, 1, NbDim);
}
myCurve->GetElement(el, Coeff->ChangeArray2());
UFirst = Knots(el); ULast = Knots(el + 1);
myCriteria[0]->Set(Coeff);
myCriteria[0]->Set(UFirst, ULast);
J1 = J1 + myCriteria[0]->Value();
myCriteria[1]->Set(Coeff);
myCriteria[1]->Set(UFirst, ULast);
J2 = J2 + myCriteria[1]->Value();
myCriteria[2]->Set(Coeff);
myCriteria[2]->Set(UFirst, ULast);
J3 = J3 + myCriteria[2]->Value();
}
// Calculation of ICDANA - see MOTEST.f
// Standard_Real JEsMin[3] = {.01, .001, .001}; // from MOTLIS.f
Standard_Real JEsMin[3]; JEsMin[0] = J1min; JEsMin[1] = J2min; JEsMin[2] = J3min;
Standard_Real ValCri[3]; ValCri[0] = J1; ValCri[1] = J2; ValCri[2] = J3;
Standard_Integer ICDANA = 0;
// (2) Test l'amelioration des estimations
// (critere sureleve => Non minimisation )
for(i = 0; i <= 2; i++)
if((ValCri[i] < 0.8 * myEstimation[i]) && (myEstimation[i] > JEsMin[i])) {
if(ICDANA < 1) ICDANA = 1;
if(ValCri[i] < 0.1 * myEstimation[i]) ICDANA = 2;
myEstimation[i] = Max(1.05*ValCri[i], JEsMin[i]);
}
// (3) Mise a jours des Estimation
// (critere sous-estimer => mauvais conditionement)
if (ValCri[0] > myEstimation[0] * 2) {
myEstimation[0] += ValCri[0] * .1;
if (ICDANA == 0) {
if (ValCri[0] > myEstimation[0] * 10) {
ICDANA = 2;
}
else ICDANA = 1;
}
else {
ICDANA = 2;
}
}
if (ValCri[1] > myEstimation[1] * 20) {
myEstimation[1] += ValCri[1] * .1;
if (ICDANA == 0) {
if (ValCri[1] > myEstimation[1] * 100) {
ICDANA = 2;
}
else ICDANA = 1;
}
else {
ICDANA = 2;
}
}
if (ValCri[2] > myEstimation[2] * 20) {
myEstimation[2] += ValCri[2] * .05;
if (ICDANA == 0) {
if (ValCri[2] > myEstimation[2] * 100) {
ICDANA = 2;
}
else ICDANA = 1;
}
else {
ICDANA = 2;
}
}
return ICDANA;
}
//=======================================================================
//function : ErrorValues
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::ErrorValues(Standard_Real& MaxError,
Standard_Real& QuadraticError,
Standard_Real& AverageError)
{
if(myCurve.IsNull()) Standard_DomainError::Raise("AppParCurves_LinearCriteria::ErrorValues");
Standard_Integer NbDim = myCurve->Dimension();
Standard_Integer myNbP2d = ToolLine::NbP2d(mySSP), myNbP3d = ToolLine::NbP3d(mySSP);
if(NbDim != (2*myNbP2d + 3*myNbP3d)) Standard_DomainError::Raise("AppParCurves_LinearCriteria::ErrorValues");
TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
TColStd_Array1OfReal BasePoint(1,NbDim);
gp_Pnt2d P2d;
gp_Pnt P3d;
Standard_Integer i, ipnt, c0 = 0;
Standard_Real SqrDist, Dist;
MaxError = QuadraticError = AverageError = 0.;
for(i = myParameters->Lower(); i <= myParameters->Upper(); i++) {
myCurve->D0(myParameters->Value(i), BasePoint);
c0 = 0;
ToolLine::Value(mySSP, i, TabP3d);
for(ipnt = 1; ipnt <= myNbP3d; ipnt++) {
P3d.SetCoord(BasePoint(c0+1), BasePoint(c0+2), BasePoint(c0+3));
SqrDist = P3d.SquareDistance(TabP3d(ipnt)); Dist = Sqrt(SqrDist);
MaxError = Max(MaxError, Dist);
QuadraticError += SqrDist;
AverageError += Dist;
c0 += 3;
}
if(myNbP3d == 0) ToolLine::Value(mySSP, i, TabP2d);
else ToolLine::Value(mySSP, i, TabP3d, TabP2d);
for(ipnt = 1; ipnt <= myNbP2d; ipnt++) {
P2d.SetCoord(BasePoint(c0+1), BasePoint(c0+2));
SqrDist = P2d.SquareDistance(TabP2d(ipnt)); Dist = Sqrt(SqrDist);
MaxError = Max(MaxError, Dist);
QuadraticError += SqrDist;
AverageError += Dist;
c0 += 2;
}
}
}
//=======================================================================
//function : Hessian
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::Hessian(const Standard_Integer Element,
const Standard_Integer Dimension1,
const Standard_Integer Dimension2,
math_Matrix& H)
{
if(myCurve.IsNull()) Standard_DomainError::Raise("AppParCurves_LinearCriteria::Hessian");
if(DependenceTable()->Value(Dimension1, Dimension2) == 0)
Standard_DomainError::Raise("AppParCurves_LinearCriteria::Hessian");
Standard_Integer //NbDim = myCurve->Dimension(),
MxDeg = myCurve->Base()->WorkDegree(),
// Deg = myCurve->Degree(Element),
Order = order(myCurve->Base());
math_Matrix AuxH(0, H.RowNumber()-1, 0, H.ColNumber()-1, 0.);
TColStd_Array1OfReal& Knots = myCurve->Knots();
Standard_Real UFirst, ULast;
UFirst = Knots(Element); ULast = Knots(Element + 1);
Standard_Integer icrit;
// Quality criterion part of Hessian
H.Init(0);
for(icrit = 0; icrit <= 2; icrit++) {
myCriteria[icrit]->Set(UFirst, ULast);
myCriteria[icrit]->Hessian(Dimension1, Dimension2, AuxH);
H += (myQualityWeight*myPercent[icrit]/myEstimation[icrit]) * AuxH;
}
// Least square part of Hessian
AuxH.Init(0.);
Standard_Real coeff = (ULast - UFirst)/2., curcoeff, poid;
Standard_Integer ipnt, ii, degH = 2 * Order+1;
Handle(PLib_Base) myBase = myCurve->Base();
Standard_Integer k1, k2, i, j, i0 = H.LowerRow(), j0 = H.LowerCol(), i1, j1,
di = myPntWeight.Lower() - myParameters->Lower();
//BuilCache
if (myE != Element) BuildCache(Element);
// Compute the least square Hessian
for(ii=1, ipnt = IF; ipnt <= IL; ipnt++, ii+=(MxDeg+1)) {
poid = myPntWeight(di + ipnt) * 2.;
const Standard_Real * BV = &myCache->Value(ii);
// Hermite*Hermite part of matrix
for(i = 0; i <= degH; i++) {
k1 = (i <= Order)? i : i - Order - 1;
curcoeff = Pow(coeff, k1) * poid * BV[i];
// Hermite*Hermite part of matrix
for(j = i; j <= degH; j++) {
k2 = (j <= Order)? j : j - Order - 1;
AuxH(i, j) += curcoeff * Pow(coeff, k2) * BV[j];
}
// Hermite*Jacobi part of matrix
for(j = degH + 1; j <= MxDeg; j++) {
AuxH(i, j) += curcoeff * BV[j];
}
}
// Jacoby*Jacobi part of matrix
for(i = degH+1; i <= MxDeg; i++) {
curcoeff = BV[i] * poid;
for(j = i; j <= MxDeg; j++) {
AuxH(i, j) += curcoeff * BV[j];
}
}
}
i1 = i0;
for(i = 0; i <= MxDeg; i++) {
j1 = j0 + i;
for(j = i; j <= MxDeg; j++) {
H(i1, j1) += myQuadraticWeight * AuxH(i, j);
H(j1, i1) = H(i1, j1);
j1++;
}
i1++;
}
}
//=======================================================================
//function : Gradient
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::Gradient(const Standard_Integer Element,
const Standard_Integer Dimension,
math_Vector& G)
{
if(myCurve.IsNull())
Standard_DomainError::Raise("AppParCurves_LinearCriteria::ErrorValues");
Standard_Integer myNbP2d = ToolLine::NbP2d(mySSP), myNbP3d = ToolLine::NbP3d(mySSP);
if(Dimension > (2*myNbP2d + 3*myNbP3d))
Standard_DomainError::Raise("AppParCurves_LinearCriteria::ErrorValues");
TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
Standard_Boolean In3d;
Standard_Integer IndPnt, IndCrd;
if(Dimension <= 3*myNbP3d) {
In3d = Standard_True;
IndCrd = Dimension % 3;
IndPnt = Dimension / 3;
if(IndCrd == 0) IndCrd = 3;
else IndPnt++;
}
else {
In3d = Standard_False;
IndCrd = (Dimension - 3*myNbP3d) % 2;
IndPnt = (Dimension - 3*myNbP3d) / 2;
if(IndCrd == 0) IndCrd = 2;
else IndPnt++;
}
TColStd_Array1OfReal& Knots = myCurve->Knots();
Standard_Real UFirst, ULast, Pnt;
UFirst = Knots(Element); ULast = Knots(Element + 1);
Standard_Real coeff = (ULast-UFirst)/2;
Standard_Integer //Deg = myCurve->Degree(Element),
Order = order(myCurve->Base());
Handle(PLib_Base) myBase = myCurve->Base();
Standard_Integer MxDeg = myBase->WorkDegree();
Standard_Real curcoeff;
Standard_Integer degH = 2 * Order + 1;
Standard_Integer ipnt, k, i, ii, i0 = G.Lower(),
di = myPntWeight.Lower() - myParameters->Lower();
if (myE != Element) BuildCache(Element);
const Standard_Real * BV = &myCache->Value(1);
BV--;
G.Init(0.);
for(ii=1,ipnt = IF; ipnt <= IL; ipnt++) {
if(In3d) {
ToolLine::Value(mySSP, ipnt, TabP3d);
Pnt = TabP3d(IndPnt).Coord(IndCrd);
}
else {
if(myNbP3d == 0) ToolLine::Value(mySSP, ipnt, TabP2d);
else ToolLine::Value(mySSP, ipnt, TabP3d, TabP2d);
Pnt = TabP2d(IndPnt).Coord(IndCrd);
}
curcoeff = Pnt * myPntWeight(di + ipnt);
for(i = 0; i <= MxDeg; i++,ii++)
G(i0 + i) += BV[ii] * curcoeff;
}
G *= 2. * myQuadraticWeight;
for(i = 0; i <= degH; i++) {
k = (i <= Order)? i : i - Order - 1;
curcoeff = Pow(coeff, k);
G(i0 + i) *= curcoeff;
}
}
//=======================================================================
//function : InputVector
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::InputVector(const math_Vector& X,
const Handle(FEmTool_HAssemblyTable)& AssTable)
{
Standard_Integer NbDim = myCurve->Dimension(),
NbElm = myCurve->NbElements() ;
Standard_Integer MxDeg = 0 ;
MxDeg = myCurve->Base()->WorkDegree();
TColStd_Array2OfReal CoeffEl(0, MxDeg, 1, NbDim);
Handle(TColStd_HArray1OfInteger) GlobIndex;
Standard_Integer el, dim, i, i0 = X.Lower() - 1;
for(el = 1; el <= NbElm; el++) {
for(dim = 1; dim <= NbDim; dim++) {
GlobIndex = AssTable->Value(dim, el);
for(i = 0; i <= MxDeg; i++) CoeffEl(i, dim) = X(i0 + GlobIndex->Value(i));
}
myCurve->SetDegree(el, MxDeg);
myCurve->SetElement(el, CoeffEl);
}
}
//=======================================================================
//function : SetWeight
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::SetWeight(const Standard_Real QuadraticWeight,
const Standard_Real QualityWeight,
const Standard_Real percentJ1,
const Standard_Real percentJ2,
const Standard_Real percentJ3)
{
if (QuadraticWeight < 0. || QualityWeight < 0.)
Standard_DomainError::Raise("AppParCurves_LinearCriteria::SetWeight");
if (percentJ1 < 0. || percentJ2 < 0. || percentJ3 < 0.)
Standard_DomainError::Raise("AppParCurves_LinearCriteria::SetWeight");
myQuadraticWeight = QuadraticWeight; myQualityWeight = QualityWeight;
Standard_Real Total = percentJ1 + percentJ2 + percentJ3;
myPercent[0] = percentJ1 / Total;
myPercent[1] = percentJ2 / Total;
myPercent[2] = percentJ3 / Total;
}
//=======================================================================
//function : GetWeight
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::GetWeight(Standard_Real& QuadraticWeight,
Standard_Real& QualityWeight) const
{
QuadraticWeight = myQuadraticWeight; QualityWeight = myQualityWeight;
}
//=======================================================================
//function : SetWeight
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::SetWeight(const TColStd_Array1OfReal& Weight)
{
myPntWeight = Weight;
}
//=======================================================================
//function : BuildCache
//purpose :
//=======================================================================
void AppParCurves_LinearCriteria::BuildCache(const Standard_Integer Element)
{
Standard_Real t;
Standard_Real UFirst, ULast;
Standard_Integer ipnt;
UFirst = myCurve->Knots()(Element);
ULast = myCurve->Knots()(Element + 1);
IF = 0;
for(ipnt = myParameters->Lower(); ipnt <= myParameters->Upper(); ipnt++) {
t = myParameters->Value(ipnt);
if((t > UFirst && t <= ULast) || (Element == 1 && t == UFirst)) {
if (IF == 0) IF=ipnt;
IL = ipnt;
}
else if (t>ULast) break;
}
if (IF != 0) {
Handle(PLib_Base) myBase = myCurve->Base();
Standard_Integer order = myBase->WorkDegree()+1, ii;
myCache = new TColStd_HArray1OfReal (1, (IL-IF+1)*(order));
ii =1;
for(ipnt = IF, ii=1; ipnt <= IL; ipnt++, ii+=order) {
Standard_Real * cache = &myCache->ChangeValue(ii);
TColStd_Array1OfReal BasicValue(cache[0], 0, order-1);
t = myParameters->Value(ipnt);
Standard_Real coeff = 2./(ULast - UFirst), c0 = -(ULast + UFirst)/2., s;
s = (t + c0) * coeff;
myBase->D0(s, BasicValue);
}
}
else { //pas de points dans l'interval.
IF = IL;
IL--;
}
myE = Element;
}

119
src/AppParCurves/AppParCurves_SmoothCriterion.cdl
View File

@@ -1,119 +0,0 @@
-- Created on: 1997-09-11
-- Created by: Philippe MANGIN
-- Copyright (c) 1997-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.
deferred class SmoothCriterion from AppParCurves
inherits TShared from MMgt
---Purpose: defined criterion to smooth points in curve
uses
Vector from math,
Matrix from math,
Curve from FEmTool,
HAssemblyTable from FEmTool,
HArray2OfInteger from TColStd,
HArray1OfReal from TColStd,
Array1OfReal from TColStd
raises
NotImplemented,
DomainError
is
SetParameters(me : mutable; Parameters : HArray1OfReal)
is deferred;
SetCurve(me : mutable; C :Curve from FEmTool)
is deferred;
GetCurve(me; C : out Curve from FEmTool)
is deferred;
SetEstimation(me : mutable; E1, E2, E3 : Real)
is deferred;
EstLength(me : mutable)
---C++: return &
returns Real is deferred;
GetEstimation(me; E1, E2, E3 : out Real)
is deferred;
AssemblyTable(me)
returns HAssemblyTable from FEmTool
is deferred;
DependenceTable(me)
returns HArray2OfInteger from TColStd
is deferred;
QualityValues (me : mutable; J1min, J2min, J3min : Real;
J1, J2, J3 : out Real)
returns Integer is deferred;
ErrorValues(me : mutable;
MaxError, QuadraticError, AverageError : out Real)
is deferred;
Hessian(me : mutable ;
Element : Integer;
Dimension1 : Integer;
Dimension2 : Integer;
H : out Matrix from math)
raises DomainError -- If DependenceTable(Dimension1,Dimension2) is False
is deferred;
Gradient(me : mutable;
Element : Integer;
Dimension : Integer;
G : out Vector from math)
is deferred;
InputVector(me : mutable; X : Vector from math;
AssTable : HAssemblyTable from FEmTool)
---Purpose: Convert the assembly Vector in an Curve;
--
raises DomainError is deferred;
SetWeight(me: mutable;
QuadraticWeight, QualityWeight : Real;
percentJ1, percentJ2, percentJ3 : Real)
is deferred;
GetWeight(me; QuadraticWeight, QualityWeight : out Real)
is deferred;
SetWeight(me: mutable;
Weight : Array1OfReal)
is deferred;
end SmoothCriterion;

17
src/AppParCurves/AppParCurves_SmoothCriterion.cxx
View File

@@ -1,17 +0,0 @@
// Created on: 1998-11-30
// Created by: Igor FEOKTISTOV
// Copyright (c) 1998-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 <AppParCurves_SmoothCriterion.ixx>

443
src/AppParCurves/AppParCurves_Variational.cdl
View File

@@ -1,443 +0,0 @@
-- Created on: 1996-05-14
-- Created by: Philippe MANGIN / Jeannine PANCIATICI
-- 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.
-- Igor FEOKTISTOV - correction 14/12/98
generic class Variational from AppParCurves
(MultiLine as any;
ToolLine as any) -- as ToolLine(MultiLine)
---Purpose: This class is used to smooth N points with constraints
-- by minimization of quadratic criterium but also
-- variational criterium in order to obtain " fair Curve "
--
uses Matrix from math,
Vector from math,
HArray1OfReal from TColStd,
Array1OfReal from TColStd,
HArray1OfInteger from TColStd,
Shape from GeomAbs,
HArray1OfConstraintCouple from AppParCurves,
MultiBSpCurve from AppParCurves,
SmoothCriterion from AppParCurves,
Curve from FEmTool,
Assembly from FEmTool,
Base from PLib,
Constraint from AppParCurves
raises OutOfRange from Standard,
DimensionError from Standard,
DomainError from Standard,
ConstructionError from Standard,
NotDone from StdFail,
VectorWithNullMagnitude from gp
class MyCriterion instantiates LinearCriteria from AppParCurves
(MultiLine, ToolLine);
is
Create(SSP: MultiLine;
FirstPoint, LastPoint: Integer;
TheConstraints: HArray1OfConstraintCouple;
MaxDegree: Integer = 14;
MaxSegment: Integer = 100;
Continuity : Shape from GeomAbs = GeomAbs_C2;
WithMinMax : Boolean = Standard_False;
WithCutting: Boolean = Standard_True;
Tolerance : Real = 1.0;
NbIterations: Integer = 2)
---Purpose: Constructor.
-- Initialization of the fields.
-- warning : Nc0 : number of PassagePoint consraints
-- Nc2 : number of TangencyPoint constraints
-- Nc3 : number of CurvaturePoint constraints
-- if
-- ((MaxDegree-Continuity)*MaxSegment -Nc0 - 2*Nc1
-- -3*Nc2)
-- is negative
-- The problem is over-constrained.
--
-- Limitation : The MultiLine has to be composed by
-- only one Line ( Dimension 2 or 3).
returns Variational from AppParCurves;
Approximate(me : in out)
---Purpose: Makes the approximation with the current fields.
raises NotDone from StdFail
is static;
-- ==================== The Selectors ===========================
IsCreated(me)
---Purpose: returns True if the creation is done
-- and correspond to the current fields.
returns Boolean
is static;
IsDone(me)
---Purpose: returns True if the approximation is ok
-- and correspond to the current fields.
returns Boolean
is static;
IsOverConstrained(me)
---Purpose: returns True if the problem is overconstrained
-- in this case, approximation cannot be done.
returns Boolean
is static;
Value(me)
---Purpose: returns all the BSpline curves approximating the
-- MultiLine SSP after minimization of the parameter.
returns MultiBSpCurve from AppParCurves
raises NotDone from StdFail
is static;
MaxError(me)
---Purpose: returns the maximum of the distances between
-- the points of the multiline and the approximation
-- curves.
returns Real
raises NotDone from StdFail
is static;
MaxErrorIndex(me)
---Purpose: returns the index of the MultiPoint of ErrorMax
returns Integer
raises NotDone from StdFail
is static;
QuadraticError(me)
---Purpose: returns the quadratic average of the distances between
-- the points of the multiline and the approximation
-- curves.
returns Real
raises NotDone from StdFail
is static;
Distance(me : in out ; mat : out Matrix from math)
---Purpose: returns the distances between the points of the
-- multiline and the approximation curves.
raises NotDone from StdFail
is static;
AverageError(me)
---Purpose: returns the average error between
-- the MultiLine and the approximation.
returns Real
raises NotDone from StdFail
is static;
Parameters(me)
---Purpose: returns the parameters uses to the approximations
---C++: return const&
returns HArray1OfReal
raises NotDone from StdFail
is static;
Knots(me)
---Purpose: returns the knots uses to the approximations
---C++: return const&
returns HArray1OfReal
raises NotDone from StdFail
is static;
Criterium(me; VFirstOrder,
VSecondOrder,
VThirdOrder : out Real)
---Purpose: returns the values of the quality criterium.
raises NotDone from StdFail
is static;
CriteriumWeight(me ;
Percent1, Percent2, Percent3 : out Real)
---Purpose: returns the Weights (as percent) associed to the criterium used in
-- the optimization.
is static;
MaxDegree(me)
---Purpose: returns the Maximum Degree used in the approximation
returns Integer
is static;
MaxSegment(me)
---Purpose: returns the Maximum of segment used in the approximation
returns Integer
is static;
Continuity(me)
---Purpose: returns the Continuity used in the approximation
returns Shape from GeomAbs
is static;
WithMinMax(me)
---Purpose: returns if the approximation search to minimize the
-- maximum Error or not.
returns Boolean
is static;
WithCutting(me)
---Purpose: returns if the approximation can insert new Knots or not.
returns Boolean
is static;
Tolerance(me)
---Purpose: returns the tolerance used in the approximation.
returns Real
is static;
NbIterations(me)
---Purpose: returns the number of iterations used in the approximation.
returns Integer
is static;
Dump(me ; o : in out OStream)
---Purpose: Prints on the stream o information on the current state
-- of the object.
-- MaxError,MaxErrorIndex,AverageError,QuadraticError,Criterium
-- Distances,Degre,Nombre de poles, parametres, noeuds
is static;
SetConstraints(me:in out; aConstrainst:HArray1OfConstraintCouple)
---Purpose: Define the constraints to approximate
-- If this value is incompatible with the others fields
-- this method modify nothing and returns false
returns Boolean
is static;
SetParameters(me:in out; param : HArray1OfReal)
---Purpose: Defines the parameters used by the approximations.
raises DimensionError
is static;
SetKnots(me:in out; knots : HArray1OfReal)
---Purpose: Defines the knots used by the approximations
-- If this value is incompatible with the others fields
-- this method modify nothing and returns false
returns Boolean
raises DimensionError,
DomainError
is static;
SetMaxDegree(me: in out; Degree : Integer)
---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
returns Boolean
is static;
SetMaxSegment(me: in out; NbSegment : Integer)
---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
returns Boolean
is static;
SetContinuity(me: in out; C : Shape from GeomAbs)
---Purpose: Define the Continuity used in the approximation
-- If this value is incompatible with the others fields
-- this method modify nothing and returns false
returns Boolean
raises ConstructionError from Standard
is static;
SetWithMinMax(me: in out; MinMax : Boolean)
---Purpose: Define if the approximation search to minimize the
-- maximum Error or not.
is static;
SetWithCutting(me : in out; Cutting : Boolean )
---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
returns Boolean
is static;
SetCriteriumWeight(me : in out;
Percent1, Percent2, Percent3 : Real)
---Purpose: define the Weights (as percent) associed to the criterium used in
-- the optimization.
--
raises DomainError -- if Percent <= 0
is static;
SetCriteriumWeight(me : in out;
Order : Integer;
Percent : Real)
---Purpose: define the Weight (as percent) associed to the
-- criterium Order used in the optimization : Others
-- weights are updated.
raises DomainError, -- if Percent < 0
OutOfRange -- if Order < 1 or Order > 3
is static;
SetTolerance(me:in out; Tol : Real)
---Purpose: define the tolerance used in the approximation.
is static;
SetNbIterations(me:in out; Iter : Integer)
---Purpose: define the number of iterations used in the approximation.
raises DomainError -- if Iter < 1
is static;
-- ====================== The Private methods ======================
TheMotor(me : in out;
J : in out SmoothCriterion from AppParCurves;
WQuadratic, WQuality : Real;
TheCurve : in out Curve from FEmTool;
Ecarts : out Array1OfReal from TColStd) is private;
Adjusting(me : in out;
J : in out SmoothCriterion from AppParCurves;
WQuadratic, WQuality : in out Real;
TheCurve : in out Curve from FEmTool;
Ecarts : out Array1OfReal from TColStd) is private;
Optimization(me;
J : in out SmoothCriterion from AppParCurves;
A : in out Assembly from FEmTool;
ToAssemble : in Boolean;
EpsDeg : Real;
Curve : out Curve from FEmTool;
Parameters : Array1OfReal from TColStd) is private;
Project(me; C : Curve from FEmTool;
Ti : Array1OfReal from TColStd;
ProjTi : out Array1OfReal from TColStd;
Distance : out Array1OfReal from TColStd;
NumPoints : out Integer;
MaxErr, QuaErr, AveErr : out Real;
NbIterations: Integer=2) is private;
ACR(me; Curve : in out Curve from FEmTool;
Ti : in out Array1OfReal from TColStd;
Decima: Integer) is private;
SplitCurve(me; InCurve : Curve from FEmTool;
Ti : Array1OfReal from TColStd;
CurveTol: Real;
OutCurve: out Curve from FEmTool;
iscut : out Boolean) is private;
Init(me : in out)
raises NotDone from StdFail,
ConstructionError from Standard,
DimensionError from Standard
is private;
InitSmoothCriterion(me : in out)
is private;
InitParameters(me : in out; Length : out Real)
raises ConstructionError from Standard
is private;
InitCriterionEstimations(me; Length : Real; J1, J2, J3 : out Real)
is private;
EstTangent(me; ipnt : Integer; VTang : out Vector from math)
is private;
EstSecnd(me; ipnt : Integer; VTang1, VTang2 : Vector from math;
Length : Real; VScnd : out Vector from math)
is private;
InitCutting(me; aBase : Base from PLib; CurvTol : Real;
aCurve : out Curve from FEmTool)
raises ConstructionError from Standard
is private;
AssemblingConstraints(me; Curve : Curve from FEmTool;
Parameters : Array1OfReal from TColStd;
CBLONG : Real from Standard;
A : out Assembly from FEmTool)
is private;
InitTthetaF(me : in out; ndimen : Integer from Standard;
typcon : Constraint from AppParCurves;
begin : Integer from Standard;
jndex : Integer from Standard)
returns Boolean
is private;
fields
-- Description of the points to smooth and the constraints
mySSP : MultiLine;
myNbP3d : Integer;
myNbP2d : Integer;
myDimension : Integer;
myFirstPoint : Integer;
myLastPoint : Integer;
myNbPoints : Integer;
myTabPoints : HArray1OfReal;
myConstraints : HArray1OfConstraintCouple;
myNbConstraints : Integer;
myTabConstraints : HArray1OfReal;
myNbPassPoints : Integer;
myNbTangPoints : Integer;
myNbCurvPoints : Integer;
myTypConstraints : HArray1OfInteger;
myTtheta : HArray1OfReal;
myTfthet : HArray1OfReal;
-- Context parameters
myMaxDegree : Integer;
myMaxSegment : Integer;
myNbIterations: Integer;
myTolerance : Real;
-- Options
myContinuity : Shape from GeomAbs;
myNivCont : Integer;
myWithMinMax : Boolean;
myWithCutting: Boolean;
myPercent : Real[3];
myCriterium : Real[4];
mySmoothCriterion : SmoothCriterion from AppParCurves;
-- output
myParameters : HArray1OfReal;
myKnots : HArray1OfReal;
myMBSpCurve : MultiBSpCurve;
myMaxError : Real;
myMaxErrorIndex: Integer;
myAverageError : Real;
myIsCreated : Boolean;
myIsDone : Boolean;
myIsOverConstr : Boolean;
end Variational;

1133
src/AppParCurves/AppParCurves_Variational.gxx
View File

File diff suppressed because it is too large Load Diff

362
src/AppParCurves/AppParCurves_Variational_1.gxx
View File

@@ -1,362 +0,0 @@
// Created on: 1997-09-17
// Created by: Philippe MANGIN /Igor Feoktistov (1998)
// Copyright (c) 1997-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 <SortTools_StraightInsertionSortOfReal.hxx>
#include <TCollection_CompareOfReal.hxx>
#include <TColStd_HArray2OfInteger.hxx>
#include <TColStd_HArray1OfReal.hxx>
#include <TColStd_Array2OfInteger.hxx>
#include <FEmTool_Assembly.hxx>
#include <FEmTool_AssemblyTable.hxx>
#include <FEmTool_Curve.hxx>
//====================== Private Methodes =============================//
//=======================================================================
//function : TheMotor
//purpose : Smoothing's motor like STRIM routine "MOTLIS"
//=======================================================================
void AppParCurves_Variational::TheMotor(
Handle(AppParCurves_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;
// SortTools_StraightInsertionSortOfReal Sort;
TCollection_CompareOfReal CompReal;
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.
// Sort.Sort(CurrentTi->ChangeArray1(), CompReal);
SortTools_StraightInsertionSortOfReal::Sort(CurrentTi->ChangeArray1(), CompReal);
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.
// Sort.Sort(CurrentTi->ChangeArray1(), CompReal);
SortTools_StraightInsertionSortOfReal::Sort(CurrentTi->ChangeArray1(), CompReal);
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.
// Sort.Sort(CurrentTi->ChangeArray1(), CompReal);
SortTools_StraightInsertionSortOfReal::Sort(CurrentTi->ChangeArray1(), CompReal);
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;
}

113
src/AppParCurves/AppParCurves_Variational_2.gxx
View File

@@ -1,113 +0,0 @@
// Created on: 1998-12-07
// Created by: Igor FEOKTISTOV
// Copyright (c) 1998-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 <math_Matrix.hxx>
#include <math_Vector.hxx>
#include <TColStd_Array2OfReal.hxx>
#include <PLib_Base.hxx>
#include <PLib_JacobiPolynomial.hxx>
#include <PLib_HermitJacobi.hxx>
#include <FEmTool_HAssemblyTable.hxx>
//=======================================================================
//function : Optimization
//purpose : (like FORTRAN subroutine MINIMI)
//=======================================================================
void AppParCurves_Variational::Optimization(Handle(AppParCurves_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);
}
}
}

132
src/AppParCurves/AppParCurves_Variational_3.gxx
View File

@@ -1,132 +0,0 @@
// Created on: 1998-12-08
// Created by: Igor FEOKTISTOV
// Copyright (c) 1998-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.
void AppParCurves_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]
}

133
src/AppParCurves/AppParCurves_Variational_4.gxx
View File

@@ -1,133 +0,0 @@
// Created on: 1998-12-08
// Created by: Igor FEOKTISTOV
// Copyright (c) 1998-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.
void AppParCurves_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.;
}

152
src/AppParCurves/AppParCurves_Variational_5.gxx
View File

@@ -1,152 +0,0 @@
// Created on: 1998-12-10
// Created by: Igor FEOKTISTOV
// Copyright (c) 1998-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 <SortTools_ShellSortOfReal.hxx>
#include <TCollection_CompareOfReal.hxx>
#include <PLib_Base.hxx>
//----------------------------------------------------------//
// 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 AppParCurves_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 DEB
Standard_Integer MaxDegree =
#endif
InCurve->Base()->WorkDegree();
Standard_Integer NbElm = NbElmOld;
TColStd_Array1OfReal NewKnots(NbElm + 1, myMaxSegment);
#ifndef DEB
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);
// SortTools_ShellSortOfReal Sort;
TCollection_CompareOfReal CompReal;
// Sort.Sort(OutKnots, CompReal);
SortTools_ShellSortOfReal::Sort(OutKnots, CompReal);
}
else
iscut = Standard_False;
}

622
src/AppParCurves/AppParCurves_Variational_6.gxx
View File

@@ -1,622 +0,0 @@
// Created on: 1998-12-21
// Created by: Igor FEOKTISTOV
// Copyright (c) 1998-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.
// Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559
#include <PLib_Base.hxx>
#include <PLib_JacobiPolynomial.hxx>
#include <PLib_HermitJacobi.hxx>
#include <FEmTool_Curve.hxx>
#include <math_Vector.hxx>
#include <TColStd_Array1OfReal.hxx>
//=======================================================================
//function : InitSmoothCriterion
//purpose : Initializes the SmoothCriterion
//=======================================================================
void AppParCurves_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 AppParCurves_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)
Standard_ConstructionError::Raise("AppParCurves_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 AppParCurves_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 AppParCurves_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 AppParCurves_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 AppParCurves_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)
Standard_ConstructionError::Raise("AppParCurves_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)
Standard_ConstructionError::Raise("AppParCurves_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;
}
}
}

178
src/AppParCurves/AppParCurves_Variational_7.gxx
View File

@@ -1,178 +0,0 @@
// Created on: 1997-09-17
// Created by: Philippe MANGIN
// Copyright (c) 1997-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.
//====================== Private Methodes =============================//
//=======================================================================
//function : Adjusting
//purpose : Smoothing's adjusting like STRIM routine "MAJLIS"
//=======================================================================
void AppParCurves_Variational::Adjusting(
Handle(AppParCurves_SmoothCriterion)& J,
Standard_Real& WQuadratic,
Standard_Real& WQuality,
Handle(FEmTool_Curve)& TheCurve,
TColStd_Array1OfReal& Ecarts)
{
// cout << "=========== Adjusting =============" << 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(AppParCurves_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 AppParCurves_MyCriterion(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;
}
}
}

292
src/AppParCurves/AppParCurves_Variational_8.gxx
View File

@@ -1,292 +0,0 @@
// Created on: 1999-02-15
// Created by: Igor FEOKTISTOV
// Copyright (c) 1999-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 <PLib_Base.hxx>
#include <PLib_JacobiPolynomial.hxx>
#include <PLib_HermitJacobi.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <math_Vector.hxx>
void AppParCurves_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)*)&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;
}
}
}

110
src/AppParCurves/AppParCurves_Variational_9.gxx
View File

@@ -1,110 +0,0 @@
// Created on: 1999-02-19
// Created by: Sergey KHROMOV
// Copyright (c) 1999-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.
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;
}
Standard_Boolean AppParCurves_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;
}

9
src/AppParCurves/FILES
View File

@@ -1,9 +0,0 @@
AppParCurves_Variational_1.gxx
AppParCurves_Variational_2.gxx
AppParCurves_Variational_3.gxx
AppParCurves_Variational_4.gxx
AppParCurves_Variational_5.gxx
AppParCurves_Variational_6.gxx
AppParCurves_Variational_7.gxx
AppParCurves_Variational_8.gxx
AppParCurves_Variational_9.gxx