[occt.git] / src / FEmTool / FEmTool_Curve.cxx

// File:	FEmTool_Curve.cxx
// Created:	Tue Nov  4 10:45:06 1997
// Author:	Roman BORISOV / Igor FEOKTISTOV
//		<rbv@redfox.nnov.matra-dtv.fr>


#define No_Standard_RangeError
#define No_Standard_OutOfRange

#include <FEmTool_Curve.ixx>
#include <PLib.hxx>
#include <PLib_JacobiPolynomial.hxx>
#include <PLib_HermitJacobi.hxx>

//=======================================================================
//function : FEmTool_Curve
//purpose  :
//=======================================================================
FEmTool_Curve::FEmTool_Curve(const Standard_Integer Dimension,
			     const Standard_Integer NbElements,
			     const Handle(PLib_Base)& TheBase, 
			     const Standard_Real Tolerance) : 
       myNbElements(NbElements), myDimension(Dimension), 
       myTolerance(Tolerance), 
       myBase(TheBase), myDegree(1, myNbElements), 
       myCoeff(1, myDimension*myNbElements*(myBase->WorkDegree() + 1)), 
       myPoly(1, myDimension*myNbElements*(myBase->WorkDegree() + 1)), 
       myDeri(1, myDimension*myNbElements*(myBase->WorkDegree())), 
       myDsecn(1, myDimension*myNbElements*(myBase->WorkDegree() - 1)), 
       HasPoly(1, myNbElements), HasDeri(1, myNbElements), 
       HasSecn(1, myNbElements), myLength(1, myNbElements),
       myIndex(0)
{
  myKnots = new TColStd_HArray1OfReal(1, myNbElements + 1);
  myDegree.Init(myBase->WorkDegree());
  HasPoly.Init(0);
  HasDeri.Init(0);
  HasSecn.Init(0);
  myLength.Init(-1);
}

 TColStd_Array1OfReal& FEmTool_Curve::Knots() const
{
  return myKnots->ChangeArray1();
}


//=======================================================================
//function : SetElement
//purpose  :
//=======================================================================
 void FEmTool_Curve::SetElement(const Standard_Integer IndexOfElement,
				const TColStd_Array2OfReal& Coeffs) 
{
  Standard_Integer i, j, degBase, deg;
  if (IndexOfElement > myNbElements || IndexOfElement < 1) Standard_OutOfRange::Raise(); 
  degBase = myBase->WorkDegree();
  deg = myDegree(IndexOfElement);
  Standard_Integer iBase = (IndexOfElement - 1)*(degBase + 1)*myDimension, 
                   i1 = iBase-myDimension, i2 = Coeffs.LowerRow() - 1, 
                   j1 = Coeffs.LowerCol() - 1;
  for(i = 1; i <= deg + 1; i++) {
    i1 += myDimension; i2++; 
    for(j = 1; j <= myDimension; j++)
      myCoeff(i1 + j) = Coeffs(i2, j1 + j);
  }

  Standard_Real stenor = (myKnots->Value(IndexOfElement + 1) - myKnots->Value(IndexOfElement)) / 2.,
                mfact;
  Handle(PLib_HermitJacobi) myHermitJacobi = (*((Handle(PLib_HermitJacobi)*)&myBase));

  i1 = iBase;
  i2 = iBase + (myHermitJacobi->NivConstr() + 1) * myDimension;
  for(i = 1; i <= myHermitJacobi->NivConstr(); i++) {
    i1 += myDimension; i2 += myDimension; 
    mfact = Pow(stenor, i);
    for(j = 1; j <= myDimension; j++) {
      myCoeff(i1 + j) *= mfact;
      myCoeff(i2 + j) *= mfact;
    }
  }

  HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
  myLength(IndexOfElement) = - 1;
}


//=======================================================================
//function : GetElement
//purpose  :
//=======================================================================
 void FEmTool_Curve::GetElement(const Standard_Integer IndexOfElement, TColStd_Array2OfReal& Coeffs) 
{
  Standard_Integer i, j, degBase, deg;
  if (IndexOfElement > myNbElements || IndexOfElement < 1) Standard_OutOfRange::Raise(); 
  degBase = myBase->WorkDegree();
  deg = myDegree(IndexOfElement);
  Standard_Integer iBase = (IndexOfElement - 1)*(degBase + 1)*myDimension, 
                   i1 = iBase-myDimension, i2 = Coeffs.LowerRow() - 1, 
                   j1 = Coeffs.LowerCol() - 1;
  for(i = 1; i <= deg + 1; i++) {
    i1 += myDimension; i2++; 
    for(j = 1; j <= myDimension; j++)
      Coeffs(i2, j1 + j) = myCoeff.Value(i1 + j);
  }

  Standard_Real stenor = 2. / (myKnots->Value(IndexOfElement + 1) - myKnots->Value(IndexOfElement)),
                mfact;

  Handle(PLib_HermitJacobi) myHermitJacobi = (*((Handle(PLib_HermitJacobi)*)&myBase));

  i2 = Coeffs.LowerRow();
  Standard_Integer i3 = i2 + myHermitJacobi->NivConstr() + 1;

  for(i = 1; i <= myHermitJacobi->NivConstr(); i++) {
    mfact = Pow(stenor, i);
    for(j = j1+1; j <= myDimension; j++) {
      Coeffs(i2 + i, j) *= mfact;
      Coeffs(i3 + i, j) *= mfact;
    }
  }
}


//=======================================================================
//function : GetPolynom
//purpose  :
//=======================================================================
 void FEmTool_Curve::GetPolynom(TColStd_Array1OfReal& Coeffs) 

{
  Standard_Integer IndexOfElement, i, di = Coeffs.Lower() - myPoly.Lower();


  for(IndexOfElement = 1; IndexOfElement <= myNbElements; IndexOfElement++)
    if (!HasPoly.Value(IndexOfElement)) Update(IndexOfElement, 0);

  for(i = myPoly.Lower(); i <= myPoly.Upper(); i++)
    Coeffs(di + i) = myPoly(i);

}


//=======================================================================
//function : D0
//purpose  :
//=======================================================================
 void FEmTool_Curve::D0(const Standard_Real U, TColStd_Array1OfReal& Pnt) 
{
  Standard_Integer deg;
  Standard_Real S;

  if (!myIndex || (U<Uf) || (U>Ul) ||
      (myKnots->Value(myIndex)!=Uf) ||
      (myKnots->Value(myIndex+1)!=Ul)) { 
    // Search the span
    if(U <= myKnots->Value(2)) myIndex = 1;
    else {
      for(myIndex = 2; myIndex <= myNbElements; myIndex++)
	if(U >= myKnots->Value(myIndex) && U <= myKnots->Value(myIndex+1)) break;
      if (myIndex > myNbElements) myIndex = myNbElements;
    }
    Uf = myKnots->Value(myIndex);
    Ul = myKnots->Value(myIndex+1);
    Denom = 1. /(Ul-Uf);
    USum = Uf+Ul;
    myPtr = (myIndex - 1)*(myBase->WorkDegree()+ 1)*myDimension + 1;
  }

  deg = myDegree(myIndex);
  if (!HasPoly.Value(myIndex)) Update(myIndex, 0);

//Parameter normalization: S [-1, 1]
  S = (2*U - USum)*Denom;
  PLib::NoDerivativeEvalPolynomial(S, deg, myDimension, deg*myDimension, 
				   myPoly(myPtr), Pnt(Pnt.Lower()));
}


//=======================================================================
//function : D1
//purpose  :
//=======================================================================
 void FEmTool_Curve::D1(const Standard_Real U, TColStd_Array1OfReal& Vec) 
{
  Standard_Integer deg, i;
  Standard_Real S;

  if (!myIndex || (U<Uf) || (U>Ul) ||
      (myKnots->Value(myIndex)!=Uf) ||
      (myKnots->Value(myIndex+1)!=Ul)) { 
    // Search the span
    if(U <= myKnots->Value(2)) myIndex = 1;
    else {
      for(myIndex = 2; myIndex <= myNbElements; myIndex++)
	if(U >= myKnots->Value(myIndex) && U <= myKnots->Value(myIndex+1)) break;
      if (myIndex > myNbElements) myIndex = myNbElements;
    }
    Uf = myKnots->Value(myIndex);
    Ul = myKnots->Value(myIndex+1);
    Denom = 1. /(Ul-Uf);
    USum = Uf+Ul;
    myPtr = (myIndex - 1)*(myBase->WorkDegree()+ 1)*myDimension + 1;
  }

  deg = myDegree(myIndex);
  if (!HasDeri.Value(myIndex)) Update(myIndex, 1);
  
//Parameter normalization: S [-1, 1]
  S = (2*U - USum) * Denom;
  PLib::NoDerivativeEvalPolynomial(S, deg-1, myDimension, (deg-1)*myDimension, 
				   myDeri(1+(myIndex - 1)*myBase->WorkDegree()*
					  myDimension), 
				   Vec(Vec.Lower()));
  
  S = 2*Denom;
  for(i=Vec.Lower(); i<=Vec.Upper(); i++) Vec(i) *= S;
}


//=======================================================================
//function : D2
//purpose  :
//=======================================================================
 void FEmTool_Curve::D2(const Standard_Real U, TColStd_Array1OfReal& Vec) 
{
  Standard_Integer deg, i;
  Standard_Real S;

  if (!myIndex || (U<Uf) || (U>Ul) ||
      (myKnots->Value(myIndex)!=Uf) ||
      (myKnots->Value(myIndex+1)!=Ul)) { 
    // Search the span
    if(U <= myKnots->Value(2)) myIndex = 1;
    else {
      for(myIndex = 2; myIndex <= myNbElements; myIndex++)
	if(U >= myKnots->Value(myIndex) && U <= myKnots->Value(myIndex+1)) break;
      if (myIndex > myNbElements) myIndex = myNbElements;
    }
    Uf = myKnots->Value(myIndex);
    Ul = myKnots->Value(myIndex+1);
    Denom = 1. /(Ul-Uf);
    USum = Uf+Ul;
    myPtr = (myIndex - 1)*(myBase->WorkDegree()+ 1)*myDimension + 1;
  }

  deg = myDegree(myIndex);
  if (!HasSecn.Value(myIndex)) Update(myIndex, 2); 

//Parameter normalization: S [-1, 1]
  S = (2*U - USum)*Denom;
  PLib::NoDerivativeEvalPolynomial(S, deg-2, myDimension, (deg-2)*myDimension, 
				   myDsecn(1+(myIndex - 1)*
					   (myBase->WorkDegree() - 1)*myDimension), 
				   Vec(Vec.Lower()));

  S = 4*Denom*Denom;
  for(i=Vec.Lower(); i<=Vec.Upper(); i++) Vec(i) *= S;
}


//=======================================================================
//function : Length
//purpose  :
//=======================================================================
 void FEmTool_Curve::Length(const Standard_Real FirstU,
			    const Standard_Real LastU,
			    Standard_Real& Length) 
{
  Standard_Integer Low, High, deg, degBase, i, Ptr;
  if(FirstU > LastU) Standard_OutOfRange::Raise("FEmTool_Curve::Length");  
  
  if(myKnots->Value(1) > FirstU) Low = 1;
  else
    for(Low = 1; Low <= myNbElements; Low++)
      if(FirstU >= myKnots->Value(Low) && FirstU <= myKnots->Value(Low+1)) break;
  if(Low > myNbElements) Low = myNbElements;
  
  if(myKnots->Value(1) > LastU) High = 1;
  else
    for(High = Low; High <= myNbElements; High++)
      if(LastU >= myKnots->Value(High) && LastU <= myKnots->Value(High+1)) break;
  if(myKnots->Value(myNbElements + 1) < LastU) High = myNbElements;
  
  Standard_Real Li;
  degBase = myBase->WorkDegree();
  Length = 0;

  Standard_Real FirstS, LastS;
  FirstS = (2*FirstU - myKnots->Value(Low) - myKnots->Value(Low + 1))/
    (myKnots->Value(Low + 1) - myKnots->Value(Low));
  LastS = (2*LastU - myKnots->Value(High) - myKnots->Value(High + 1))/
    (myKnots->Value(High + 1) - myKnots->Value(High));
  
  if(Low == High) {
    
    Ptr = (Low - 1)*(degBase + 1)*myDimension + 1;
    deg = myDegree(Low);
    
    if(!HasPoly(Low)) Update(Low, 0);
    PLib::EvalLength(deg, myDimension, myPoly(Ptr), 
		     FirstS, LastS, Length);
    return;
  }
  
  deg = myDegree(Low);
  Ptr = (Low - 1)*(degBase + 1)*myDimension + 1;

  if(!HasPoly(Low))  Update(Low, 0);
  if(FirstS < -1.) {
    PLib::EvalLength(deg, myDimension, myPoly(Ptr), FirstS, -1., Li);
    Length += Li;
    if(myLength(Low) < 0.) { 
      PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., 1.,Li);
      myLength(Low) = Li;
    }
    Length += myLength(Low);
  }
  else {
    PLib::EvalLength(deg, myDimension, myPoly(Ptr), FirstS, 1., Li);
    Length += Li;
  }

  deg = myDegree(High);
  Ptr = (High - 1)*(degBase + 1)*myDimension + 1;

  if(!HasPoly(High)) Update(High, 0);
  if(LastS > 1.) {
    PLib::EvalLength(deg, myDimension, myPoly(Ptr), 1., LastS, Li);
    Length += Li;
    if(myLength(High) < 0.) { 
      PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., 1., Li);
      myLength(High) = Li;
    }
    Length += myLength(High);
  }
  else {
    PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., LastS, Li);
    Length += Li;
  }

  
  for(i = Low + 1; i < High; i++) {
    if (myLength.Value(i) < 0) {
      Ptr = (i - 1)*(degBase + 1)*myDimension + 1;
      deg = myDegree(i);
      if(!HasPoly(i)) Update(i, 0);
      PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., 1., Li);
      myLength(i) = Li; 
    }
    Length += myLength(i);
  }
}

 Standard_Integer FEmTool_Curve::NbElements() const
{
  return myNbElements;
}

 Standard_Integer FEmTool_Curve::Dimension() const
{
  return myDimension;
}

 Handle(PLib_Base) FEmTool_Curve::Base() const
{
  return myBase;
}

 Standard_Integer FEmTool_Curve::Degree(const Standard_Integer IndexOfElement) const
{
  return myDegree.Value(IndexOfElement);
}


//=======================================================================
//function : SetDegree
//purpose  :
//=======================================================================
 void FEmTool_Curve::SetDegree(const Standard_Integer IndexOfElement,
			       const Standard_Integer Degree)
{
  if (Degree <= myBase->WorkDegree()) {
    myDegree(IndexOfElement) = Degree;
    HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
    myLength(IndexOfElement) = -1;
  }
  else if(Degree > myBase->WorkDegree()) Standard_OutOfRange::Raise("FEmTool_Curve::SetDegree");
}


//=======================================================================
//function : ReduceDegree
//purpose  :
//=======================================================================
 void FEmTool_Curve::ReduceDegree(const Standard_Integer IndexOfElement,
				  const Standard_Real Tol,
				  Standard_Integer& NewDegree,
				  Standard_Real& MaxError)
{
  Standard_Integer deg = myDegree(IndexOfElement);
  
  Standard_Integer Ptr = (IndexOfElement - 1)*
    (myBase->WorkDegree() + 1)*myDimension + 1;
  
  myBase->ReduceDegree(myDimension, deg, Tol, myCoeff.ChangeValue(Ptr), NewDegree, MaxError);
  Handle(PLib_HermitJacobi) myHermitJacobi = (*((Handle(PLib_HermitJacobi)*)&myBase));
  
  NewDegree = Max(NewDegree, 2 * myHermitJacobi->NivConstr() + 1);
  
  if (NewDegree < deg) {
    myDegree(IndexOfElement) = NewDegree;
    HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
    myLength(IndexOfElement) = -1;
  }
}

//=======================================================================
//function : Update
//purpose  :
//=======================================================================
 void FEmTool_Curve::Update(const Standard_Integer Index,
			    const Standard_Integer Order)
{
  Standard_Integer  degBase = myBase->WorkDegree(), deg = myDegree(Index);

  if (!HasPoly(Index)) {
    Standard_Integer Ptr = (Index - 1)*(degBase + 1)*myDimension + 1;

    TColStd_Array1OfReal Coeff(myPoly.ChangeValue(Ptr), 0, myDimension*(deg + 1) - 1);
    TColStd_Array1OfReal BaseCoeff(myCoeff.ChangeValue(Ptr), 0, myDimension*(deg + 1) - 1);

    myBase->ToCoefficients(myDimension, deg, BaseCoeff, Coeff);
    HasPoly(Index) = 1;
  }
 
  if (Order >=1) {
    Standard_Integer i1 = (Index - 1)*(degBase)*myDimension - myDimension, 
                     i2 = (Index - 1)*(degBase + 1)*myDimension;
    Standard_Integer i,j;
    if (!HasDeri.Value(Index)) {
      for(i = 1; i <= deg; i++) {
	i1 += myDimension; i2 += myDimension;
	for(j = 1; j<= myDimension; j++)
	  myDeri(i1 + j) = i*myPoly.Value(i2 + j);
      }
      HasDeri(Index) = 1;
    }
    if ((Order>=2) &&!HasSecn.Value(Index)) {
      i1 = (Index - 1)*(degBase-1)*myDimension - myDimension, 
      i2 = (Index - 1)*(degBase)*myDimension;
      for(i = 1; i < deg; i++) {
	i1 += myDimension; i2 += myDimension;
	for(j = 1; j<= myDimension; j++)
	  myDsecn(i1 + j) = i*myDeri.Value(i2 + j);
      }
      HasSecn(Index) = 1;
    }  
  }
}
Commit	Line	Data
7fd59977	1	// File: FEmTool_Curve.cxx
	2	// Created: Tue Nov 4 10:45:06 1997
	3	// Author: Roman BORISOV / Igor FEOKTISTOV
	4	// <rbv@redfox.nnov.matra-dtv.fr>
	5
	6
	7	#define No_Standard_RangeError
	8	#define No_Standard_OutOfRange
	9
	10	#include <FEmTool_Curve.ixx>
	11	#include <PLib.hxx>
	12	#include <PLib_JacobiPolynomial.hxx>
	13	#include <PLib_HermitJacobi.hxx>
	14
	15	//=======================================================================
	16	//function : FEmTool_Curve
	17	//purpose :
	18	//=======================================================================
	19	FEmTool_Curve::FEmTool_Curve(const Standard_Integer Dimension,
	20	const Standard_Integer NbElements,
	21	const Handle(PLib_Base)& TheBase,
	22	const Standard_Real Tolerance) :
	23	myNbElements(NbElements), myDimension(Dimension),
	24	myTolerance(Tolerance),
	25	myBase(TheBase), myDegree(1, myNbElements),
	26	myCoeff(1, myDimensionmyNbElements(myBase->WorkDegree() + 1)),
	27	myPoly(1, myDimensionmyNbElements(myBase->WorkDegree() + 1)),
	28	myDeri(1, myDimensionmyNbElements(myBase->WorkDegree())),
	29	myDsecn(1, myDimensionmyNbElements(myBase->WorkDegree() - 1)),
	30	HasPoly(1, myNbElements), HasDeri(1, myNbElements),
	31	HasSecn(1, myNbElements), myLength(1, myNbElements),
	32	myIndex(0)
	33	{
	34	myKnots = new TColStd_HArray1OfReal(1, myNbElements + 1);
	35	myDegree.Init(myBase->WorkDegree());
	36	HasPoly.Init(0);
	37	HasDeri.Init(0);
	38	HasSecn.Init(0);
	39	myLength.Init(-1);
	40	}
	41
	42	TColStd_Array1OfReal& FEmTool_Curve::Knots() const
	43	{
	44	return myKnots->ChangeArray1();
	45	}
	46
	47
	48
	49	//=======================================================================
	50	//function : SetElement
	51	//purpose :
	52	//=======================================================================
	53	void FEmTool_Curve::SetElement(const Standard_Integer IndexOfElement,
	54	const TColStd_Array2OfReal& Coeffs)
	55	{
	56	Standard_Integer i, j, degBase, deg;
	57	if (IndexOfElement > myNbElements \|\| IndexOfElement < 1) Standard_OutOfRange::Raise();
	58	degBase = myBase->WorkDegree();
	59	deg = myDegree(IndexOfElement);
	60	Standard_Integer iBase = (IndexOfElement - 1)(degBase + 1)myDimension,
	61	i1 = iBase-myDimension, i2 = Coeffs.LowerRow() - 1,
	62	j1 = Coeffs.LowerCol() - 1;
	63	for(i = 1; i <= deg + 1; i++) {
	64	i1 += myDimension; i2++;
65	for(j = 1; j <= myDimension; j++)
66	myCoeff(i1 + j) = Coeffs(i2, j1 + j);
67	}
68
69	Standard_Real stenor = (myKnots->Value(IndexOfElement + 1) - myKnots->Value(IndexOfElement)) / 2.,
70	mfact;
71	Handle(PLib_HermitJacobi) myHermitJacobi = (((Handle(PLib_HermitJacobi))&myBase));
72
73	i1 = iBase;
74	i2 = iBase + (myHermitJacobi->NivConstr() + 1) * myDimension;
75	for(i = 1; i <= myHermitJacobi->NivConstr(); i++) {
76	i1 += myDimension; i2 += myDimension;
77	mfact = Pow(stenor, i);
78	for(j = 1; j <= myDimension; j++) {
79	myCoeff(i1 + j) *= mfact;
80	myCoeff(i2 + j) *= mfact;
81	}
82	}
83
84	HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
85	myLength(IndexOfElement) = - 1;
86	}
87
88
89	//=======================================================================
90	//function : GetElement
91	//purpose :
92	//=======================================================================
93	void FEmTool_Curve::GetElement(const Standard_Integer IndexOfElement, TColStd_Array2OfReal& Coeffs)
94	{
95	Standard_Integer i, j, degBase, deg;
96	if (IndexOfElement > myNbElements \|\| IndexOfElement < 1) Standard_OutOfRange::Raise();
97	degBase = myBase->WorkDegree();
98	deg = myDegree(IndexOfElement);
99	Standard_Integer iBase = (IndexOfElement - 1)(degBase + 1)myDimension,
100	i1 = iBase-myDimension, i2 = Coeffs.LowerRow() - 1,
101	j1 = Coeffs.LowerCol() - 1;
102	for(i = 1; i <= deg + 1; i++) {
103	i1 += myDimension; i2++;
104	for(j = 1; j <= myDimension; j++)
105	Coeffs(i2, j1 + j) = myCoeff.Value(i1 + j);
106	}
107
108	Standard_Real stenor = 2. / (myKnots->Value(IndexOfElement + 1) - myKnots->Value(IndexOfElement)),
109	mfact;
110
111	Handle(PLib_HermitJacobi) myHermitJacobi = (((Handle(PLib_HermitJacobi))&myBase));
112
113	i2 = Coeffs.LowerRow();
114	Standard_Integer i3 = i2 + myHermitJacobi->NivConstr() + 1;
115
116	for(i = 1; i <= myHermitJacobi->NivConstr(); i++) {
117	mfact = Pow(stenor, i);
118	for(j = j1+1; j <= myDimension; j++) {
119	Coeffs(i2 + i, j) *= mfact;
120	Coeffs(i3 + i, j) *= mfact;
121	}
122	}
123	}
124
125
126	//=======================================================================
127	//function : GetPolynom
128	//purpose :
129	//=======================================================================
130	void FEmTool_Curve::GetPolynom(TColStd_Array1OfReal& Coeffs)
131
132	{
133	Standard_Integer IndexOfElement, i, di = Coeffs.Lower() - myPoly.Lower();
134
135
136	for(IndexOfElement = 1; IndexOfElement <= myNbElements; IndexOfElement++)
137	if (!HasPoly.Value(IndexOfElement)) Update(IndexOfElement, 0);
138
139	for(i = myPoly.Lower(); i <= myPoly.Upper(); i++)
140	Coeffs(di + i) = myPoly(i);
141
142	}
143
144
145	//=======================================================================
146	//function : D0
147	//purpose :
148	//=======================================================================
149	void FEmTool_Curve::D0(const Standard_Real U, TColStd_Array1OfReal& Pnt)
150	{
151	Standard_Integer deg;
152	Standard_Real S;
153
154	if (!myIndex \|\| (U<Uf) \|\| (U>Ul) \|\|
155	(myKnots->Value(myIndex)!=Uf) \|\|
156	(myKnots->Value(myIndex+1)!=Ul)) {
157	// Search the span
158	if(U <= myKnots->Value(2)) myIndex = 1;
159	else {
160	for(myIndex = 2; myIndex <= myNbElements; myIndex++)
161	if(U >= myKnots->Value(myIndex) && U <= myKnots->Value(myIndex+1)) break;
162	if (myIndex > myNbElements) myIndex = myNbElements;
163	}
164	Uf = myKnots->Value(myIndex);
165	Ul = myKnots->Value(myIndex+1);
166	Denom = 1. /(Ul-Uf);
167	USum = Uf+Ul;
168	myPtr = (myIndex - 1)(myBase->WorkDegree()+ 1)myDimension + 1;
169	}
170
171	deg = myDegree(myIndex);
172	if (!HasPoly.Value(myIndex)) Update(myIndex, 0);
173
174	//Parameter normalization: S [-1, 1]
175	S = (2U - USum)Denom;
176	PLib::NoDerivativeEvalPolynomial(S, deg, myDimension, deg*myDimension,
177	myPoly(myPtr), Pnt(Pnt.Lower()));
178	}
179
180
181	//=======================================================================
182	//function : D1
183	//purpose :
184	//=======================================================================
185	void FEmTool_Curve::D1(const Standard_Real U, TColStd_Array1OfReal& Vec)
186	{
187	Standard_Integer deg, i;
188	Standard_Real S;
189
190	if (!myIndex \|\| (U<Uf) \|\| (U>Ul) \|\|
191	(myKnots->Value(myIndex)!=Uf) \|\|
192	(myKnots->Value(myIndex+1)!=Ul)) {
193	// Search the span
194	if(U <= myKnots->Value(2)) myIndex = 1;
195	else {
196	for(myIndex = 2; myIndex <= myNbElements; myIndex++)
197	if(U >= myKnots->Value(myIndex) && U <= myKnots->Value(myIndex+1)) break;
198	if (myIndex > myNbElements) myIndex = myNbElements;
199	}
200	Uf = myKnots->Value(myIndex);
201	Ul = myKnots->Value(myIndex+1);
202	Denom = 1. /(Ul-Uf);
203	USum = Uf+Ul;
204	myPtr = (myIndex - 1)(myBase->WorkDegree()+ 1)myDimension + 1;
205	}
206
207	deg = myDegree(myIndex);
208	if (!HasDeri.Value(myIndex)) Update(myIndex, 1);
209
210	//Parameter normalization: S [-1, 1]
211	S = (2U - USum) Denom;
212	PLib::NoDerivativeEvalPolynomial(S, deg-1, myDimension, (deg-1)*myDimension,
213	myDeri(1+(myIndex - 1)myBase->WorkDegree()
214	myDimension),
215	Vec(Vec.Lower()));
216
217	S = 2*Denom;
218	for(i=Vec.Lower(); i<=Vec.Upper(); i++) Vec(i) *= S;
219	}
220
221
222
223	//=======================================================================
224	//function : D2
225	//purpose :
226	//=======================================================================
227	void FEmTool_Curve::D2(const Standard_Real U, TColStd_Array1OfReal& Vec)
228	{
229	Standard_Integer deg, i;
230	Standard_Real S;
231
232	if (!myIndex \|\| (U<Uf) \|\| (U>Ul) \|\|
233	(myKnots->Value(myIndex)!=Uf) \|\|
234	(myKnots->Value(myIndex+1)!=Ul)) {
235	// Search the span
236	if(U <= myKnots->Value(2)) myIndex = 1;
237	else {
238	for(myIndex = 2; myIndex <= myNbElements; myIndex++)
239	if(U >= myKnots->Value(myIndex) && U <= myKnots->Value(myIndex+1)) break;
240	if (myIndex > myNbElements) myIndex = myNbElements;
241	}
242	Uf = myKnots->Value(myIndex);
243	Ul = myKnots->Value(myIndex+1);
244	Denom = 1. /(Ul-Uf);
245	USum = Uf+Ul;
246	myPtr = (myIndex - 1)(myBase->WorkDegree()+ 1)myDimension + 1;
247	}
248
249	deg = myDegree(myIndex);
250	if (!HasSecn.Value(myIndex)) Update(myIndex, 2);
251
252	//Parameter normalization: S [-1, 1]
253	S = (2U - USum)Denom;
254	PLib::NoDerivativeEvalPolynomial(S, deg-2, myDimension, (deg-2)*myDimension,
255	myDsecn(1+(myIndex - 1)*
256	(myBase->WorkDegree() - 1)*myDimension),
257	Vec(Vec.Lower()));
258
259	S = 4DenomDenom;
260	for(i=Vec.Lower(); i<=Vec.Upper(); i++) Vec(i) *= S;
261	}
262
263
264
265	//=======================================================================
266	//function : Length
267	//purpose :
268	//=======================================================================
269	void FEmTool_Curve::Length(const Standard_Real FirstU,
270	const Standard_Real LastU,
271	Standard_Real& Length)
272	{
273	Standard_Integer Low, High, deg, degBase, i, Ptr;
274	if(FirstU > LastU) Standard_OutOfRange::Raise("FEmTool_Curve::Length");
275
276	if(myKnots->Value(1) > FirstU) Low = 1;
277	else
278	for(Low = 1; Low <= myNbElements; Low++)
279	if(FirstU >= myKnots->Value(Low) && FirstU <= myKnots->Value(Low+1)) break;
280	if(Low > myNbElements) Low = myNbElements;
281
282	if(myKnots->Value(1) > LastU) High = 1;
283	else
284	for(High = Low; High <= myNbElements; High++)
285	if(LastU >= myKnots->Value(High) && LastU <= myKnots->Value(High+1)) break;
286	if(myKnots->Value(myNbElements + 1) < LastU) High = myNbElements;
287
288	Standard_Real Li;
289	degBase = myBase->WorkDegree();
290	Length = 0;
291
292	Standard_Real FirstS, LastS;
293	FirstS = (2*FirstU - myKnots->Value(Low) - myKnots->Value(Low + 1))/
294	(myKnots->Value(Low + 1) - myKnots->Value(Low));
295	LastS = (2*LastU - myKnots->Value(High) - myKnots->Value(High + 1))/
296	(myKnots->Value(High + 1) - myKnots->Value(High));
297
298	if(Low == High) {
299
300	Ptr = (Low - 1)(degBase + 1)myDimension + 1;
301	deg = myDegree(Low);
302
303	if(!HasPoly(Low)) Update(Low, 0);
304	PLib::EvalLength(deg, myDimension, myPoly(Ptr),
305	FirstS, LastS, Length);
306	return;
307	}
308
309	deg = myDegree(Low);
310	Ptr = (Low - 1)(degBase + 1)myDimension + 1;
311
312	if(!HasPoly(Low)) Update(Low, 0);
313	if(FirstS < -1.) {
314	PLib::EvalLength(deg, myDimension, myPoly(Ptr), FirstS, -1., Li);
315	Length += Li;
316	if(myLength(Low) < 0.) {
317	PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., 1.,Li);
318	myLength(Low) = Li;
319	}
320	Length += myLength(Low);
321	}
322	else {
323	PLib::EvalLength(deg, myDimension, myPoly(Ptr), FirstS, 1., Li);
324	Length += Li;
325	}
326
327	deg = myDegree(High);
328	Ptr = (High - 1)(degBase + 1)myDimension + 1;
329
330	if(!HasPoly(High)) Update(High, 0);
331	if(LastS > 1.) {
332	PLib::EvalLength(deg, myDimension, myPoly(Ptr), 1., LastS, Li);
333	Length += Li;
334	if(myLength(High) < 0.) {
335	PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., 1., Li);
336	myLength(High) = Li;
337	}
338	Length += myLength(High);
339	}
340	else {
341	PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., LastS, Li);
342	Length += Li;
343	}
344
345
346	for(i = Low + 1; i < High; i++) {
347	if (myLength.Value(i) < 0) {
348	Ptr = (i - 1)(degBase + 1)myDimension + 1;
349	deg = myDegree(i);
350	if(!HasPoly(i)) Update(i, 0);
351	PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., 1., Li);
352	myLength(i) = Li;
353	}
354	Length += myLength(i);
355	}
356	}
357
358	Standard_Integer FEmTool_Curve::NbElements() const
359	{
360	return myNbElements;
361	}
362
363	Standard_Integer FEmTool_Curve::Dimension() const
364	{
365	return myDimension;
366	}
367
368	Handle(PLib_Base) FEmTool_Curve::Base() const
369	{
370	return myBase;
371	}
372
373	Standard_Integer FEmTool_Curve::Degree(const Standard_Integer IndexOfElement) const
374	{
375	return myDegree.Value(IndexOfElement);
376	}
377
378
379	//=======================================================================
380	//function : SetDegree
381	//purpose :
382	//=======================================================================
383	void FEmTool_Curve::SetDegree(const Standard_Integer IndexOfElement,
384	const Standard_Integer Degree)
385	{
386	if (Degree <= myBase->WorkDegree()) {
387	myDegree(IndexOfElement) = Degree;
388	HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
389	myLength(IndexOfElement) = -1;
390	}
391	else if(Degree > myBase->WorkDegree()) Standard_OutOfRange::Raise("FEmTool_Curve::SetDegree");
392	}
393
394
395
396	//=======================================================================
397	//function : ReduceDegree
398	//purpose :
399	//=======================================================================
400	void FEmTool_Curve::ReduceDegree(const Standard_Integer IndexOfElement,
401	const Standard_Real Tol,
402	Standard_Integer& NewDegree,
403	Standard_Real& MaxError)
404	{
405	Standard_Integer deg = myDegree(IndexOfElement);
406
407	Standard_Integer Ptr = (IndexOfElement - 1)*
408	(myBase->WorkDegree() + 1)*myDimension + 1;
409
410	myBase->ReduceDegree(myDimension, deg, Tol, myCoeff.ChangeValue(Ptr), NewDegree, MaxError);
411	Handle(PLib_HermitJacobi) myHermitJacobi = (((Handle(PLib_HermitJacobi))&myBase));
412
413	NewDegree = Max(NewDegree, 2 * myHermitJacobi->NivConstr() + 1);
414
415	if (NewDegree < deg) {
416	myDegree(IndexOfElement) = NewDegree;
417	HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
418	myLength(IndexOfElement) = -1;
419	}
420	}
421
422	//=======================================================================
423	//function : Update
424	//purpose :
425	//=======================================================================
426	void FEmTool_Curve::Update(const Standard_Integer Index,
427	const Standard_Integer Order)
428	{
429	Standard_Integer degBase = myBase->WorkDegree(), deg = myDegree(Index);
430
431	if (!HasPoly(Index)) {
432	Standard_Integer Ptr = (Index - 1)(degBase + 1)myDimension + 1;
433
434	TColStd_Array1OfReal Coeff(myPoly.ChangeValue(Ptr), 0, myDimension*(deg + 1) - 1);
435	TColStd_Array1OfReal BaseCoeff(myCoeff.ChangeValue(Ptr), 0, myDimension*(deg + 1) - 1);
436
437	myBase->ToCoefficients(myDimension, deg, BaseCoeff, Coeff);
438	HasPoly(Index) = 1;
439	}
440
441	if (Order >=1) {
442	Standard_Integer i1 = (Index - 1)(degBase)myDimension - myDimension,
443	i2 = (Index - 1)(degBase + 1)myDimension;
444	Standard_Integer i,j;
445	if (!HasDeri.Value(Index)) {
446	for(i = 1; i <= deg; i++) {
447	i1 += myDimension; i2 += myDimension;
448	for(j = 1; j<= myDimension; j++)
449	myDeri(i1 + j) = i*myPoly.Value(i2 + j);
450	}
451	HasDeri(Index) = 1;
452	}
453	if ((Order>=2) &&!HasSecn.Value(Index)) {
454	i1 = (Index - 1)(degBase-1)myDimension - myDimension,
455	i2 = (Index - 1)(degBase)myDimension;
456	for(i = 1; i < deg; i++) {
457	i1 += myDimension; i2 += myDimension;
458	for(j = 1; j<= myDimension; j++)
459	myDsecn(i1 + j) = i*myDeri.Value(i2 + j);
460	}
461	HasSecn(Index) = 1;
462	}
463	}
464	}