0024624: Lost word in license statement in source files

[occt.git] / src / AppParCurves / AppParCurves_Variational_2.gxx

// Created on: 1998-12-07
// Created by: Igor FEOKTISTOV
// Copyright (c) 1998-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#include <math_Matrix.hxx>
#include <math_Vector.hxx>
#include <TColStd_Array2OfReal.hxx>
#include <PLib_Base.hxx>
#include <PLib_JacobiPolynomial.hxx>
#include <PLib_HermitJacobi.hxx>
#include <FEmTool_HAssemblyTable.hxx>



//=======================================================================
//function : Optimization
//purpose  :  (like FORTRAN subroutine MINIMI)
//=======================================================================
void AppParCurves_Variational::Optimization(Handle(AppParCurves_SmoothCriterion)& J,
					    FEmTool_Assembly& A,
					    const Standard_Boolean ToAssemble,
					    const Standard_Real EpsDeg,
					    Handle(FEmTool_Curve)& Curve,
					    const TColStd_Array1OfReal& Parameters) const
{
  Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
                   NbElm = Curve->NbElements(),
                   NbDim = Curve->Dimension();

  Handle(FEmTool_HAssemblyTable) AssTable;

  math_Matrix H(0, MxDeg, 0, MxDeg);
  math_Vector G(0, MxDeg), Sol(1, A.NbGlobVar());

  Standard_Integer el, dim;

  A.GetAssemblyTable(AssTable);
  Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;

  Standard_Real CBLONG = J->EstLength();

  // Updating Assembly
  if (ToAssemble) 
    A.NullifyMatrix();
  A.NullifyVector();

  
  for (el = 1; el <= NbElm; el++) {
    if (ToAssemble) {
      J->Hessian(el, 1, 1, H);
      for(dim = 1; dim <= NbDim; dim++)
	A.AddMatrix(el, dim, dim, H);
    }

  for(dim = 1; dim <= NbDim; dim++) {
      J->Gradient(el, dim, G);
      A.AddVector(el, dim, G);
    }
  }


  // Solution of system 
  if (ToAssemble) {
    if(NbConstr != 0) { //Treatment of constraints
      AssemblingConstraints(Curve, Parameters, CBLONG, A);
    }      
    A.Solve();
  }
  A.Solution(Sol);


  // Updating J
  J->SetCurve(Curve);
  J->InputVector(Sol, AssTable);

  // Updating Curve and reduction of degree
  
  Standard_Integer Newdeg;
  Standard_Real MaxError;
  
  if(NbConstr == 0) {
    for(el = 1; el <= NbElm; el++) {
      Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError);
    }
  }
  else {
    
    TColStd_Array1OfReal& TabInt = Curve->Knots();
    Standard_Integer Icnt = 1, p0 = Parameters.Lower() - myFirstPoint, point;
    for(el = 1; el <= NbElm; el++) {     
      while((Icnt < NbConstr) && 
	    (Parameters(p0 + myTypConstraints->Value(2 * Icnt - 1)) <= TabInt(el))) Icnt++;
      point = p0 + myTypConstraints->Value(2 * Icnt - 1);
      if(Parameters(point) <= TabInt(el) || Parameters(point) >= TabInt(el + 1))
	Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError);
      else
	if(Curve->Degree(el) < MxDeg) Curve->SetDegree(el, MxDeg);
    }
  }   
}