[occt.git] / src / Extrema / Extrema_FuncExtCS.cxx

// Created on: 1996-01-09
// Created by: Laurent PAINNOT
// Copyright (c) 1996-1999 Matra Datavision
// Copyright (c) 1999-2012 OPEN CASCADE SAS
//
// The content of this file is subject to the Open CASCADE Technology Public
// License Version 6.5 (the "License"). You may not use the content of this file
// except in compliance with the License. Please obtain a copy of the License
// at http://www.opencascade.org and read it completely before using this file.
//
// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
//
// The Original Code and all software distributed under the License is
// distributed on an "AS IS" basis, without warranty of any kind, and the
// Initial Developer hereby disclaims all such warranties, including without
// limitation, any warranties of merchantability, fitness for a particular
// purpose or non-infringement. Please see the License for the specific terms
// and conditions governing the rights and limitations under the License.


#include <Extrema_FuncExtCS.ixx>
#include <gp_Vec.hxx>
#include <Standard_TypeMismatch.hxx>
#include <Precision.hxx>

/*-----------------------------------------------------------------------------
 Fonction permettant de rechercher une distance extremale entre une courbe C 
et une surface S.
 Cette classe herite de math_FunctionWithDerivative et est utilisee par
les algorithmes math_FunctionRoot et math_FunctionRoots.

{ F1(t,u,v) = (C(t)-S(u,v)).Dtc(t) }
{ F2(t,u,v) = (C(t)-S(u,v)).Dus(u,v) }
{ F3(t,u,v) = (C(t)-S(u,v)).Dvs(u,v) }

{ Dtf1(t,u,v) = Dtc(t).Dtc(t)+(C(t)-S(u,v)).Dttc(t) 
              = ||Dtc(t)||**2+(C(t)-S(u,v)).Dttc(t) }
{ Duf1(t,u,v) = -Dus(u,v).Dtc(t) }
{ Dvf1(t,u,v) = -Dvs(u,v).Dtc(t) }

{ Dtf2(t,u,v) = Dtc(t).Dus(u,v) }
{ Duf2(t,u,v) = -Dus(u,v).Dus(u,v)+(C(t)-S(u,v)).Duus(u,v)
              = -||Dus(u,v)||**2+(C(t)-S(u,v)).Duus(u,v) }
{ Dvf2(t,u,v) = -Dvs(u,v).Dus(u,v)+(C(t)-S(u,v)).Duvs(u,v) }

{ Dtf3(t,u,v) = Dtc(t).Dvs(u,v) }
{ Duf3(t,u,v) = -Dus(u,v).Dvs(u,v)+(C(t)-S(u,v)).Duvs(u,v) }
{ Dvf3(t,u,v) = -Dvs(u,v).Dvs(u,v)+(C(t)-S(u,v)).Dvvs(u,v) }

----------------------------------------------------------------------------*/


//=======================================================================
//function : Extrema_FuncExtCS
//purpose  : 
//=======================================================================

 Extrema_FuncExtCS::Extrema_FuncExtCS()
{
  myCinit = Standard_False;
  mySinit = Standard_False;
}

//=======================================================================
//function : Extrema_FuncExtCS
//purpose  : 
//=======================================================================

 Extrema_FuncExtCS::Extrema_FuncExtCS(const Adaptor3d_Curve& C, 
				      const Adaptor3d_Surface& S)
{
  Initialize(C, S);
}

//=======================================================================
//function : Initialize
//purpose  : 
//=======================================================================

void Extrema_FuncExtCS::Initialize(const Adaptor3d_Curve& C, 
				   const Adaptor3d_Surface& S)
{
  myC = (Adaptor3d_CurvePtr)&C;
  myS = (Adaptor3d_SurfacePtr)&S;
  myCinit = Standard_True;
  mySinit = Standard_True;
  myPoint1.Clear();
  myPoint2.Clear();
  mySqDist.Clear();
}

//=======================================================================
//function : NbVariables
//purpose  : 
//=======================================================================

Standard_Integer Extrema_FuncExtCS::NbVariables() const 
{
  return (3);
}

//=======================================================================
//function : NbEquations
//purpose  : 
//=======================================================================

Standard_Integer Extrema_FuncExtCS::NbEquations() const 
{
  return (3);
}

//=======================================================================
//function : Value
//purpose  : 
//=======================================================================

Standard_Boolean Extrema_FuncExtCS::Value(const math_Vector& UV, 
					  math_Vector& F)
{
  if (!myCinit || !mySinit) Standard_TypeMismatch::Raise();

  myt = UV(1);
  myU = UV(2);
  myV = UV(3);

//  gp_Vec Dtc, Dttc;
  gp_Vec Dtc;
///  gp_Vec Dus, Dvs, Duvs, Duus, Dvvs;
  gp_Vec Dus, Dvs;
  myC->D1(myt, myP1, Dtc);
  myS->D1(myU,myV,myP2,Dus,Dvs);

  gp_Vec P1P2 (myP2,myP1);

  F(1) = P1P2.Dot(Dtc);
  F(2) = P1P2.Dot(Dus);
  F(3) = P1P2.Dot(Dvs);
  
  return  Standard_True;
}

//=======================================================================
//function : Derivatives
//purpose  : 
//=======================================================================

Standard_Boolean Extrema_FuncExtCS::Derivatives(const math_Vector& UV, 
						math_Matrix& DF)
{
  math_Vector F(1,3);
  return Values(UV,F,DF);
}

//=======================================================================
//function : Values
//purpose  : 
//=======================================================================

Standard_Boolean Extrema_FuncExtCS::Values(const math_Vector& UV, 
					   math_Vector& F, 
					   math_Matrix& Df)
{
  if (!myCinit || !mySinit) Standard_TypeMismatch::Raise();

  myt = UV(1);
  myU = UV(2);
  myV = UV(3);

  gp_Vec Dtc, Dttc;
  gp_Vec Dus, Dvs, Duvs, Duus, Dvvs;
  myC->D2(myt, myP1, Dtc, Dttc);
  myS->D2(myU,myV,myP2,Dus,Dvs,Duus,Dvvs,Duvs);

  gp_Vec P1P2 (myP2,myP1);

  F(1) = P1P2.Dot(Dtc);
  F(2) = P1P2.Dot(Dus);
  F(3) = P1P2.Dot(Dvs);

  Df(1,1) = Dtc.SquareMagnitude() + P1P2.Dot(Dttc);
  Df(1,2) = -Dus.Dot(Dtc);
  Df(1,3) = -Dvs.Dot(Dtc);

  Df(2,1) = -Df(1, 2);   // Dtc.Dot(Dus);
  Df(2,2) = -Dus.SquareMagnitude()+P1P2.Dot(Duus);
  Df(2,3) = -Dvs.Dot(Dus)+P1P2.Dot(Duvs);

  Df(3,1) = -Df(1,3);    // Dtc.Dot(Dvs);
  Df(3,2) = Df(2,3);     // -Dus.Dot(Dvs)+P1P2.Dot(Duvs);
  Df(3,3) = -Dvs.SquareMagnitude()+P1P2.Dot(Dvvs);

  return Standard_True;

}

//=======================================================================
//function : GetStateNumber
//purpose  : 
//=======================================================================

Standard_Integer Extrema_FuncExtCS::GetStateNumber()
{
  if (!myCinit || !mySinit) Standard_TypeMismatch::Raise();
#if 0
  math_Vector Sol(1, 3), UVSol(1, 3);
  UVSol(1) = myt; UVSol(2) = myU; UVSol(3) = myV;
  Value(UVSol, Sol);
  cout <<"F(1)= "<<Sol(1)<<" F(2)= "<<Sol(2)<<" F(3)= "<<Sol(3)<<endl;
#endif
  //comparison of solution with previous solutions
  Standard_Real tol2d = Precision::PConfusion() * Precision::PConfusion();
  Standard_Integer i = 1, nbSol = mySqDist.Length();
  for( ; i <=  nbSol; i++)
  {
    Standard_Real aU = myPoint1(i).Parameter();
    if( (myU - aU) * (myU - aU) <= tol2d )
      break;
  }
  if (i <= nbSol)
    return 0;
  mySqDist.Append(myP1.SquareDistance(myP2));
  myPoint1.Append(Extrema_POnCurv(myt,myP1));
  myPoint2.Append(Extrema_POnSurf(myU,myV,myP2));
  return 0;
}

//=======================================================================
//function : NbExt
//purpose  : 
//=======================================================================

Standard_Integer Extrema_FuncExtCS::NbExt() const 
{
  return mySqDist.Length();
}

//=======================================================================
//function : SquareDistance
//purpose  : 
//=======================================================================

Standard_Real Extrema_FuncExtCS::SquareDistance(const Standard_Integer N) const 
{
  if (!myCinit || !mySinit) Standard_TypeMismatch::Raise();
  return mySqDist.Value(N);
}

//=======================================================================
//function : PointOnCurve
//purpose  : 
//=======================================================================

const Extrema_POnCurv& Extrema_FuncExtCS::PointOnCurve(const Standard_Integer N) const 
{
  if (!myCinit || !mySinit) Standard_TypeMismatch::Raise();
  return myPoint1.Value(N);
}

//=======================================================================
//function : PointOnSurface
//purpose  : 
//=======================================================================

const Extrema_POnSurf& Extrema_FuncExtCS::PointOnSurface(const Standard_Integer N) const 
{
  if (!myCinit || !mySinit) Standard_TypeMismatch::Raise();
  return myPoint2.Value(N);
}

//=======================================================================
//function : Bidon1
//purpose  : 
//=======================================================================

Adaptor3d_SurfacePtr Extrema_FuncExtCS::Bidon1() const 
{
  return (Adaptor3d_SurfacePtr)0L;
}

//=======================================================================
//function : Bidon2
//purpose  : 
//=======================================================================

Adaptor3d_CurvePtr Extrema_FuncExtCS::Bidon2() const 
{
  return (Adaptor3d_CurvePtr)0L;
}
Commit	Line	Data
b311480e	1	// Created on: 1996-01-09
	2	// Created by: Laurent PAINNOT
	3	// Copyright (c) 1996-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
7fd59977	21
	22
	23	#include <Extrema_FuncExtCS.ixx>
	24	#include <gp_Vec.hxx>
	25	#include <Standard_TypeMismatch.hxx>
5368adff	26	#include <Precision.hxx>
7fd59977	27
	28	/*-----------------------------------------------------------------------------
	29	Fonction permettant de rechercher une distance extremale entre une courbe C
	30	et une surface S.
	31	Cette classe herite de math_FunctionWithDerivative et est utilisee par
	32	les algorithmes math_FunctionRoot et math_FunctionRoots.
	33
	34	{ F1(t,u,v) = (C(t)-S(u,v)).Dtc(t) }
	35	{ F2(t,u,v) = (C(t)-S(u,v)).Dus(u,v) }
	36	{ F3(t,u,v) = (C(t)-S(u,v)).Dvs(u,v) }
	37
	38	{ Dtf1(t,u,v) = Dtc(t).Dtc(t)+(C(t)-S(u,v)).Dttc(t)
	39	= \|\|Dtc(t)\|\|**2+(C(t)-S(u,v)).Dttc(t) }
	40	{ Duf1(t,u,v) = -Dus(u,v).Dtc(t) }
	41	{ Dvf1(t,u,v) = -Dvs(u,v).Dtc(t) }
	42
	43	{ Dtf2(t,u,v) = Dtc(t).Dus(u,v) }
	44	{ Duf2(t,u,v) = -Dus(u,v).Dus(u,v)+(C(t)-S(u,v)).Duus(u,v)
	45	= -\|\|Dus(u,v)\|\|**2+(C(t)-S(u,v)).Duus(u,v) }
	46	{ Dvf2(t,u,v) = -Dvs(u,v).Dus(u,v)+(C(t)-S(u,v)).Duvs(u,v) }
	47
	48	{ Dtf3(t,u,v) = Dtc(t).Dvs(u,v) }
	49	{ Duf3(t,u,v) = -Dus(u,v).Dvs(u,v)+(C(t)-S(u,v)).Duvs(u,v) }
	50	{ Dvf3(t,u,v) = -Dvs(u,v).Dvs(u,v)+(C(t)-S(u,v)).Dvvs(u,v) }
	51
	52	----------------------------------------------------------------------------*/
	53
	54
	55
	56	//=======================================================================
	57	//function : Extrema_FuncExtCS
	58	//purpose :
	59	//=======================================================================
	60
	61	Extrema_FuncExtCS::Extrema_FuncExtCS()
	62	{
	63	myCinit = Standard_False;
	64	mySinit = Standard_False;
	65	}
	66
	67	//=======================================================================
	68	//function : Extrema_FuncExtCS
	69	//purpose :
	70	//=======================================================================
	71
	72	Extrema_FuncExtCS::Extrema_FuncExtCS(const Adaptor3d_Curve& C,
	73	const Adaptor3d_Surface& S)
	74	{
	75	Initialize(C, S);
	76	}
	77
	78	//=======================================================================
	79	//function : Initialize
	80	//purpose :
	81	//=======================================================================
	82
	83	void Extrema_FuncExtCS::Initialize(const Adaptor3d_Curve& C,
	84	const Adaptor3d_Surface& S)
	85	{
	86	myC = (Adaptor3d_CurvePtr)&C;
	87	myS = (Adaptor3d_SurfacePtr)&S;
	88	myCinit = Standard_True;
	89	mySinit = Standard_True;
	90	myPoint1.Clear();
91	myPoint2.Clear();
92	mySqDist.Clear();
93	}
94
95	//=======================================================================
96	//function : NbVariables
97	//purpose :
98	//=======================================================================
99
100	Standard_Integer Extrema_FuncExtCS::NbVariables() const
101	{
102	return (3);
103	}
104
105	//=======================================================================
106	//function : NbEquations
107	//purpose :
108	//=======================================================================
109
110	Standard_Integer Extrema_FuncExtCS::NbEquations() const
111	{
112	return (3);
113	}
114
115	//=======================================================================
116	//function : Value
117	//purpose :
118	//=======================================================================
119
120	Standard_Boolean Extrema_FuncExtCS::Value(const math_Vector& UV,
121	math_Vector& F)
122	{
123	if (!myCinit \|\| !mySinit) Standard_TypeMismatch::Raise();
124
125	myt = UV(1);
126	myU = UV(2);
127	myV = UV(3);
128
129	// gp_Vec Dtc, Dttc;
130	gp_Vec Dtc;
131	/// gp_Vec Dus, Dvs, Duvs, Duus, Dvvs;
132	gp_Vec Dus, Dvs;
133	myC->D1(myt, myP1, Dtc);
134	myS->D1(myU,myV,myP2,Dus,Dvs);
135
136	gp_Vec P1P2 (myP2,myP1);
137
138	F(1) = P1P2.Dot(Dtc);
139	F(2) = P1P2.Dot(Dus);
140	F(3) = P1P2.Dot(Dvs);
141
142	return Standard_True;
143	}
144
145	//=======================================================================
146	//function : Derivatives
147	//purpose :
148	//=======================================================================
149
150	Standard_Boolean Extrema_FuncExtCS::Derivatives(const math_Vector& UV,
151	math_Matrix& DF)
152	{
153	math_Vector F(1,3);
154	return Values(UV,F,DF);
155	}
156
157	//=======================================================================
158	//function : Values
159	//purpose :
160	//=======================================================================
161
162	Standard_Boolean Extrema_FuncExtCS::Values(const math_Vector& UV,
163	math_Vector& F,
164	math_Matrix& Df)
165	{
166	if (!myCinit \|\| !mySinit) Standard_TypeMismatch::Raise();
167
168	myt = UV(1);
169	myU = UV(2);
170	myV = UV(3);
171
172	gp_Vec Dtc, Dttc;
173	gp_Vec Dus, Dvs, Duvs, Duus, Dvvs;
174	myC->D2(myt, myP1, Dtc, Dttc);
175	myS->D2(myU,myV,myP2,Dus,Dvs,Duus,Dvvs,Duvs);
176
177	gp_Vec P1P2 (myP2,myP1);
178
179	F(1) = P1P2.Dot(Dtc);
180	F(2) = P1P2.Dot(Dus);
181	F(3) = P1P2.Dot(Dvs);
182
183	Df(1,1) = Dtc.SquareMagnitude() + P1P2.Dot(Dttc);
184	Df(1,2) = -Dus.Dot(Dtc);
185	Df(1,3) = -Dvs.Dot(Dtc);
186
187	Df(2,1) = -Df(1, 2); // Dtc.Dot(Dus);
188	Df(2,2) = -Dus.SquareMagnitude()+P1P2.Dot(Duus);
189	Df(2,3) = -Dvs.Dot(Dus)+P1P2.Dot(Duvs);
190
191	Df(3,1) = -Df(1,3); // Dtc.Dot(Dvs);
192	Df(3,2) = Df(2,3); // -Dus.Dot(Dvs)+P1P2.Dot(Duvs);
193	Df(3,3) = -Dvs.SquareMagnitude()+P1P2.Dot(Dvvs);
194
195	return Standard_True;
196
197	}
198
199	//=======================================================================
200	//function : GetStateNumber
201	//purpose :
202	//=======================================================================
203
204	Standard_Integer Extrema_FuncExtCS::GetStateNumber()
205	{
206	if (!myCinit \|\| !mySinit) Standard_TypeMismatch::Raise();
207	#if 0
208	math_Vector Sol(1, 3), UVSol(1, 3);
209	UVSol(1) = myt; UVSol(2) = myU; UVSol(3) = myV;
210	Value(UVSol, Sol);
211	cout <<"F(1)= "<<Sol(1)<<" F(2)= "<<Sol(2)<<" F(3)= "<<Sol(3)<<endl;
212	#endif
5368adff	213	//comparison of solution with previous solutions
	214	Standard_Real tol2d = Precision::PConfusion() * Precision::PConfusion();
	215	Standard_Integer i = 1, nbSol = mySqDist.Length();
	216	for( ; i <= nbSol; i++)
	217	{
	218	Standard_Real aU = myPoint1(i).Parameter();
	219	if( (myU - aU) * (myU - aU) <= tol2d )
	220	break;
	221	}
	222	if (i <= nbSol)
	223	return 0;
7fd59977	224	mySqDist.Append(myP1.SquareDistance(myP2));
	225	myPoint1.Append(Extrema_POnCurv(myt,myP1));
	226	myPoint2.Append(Extrema_POnSurf(myU,myV,myP2));
	227	return 0;
	228	}
	229
	230	//=======================================================================
	231	//function : NbExt
	232	//purpose :
	233	//=======================================================================
	234
	235	Standard_Integer Extrema_FuncExtCS::NbExt() const
	236	{
	237	return mySqDist.Length();
	238	}
	239
	240	//=======================================================================
	241	//function : SquareDistance
	242	//purpose :
	243	//=======================================================================
	244
	245	Standard_Real Extrema_FuncExtCS::SquareDistance(const Standard_Integer N) const
	246	{
	247	if (!myCinit \|\| !mySinit) Standard_TypeMismatch::Raise();
	248	return mySqDist.Value(N);
	249	}
	250
	251	//=======================================================================
	252	//function : PointOnCurve
	253	//purpose :
	254	//=======================================================================
	255
	256	const Extrema_POnCurv& Extrema_FuncExtCS::PointOnCurve(const Standard_Integer N) const
	257	{
	258	if (!myCinit \|\| !mySinit) Standard_TypeMismatch::Raise();
	259	return myPoint1.Value(N);
	260	}
	261
	262	//=======================================================================
	263	//function : PointOnSurface
	264	//purpose :
	265	//=======================================================================
	266
	267	const Extrema_POnSurf& Extrema_FuncExtCS::PointOnSurface(const Standard_Integer N) const
	268	{
	269	if (!myCinit \|\| !mySinit) Standard_TypeMismatch::Raise();
	270	return myPoint2.Value(N);
	271	}
	272
	273	//=======================================================================
	274	//function : Bidon1
	275	//purpose :
	276	//=======================================================================
	277
	278	Adaptor3d_SurfacePtr Extrema_FuncExtCS::Bidon1() const
	279	{
	280	return (Adaptor3d_SurfacePtr)0L;
	281	}
	282
	283	//=======================================================================
	284	//function : Bidon2
	285	//purpose :
	286	//=======================================================================
	287
288	Adaptor3d_CurvePtr Extrema_FuncExtCS::Bidon2() const
289	{
290	return (Adaptor3d_CurvePtr)0L;
291	}
292