[occt.git] / src / FairCurve / FairCurve_DistributionOfTension.cxx

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

// 30-01-1996 : PMN Version originale

#ifndef DEB
#define No_Standard_RangeError
#define No_Standard_OutOfRange
#endif

#include <FairCurve_DistributionOfTension.ixx>


#include <gp_XY.hxx>
#include <gp_Pnt2d.hxx>
#include <math_Matrix.hxx>
#include <math_Vector.hxx>
#include <BSplCLib.hxx>


 FairCurve_DistributionOfTension::FairCurve_DistributionOfTension(const Standard_Integer BSplOrder,
								  const Handle(TColStd_HArray1OfReal)& FlatKnots, 
								  const Handle(TColgp_HArray1OfPnt2d)& Poles,
								  const Standard_Integer DerivativeOrder,
								  const Standard_Real LengthSliding, 
								  const FairCurve_BattenLaw& Law, 
								  const Standard_Integer NbValAux,
								  const Standard_Boolean Uniform ) :
                                   FairCurve_DistributionOfEnergy(BSplOrder,  
								  FlatKnots, 
								  Poles, 
								  DerivativeOrder,
								  NbValAux) ,
								  MyLengthSliding (LengthSliding),  
								  MyLaw (Law)
{
 if (Uniform) {MyLaw.Value(0.5, MyHeight);} // it used in MVC to avoid Parametrization Problemes
 else         {MyHeight = 0;}
}


Standard_Boolean FairCurve_DistributionOfTension::Value(const math_Vector& TParam, math_Vector& FTension)
{
  Standard_Boolean Ok = Standard_True;
  Standard_Integer ier, ii, jj, kk;
  gp_XY CPrim  (0., 0.);
  Standard_Integer  LastGradientIndex, FirstNonZero, LastZero; 

  // (0.0) initialisations generales
  FTension.Init(0.0);
  math_Matrix Base(1, 3, 1, MyBSplOrder ); // On shouhaite utiliser la derive premieres
                                           // Dans EvalBsplineBasis C' <=> DerivOrder = 2
                                           // et il faut ajouter 1 rang dans la matrice Base => 3 rang
   
  ier  =  BSplCLib::EvalBsplineBasis( 1, 1,  MyBSplOrder, 
                                    MyFlatKnots->Array1(), TParam(TParam.Lower()),
                                    FirstNonZero, Base );
  if (ier != 0) return Standard_False;
  LastZero = FirstNonZero - 1;
  FirstNonZero = 2*LastZero+1;

  // (0.1) evaluation de CPrim
  for (ii= 1; ii<= MyBSplOrder; ii++) {
      CPrim += Base(2, ii) * MyPoles->Value(ii+LastZero).Coord();
      }

  // (1) Evaluation de la tension locale --------------------------------
  Standard_Real NormeCPrim = CPrim.Modulus();
  Standard_Real Hauteur, Difference;

  if (MyHeight > 0) {Hauteur = MyHeight;} // it used in MVC to avoid Parametrization Problemes
  else { 
    Ok = MyLaw.Value (TParam(TParam.Lower()), Hauteur); 
    if (!Ok) return Ok;
  }
  Difference = NormeCPrim - MyLengthSliding;
 
  FTension(FTension.Lower()) = Hauteur * pow(Difference, 2) / MyLengthSliding ;

  if (MyDerivativeOrder >= 1) {
  // (2) Evaluation du gradient de la tension locale ----------------------
      math_Vector GradDifference (1, 2*MyBSplOrder+MyNbValAux);
      Standard_Real Xaux, Yaux, Facteur;

      Xaux = CPrim.X() / NormeCPrim;
      Yaux = CPrim.Y() / NormeCPrim;
      Facteur = 2 * Hauteur * Difference / MyLengthSliding;

      kk = FTension.Lower() + FirstNonZero;
      jj = 1;
      for (ii=1; ii<= MyBSplOrder; ii++) {
	     GradDifference(jj)   = Base(2, ii) * Xaux;
             FTension(kk) = Facteur *  GradDifference(jj);
             jj +=1;
	     GradDifference(jj) = Base(2, ii) * Yaux;
             FTension(kk+1) = Facteur *  GradDifference(jj);
             jj += 1;
             kk += 2;
      }
      if (MyNbValAux == 1) {
         LastGradientIndex = FTension.Lower() + 2*MyPoles->Length() + 1;
         GradDifference( GradDifference.Upper()) =  (1 - pow( NormeCPrim/MyLengthSliding, 2));
         FTension(LastGradientIndex) = Hauteur * GradDifference(GradDifference.Upper());        
       }

      else { LastGradientIndex = FTension.Lower() + 2*MyPoles->Length(); }


      if (MyDerivativeOrder >= 2) { 
   
// (3) Evaluation du Hessien de la tension locale ----------------------

         Standard_Real FacteurX =  Difference * (1-pow(Xaux,2)) / NormeCPrim;
         Standard_Real FacteurY =  Difference * (1-pow(Yaux,2)) / NormeCPrim;
         Standard_Real FacteurXY = - Difference * Xaux*Yaux / NormeCPrim;
         Standard_Real Produit;
         Standard_Integer k1, k2;
   
         Facteur = 2 * Hauteur / MyLengthSliding;

         kk = FirstNonZero;
         k2 = LastGradientIndex + (kk-1)*kk/2;

         for (ii=2; ii<= 2*MyBSplOrder; ii+=2) {
	   k1 = k2+FirstNonZero;
           k2 = k1+kk;
           kk += 2;              
           for (jj=2; jj< ii; jj+=2) {
             Produit =  Base(2, ii/2) *  Base(2, jj/2);

	     FTension(k1) = Facteur * ( GradDifference(ii-1)*GradDifference(jj-1)
                                      + FacteurX * Produit) ;  // derivation en XiXj
             k1++;
	     FTension(k1) = Facteur * ( GradDifference(ii)*GradDifference(jj-1)	
                                      + FacteurXY * Produit); // derivation en YiXj
	     k1++;
	     FTension(k2) = Facteur * ( GradDifference(ii-1)*GradDifference(jj)
                                    + FacteurXY * Produit); // derivation en XiYj
             k2++;
	     FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(jj)
                                      + FacteurY * Produit); // derivation en YiYj
	     k2++;
	     }
         // cas ou jj = ii : remplisage en triangle
             Produit =  pow (Base(2, ii/2), 2);

	     FTension(k1) = Facteur * ( GradDifference(ii-1)*GradDifference(ii-1)
                                      + FacteurX * Produit) ;  // derivation en XiXi
	     FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(ii-1)
                                      + FacteurXY * Produit); // derivation en XiYi
             k2++;
	     FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(ii)
                                      + FacteurY * Produit); // derivation en YiYi
	 }
         if (MyNbValAux == 1) {
           FacteurX = -2*CPrim.X()*Hauteur / pow (MyLengthSliding, 2);
           FacteurY = -2*CPrim.Y()*Hauteur / pow (MyLengthSliding, 2);

           ii = LastGradientIndex-FTension.Lower();
	   kk = LastGradientIndex + (ii-1)*ii/2 +  FirstNonZero;
	   for (ii=1; ii<= MyBSplOrder; ii++) {
             FTension(kk) = FacteurX * Base(2, ii);
             kk ++;
             FTension(kk) = FacteurY * Base(2, ii);
             kk ++;
	   }
           FTension(FTension.Upper()) = 2 * Hauteur * pow (NormeCPrim/MyLengthSliding, 2)
	                              / MyLengthSliding;
	 }
       }
    }
          
// sortie standard           
  return Ok;
}
Commit	Line	Data
973c2be1	1	// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e	2	//
973c2be1	3	// This file is part of Open CASCADE Technology software library.
b311480e	4	//
d5f74e42	5	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	6	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	7	// by the Free Software Foundation, with special exception defined in the file
	8	// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
	9	// distribution for complete text of the license and disclaimer of any warranty.
b311480e	10	//
973c2be1	11	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	12	// commercial license or contractual agreement.
b311480e	13
7fd59977	14	// 30-01-1996 : PMN Version originale
	15
	16	#ifndef DEB
	17	#define No_Standard_RangeError
	18	#define No_Standard_OutOfRange
	19	#endif
	20
	21	#include <FairCurve_DistributionOfTension.ixx>
	22
	23
	24	#include <gp_XY.hxx>
	25	#include <gp_Pnt2d.hxx>
	26	#include <math_Matrix.hxx>
	27	#include <math_Vector.hxx>
	28	#include <BSplCLib.hxx>
	29
	30
	31
	32	FairCurve_DistributionOfTension::FairCurve_DistributionOfTension(const Standard_Integer BSplOrder,
	33	const Handle(TColStd_HArray1OfReal)& FlatKnots,
	34	const Handle(TColgp_HArray1OfPnt2d)& Poles,
	35	const Standard_Integer DerivativeOrder,
	36	const Standard_Real LengthSliding,
	37	const FairCurve_BattenLaw& Law,
	38	const Standard_Integer NbValAux,
	39	const Standard_Boolean Uniform ) :
	40	FairCurve_DistributionOfEnergy(BSplOrder,
	41	FlatKnots,
	42	Poles,
	43	DerivativeOrder,
	44	NbValAux) ,
	45	MyLengthSliding (LengthSliding),
	46	MyLaw (Law)
	47	{
	48	if (Uniform) {MyLaw.Value(0.5, MyHeight);} // it used in MVC to avoid Parametrization Problemes
	49	else {MyHeight = 0;}
	50	}
	51
	52
	53
	54	Standard_Boolean FairCurve_DistributionOfTension::Value(const math_Vector& TParam, math_Vector& FTension)
	55	{
	56	Standard_Boolean Ok = Standard_True;
	57	Standard_Integer ier, ii, jj, kk;
	58	gp_XY CPrim (0., 0.);
	59	Standard_Integer LastGradientIndex, FirstNonZero, LastZero;
	60
	61	// (0.0) initialisations generales
	62	FTension.Init(0.0);
	63	math_Matrix Base(1, 3, 1, MyBSplOrder ); // On shouhaite utiliser la derive premieres
	64	// Dans EvalBsplineBasis C' <=> DerivOrder = 2
	65	// et il faut ajouter 1 rang dans la matrice Base => 3 rang
	66
	67	ier = BSplCLib::EvalBsplineBasis( 1, 1, MyBSplOrder,
	68	MyFlatKnots->Array1(), TParam(TParam.Lower()),
	69	FirstNonZero, Base );
	70	if (ier != 0) return Standard_False;
	71	LastZero = FirstNonZero - 1;
	72	FirstNonZero = 2*LastZero+1;
	73
	74	// (0.1) evaluation de CPrim
	75	for (ii= 1; ii<= MyBSplOrder; ii++) {
	76	CPrim += Base(2, ii) * MyPoles->Value(ii+LastZero).Coord();
	77	}
78
79	// (1) Evaluation de la tension locale --------------------------------
80	Standard_Real NormeCPrim = CPrim.Modulus();
81	Standard_Real Hauteur, Difference;
82
83	if (MyHeight > 0) {Hauteur = MyHeight;} // it used in MVC to avoid Parametrization Problemes
84	else {
85	Ok = MyLaw.Value (TParam(TParam.Lower()), Hauteur);
86	if (!Ok) return Ok;
87	}
88	Difference = NormeCPrim - MyLengthSliding;
89
90	FTension(FTension.Lower()) = Hauteur * pow(Difference, 2) / MyLengthSliding ;
91
92	if (MyDerivativeOrder >= 1) {
93	// (2) Evaluation du gradient de la tension locale ----------------------
94	math_Vector GradDifference (1, 2*MyBSplOrder+MyNbValAux);
95	Standard_Real Xaux, Yaux, Facteur;
96
97	Xaux = CPrim.X() / NormeCPrim;
98	Yaux = CPrim.Y() / NormeCPrim;
99	Facteur = 2 * Hauteur * Difference / MyLengthSliding;
100
101	kk = FTension.Lower() + FirstNonZero;
102	jj = 1;
103	for (ii=1; ii<= MyBSplOrder; ii++) {
104	GradDifference(jj) = Base(2, ii) * Xaux;
105	FTension(kk) = Facteur * GradDifference(jj);
106	jj +=1;
107	GradDifference(jj) = Base(2, ii) * Yaux;
108	FTension(kk+1) = Facteur * GradDifference(jj);
109	jj += 1;
110	kk += 2;
111	}
112	if (MyNbValAux == 1) {
113	LastGradientIndex = FTension.Lower() + 2*MyPoles->Length() + 1;
114	GradDifference( GradDifference.Upper()) = (1 - pow( NormeCPrim/MyLengthSliding, 2));
115	FTension(LastGradientIndex) = Hauteur * GradDifference(GradDifference.Upper());
116	}
117
118	else { LastGradientIndex = FTension.Lower() + 2*MyPoles->Length(); }
119
120
121	if (MyDerivativeOrder >= 2) {
122
123	// (3) Evaluation du Hessien de la tension locale ----------------------
124
125	Standard_Real FacteurX = Difference * (1-pow(Xaux,2)) / NormeCPrim;
126	Standard_Real FacteurY = Difference * (1-pow(Yaux,2)) / NormeCPrim;
127	Standard_Real FacteurXY = - Difference * Xaux*Yaux / NormeCPrim;
128	Standard_Real Produit;
129	Standard_Integer k1, k2;
130
131	Facteur = 2 * Hauteur / MyLengthSliding;
132
133	kk = FirstNonZero;
134	k2 = LastGradientIndex + (kk-1)*kk/2;
135
136	for (ii=2; ii<= 2*MyBSplOrder; ii+=2) {
137	k1 = k2+FirstNonZero;
138	k2 = k1+kk;
139	kk += 2;
140	for (jj=2; jj< ii; jj+=2) {
141	Produit = Base(2, ii/2) * Base(2, jj/2);
142
143	FTension(k1) = Facteur * ( GradDifference(ii-1)*GradDifference(jj-1)
144	+ FacteurX * Produit) ; // derivation en XiXj
145	k1++;
146	FTension(k1) = Facteur * ( GradDifference(ii)*GradDifference(jj-1)
147	+ FacteurXY * Produit); // derivation en YiXj
148	k1++;
149	FTension(k2) = Facteur * ( GradDifference(ii-1)*GradDifference(jj)
150	+ FacteurXY * Produit); // derivation en XiYj
151	k2++;
152	FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(jj)
153	+ FacteurY * Produit); // derivation en YiYj
154	k2++;
155	}
156	// cas ou jj = ii : remplisage en triangle
157	Produit = pow (Base(2, ii/2), 2);
158
159	FTension(k1) = Facteur * ( GradDifference(ii-1)*GradDifference(ii-1)
160	+ FacteurX * Produit) ; // derivation en XiXi
161	FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(ii-1)
162	+ FacteurXY * Produit); // derivation en XiYi
163	k2++;
164	FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(ii)
165	+ FacteurY * Produit); // derivation en YiYi
166	}
167	if (MyNbValAux == 1) {
168	FacteurX = -2CPrim.X()Hauteur / pow (MyLengthSliding, 2);
169	FacteurY = -2CPrim.Y()Hauteur / pow (MyLengthSliding, 2);
170
171	ii = LastGradientIndex-FTension.Lower();
172	kk = LastGradientIndex + (ii-1)*ii/2 + FirstNonZero;
173	for (ii=1; ii<= MyBSplOrder; ii++) {
174	FTension(kk) = FacteurX * Base(2, ii);
175	kk ++;
176	FTension(kk) = FacteurY * Base(2, ii);
177	kk ++;
178	}
179	FTension(FTension.Upper()) = 2 * Hauteur * pow (NormeCPrim/MyLengthSliding, 2)
180	/ MyLengthSliding;
181	}
182	}
183	}
184
185	// sortie standard
186	return Ok;
187	}
188