0024428: Implementation of LGPL license

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

// Created on: 1997-09-17
// Created by: Philippe MANGIN
// 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.

//====================== Private Methodes =============================//
//=======================================================================
//function : Adjusting
//purpose  : Smoothing's  adjusting like STRIM routine "MAJLIS"
//=======================================================================
void AppParCurves_Variational::Adjusting(
			 Handle(AppParCurves_SmoothCriterion)& J,
			 Standard_Real& WQuadratic,
			 Standard_Real& WQuality,
			 Handle(FEmTool_Curve)& TheCurve,
			 TColStd_Array1OfReal& Ecarts) 
{

//  cout << "=========== Adjusting =============" << endl;

    /* Initialized data */

  const Standard_Integer mxiter = 2;
  const Standard_Real eps1 = 1e-6;
  Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1;
  Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
  Standard_Real CurvTol = eps1 * J->EstLength() / NbrPnt;


    /* Local variables */
  Standard_Integer iter, ipnt;
  Standard_Real ecart, emold, erold, tpara;
  Standard_Real vocri[4], j1cibl, vtest, vseuil;
  Standard_Integer i, numint, flag;
  TColStd_Array1OfReal tbpoid(myFirstPoint, myLastPoint);
  Standard_Boolean loptim, lrejet;
  Handle(AppParCurves_SmoothCriterion) JNew;
  Handle(FEmTool_Curve) CNew;
  Standard_Real E1, E2, E3;
  

/* (0.b) Initialisations */

  loptim = Standard_True;
  iter = 0;
  tbpoid.Init(1.);


/* ============   boucle sur le moteur de lissage  ============== */

  vtest = WQuality * .9;
  j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint));
  
  while(loptim) {
    
    ++iter;
    
/*     (1) Sauvegarde de l'etat precedents */

    vocri[0] = myCriterium[0];
    vocri[1] = myCriterium[1];
    vocri[2] = myCriterium[2];
    vocri[3] = myCriterium[3];
    erold = myMaxError; 
    emold = myAverageError; 

/*     (2) Augmentation du poids des moindre carre */

    if (j1cibl > vtest) {
      WQuadratic = j1cibl / vtest * WQuadratic;
    }

/*     (3) Augmentation du poid associe aux points a problemes */

    vseuil = WQuality * .88;
    
    for (ipnt = myFirstPoint; ipnt <= myLastPoint; ++ipnt) {
      if (Ecarts(ipnt) > vtest) {
	ecart = (Ecarts(ipnt) - vseuil) / WQuality;
	tbpoid(ipnt) = (ecart * 3 + 1.) * tbpoid(ipnt);
      }
    }

/*     (4) Decoupe force */

    if (TheCurve->NbElements() < myMaxSegment && myWithCutting) {

      numint = NearIndex(myParameters->Value(myMaxErrorIndex), TheCurve->Knots(), 0, flag);

      tpara = (TheCurve->Knots()(numint) + TheCurve->Knots()(numint + 1) + 
	       myParameters->Value(myMaxErrorIndex) * 2) / 4;
      
      CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements() + 1, 
			       TheCurve->Base(), CurvTol);

      for(i = 1; i <= numint; i++) CNew->Knots()(i) = TheCurve->Knots()(i);
      for(i = numint + 1; i <= TheCurve->Knots().Length(); i++) 
	CNew->Knots()(i + 1) = TheCurve->Knots()(i);

      CNew->Knots()(numint + 1) = tpara;

    } else {

      CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements(), TheCurve->Base(), CurvTol);

      CNew->Knots() = TheCurve->Knots();
    }


    JNew = new AppParCurves_MyCriterion(mySSP, myFirstPoint, myLastPoint);

    JNew->EstLength() = J->EstLength();

    J->GetEstimation(E1, E2, E3);

    JNew->SetEstimation(E1, E2, E3);

    JNew->SetParameters(myParameters);

    JNew->SetWeight(WQuadratic, WQuality, myPercent[0], myPercent[1], myPercent[2]);

    JNew->SetWeight(tbpoid);

    JNew->SetCurve(CNew);

/*     (5) Relissage */

    TheMotor(JNew, WQuadratic, WQuality, CNew, Ecarts);

/*     (6) Tests de rejet */

    j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint));
    vseuil = Sqrt(vocri[1]) + (erold - myMaxError) * 4;
    
    lrejet = ((myMaxError > WQuality) && (myMaxError > erold * 1.01)) || (Sqrt(myCriterium[1]) > vseuil * 1.05);

    if (lrejet) {
      myCriterium[0] = vocri[0];
      myCriterium[1] = vocri[1];
      myCriterium[2] = vocri[2];
      myCriterium[3] = vocri[3];
      myMaxError = erold;
      myAverageError = emold;
      
      loptim = Standard_False;
    }
    else {
      J = JNew;
      TheCurve = CNew;
      J->SetCurve(TheCurve);
    }

/*     (7) Test de convergence */

    if (((iter >= mxiter) && (myMaxSegment == CNew->NbElements())) || myMaxError < WQuality) {
      loptim = Standard_False;
    }
    
  }
  
  
}