OCCT

// 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 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.

#define No_Standard_RangeError
#define No_Standard_OutOfRange

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

IMPLEMENT_STANDARD_RTTIEXT(FEmTool_Curve, Standard_Transient)

//=================================================================================================

FEmTool_Curve::FEmTool_Curve(const int                Dimension,
                             const int                NbElements,
                             const PLib_HermitJacobi& TheBase,
                             const double)
    : 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 NCollection_HArray1<double>(1, myNbElements + 1);
  myDegree.Init(myBase.WorkDegree());
  HasPoly.Init(0);
  HasDeri.Init(0);
  HasSecn.Init(0);
  myLength.Init(-1);
}

NCollection_Array1<double>& FEmTool_Curve::Knots() const
{
  return myKnots->ChangeArray1();
}

//=================================================================================================

void FEmTool_Curve::SetElement(const int IndexOfElement, const NCollection_Array2<double>& Coeffs)
{
  int i, j, degBase, deg;
  if (IndexOfElement > myNbElements || IndexOfElement < 1)
  {
    throw Standard_OutOfRange();
  }
  degBase   = myBase.WorkDegree();
  deg       = myDegree(IndexOfElement);
  int 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);
    }
  }

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

  i1 = iBase;
  i2 = iBase + (myBase.NivConstr() + 1) * myDimension;
  for (i = 1; i <= myBase.NivConstr(); i++)
  {
    i1 += myDimension;
    i2 += myDimension;
    mfact = std::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;
}

//=================================================================================================

void FEmTool_Curve::GetElement(const int IndexOfElement, NCollection_Array2<double>& Coeffs)
{
  int i, j, degBase, deg;
  if (IndexOfElement > myNbElements || IndexOfElement < 1)
  {
    throw Standard_OutOfRange();
  }
  degBase   = myBase.WorkDegree();
  deg       = myDegree(IndexOfElement);
  int 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);
    }
  }

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

  i2     = Coeffs.LowerRow();
  int i3 = i2 + myBase.NivConstr() + 1;

  for (i = 1; i <= myBase.NivConstr(); i++)
  {
    mfact = std::pow(stenor, i);
    for (j = j1 + 1; j <= myDimension; j++)
    {
      Coeffs(i2 + i, j) *= mfact;
      Coeffs(i3 + i, j) *= mfact;
    }
  }
}

//=================================================================================================

void FEmTool_Curve::GetPolynom(NCollection_Array1<double>& Coeffs)

{
  int 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);
  }
}

//=================================================================================================

void FEmTool_Curve::D0(const double U, NCollection_Array1<double>& Pnt)
{
  int    deg;
  double 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()));
}

//=================================================================================================

void FEmTool_Curve::D1(const double U, NCollection_Array1<double>& Vec)
{
  int    deg, i;
  double 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;
  }
}

//=================================================================================================

void FEmTool_Curve::D2(const double U, NCollection_Array1<double>& Vec)
{
  int    deg, i;
  double 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;
  }
}

//=================================================================================================

void FEmTool_Curve::Length(const double FirstU, const double LastU, double& Length)
{
  int Low, High, deg, degBase, i, Ptr;
  if (FirstU > LastU)
  {
    throw Standard_OutOfRange("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;
  }

  double Li;
  degBase = myBase.WorkDegree();
  Length  = 0;

  double 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);
  }
}

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

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

const PLib_HermitJacobi& FEmTool_Curve::Base() const
{
  return myBase;
}

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

//=================================================================================================

void FEmTool_Curve::SetDegree(const int IndexOfElement, const int Degree)
{
  if (Degree <= myBase.WorkDegree())
  {
    myDegree(IndexOfElement) = Degree;
    HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
    myLength(IndexOfElement)                                                    = -1;
  }
  else if (Degree > myBase.WorkDegree())
  {
    throw Standard_OutOfRange("FEmTool_Curve::SetDegree");
  }
}

//=================================================================================================

void FEmTool_Curve::ReduceDegree(const int    IndexOfElement,
                                 const double Tol,
                                 int&         NewDegree,
                                 double&      MaxError)
{
  int deg = myDegree(IndexOfElement);

  int Ptr = (IndexOfElement - 1) * (myBase.WorkDegree() + 1) * myDimension + 1;

  myBase.ReduceDegree(myDimension, deg, Tol, myCoeff.ChangeValue(Ptr), NewDegree, MaxError);

  NewDegree = std::max(NewDegree, 2 * myBase.NivConstr() + 1);

  if (NewDegree < deg)
  {
    myDegree(IndexOfElement) = NewDegree;
    HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
    myLength(IndexOfElement)                                                    = -1;
  }
}

//=================================================================================================

void FEmTool_Curve::Update(const int Index, const int Order)
{
  int degBase = myBase.WorkDegree(), deg = myDegree(Index);

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

    NCollection_Array1<double> Coeff(myPoly.ChangeValue(Ptr), 0, myDimension * (deg + 1) - 1);
    NCollection_Array1<double> BaseCoeff(myCoeff.ChangeValue(Ptr), 0, myDimension * (deg + 1) - 1);

    myBase.ToCoefficients(myDimension, deg, BaseCoeff, Coeff);
    HasPoly(Index) = 1;
  }

  if (Order >= 1)
  {
    int i1 = (Index - 1) * (degBase)*myDimension - myDimension,
        i2 = (Index - 1) * (degBase + 1) * myDimension;
    int 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;
    }
  }
}