src/AppParCurves/AppParCurves_Variational_2.gxx

   1 // Created on: 1998-12-07
   2 // Created by: Igor FEOKTISTOV
   3 // Copyright (c) 1998-1999 Matra Datavision
   4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
   5 //
   6 // The content of this file is subject to the Open CASCADE Technology Public
   7 // License Version 6.5 (the "License"). You may not use the content of this file
   8 // except in compliance with the License. Please obtain a copy of the License
   9 // at http://www.opencascade.org and read it completely before using this file.
  10 //
  11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  13 //
  14 // The Original Code and all software distributed under the License is
  15 // distributed on an "AS IS" basis, without warranty of any kind, and the
  16 // Initial Developer hereby disclaims all such warranties, including without
  17 // limitation, any warranties of merchantability, fitness for a particular
  18 // purpose or non-infringement. Please see the License for the specific terms
  19 // and conditions governing the rights and limitations under the License.
  20
  21
  22 #include <math_Matrix.hxx>
  23 #include <math_Vector.hxx>
  24 #include <TColStd_Array2OfReal.hxx>
  25 #include <PLib_Base.hxx>
  26 #include <PLib_JacobiPolynomial.hxx>
  27 #include <PLib_HermitJacobi.hxx>
  28 #include <FEmTool_HAssemblyTable.hxx>
  29
  30
  31
  32 //=======================================================================
  33 //function : Optimization
  34 //purpose  :  (like FORTRAN subroutine MINIMI)
  35 //=======================================================================
  36 void AppParCurves_Variational::Optimization(Handle(AppParCurves_SmoothCriterion)& J,
  37                                             FEmTool_Assembly& A,
  38                                             const Standard_Boolean ToAssemble,
  39                                             const Standard_Real EpsDeg,
  40                                             Handle(FEmTool_Curve)& Curve,
  41                                             const TColStd_Array1OfReal& Parameters) const
  42 {
  43   Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
  44                    NbElm = Curve->NbElements(),
  45                    NbDim = Curve->Dimension();
  46
  47   Handle(FEmTool_HAssemblyTable) AssTable;
  48
  49   math_Matrix H(0, MxDeg, 0, MxDeg);
  50   math_Vector G(0, MxDeg), Sol(1, A.NbGlobVar());
  51
  52   Standard_Integer el, dim;
  53
  54   A.GetAssemblyTable(AssTable);
  55   Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
  56
  57   Standard_Real CBLONG = J->EstLength();
  58
  59   // Updating Assembly
  60   if (ToAssemble)
  61     A.NullifyMatrix();
  62   A.NullifyVector();
  63
  64
  65   for (el = 1; el <= NbElm; el++) {
  66     if (ToAssemble) {
  67       J->Hessian(el, 1, 1, H);
  68       for(dim = 1; dim <= NbDim; dim++)
  69         A.AddMatrix(el, dim, dim, H);
  70     }
  71
  72   for(dim = 1; dim <= NbDim; dim++) {
  73       J->Gradient(el, dim, G);
  74       A.AddVector(el, dim, G);
  75     }
  76   }
  77
  78
  79   // Solution of system
  80   if (ToAssemble) {
  81     if(NbConstr != 0) { //Treatment of constraints
  82       AssemblingConstraints(Curve, Parameters, CBLONG, A);
  83     }
  84     A.Solve();
  85   }
  86   A.Solution(Sol);
  87
  88
  89   // Updating J
  90   J->SetCurve(Curve);
  91   J->InputVector(Sol, AssTable);
  92
  93   // Updating Curve and reduction of degree
  94
  95   Standard_Integer Newdeg;
  96   Standard_Real MaxError;
  97
  98   if(NbConstr == 0) {
  99     for(el = 1; el <= NbElm; el++) {
 100       Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError);
 101     }
 102   }
 103   else {
 104
 105     TColStd_Array1OfReal& TabInt = Curve->Knots();
 106     Standard_Integer Icnt = 1, p0 = Parameters.Lower() - myFirstPoint, point;
 107     for(el = 1; el <= NbElm; el++) {
 108       while((Icnt < NbConstr) &&
 109             (Parameters(p0 + myTypConstraints->Value(2 * Icnt - 1)) <= TabInt(el))) Icnt++;
 110       point = p0 + myTypConstraints->Value(2 * Icnt - 1);
 111       if(Parameters(point) <= TabInt(el) || Parameters(point) >= TabInt(el + 1))
 112         Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError);
 113       else
 114         if(Curve->Degree(el) < MxDeg) Curve->SetDegree(el, MxDeg);
 115     }
 116   }
 117 }
 118