mirror of
https://git.dev.opencascade.org/repos/occt.git
synced 2025-07-30 13:05:50 +03:00
d5f74e42d6
License statement text corrected; compiler warnings caused by Bison 2.41 disabled for MSVC; a few other compiler warnings on 54-bit Windows eliminated by appropriate type cast Wrong license statements corrected in several files. Copyright and license statements added in XSD and GLSL files. Copyright year updated in some files. Obsolete documentation files removed from DrawResources.
868 lines
23 KiB
Plaintext
868 lines
23 KiB
Plaintext
// Copyright (c) 1995-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.
|
|
|
|
// lpa, le 11/09/91
|
|
|
|
|
|
// Approximation d un ensemble de points contraints (MultiLine) avec une
|
|
// solution approchee (MultiCurve). L algorithme utilise est l algorithme
|
|
// d Uzawa du package mathematique.
|
|
|
|
#define No_Standard_RangeError
|
|
#define No_Standard_OutOfRange
|
|
|
|
#include <math_Vector.hxx>
|
|
#include <math_Matrix.hxx>
|
|
#include <AppParCurves_Constraint.hxx>
|
|
#include <AppParCurves_ConstraintCouple.hxx>
|
|
#include <AppParCurves_MultiPoint.hxx>
|
|
#include <AppParCurves.hxx>
|
|
#include <Standard_DimensionError.hxx>
|
|
#include <math_Uzawa.hxx>
|
|
#include <TColStd_Array1OfInteger.hxx>
|
|
#include <TColStd_Array2OfInteger.hxx>
|
|
#include <gp_Pnt.hxx>
|
|
#include <gp_Pnt2d.hxx>
|
|
#include <gp_Vec.hxx>
|
|
#include <gp_Vec2d.hxx>
|
|
#include <TColgp_Array1OfPnt.hxx>
|
|
#include <TColgp_Array1OfPnt2d.hxx>
|
|
#include <TColgp_Array1OfVec.hxx>
|
|
#include <TColgp_Array1OfVec2d.hxx>
|
|
|
|
|
|
AppParCurves_ResolConstraint::
|
|
AppParCurves_ResolConstraint(const MultiLine& SSP,
|
|
AppParCurves_MultiCurve& SCurv,
|
|
const Standard_Integer FirstPoint,
|
|
const Standard_Integer LastPoint,
|
|
const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
|
|
const math_Matrix& Bern,
|
|
const math_Matrix& DerivativeBern,
|
|
const Standard_Real Tolerance):
|
|
Cont(1, NbConstraints(SSP, FirstPoint,
|
|
LastPoint, TheConstraints),
|
|
1, NbColumns(SSP, SCurv.NbPoles()-1), 0.0),
|
|
DeCont(1, NbConstraints(SSP, FirstPoint,
|
|
LastPoint, TheConstraints),
|
|
1, NbColumns(SSP, SCurv.NbPoles()-1), 0.0),
|
|
Secont(1, NbConstraints(SSP, FirstPoint,
|
|
LastPoint, TheConstraints), 0.0),
|
|
CTCinv(1, NbConstraints(SSP, FirstPoint,
|
|
LastPoint, TheConstraints),
|
|
1, NbConstraints(SSP, FirstPoint,
|
|
LastPoint, TheConstraints)),
|
|
Vardua(1, NbConstraints(SSP, FirstPoint,
|
|
LastPoint, TheConstraints)),
|
|
IPas(1, LastPoint-FirstPoint+1),
|
|
ITan(1, LastPoint-FirstPoint+1),
|
|
ICurv(1, LastPoint-FirstPoint+1)
|
|
{
|
|
Standard_Integer i, j, k, NbCu= SCurv.NbCurves();
|
|
Standard_Integer Npt;
|
|
Standard_Integer Inc3, IncSec, IncCol, IP = 0, CCol;
|
|
Standard_Integer myindex, Def = SCurv.NbPoles()-1;
|
|
Standard_Integer nb3d, nb2d, Npol= Def+1, Npol2 = 2*Npol;
|
|
Standard_Real T1 = 0., T2 = 0., T3, Tmax, Daij;
|
|
Standard_Boolean Ok;
|
|
gp_Vec V;
|
|
gp_Vec2d V2d;
|
|
gp_Pnt Poi;
|
|
gp_Pnt2d Poi2d;
|
|
AppParCurves_ConstraintCouple mycouple;
|
|
AppParCurves_Constraint FC = AppParCurves_NoConstraint,
|
|
LC = AppParCurves_NoConstraint, Cons;
|
|
|
|
|
|
|
|
// Boucle de calcul du nombre de points de passage afin de dimensionner
|
|
// les matrices.
|
|
IncPass = 0;
|
|
IncTan = 0;
|
|
IncCurv = 0;
|
|
for (i = TheConstraints->Lower(); i <= TheConstraints->Upper(); i++) {
|
|
mycouple = TheConstraints->Value(i);
|
|
myindex = mycouple.Index();
|
|
Cons = mycouple.Constraint();
|
|
if (myindex == FirstPoint) FC = Cons;
|
|
if (myindex == LastPoint) LC = Cons;
|
|
if (Cons >= 1) {
|
|
IncPass++; // IncPass = nbre de points de passage.
|
|
IPas(IncPass) = myindex;
|
|
}
|
|
if (Cons >= 2) {
|
|
IncTan++; // IncTan= nbre de points de tangence.
|
|
ITan(IncTan) = myindex;
|
|
}
|
|
if (Cons == 3) {
|
|
IncCurv++; // IncCurv = nbre de pts de courbure.
|
|
ICurv(IncCurv) = myindex;
|
|
}
|
|
}
|
|
|
|
|
|
if (IncPass == 0) {
|
|
Done = Standard_True;
|
|
return;
|
|
}
|
|
|
|
nb3d = ToolLine::NbP3d(SSP);
|
|
nb2d = ToolLine::NbP2d(SSP);
|
|
Standard_Integer mynb3d=nb3d, mynb2d=nb2d;
|
|
if (nb3d == 0) mynb3d = 1;
|
|
if (nb2d == 0) mynb2d = 1;
|
|
CCol = nb3d* 3 + nb2d* 2;
|
|
|
|
|
|
// Declaration et initialisation des matrices et vecteurs de contraintes:
|
|
math_Matrix ContInit(1, IncPass, 1, Npol);
|
|
math_Vector Start(1, CCol*Npol);
|
|
TColStd_Array2OfInteger Ibont(1, NbCu, 1, IncTan);
|
|
|
|
|
|
// Remplissage de Cont pour les points de passage:
|
|
// =================================================
|
|
for (i =1; i <= IncPass; i++) { // Cette partie ne depend que de Bernstein
|
|
Npt = IPas(i);
|
|
for (j = 1; j <= Npol; j++) {
|
|
ContInit(i, j) = Bern(Npt, j);
|
|
}
|
|
}
|
|
for (i = 1; i <= CCol; i++) {
|
|
Cont.Set(IncPass*(i-1)+1, IncPass*i, Npol*(i-1)+1, Npol*i, ContInit);
|
|
}
|
|
|
|
|
|
// recuperation des vecteurs de depart pour Uzawa. Ce vecteur represente les
|
|
// poles de SCurv.
|
|
// Remplissage de secont et resolution.
|
|
|
|
TColgp_Array1OfVec tabV(1, mynb3d);
|
|
TColgp_Array1OfVec2d tabV2d(1, mynb2d);
|
|
TColgp_Array1OfPnt tabP(1, mynb3d);
|
|
TColgp_Array1OfPnt2d tabP2d(1, mynb2d);
|
|
|
|
Inc3 = CCol*IncPass +1;
|
|
IncCol = 0;
|
|
IncSec = 0;
|
|
for (k = 1; k <= NbCu; k++) {
|
|
if (k <= nb3d) {
|
|
for (i = 1; i <= IncTan; i++) {
|
|
Npt = ITan(i);
|
|
// choix du maximum de tangence pour exprimer la colinearite:
|
|
ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k);
|
|
V.Coord(T1, T2, T3);
|
|
Tmax = Abs(T1);
|
|
Ibont(k, i) = 1;
|
|
if (Abs(T2) > Tmax) {
|
|
Tmax = Abs(T2);
|
|
Ibont(k, i) = 2;
|
|
}
|
|
if (Abs(T3) > Tmax) {
|
|
Tmax = Abs(T3);
|
|
Ibont(k, i) = 3;
|
|
}
|
|
if (Ibont(k, i) == 3) {
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DerivativeBern(Npt, j);
|
|
Cont(Inc3, j+ Npol + IncCol) = Daij*T3/Tmax;
|
|
Cont(Inc3, j + Npol2 + IncCol) = -Daij *T2/Tmax;
|
|
Cont(Inc3+1, j+ IncCol) = Daij* T3/Tmax;
|
|
Cont(Inc3+1, j+ Npol2 + IncCol) = -Daij*T1/Tmax;
|
|
}
|
|
}
|
|
else if (Ibont(k, i) == 1) {
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DerivativeBern(Npt, j);
|
|
Cont(Inc3, j + IncCol) = Daij*T3/Tmax;
|
|
Cont(Inc3, j + Npol2 + IncCol) = -Daij *T1/Tmax;
|
|
Cont(Inc3+1, j+ IncCol) = Daij* T2/Tmax;
|
|
Cont(Inc3+1, j+ Npol +IncCol) = -Daij*T1/Tmax;
|
|
}
|
|
}
|
|
else if (Ibont(k, i) == 2) {
|
|
for (j = 1; j <= Def+1; j++) {
|
|
Daij = DerivativeBern(Npt, j);
|
|
Cont(Inc3, j+ Npol + IncCol) = Daij*T3/Tmax;
|
|
Cont(Inc3, j + Npol2 + IncCol) = -Daij *T2/Tmax;
|
|
Cont(Inc3+1, j+ IncCol) = Daij* T2/Tmax;
|
|
Cont(Inc3+1, j+ Npol + IncCol) = -Daij*T1/Tmax;
|
|
}
|
|
}
|
|
Inc3 = Inc3 + 2;
|
|
}
|
|
|
|
// Remplissage du second membre:
|
|
for (i = 1; i <= IncPass; i++) {
|
|
ToolLine::Value(SSP, IPas(i), tabP);
|
|
Poi = tabP(k);
|
|
Secont(i + IncSec) = Poi.X();
|
|
Secont(i + IncPass + IncSec) = Poi.Y();
|
|
Secont(i + 2*IncPass + IncSec) = Poi.Z();
|
|
}
|
|
IncSec = IncSec + 3*IncPass;
|
|
|
|
// Vecteur de depart:
|
|
for (j =1; j <= Npol; j++) {
|
|
Poi = SCurv.Value(j).Point(k);
|
|
Start(j + IncCol) = Poi.X();
|
|
Start(j+ Npol + IncCol) = Poi.Y();
|
|
Start(j + Npol2 + IncCol) = Poi.Z();
|
|
}
|
|
IncCol = IncCol + 3*Npol;
|
|
}
|
|
|
|
|
|
else {
|
|
for (i = 1; i <= IncTan; i++) {
|
|
Npt = ITan(i);
|
|
ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k-nb3d);
|
|
T1 = V2d.X();
|
|
T2 = V2d.Y();
|
|
Tmax = Abs(T1);
|
|
Ibont(k, i) = 1;
|
|
if (Abs(T2) > Tmax) {
|
|
Tmax = Abs(T2);
|
|
Ibont(k, i) = 2;
|
|
}
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DerivativeBern(Npt, j);
|
|
Cont(Inc3, j + IncCol) = Daij*T2;
|
|
Cont(Inc3, j + Npol + IncCol) = -Daij*T1;
|
|
}
|
|
Inc3 = Inc3 +1;
|
|
}
|
|
|
|
// Remplissage du second membre:
|
|
for (i = 1; i <= IncPass; i++) {
|
|
ToolLine::Value(SSP, IPas(i), tabP2d);
|
|
Poi2d = tabP2d(i-nb3d);
|
|
Secont(i + IncSec) = Poi2d.X();
|
|
Secont(i + IncPass + IncSec) = Poi2d.Y();
|
|
}
|
|
IncSec = IncSec + 2*IncPass;
|
|
|
|
// Remplissage du vecteur de depart:
|
|
for (j = 1; j <= Npol; j++) {
|
|
Poi2d = SCurv.Value(j).Point2d(k);
|
|
Start(j + IncCol) = Poi2d.X();
|
|
Start(j + Npol + IncCol) = Poi2d.Y();
|
|
}
|
|
IncCol = IncCol + Npol2;
|
|
}
|
|
}
|
|
|
|
// Equations exprimant le meme rapport de tangence sur chaque courbe:
|
|
// On prend les coordonnees les plus significatives.
|
|
|
|
--Inc3;
|
|
for (i = 1; i <= IncTan; ++i) {
|
|
IncCol = 0;
|
|
Npt = ITan(i);
|
|
for (k = 1; k <= NbCu-1; ++k) {
|
|
++Inc3;
|
|
// Initialize first relation variable (T1)
|
|
Standard_Integer addIndex_1 = 0, aVal = Ibont(k, i);
|
|
switch (aVal)
|
|
{
|
|
case 1: //T1 ~ T1x
|
|
{
|
|
if (k <= nb3d) {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k);
|
|
T1 = V.X();
|
|
IP = 3 * Npol;
|
|
}
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k-nb3d);
|
|
T1 = V2d.X();
|
|
IP = 2 * Npol;
|
|
}
|
|
addIndex_1 = 0;
|
|
break;
|
|
}
|
|
case 2: //T1 ~ T1y
|
|
{
|
|
if (k <= nb3d) {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k);
|
|
IP = 3 * Npol;
|
|
}
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k-nb3d);
|
|
T1 = V2d.Y();
|
|
IP = 2 * Npol;
|
|
}
|
|
addIndex_1 = Npol;
|
|
break;
|
|
}
|
|
case 3: //T1 ~ T1z
|
|
{
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k);
|
|
T1 = V.Z();
|
|
IP = 3 * Npol;
|
|
addIndex_1 = 2 * Npol;
|
|
break;
|
|
}
|
|
}
|
|
// Initialize second relation variable (T2)
|
|
Standard_Integer addIndex_2 = 0, aNextVal = Ibont(k + 1, i);
|
|
switch (aNextVal)
|
|
{
|
|
case 1: //T2 ~ T2x
|
|
{
|
|
if ((k+1) <= nb3d) {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k+1);
|
|
T2 = V.X();
|
|
}
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k+1-nb3d);
|
|
T2 = V2d.X();
|
|
}
|
|
addIndex_2 = 0;
|
|
break;
|
|
}
|
|
case 2: //T2 ~ T2y
|
|
{
|
|
if ((k+1) <= nb3d) {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k+1);
|
|
T2 = V.Y();
|
|
}
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k+1-nb3d);
|
|
T2 = V2d.Y();
|
|
}
|
|
addIndex_2 = Npol;
|
|
break;
|
|
}
|
|
case 3: //T2 ~ T2z
|
|
{
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k+1);
|
|
T2 = V.Z();
|
|
addIndex_2 = 2 * Npol;
|
|
break;
|
|
}
|
|
}
|
|
|
|
// Relations between T1 and T2:
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DerivativeBern(Npt, j);
|
|
Cont(Inc3, j + IncCol + addIndex_1) = Daij*T2;
|
|
Cont(Inc3, j + IP + IncCol + addIndex_2) = -Daij*T1;
|
|
}
|
|
IncCol += IP;
|
|
}
|
|
}
|
|
|
|
|
|
|
|
// Equations concernant la courbure:
|
|
|
|
|
|
/* Inc3 = Inc3 +1;
|
|
IncCol = 0;
|
|
for (k = 1; k <= NbCu; k++) {
|
|
for (i = 1; i <= IncCurv; i++) {
|
|
Npt = ICurv(i);
|
|
AppDef_MultiPointConstraint MP = SSP.Value(Npt);
|
|
DDA = SecondDerivativeBern(Parameters(Npt));
|
|
if (SSP.Value(1).Dimension(k) == 3) {
|
|
C1 = MP.Curv(k).X();
|
|
C2 = MP.Curv(k).Y();
|
|
C3 = MP.Curv(k).Z();
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DerivativeBern(Npt, j);
|
|
D2Aij = DDA(j);
|
|
Cont(Inc3, j + IncCol) = D2Aij;
|
|
Cont(Inc3, j + Npol2 + IncCol) = -C2*Daij;
|
|
Cont(Inc3, j + Npol + IncCol) = C3*Daij;
|
|
|
|
Cont(Inc3+1, j + Npol + IncCol) = D2Aij;
|
|
Cont(Inc3+1, j +IncCol) = -C3*Daij;
|
|
Cont(Inc3+1, j + Npol2 + IncCol) = C1*Daij;
|
|
|
|
Cont(Inc3+2, j + Npol2+IncCol) = D2Aij;
|
|
Cont(Inc3+2, j + Npol+IncCol) = -C1*Daij;
|
|
Cont(Inc3+2, j + IncCol) = C2*Daij;
|
|
}
|
|
Inc3 = Inc3 + 3;
|
|
}
|
|
else { // Dimension 2:
|
|
C1 = MP.Curv2d(k).X();
|
|
C2 = MP.Curv2d(k).Y();
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DerivativeBern(Npt, j);
|
|
D2Aij = DDA(j);
|
|
Cont(Inc3, j + IncCol) = D2Aij*C1;
|
|
Cont(Inc3+1, j + Npol + IncCol) = D2Aij*C2;
|
|
}
|
|
Inc3 = Inc3 + 2;
|
|
}
|
|
}
|
|
}
|
|
|
|
*/
|
|
// Resolution par Uzawa:
|
|
|
|
math_Uzawa UzaResol(Cont, Secont, Start, Tolerance);
|
|
if (!(UzaResol.IsDone())) {
|
|
Done = Standard_False;
|
|
return;
|
|
}
|
|
CTCinv = UzaResol.InverseCont();
|
|
UzaResol.Duale(Vardua);
|
|
for (i = 1; i <= CTCinv.RowNumber(); i++) {
|
|
for (j = i; j <= CTCinv.RowNumber(); j++) {
|
|
CTCinv(i, j) = CTCinv(j, i);
|
|
}
|
|
}
|
|
Done = Standard_True;
|
|
math_Vector VecPoles (1, CCol*Npol);
|
|
VecPoles = UzaResol.Value();
|
|
|
|
Standard_Integer polinit = 1, polfin=Npol;
|
|
if (FC >= 1) polinit = 2;
|
|
if (LC >= 1) polfin = Npol-1;
|
|
|
|
for (i = polinit; i <= polfin; i++) {
|
|
IncCol = 0;
|
|
AppParCurves_MultiPoint MPol(nb3d, nb2d);
|
|
for (k = 1; k <= NbCu; k++) {
|
|
if (k <= nb3d) {
|
|
gp_Pnt Pol(VecPoles(IncCol + i),
|
|
VecPoles(IncCol + Npol +i),
|
|
VecPoles(IncCol + 2*Npol + i));
|
|
MPol.SetPoint(k, Pol);
|
|
IncCol = IncCol + 3*Npol;
|
|
}
|
|
else {
|
|
gp_Pnt2d Pol2d(VecPoles(IncCol + i),
|
|
VecPoles(IncCol + Npol + i));
|
|
MPol.SetPoint2d(k, Pol2d);
|
|
IncCol = IncCol + 2*Npol;
|
|
}
|
|
}
|
|
SCurv.SetValue(i, MPol);
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
Standard_Boolean AppParCurves_ResolConstraint::IsDone () const {
|
|
return Done;
|
|
}
|
|
|
|
|
|
Standard_Integer AppParCurves_ResolConstraint::
|
|
NbConstraints(const MultiLine& SSP,
|
|
const Standard_Integer,
|
|
const Standard_Integer,
|
|
const Handle(AppParCurves_HArray1OfConstraintCouple)&
|
|
TheConstraints)
|
|
const
|
|
{
|
|
// Boucle de calcul du nombre de points de passage afin de dimensionner
|
|
// les matrices.
|
|
Standard_Integer aIncPass, aIncTan, aIncCurv, aCCol;
|
|
Standard_Integer i;
|
|
AppParCurves_Constraint Cons;
|
|
|
|
aIncPass = 0;
|
|
aIncTan = 0;
|
|
aIncCurv = 0;
|
|
|
|
for (i = TheConstraints->Lower(); i <= TheConstraints->Upper(); i++) {
|
|
Cons = (TheConstraints->Value(i)).Constraint();
|
|
if (Cons >= 1) {
|
|
aIncPass++; // IncPass = nbre de points de passage.
|
|
}
|
|
if (Cons >= 2) {
|
|
aIncTan++; // IncTan= nbre de points de tangence.
|
|
}
|
|
if (Cons == 3) {
|
|
aIncCurv++; // IncCurv = nbre de pts de courbure.
|
|
}
|
|
}
|
|
Standard_Integer nb3d = ToolLine::NbP3d(SSP);
|
|
Standard_Integer nb2d = ToolLine::NbP2d(SSP);
|
|
aCCol = nb3d* 3 + nb2d* 2;
|
|
|
|
return aCCol*aIncPass + aIncTan*(aCCol-1) + 3*aIncCurv;
|
|
|
|
}
|
|
|
|
|
|
Standard_Integer AppParCurves_ResolConstraint::NbColumns(const MultiLine& SSP,
|
|
const Standard_Integer Deg)
|
|
const
|
|
{
|
|
Standard_Integer nb3d = ToolLine::NbP3d(SSP);
|
|
Standard_Integer nb2d = ToolLine::NbP2d(SSP);
|
|
Standard_Integer CCol = nb3d* 3 + nb2d* 2;
|
|
|
|
return CCol*(Deg+1);
|
|
}
|
|
|
|
const math_Matrix& AppParCurves_ResolConstraint::ConstraintMatrix() const
|
|
{
|
|
return Cont;
|
|
}
|
|
|
|
const math_Matrix& AppParCurves_ResolConstraint::InverseMatrix() const
|
|
{
|
|
return CTCinv;
|
|
}
|
|
|
|
|
|
const math_Vector& AppParCurves_ResolConstraint::Duale() const
|
|
{
|
|
return Vardua;
|
|
}
|
|
|
|
|
|
|
|
const math_Matrix& AppParCurves_ResolConstraint::
|
|
ConstraintDerivative(const MultiLine& SSP,
|
|
const math_Vector& Parameters,
|
|
const Standard_Integer Deg,
|
|
const math_Matrix& DA)
|
|
{
|
|
Standard_Integer i, j, k, nb2d, nb3d, CCol, Inc3, IncCol, Npt;
|
|
Standard_Integer Npol, Npol2, NbCu =ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP);
|
|
Standard_Integer IP;
|
|
Standard_Real Daij;
|
|
Npol = Deg+1; Npol2 = 2*Npol;
|
|
Standard_Boolean Ok;
|
|
TColStd_Array2OfInteger Ibont(1, NbCu, 1, IncTan);
|
|
math_Matrix DeCInit(1, IncPass, 1, Npol);
|
|
math_Vector DDA(1, Npol);
|
|
gp_Vec V;
|
|
gp_Vec2d V2d;
|
|
Standard_Real T1, T2, T3, Tmax, DDaij;
|
|
// Standard_Integer FirstP = IPas(1);
|
|
nb3d = ToolLine::NbP3d(SSP);
|
|
nb2d = ToolLine::NbP2d(SSP);
|
|
Standard_Integer mynb3d = nb3d, mynb2d = nb2d;
|
|
if (nb3d == 0) mynb3d = 1;
|
|
if (nb2d == 0) mynb2d = 1;
|
|
CCol = nb3d* 3 + nb2d* 2;
|
|
|
|
TColgp_Array1OfVec tabV(1, mynb3d);
|
|
TColgp_Array1OfVec2d tabV2d(1, mynb2d);
|
|
TColgp_Array1OfPnt tabP(1, mynb3d);
|
|
TColgp_Array1OfPnt2d tabP2d(1, mynb2d);
|
|
|
|
for (i = 1; i <= DeCont.RowNumber(); i++)
|
|
for (j = 1; j <= DeCont.ColNumber(); j++)
|
|
DeCont(i, j) = 0.0;
|
|
|
|
|
|
// Remplissage de DK pour les points de passages:
|
|
|
|
for (i =1; i <= IncPass; i++) {
|
|
Npt = IPas(i);
|
|
for (j = 1; j <= Npol; j++) DeCInit(i, j) = DA(Npt, j);
|
|
}
|
|
for (i = 1; i <= CCol; i++) {
|
|
DeCont.Set(IncPass*(i-1)+1, IncPass*i, Npol*(i-1)+1, Npol*i, DeCInit);
|
|
}
|
|
|
|
|
|
// Pour les points de tangence:
|
|
|
|
Inc3 = CCol*IncPass +1;
|
|
IncCol = 0;
|
|
|
|
for (k = 1; k <= NbCu; k++) {
|
|
if (k <= nb3d) {
|
|
for (i = 1; i <= IncTan; i++) {
|
|
Npt = ITan(i);
|
|
// MultiPoint MPoint = ToolLine::Value(SSP, Npt);
|
|
// choix du maximum de tangence pour exprimer la colinearite:
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k);
|
|
V.Coord(T1, T2, T3);
|
|
Tmax = Abs(T1);
|
|
Ibont(k, i) = 1;
|
|
if (Abs(T2) > Tmax) {
|
|
Tmax = Abs(T2);
|
|
Ibont(k, i) = 2;
|
|
}
|
|
if (Abs(T3) > Tmax) {
|
|
Tmax = Abs(T3);
|
|
Ibont(k, i) = 3;
|
|
}
|
|
AppParCurves::SecondDerivativeBernstein(Parameters(Npt), DDA);
|
|
if (Ibont(k, i) == 3) {
|
|
for (j = 1; j <= Npol; j++) {
|
|
DDaij = DDA(j);
|
|
DeCont(Inc3, j+ Npol + IncCol) = DDaij*T3/Tmax;
|
|
DeCont(Inc3, j + Npol2 + IncCol) = -DDaij *T2/Tmax;
|
|
DeCont(Inc3+1, j+ IncCol) = DDaij* T3/Tmax;
|
|
DeCont(Inc3+1, j+ Npol2 + IncCol) = -DDaij*T1/Tmax;
|
|
}
|
|
}
|
|
else if (Ibont(k, i) == 1) {
|
|
for (j = 1; j <= Npol; j++) {
|
|
DDaij = DDA(j);
|
|
DeCont(Inc3, j + IncCol) = DDaij*T3/Tmax;
|
|
DeCont(Inc3, j + Npol2 + IncCol) = -DDaij *T1/Tmax;
|
|
DeCont(Inc3+1, j+ IncCol) = DDaij* T2/Tmax;
|
|
DeCont(Inc3+1, j+ Npol +IncCol) = -DDaij*T1/Tmax;
|
|
}
|
|
}
|
|
else if (Ibont(k, i) == 2) {
|
|
for (j = 1; j <= Npol; j++) {
|
|
DDaij = DDA(j);
|
|
DeCont(Inc3, j+ Npol + IncCol) = DDaij*T3/Tmax;
|
|
DeCont(Inc3, j + Npol2 + IncCol) = -DDaij *T2/Tmax;
|
|
DeCont(Inc3+1, j+ IncCol) = DDaij* T2/Tmax;
|
|
DeCont(Inc3+1, j+ Npol + IncCol) = -DDaij*T1/Tmax;
|
|
}
|
|
}
|
|
Inc3 = Inc3 + 2;
|
|
}
|
|
IncCol = IncCol + 3*Npol;
|
|
}
|
|
else {
|
|
for (i = 1; i <= IncTan; i++) {
|
|
Npt = ITan(i);
|
|
AppParCurves::SecondDerivativeBernstein(Parameters(Npt), DDA);
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k);
|
|
V2d.Coord(T1, T2);
|
|
Tmax = Abs(T1);
|
|
Ibont(k, i) = 1;
|
|
if (Abs(T2) > Tmax) {
|
|
Tmax = Abs(T2);
|
|
Ibont(k, i) = 2;
|
|
}
|
|
for (j = 1; j <= Npol; j++) {
|
|
DDaij = DDA(j);
|
|
DeCont(Inc3, j + IncCol) = DDaij*T2;
|
|
DeCont(Inc3, j + Npol + IncCol) = -DDaij*T1;
|
|
}
|
|
Inc3 = Inc3 +1;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Equations exprimant le meme rapport de tangence sur chaque courbe:
|
|
// On prend les coordonnees les plus significatives.
|
|
|
|
Inc3 = Inc3 -1;
|
|
for (i =1; i <= IncTan; i++) {
|
|
IncCol = 0;
|
|
Npt = ITan(i);
|
|
AppParCurves::SecondDerivativeBernstein(Parameters(Npt), DDA);
|
|
// MultiPoint MP = ToolLine::Value(SSP, Npt);
|
|
for (k = 1; k <= NbCu-1; k++) {
|
|
Inc3 = Inc3 + 1;
|
|
if (Ibont(k, i) == 1) {
|
|
if (k <= nb3d) {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k);
|
|
T1 = V.X();
|
|
IP = 3*Npol;
|
|
}
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k);
|
|
T1 = V2d.X();
|
|
IP = 2*Npol;
|
|
}
|
|
if (Ibont(k+1, i) == 1) { // Relations entre T1x et T2x:
|
|
if ((k+1) <= nb3d) {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k+1);
|
|
T2 = V.X();
|
|
}
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k+1);
|
|
T2 = V2d.X();
|
|
}
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DDA(j);
|
|
Cont(Inc3, j + IncCol) = Daij*T2;
|
|
Cont(Inc3, j + IP + IncCol) = -Daij*T1;
|
|
}
|
|
IncCol = IncCol + IP;
|
|
}
|
|
else if (Ibont(k+1, i) == 2) { // Relations entre T1x et T2y:
|
|
if ((k+1) <= nb3d) {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k+1);
|
|
T2 = V.Y();
|
|
}
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k+1);
|
|
T2 = V2d.Y();
|
|
}
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DDA(j);
|
|
Cont(Inc3, j + IncCol) = Daij*T2;
|
|
Cont(Inc3, j + IP + Npol + IncCol) = -Daij*T1;
|
|
}
|
|
IncCol = IncCol + IP;
|
|
}
|
|
else if (Ibont(k+1,i) == 3) { // Relations entre T1x et T2z:
|
|
ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k+1);
|
|
T2 = V.Z();
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DDA(j);
|
|
Cont(Inc3, j + IncCol) = Daij*T2;
|
|
Cont(Inc3, j + IP + 2*Npol + IncCol) = -Daij*T1;
|
|
}
|
|
IncCol = IncCol + IP;
|
|
}
|
|
}
|
|
else if (Ibont(k,i) == 2) {
|
|
if (k <= nb3d) {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k);
|
|
T1 = V.Y();
|
|
IP = 3*Npol;
|
|
}
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k);
|
|
T1 = V2d.Y();
|
|
IP = 2*Npol;
|
|
}
|
|
if (Ibont(k+1, i) == 1) { // Relations entre T1y et T2x:
|
|
if ((k+1) <= nb3d) {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k+1);
|
|
T2 = V.X();
|
|
}
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k+1);
|
|
T2 = V2d.X();
|
|
}
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DDA( j);
|
|
Cont(Inc3, j + Npol + IncCol) = Daij*T2;
|
|
Cont(Inc3, j + IP + IncCol) = -Daij*T1;
|
|
}
|
|
IncCol = IncCol + IP;
|
|
|
|
}
|
|
else if (Ibont(k+1, i) == 2) { // Relations entre T1y et T2y:
|
|
if ((k+1) <= nb3d) {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k+1);
|
|
T2 = V.Y();
|
|
}
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k+1);
|
|
T2 = V2d.Y();
|
|
}
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DDA(j);
|
|
Cont(Inc3, j + Npol + IncCol) = Daij*T2;
|
|
Cont(Inc3, j + IP + Npol + IncCol) = -Daij*T1;
|
|
}
|
|
IncCol = IncCol + IP;
|
|
|
|
}
|
|
else if (Ibont(k+1,i)== 3) { // Relations entre T1y et T2z:
|
|
ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k+1);
|
|
T2 = V.Z();
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DDA(j);
|
|
Cont(Inc3, j + Npol +IncCol) = Daij*T2;
|
|
Cont(Inc3, j + IP + 2*Npol + IncCol) = -Daij*T1;
|
|
}
|
|
IncCol = IncCol + IP;
|
|
}
|
|
}
|
|
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k);
|
|
T1 = V.Z();
|
|
IP = 3*Npol;
|
|
if (Ibont(k+1, i) == 1) { // Relations entre T1z et T2x:
|
|
if ((k+1) <= nb3d) {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k+1);
|
|
T2 = V.X();
|
|
}
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k+1);
|
|
T2 = V2d.X();
|
|
}
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DDA(j);
|
|
Cont(Inc3, j + 2*Npol +IncCol) = Daij*T2;
|
|
Cont(Inc3, j + IP + IncCol) = -Daij*T1;
|
|
}
|
|
IncCol = IncCol + IP;
|
|
}
|
|
|
|
else if (Ibont(k+1, i) == 2) { // Relations entre T1z et T2y:
|
|
if ((k+1) <= nb3d) {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k+1);
|
|
T2 = V.Y();
|
|
}
|
|
else {
|
|
Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
|
|
V2d = tabV2d(k+1);
|
|
T2 = V2d.Y();
|
|
}
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DDA(j);
|
|
Cont(Inc3, j + 2*Npol +IncCol) = Daij*T2;
|
|
Cont(Inc3, j + IP + Npol + IncCol) = -Daij*T1;
|
|
}
|
|
IncCol = IncCol + IP;
|
|
}
|
|
|
|
else if (Ibont(k+1,i)==3) { // Relations entre T1z et T2z:
|
|
ToolLine::Tangency(SSP, Npt, tabV);
|
|
V = tabV(k+1);
|
|
T2 = V.Z();
|
|
for (j = 1; j <= Npol; j++) {
|
|
Daij = DDA(j);
|
|
Cont(Inc3, j + 2*Npol + IncCol) = Daij*T2;
|
|
Cont(Inc3, j + IP + 2*Npol + IncCol) = -Daij*T1;
|
|
}
|
|
IncCol = IncCol + IP;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return DeCont;
|
|
}
|