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

// Created on: 1997-11-04
// Created by: Roman BORISOV / Igor FEOKTISTOV
// 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 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.

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