| 1 | // Copyright (c) 1996-1999 Matra Datavision |
| 2 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
| 3 | // |
| 4 | // This file is part of Open CASCADE Technology software library. |
| 5 | // |
| 6 | // This library is free software; you can redistribute it and/or modify it under |
| 7 | // the terms of the GNU Lesser General Public License version 2.1 as published |
| 8 | // by the Free Software Foundation, with special exception defined in the file |
| 9 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
| 10 | // distribution for complete text of the license and disclaimer of any warranty. |
| 11 | // |
| 12 | // Alternatively, this file may be used under the terms of Open CASCADE |
| 13 | // commercial license or contractual agreement. |
| 14 | |
| 15 | // Jeannine PANCIATICI le 06/06/96 |
| 16 | // Igor FEOKTISTOV 14/12/98 - correction of Approximate() and Init(). |
| 17 | // Approximation d une MultiLine de points decrite par le tool MLineTool. |
| 18 | // avec criteres variationnels |
| 19 | |
| 20 | #include <AppDef_MultiLine.hxx> |
| 21 | #include <AppDef_SmoothCriterion.hxx> |
| 22 | #include <AppDef_Variational.hxx> |
| 23 | #include <AppParCurves_MultiBSpCurve.hxx> |
| 24 | #include <FEmTool_Assembly.hxx> |
| 25 | #include <FEmTool_Curve.hxx> |
| 26 | #include <gp_VectorWithNullMagnitude.hxx> |
| 27 | #include <math_Matrix.hxx> |
| 28 | #include <PLib_Base.hxx> |
| 29 | #include <Standard_ConstructionError.hxx> |
| 30 | #include <Standard_DimensionError.hxx> |
| 31 | #include <Standard_DomainError.hxx> |
| 32 | #include <Standard_OutOfRange.hxx> |
| 33 | #include <StdFail_NotDone.hxx> |
| 34 | |
| 35 | #define No_Standard_RangeError |
| 36 | #define No_Standard_OutOfRange |
| 37 | #define No_Standard_DimensionError |
| 38 | #define No_Standard_ConstructionError |
| 39 | |
| 40 | #include <Standard_Stream.hxx> |
| 41 | |
| 42 | #include <AppParCurves.hxx> |
| 43 | #include <AppParCurves_Constraint.hxx> |
| 44 | #include <AppParCurves_HArray1OfConstraintCouple.hxx> |
| 45 | #include <AppParCurves_Array1OfMultiPoint.hxx> |
| 46 | #include <AppParCurves_MultiPoint.hxx> |
| 47 | #include <AppParCurves_MultiBSpCurve.hxx> |
| 48 | #include <AppDef_LinearCriteria.hxx> |
| 49 | #include <Convert_CompPolynomialToPoles.hxx> |
| 50 | #include <gp_Pnt.hxx> |
| 51 | #include <gp_Pnt2d.hxx> |
| 52 | #include <gp_Vec.hxx> |
| 53 | #include <gp_Vec2d.hxx> |
| 54 | #include <TColgp_Array1OfPnt.hxx> |
| 55 | #include <TColgp_Array1OfPnt2d.hxx> |
| 56 | #include <TColgp_Array1OfVec.hxx> |
| 57 | #include <TColgp_Array1OfVec2d.hxx> |
| 58 | #include <TColStd_Array1OfInteger.hxx> |
| 59 | #include <TColStd_HArray1OfInteger.hxx> |
| 60 | #include <TColStd_Array1OfReal.hxx> |
| 61 | #include <TColStd_HArray1OfReal.hxx> |
| 62 | #include <TColStd_HArray2OfReal.hxx> |
| 63 | #include <StdFail_NotDone.hxx> |
| 64 | #include <Standard_SStream.hxx> |
| 65 | #include <Standard_NoSuchObject.hxx> |
| 66 | #include <Precision.hxx> |
| 67 | #include <AppDef_MyLineTool.hxx> |
| 68 | |
| 69 | #include <TColStd_HArray2OfInteger.hxx> |
| 70 | #include <TColStd_Array2OfInteger.hxx> |
| 71 | #include <TColStd_Array2OfReal.hxx> |
| 72 | #include <FEmTool_Assembly.hxx> |
| 73 | #include <FEmTool_AssemblyTable.hxx> |
| 74 | #include <FEmTool_Curve.hxx> |
| 75 | #include <math_Matrix.hxx> |
| 76 | #include <math_Vector.hxx> |
| 77 | #include <PLib_Base.hxx> |
| 78 | #include <PLib_JacobiPolynomial.hxx> |
| 79 | #include <PLib_HermitJacobi.hxx> |
| 80 | #include <FEmTool_HAssemblyTable.hxx> |
| 81 | |
| 82 | // Add this line: |
| 83 | #include <algorithm> |
| 84 | |
| 85 | #if defined(_MSC_VER) |
| 86 | # include <stdio.h> |
| 87 | # include <stdarg.h> |
| 88 | #endif /* _MSC_VER */ |
| 89 | |
| 90 | // |
| 91 | //======================================================================= |
| 92 | //function : AppDef_Variational |
| 93 | //purpose : Initialization of the fields. |
| 94 | //======================================================================= |
| 95 | // |
| 96 | AppDef_Variational::AppDef_Variational(const AppDef_MultiLine& SSP, |
| 97 | const Standard_Integer FirstPoint, |
| 98 | const Standard_Integer LastPoint, |
| 99 | const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints, |
| 100 | const Standard_Integer MaxDegree, |
| 101 | const Standard_Integer MaxSegment, |
| 102 | const GeomAbs_Shape Continuity, |
| 103 | const Standard_Boolean WithMinMax, |
| 104 | const Standard_Boolean WithCutting, |
| 105 | const Standard_Real Tolerance, |
| 106 | const Standard_Integer NbIterations): |
| 107 | mySSP(SSP), |
| 108 | myFirstPoint(FirstPoint), |
| 109 | myLastPoint(LastPoint), |
| 110 | myConstraints(TheConstraints), |
| 111 | myMaxDegree(MaxDegree), |
| 112 | myMaxSegment(MaxSegment), |
| 113 | myNbIterations(NbIterations), |
| 114 | myTolerance(Tolerance), |
| 115 | myContinuity(Continuity), |
| 116 | myWithMinMax(WithMinMax), |
| 117 | myWithCutting(WithCutting) |
| 118 | { |
| 119 | // Verifications: |
| 120 | if (myMaxDegree < 1) Standard_DomainError::Raise(); |
| 121 | myMaxDegree = Min (30, myMaxDegree); |
| 122 | // |
| 123 | if (myMaxSegment < 1) Standard_DomainError::Raise(); |
| 124 | // |
| 125 | if (myWithMinMax != 0 && myWithMinMax !=1 ) Standard_DomainError::Raise(); |
| 126 | if (myWithCutting != 0 && myWithCutting !=1 ) Standard_DomainError::Raise(); |
| 127 | // |
| 128 | myIsOverConstr = Standard_False; |
| 129 | myIsCreated = Standard_False; |
| 130 | myIsDone = Standard_False; |
| 131 | switch (myContinuity) { |
| 132 | case GeomAbs_C0: |
| 133 | myNivCont=0; |
| 134 | break ; |
| 135 | case GeomAbs_C1: |
| 136 | myNivCont=1; |
| 137 | break ; |
| 138 | case GeomAbs_C2: |
| 139 | myNivCont=2; |
| 140 | break ; |
| 141 | default: |
| 142 | Standard_ConstructionError::Raise(); |
| 143 | } |
| 144 | // |
| 145 | myNbP2d = AppDef_MyLineTool::NbP2d(SSP); |
| 146 | myNbP3d = AppDef_MyLineTool::NbP3d(SSP); |
| 147 | myDimension = 2 * myNbP2d + 3* myNbP3d ; |
| 148 | // |
| 149 | myPercent[0]=0.4; |
| 150 | myPercent[1]=0.2; |
| 151 | myPercent[2]=0.4; |
| 152 | myKnots= new TColStd_HArray1OfReal(1,2); |
| 153 | myKnots->SetValue(1,0.); |
| 154 | myKnots->SetValue(2,1.); |
| 155 | |
| 156 | // Declaration |
| 157 | // |
| 158 | mySmoothCriterion = new AppDef_LinearCriteria(mySSP, myFirstPoint, myLastPoint); |
| 159 | myParameters = new TColStd_HArray1OfReal(myFirstPoint, myLastPoint); |
| 160 | myNbPoints=myLastPoint-myFirstPoint+1; |
| 161 | if (myNbPoints <= 0) Standard_ConstructionError::Raise(); |
| 162 | // |
| 163 | myTabPoints= new TColStd_HArray1OfReal(1,myDimension*myNbPoints); |
| 164 | // |
| 165 | // Table of Points initialization |
| 166 | // |
| 167 | Standard_Integer ipoint,jp2d,jp3d,index; |
| 168 | TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d)); |
| 169 | TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d)); |
| 170 | gp_Pnt2d P2d; |
| 171 | gp_Pnt P3d; |
| 172 | index=1; |
| 173 | |
| 174 | for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++) |
| 175 | { |
| 176 | |
| 177 | if(myNbP2d !=0 && myNbP3d ==0 ) { |
| 178 | AppDef_MyLineTool::Value(mySSP,ipoint,TabP2d); |
| 179 | |
| 180 | for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++) |
| 181 | { P2d = TabP2d.Value(jp2d); |
| 182 | |
| 183 | myTabPoints->SetValue(index++,P2d.X()); |
| 184 | myTabPoints->SetValue(index++,P2d.Y()); |
| 185 | } |
| 186 | } |
| 187 | if(myNbP3d !=0 && myNbP2d == 0 ) { |
| 188 | AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d); |
| 189 | |
| 190 | for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++) |
| 191 | |
| 192 | { P3d=TabP3d.Value(jp3d); |
| 193 | |
| 194 | myTabPoints->SetValue(index++,P3d.X()); |
| 195 | myTabPoints->SetValue(index++,P3d.Y()); |
| 196 | myTabPoints->SetValue(index++,P3d.Z()); |
| 197 | |
| 198 | } |
| 199 | } |
| 200 | if(myNbP3d !=0 && myNbP2d != 0 ) { |
| 201 | AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d,TabP2d); |
| 202 | |
| 203 | for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++) |
| 204 | |
| 205 | { P3d=TabP3d.Value(jp3d); |
| 206 | |
| 207 | myTabPoints->SetValue(index++,P3d.X()); |
| 208 | myTabPoints->SetValue(index++,P3d.Y()); |
| 209 | myTabPoints->SetValue(index++,P3d.Z()); |
| 210 | |
| 211 | } |
| 212 | for ( jp2d = 1 ; jp2d <= myNbP2d ; jp2d++) |
| 213 | |
| 214 | { P2d=TabP2d.Value(jp2d); |
| 215 | |
| 216 | myTabPoints->SetValue(index++,P2d.X()); |
| 217 | myTabPoints->SetValue(index++,P2d.Y()); |
| 218 | } |
| 219 | } |
| 220 | } |
| 221 | Init(); |
| 222 | } |
| 223 | // |
| 224 | //======================================================================= |
| 225 | //function : Init |
| 226 | //purpose : Initializes the tables of constraints |
| 227 | // and verifies if the problem is not over-constrained. |
| 228 | // This method is used in the Create and the method SetConstraint. |
| 229 | //======================================================================= |
| 230 | // |
| 231 | void AppDef_Variational::Init() |
| 232 | { |
| 233 | |
| 234 | Standard_Integer ipoint,jp2d,jp3d,index,jndex; |
| 235 | Standard_Integer CurMultyPoint; |
| 236 | TColgp_Array1OfVec TabV3d(1, Max(1,myNbP3d)); |
| 237 | TColgp_Array1OfVec2d TabV2d(1, Max(1,myNbP2d)); |
| 238 | TColgp_Array1OfVec TabV3dcurv(1, Max(1,myNbP3d)); |
| 239 | TColgp_Array1OfVec2d TabV2dcurv(1, Max(1,myNbP2d)); |
| 240 | |
| 241 | gp_Vec Vt3d, Vc3d; |
| 242 | gp_Vec2d Vt2d, Vc2d; |
| 243 | |
| 244 | myNbConstraints=myConstraints->Length(); |
| 245 | if (myNbConstraints < 0) Standard_ConstructionError::Raise(); |
| 246 | |
| 247 | myTypConstraints = new TColStd_HArray1OfInteger(1,Max(1,2*myNbConstraints)); |
| 248 | myTabConstraints = new TColStd_HArray1OfReal(1,Max(1,2*myDimension*myNbConstraints)); |
| 249 | myTtheta = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints)); |
| 250 | myTfthet = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints)); |
| 251 | |
| 252 | |
| 253 | // |
| 254 | // Table of types initialization |
| 255 | Standard_Integer iconstr; |
| 256 | index=1; |
| 257 | jndex=1; |
| 258 | CurMultyPoint = 1; |
| 259 | myNbPassPoints=0; |
| 260 | myNbTangPoints=0; |
| 261 | myNbCurvPoints=0; |
| 262 | AppParCurves_Constraint valcontr; |
| 263 | |
| 264 | for ( iconstr = myConstraints->Lower() ; iconstr <= myConstraints->Upper() ; iconstr++) |
| 265 | { |
| 266 | ipoint=(myConstraints->Value(iconstr)).Index(); |
| 267 | |
| 268 | valcontr=(myConstraints->Value(iconstr)).Constraint(); |
| 269 | switch (valcontr) { |
| 270 | case AppParCurves_NoConstraint: |
| 271 | CurMultyPoint -= myNbP3d * 6 + myNbP2d * 2; |
| 272 | break ; |
| 273 | case AppParCurves_PassPoint: |
| 274 | myTypConstraints->SetValue(index++,ipoint); |
| 275 | myTypConstraints->SetValue(index++,0); |
| 276 | myNbPassPoints++; |
| 277 | if(myNbP2d !=0 ) jndex=jndex+4*myNbP2d; |
| 278 | if(myNbP3d !=0 ) jndex=jndex+6*myNbP3d; |
| 279 | break ; |
| 280 | case AppParCurves_TangencyPoint: |
| 281 | myTypConstraints->SetValue(index++,ipoint); |
| 282 | myTypConstraints->SetValue(index++,1); |
| 283 | myNbTangPoints++; |
| 284 | if(myNbP2d !=0 && myNbP3d == 0 ) |
| 285 | { |
| 286 | if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV2d) == Standard_False) |
| 287 | Standard_ConstructionError::Raise(); |
| 288 | for (jp2d=1;jp2d<=myNbP2d;jp2d++) |
| 289 | { |
| 290 | Vt2d=TabV2d.Value(jp2d); |
| 291 | Vt2d.Normalize(); |
| 292 | myTabConstraints->SetValue(jndex++,Vt2d.X()); |
| 293 | myTabConstraints->SetValue(jndex++,Vt2d.Y()); |
| 294 | jndex=jndex+2; |
| 295 | InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4); |
| 296 | } |
| 297 | } |
| 298 | if(myNbP3d !=0 && myNbP2d == 0) |
| 299 | { |
| 300 | if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d) == Standard_False) |
| 301 | Standard_ConstructionError::Raise(); |
| 302 | for (jp3d=1;jp3d<=myNbP3d;jp3d++) |
| 303 | { |
| 304 | Vt3d=TabV3d.Value(jp3d); |
| 305 | Vt3d.Normalize(); |
| 306 | myTabConstraints->SetValue(jndex++,Vt3d.X()); |
| 307 | |
| 308 | myTabConstraints->SetValue(jndex++,Vt3d.Y()); |
| 309 | |
| 310 | myTabConstraints->SetValue(jndex++,Vt3d.Z()); |
| 311 | jndex=jndex+3; |
| 312 | InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6); |
| 313 | } |
| 314 | } |
| 315 | if(myNbP3d !=0 && myNbP2d != 0) |
| 316 | { |
| 317 | if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False) |
| 318 | Standard_ConstructionError::Raise(); |
| 319 | for (jp3d=1;jp3d<=myNbP3d;jp3d++) |
| 320 | { |
| 321 | Vt3d=TabV3d.Value(jp3d); |
| 322 | Vt3d.Normalize(); |
| 323 | myTabConstraints->SetValue(jndex++,Vt3d.X()); |
| 324 | myTabConstraints->SetValue(jndex++,Vt3d.Y()); |
| 325 | myTabConstraints->SetValue(jndex++,Vt3d.Z()); |
| 326 | jndex=jndex+3; |
| 327 | InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6); |
| 328 | } |
| 329 | |
| 330 | for (jp2d=1;jp2d<=myNbP2d;jp2d++) |
| 331 | { |
| 332 | Vt2d=TabV2d.Value(jp2d); |
| 333 | Vt2d.Normalize(); |
| 334 | myTabConstraints->SetValue(jndex++,Vt2d.X()); |
| 335 | myTabConstraints->SetValue(jndex++,Vt2d.Y()); |
| 336 | jndex=jndex+2; |
| 337 | InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4); |
| 338 | } |
| 339 | } |
| 340 | |
| 341 | |
| 342 | break ; |
| 343 | |
| 344 | case AppParCurves_CurvaturePoint: |
| 345 | myTypConstraints->SetValue(index++,ipoint); |
| 346 | myTypConstraints->SetValue(index++,2); |
| 347 | myNbCurvPoints++; |
| 348 | if(myNbP2d !=0 && myNbP3d == 0) |
| 349 | { |
| 350 | if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV2d) == Standard_False ) |
| 351 | Standard_ConstructionError::Raise(); |
| 352 | if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV2dcurv) == Standard_False) |
| 353 | Standard_ConstructionError::Raise(); |
| 354 | for (jp2d=1;jp2d<=myNbP2d;jp2d++) |
| 355 | { |
| 356 | Vt2d=TabV2d.Value(jp2d); |
| 357 | Vt2d.Normalize(); |
| 358 | Vc2d=TabV2dcurv.Value(jp2d); |
| 359 | if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular()) |
| 360 | Standard_ConstructionError::Raise(); |
| 361 | myTabConstraints->SetValue(jndex++,Vt2d.X()); |
| 362 | myTabConstraints->SetValue(jndex++,Vt2d.Y()); |
| 363 | myTabConstraints->SetValue(jndex++,Vc2d.X()); |
| 364 | myTabConstraints->SetValue(jndex++,Vc2d.Y()); |
| 365 | InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4); |
| 366 | } |
| 367 | } |
| 368 | |
| 369 | if(myNbP3d !=0 && myNbP2d == 0 ) |
| 370 | { |
| 371 | if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d) == Standard_False ) |
| 372 | Standard_ConstructionError::Raise(); |
| 373 | if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV3dcurv) == Standard_False) |
| 374 | Standard_ConstructionError::Raise(); |
| 375 | for (jp3d=1;jp3d<=myNbP3d;jp3d++) |
| 376 | { |
| 377 | Vt3d=TabV3d.Value(jp3d); |
| 378 | Vt3d.Normalize(); |
| 379 | Vc3d=TabV3dcurv.Value(jp3d); |
| 380 | if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False) |
| 381 | Standard_ConstructionError::Raise(); |
| 382 | myTabConstraints->SetValue(jndex++,Vt3d.X()); |
| 383 | myTabConstraints->SetValue(jndex++,Vt3d.Y()); |
| 384 | myTabConstraints->SetValue(jndex++,Vt3d.Z()); |
| 385 | myTabConstraints->SetValue(jndex++,Vc3d.X()); |
| 386 | myTabConstraints->SetValue(jndex++,Vc3d.Y()); |
| 387 | myTabConstraints->SetValue(jndex++,Vc3d.Z()); |
| 388 | InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6); |
| 389 | } |
| 390 | } |
| 391 | if(myNbP3d !=0 && myNbP2d != 0 ) |
| 392 | { |
| 393 | if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False ) |
| 394 | Standard_ConstructionError::Raise(); |
| 395 | if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV3dcurv,TabV2dcurv) == Standard_False) |
| 396 | Standard_ConstructionError::Raise(); |
| 397 | for (jp3d=1;jp3d<=myNbP3d;jp3d++) |
| 398 | { |
| 399 | Vt3d=TabV3d.Value(jp3d); |
| 400 | Vt3d.Normalize(); |
| 401 | Vc3d=TabV3dcurv.Value(jp3d); |
| 402 | if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False) |
| 403 | Standard_ConstructionError::Raise(); |
| 404 | myTabConstraints->SetValue(jndex++,Vt3d.X()); |
| 405 | myTabConstraints->SetValue(jndex++,Vt3d.Y()); |
| 406 | myTabConstraints->SetValue(jndex++,Vt3d.Z()); |
| 407 | myTabConstraints->SetValue(jndex++,Vc3d.X()); |
| 408 | myTabConstraints->SetValue(jndex++,Vc3d.Y()); |
| 409 | myTabConstraints->SetValue(jndex++,Vc3d.Z()); |
| 410 | InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6); |
| 411 | } |
| 412 | for (jp2d=1;jp2d<=myNbP2d;jp2d++) |
| 413 | { |
| 414 | Vt2d=TabV2d.Value(jp2d); |
| 415 | Vt2d.Normalize(); |
| 416 | Vc2d=TabV2dcurv.Value(jp2d); |
| 417 | if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular()) |
| 418 | Standard_ConstructionError::Raise(); |
| 419 | myTabConstraints->SetValue(jndex++,Vt2d.X()); |
| 420 | myTabConstraints->SetValue(jndex++,Vt2d.Y()); |
| 421 | myTabConstraints->SetValue(jndex++,Vc2d.X()); |
| 422 | myTabConstraints->SetValue(jndex++,Vc2d.Y()); |
| 423 | InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4); |
| 424 | } |
| 425 | |
| 426 | } |
| 427 | break ; |
| 428 | default: |
| 429 | Standard_ConstructionError::Raise(); |
| 430 | } |
| 431 | CurMultyPoint += myNbP3d * 6 + myNbP2d * 2; |
| 432 | } |
| 433 | // OverConstraint Detection |
| 434 | Standard_Integer MaxSeg; |
| 435 | if(myWithCutting == Standard_True) MaxSeg = myMaxSegment ; |
| 436 | else MaxSeg = 1; |
| 437 | if (((myMaxDegree-myNivCont)*MaxSeg-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 ) |
| 438 | { |
| 439 | myIsOverConstr =Standard_True; |
| 440 | myIsCreated = Standard_False; |
| 441 | } |
| 442 | else { |
| 443 | InitSmoothCriterion(); |
| 444 | myIsCreated = Standard_True; |
| 445 | } |
| 446 | |
| 447 | |
| 448 | } |
| 449 | // |
| 450 | //======================================================================= |
| 451 | //function : Approximate |
| 452 | //purpose : Makes the approximation with the current fields. |
| 453 | //======================================================================= |
| 454 | // |
| 455 | void AppDef_Variational::Approximate() |
| 456 | |
| 457 | { |
| 458 | if (myIsCreated == Standard_False ) StdFail_NotDone:: Raise(); |
| 459 | |
| 460 | |
| 461 | Standard_Real WQuadratic, WQuality; |
| 462 | |
| 463 | TColStd_Array1OfReal Ecarts(myFirstPoint, myLastPoint); |
| 464 | |
| 465 | mySmoothCriterion->GetWeight(WQuadratic, WQuality); |
| 466 | |
| 467 | Handle(FEmTool_Curve) TheCurve; |
| 468 | |
| 469 | mySmoothCriterion->GetCurve(TheCurve); |
| 470 | |
| 471 | //--------------------------------------------------------------------- |
| 472 | |
| 473 | TheMotor(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts); |
| 474 | |
| 475 | |
| 476 | if(myWithMinMax && myTolerance < myMaxError) |
| 477 | Adjusting(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts); |
| 478 | |
| 479 | //--------------------------------------------------------------------- |
| 480 | |
| 481 | Standard_Integer jp2d,jp3d,index,ipole, |
| 482 | NbElem = TheCurve->NbElements(); |
| 483 | |
| 484 | TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d)); |
| 485 | TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d)); |
| 486 | Standard_Real debfin[2] = {-1., 1}; |
| 487 | |
| 488 | gp_Pnt2d P2d; |
| 489 | gp_Pnt P3d; |
| 490 | |
| 491 | index=0; |
| 492 | |
| 493 | { |
| 494 | Handle(TColStd_HArray2OfReal) PolynomialIntervalsPtr = |
| 495 | new TColStd_HArray2OfReal(1,NbElem,1,2) ; |
| 496 | |
| 497 | Handle(TColStd_HArray1OfInteger) NbCoeffPtr = |
| 498 | new TColStd_HArray1OfInteger(1, myMaxSegment); |
| 499 | |
| 500 | Standard_Integer size = myMaxSegment * (myMaxDegree + 1) * myDimension ; |
| 501 | Handle(TColStd_HArray1OfReal) CoeffPtr = new TColStd_HArray1OfReal(1,size); |
| 502 | |
| 503 | CoeffPtr->Init(0.); |
| 504 | |
| 505 | Handle(TColStd_HArray1OfReal) IntervallesPtr = |
| 506 | new TColStd_HArray1OfReal(1, NbElem + 1); |
| 507 | |
| 508 | IntervallesPtr->ChangeArray1() = TheCurve->Knots(); |
| 509 | |
| 510 | TheCurve->GetPolynom(CoeffPtr->ChangeArray1()); |
| 511 | |
| 512 | Standard_Integer ii; |
| 513 | |
| 514 | for(ii = 1; ii <= NbElem; ii++) |
| 515 | NbCoeffPtr->SetValue(ii, TheCurve->Degree(ii)+1); |
| 516 | |
| 517 | |
| 518 | for (ii = PolynomialIntervalsPtr->LowerRow() ; |
| 519 | ii <= PolynomialIntervalsPtr->UpperRow() ;ii++) |
| 520 | { |
| 521 | PolynomialIntervalsPtr->SetValue(ii,1,debfin[0]) ; |
| 522 | PolynomialIntervalsPtr->SetValue(ii,2,debfin[1]) ; |
| 523 | } |
| 524 | /* |
| 525 | printf("\n =========== Parameters for converting\n"); |
| 526 | printf("nb_courbes : %d \n", NbElem); |
| 527 | printf("tab_option[4] : %d \n", myNivCont); |
| 528 | printf("myDimension : %d \n", myDimension); |
| 529 | printf("myMaxDegree : %d \n", myMaxDegree); |
| 530 | printf("\n NbcoeffPtr :\n"); |
| 531 | for(ii = 1; ii <= NbElem; ii++) printf("NbCoeffPtr(%d) = %d \n", ii, NbCoeffPtr->Value(ii)); |
| 532 | size = NbElem*(myMaxDegree + 1) * myDimension; |
| 533 | printf("\n CoeffPtr :\n"); |
| 534 | for(ii = 1; ii <= size; ii++) printf("CoeffPtr(%d) = %.8e \n", ii, CoeffPtr->Value(ii)); |
| 535 | printf("\n PolinimialIntervalsPtr :\n"); |
| 536 | for (ii = PolynomialIntervalsPtr->LowerRow() ; |
| 537 | ii <= PolynomialIntervalsPtr->UpperRow() ;ii++) |
| 538 | { |
| 539 | printf(" %d : [%f, %f] \n", ii, PolynomialIntervalsPtr->Value(ii,1), |
| 540 | PolynomialIntervalsPtr->Value(ii,2)) ; |
| 541 | } |
| 542 | printf("\n IntervallesPtr :\n"); |
| 543 | for (ii = IntervallesPtr->Lower(); |
| 544 | ii <= IntervallesPtr->Upper() - 1; ii++) |
| 545 | { |
| 546 | printf(" %d : [%f, %f] \n", ii, IntervallesPtr->Value(ii), |
| 547 | IntervallesPtr->Value(ii+1)) ; |
| 548 | } |
| 549 | */ |
| 550 | Convert_CompPolynomialToPoles AConverter(NbElem, |
| 551 | myNivCont, |
| 552 | myDimension, |
| 553 | myMaxDegree, |
| 554 | NbCoeffPtr, |
| 555 | CoeffPtr, |
| 556 | PolynomialIntervalsPtr, |
| 557 | IntervallesPtr) ; |
| 558 | if (AConverter.IsDone()) |
| 559 | { |
| 560 | Handle(TColStd_HArray2OfReal) PolesPtr ; |
| 561 | Handle(TColStd_HArray1OfInteger) Mults; |
| 562 | Standard_Integer NbPoles=AConverter.NbPoles(); |
| 563 | // Standard_Integer Deg=AConverter.Degree(); |
| 564 | AppParCurves_Array1OfMultiPoint TabMU(1,NbPoles); |
| 565 | AConverter.Poles(PolesPtr) ; |
| 566 | AConverter.Knots(myKnots) ; |
| 567 | AConverter.Multiplicities(Mults) ; |
| 568 | |
| 569 | for (ipole=PolesPtr->LowerRow();ipole<=PolesPtr->UpperRow();ipole++) |
| 570 | { |
| 571 | index=PolesPtr->LowerCol(); |
| 572 | /* if(myNbP2d !=0 ) |
| 573 | { |
| 574 | for (jp2d=1;jp2d<=myNbP2d;jp2d++) |
| 575 | { |
| 576 | P2d.SetX(PolesPtr->Value(ipole,index++)); |
| 577 | P2d.SetY(PolesPtr->Value(ipole,index++)); |
| 578 | TabP2d.SetValue(jp2d,P2d); |
| 579 | } |
| 580 | }*/ |
| 581 | if(myNbP3d !=0 ) |
| 582 | { |
| 583 | for (jp3d=1;jp3d<=myNbP3d;jp3d++) |
| 584 | { |
| 585 | // cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index); |
| 586 | P3d.SetX(PolesPtr->Value(ipole,index++)); |
| 587 | // cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index); |
| 588 | P3d.SetY(PolesPtr->Value(ipole,index++)); |
| 589 | // cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index); |
| 590 | P3d.SetZ(PolesPtr->Value(ipole,index++)); |
| 591 | TabP3d.SetValue(jp3d,P3d); |
| 592 | } |
| 593 | } |
| 594 | if(myNbP2d !=0 ) |
| 595 | { |
| 596 | for (jp2d=1;jp2d<=myNbP2d;jp2d++) |
| 597 | { |
| 598 | P2d.SetX(PolesPtr->Value(ipole,index++)); |
| 599 | P2d.SetY(PolesPtr->Value(ipole,index++)); |
| 600 | TabP2d.SetValue(jp2d,P2d); |
| 601 | } |
| 602 | } |
| 603 | if(myNbP2d !=0 && myNbP3d !=0) |
| 604 | { |
| 605 | AppParCurves_MultiPoint aMultiPoint(TabP3d,TabP2d); |
| 606 | TabMU.SetValue(ipole,aMultiPoint); |
| 607 | } |
| 608 | else if (myNbP2d !=0) |
| 609 | { |
| 610 | AppParCurves_MultiPoint aMultiPoint(TabP2d); |
| 611 | TabMU.SetValue(ipole,aMultiPoint); |
| 612 | } |
| 613 | else |
| 614 | { |
| 615 | |
| 616 | AppParCurves_MultiPoint aMultiPoint(TabP3d); |
| 617 | TabMU.SetValue(ipole,aMultiPoint); |
| 618 | } |
| 619 | |
| 620 | } |
| 621 | AppParCurves_MultiBSpCurve aCurve(TabMU,myKnots->Array1(),Mults->Array1()); |
| 622 | myMBSpCurve=aCurve; |
| 623 | myIsDone = Standard_True; |
| 624 | } |
| 625 | } |
| 626 | |
| 627 | } |
| 628 | // |
| 629 | //======================================================================= |
| 630 | //function : IsCreated |
| 631 | //purpose : returns True if the creation is done |
| 632 | //======================================================================= |
| 633 | // |
| 634 | Standard_Boolean AppDef_Variational::IsCreated() const |
| 635 | { |
| 636 | return myIsCreated; |
| 637 | } |
| 638 | // |
| 639 | //======================================================================= |
| 640 | //function : IsDone |
| 641 | //purpose : returns True if the approximation is ok |
| 642 | //======================================================================= |
| 643 | // |
| 644 | Standard_Boolean AppDef_Variational::IsDone() const |
| 645 | { |
| 646 | return myIsDone; |
| 647 | } |
| 648 | // |
| 649 | //======================================================================= |
| 650 | //function : IsOverConstrained |
| 651 | //purpose : returns True if the problem is overconstrained |
| 652 | // in this case, approximation cannot be done. |
| 653 | //======================================================================= |
| 654 | // |
| 655 | Standard_Boolean AppDef_Variational::IsOverConstrained() const |
| 656 | { |
| 657 | return myIsOverConstr; |
| 658 | } |
| 659 | // |
| 660 | //======================================================================= |
| 661 | //function : Value |
| 662 | //purpose : returns all the BSpline curves approximating the |
| 663 | // MultiLine SSP after minimization of the parameter. |
| 664 | |
| 665 | //======================================================================= |
| 666 | // |
| 667 | AppParCurves_MultiBSpCurve AppDef_Variational::Value() const |
| 668 | { |
| 669 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
| 670 | return myMBSpCurve; |
| 671 | |
| 672 | } |
| 673 | // |
| 674 | //======================================================================= |
| 675 | //function : MaxError |
| 676 | //purpose : returns the maximum of the distances between |
| 677 | // the points of the multiline and the approximation |
| 678 | // curves. |
| 679 | //======================================================================= |
| 680 | // |
| 681 | Standard_Real AppDef_Variational::MaxError() const |
| 682 | { |
| 683 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
| 684 | return myMaxError; |
| 685 | } |
| 686 | // |
| 687 | //======================================================================= |
| 688 | //function : MaxErrorIndex |
| 689 | //purpose : returns the index of the MultiPoint of ErrorMax |
| 690 | //======================================================================= |
| 691 | // |
| 692 | Standard_Integer AppDef_Variational::MaxErrorIndex() const |
| 693 | { |
| 694 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
| 695 | return myMaxErrorIndex; |
| 696 | } |
| 697 | // |
| 698 | //======================================================================= |
| 699 | //function : QuadraticError |
| 700 | //purpose : returns the quadratic average of the distances between |
| 701 | // the points of the multiline and the approximation |
| 702 | // curves. |
| 703 | //======================================================================= |
| 704 | // |
| 705 | Standard_Real AppDef_Variational::QuadraticError() const |
| 706 | { |
| 707 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
| 708 | return myCriterium[0]; |
| 709 | } |
| 710 | // |
| 711 | //======================================================================= |
| 712 | //function : Distance |
| 713 | //purpose : returns the distances between the points of the |
| 714 | // multiline and the approximation curves. |
| 715 | //======================================================================= |
| 716 | // |
| 717 | void AppDef_Variational::Distance(math_Matrix& mat) |
| 718 | |
| 719 | { |
| 720 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
| 721 | Standard_Integer ipoint,jp2d,jp3d,index; |
| 722 | TColgp_Array1OfPnt TabP3d(1,Max(1,myNbP3d)); |
| 723 | TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d)); |
| 724 | Standard_Integer j0 = mat.LowerCol() - myFirstPoint; |
| 725 | |
| 726 | gp_Pnt2d P2d; |
| 727 | gp_Pnt P3d; |
| 728 | |
| 729 | |
| 730 | gp_Pnt Pt3d; |
| 731 | gp_Pnt2d Pt2d; |
| 732 | |
| 733 | for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++) |
| 734 | { |
| 735 | index=1; |
| 736 | if(myNbP3d !=0 ) { |
| 737 | AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d); |
| 738 | |
| 739 | for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++) |
| 740 | |
| 741 | { P3d=TabP3d.Value(jp3d); |
| 742 | myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt3d); |
| 743 | mat(index++, j0 + ipoint)=P3d.Distance(Pt3d); |
| 744 | |
| 745 | } |
| 746 | } |
| 747 | if(myNbP2d !=0 ) { |
| 748 | if(myNbP3d == 0 ) AppDef_MyLineTool::Value(mySSP,ipoint,TabP2d); |
| 749 | else AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d,TabP2d); |
| 750 | for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++) |
| 751 | |
| 752 | { P2d = TabP2d.Value(jp2d); |
| 753 | myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt2d); |
| 754 | mat(index++, j0 + ipoint)=P2d.Distance(Pt2d); |
| 755 | } |
| 756 | } |
| 757 | } |
| 758 | |
| 759 | } |
| 760 | // |
| 761 | //======================================================================= |
| 762 | //function : AverageError |
| 763 | //purpose : returns the average error between |
| 764 | // the MultiLine and the approximation. |
| 765 | //======================================================================= |
| 766 | // |
| 767 | Standard_Real AppDef_Variational::AverageError() const |
| 768 | { |
| 769 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
| 770 | return myAverageError; |
| 771 | } |
| 772 | // |
| 773 | //======================================================================= |
| 774 | //function : Parameters |
| 775 | //purpose : returns the parameters uses to the approximations |
| 776 | //======================================================================= |
| 777 | // |
| 778 | const Handle(TColStd_HArray1OfReal)& AppDef_Variational::Parameters() const |
| 779 | { |
| 780 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
| 781 | return myParameters; |
| 782 | } |
| 783 | // |
| 784 | //======================================================================= |
| 785 | //function : Knots |
| 786 | //purpose : returns the knots uses to the approximations |
| 787 | //======================================================================= |
| 788 | // |
| 789 | const Handle(TColStd_HArray1OfReal)& AppDef_Variational::Knots() const |
| 790 | { |
| 791 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
| 792 | return myKnots; |
| 793 | } |
| 794 | // |
| 795 | //======================================================================= |
| 796 | //function : Criterium |
| 797 | //purpose : returns the values of the quality criterium. |
| 798 | //======================================================================= |
| 799 | // |
| 800 | void AppDef_Variational::Criterium(Standard_Real& VFirstOrder, Standard_Real& VSecondOrder, Standard_Real& VThirdOrder) const |
| 801 | { |
| 802 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
| 803 | VFirstOrder=myCriterium[1] ; |
| 804 | VSecondOrder=myCriterium[2]; |
| 805 | VThirdOrder=myCriterium[3]; |
| 806 | } |
| 807 | // |
| 808 | //======================================================================= |
| 809 | //function : CriteriumWeight |
| 810 | //purpose : returns the Weights (as percent) associed to the criterium used in |
| 811 | // the optimization. |
| 812 | //======================================================================= |
| 813 | // |
| 814 | void AppDef_Variational::CriteriumWeight(Standard_Real& Percent1, Standard_Real& Percent2, Standard_Real& Percent3) const |
| 815 | { |
| 816 | Percent1 = myPercent[0]; |
| 817 | Percent2 = myPercent[1]; |
| 818 | Percent3 = myPercent[2] ; |
| 819 | } |
| 820 | // |
| 821 | //======================================================================= |
| 822 | //function : MaxDegree |
| 823 | //purpose : returns the Maximum Degree used in the approximation |
| 824 | //======================================================================= |
| 825 | // |
| 826 | Standard_Integer AppDef_Variational::MaxDegree() const |
| 827 | { |
| 828 | return myMaxDegree; |
| 829 | } |
| 830 | // |
| 831 | //======================================================================= |
| 832 | //function : MaxSegment |
| 833 | //purpose : returns the Maximum of segment used in the approximation |
| 834 | //======================================================================= |
| 835 | // |
| 836 | Standard_Integer AppDef_Variational::MaxSegment() const |
| 837 | { |
| 838 | return myMaxSegment; |
| 839 | } |
| 840 | // |
| 841 | //======================================================================= |
| 842 | //function : Continuity |
| 843 | //purpose : returns the Continuity used in the approximation |
| 844 | //======================================================================= |
| 845 | // |
| 846 | GeomAbs_Shape AppDef_Variational::Continuity() const |
| 847 | { |
| 848 | return myContinuity; |
| 849 | } |
| 850 | // |
| 851 | //======================================================================= |
| 852 | //function : WithMinMax |
| 853 | //purpose : returns if the approximation search to minimize the |
| 854 | // maximum Error or not. |
| 855 | //======================================================================= |
| 856 | // |
| 857 | Standard_Boolean AppDef_Variational::WithMinMax() const |
| 858 | { |
| 859 | return myWithMinMax; |
| 860 | } |
| 861 | // |
| 862 | //======================================================================= |
| 863 | //function : WithCutting |
| 864 | //purpose : returns if the approximation can insert new Knots or not. |
| 865 | //======================================================================= |
| 866 | // |
| 867 | Standard_Boolean AppDef_Variational::WithCutting() const |
| 868 | { |
| 869 | return myWithCutting; |
| 870 | } |
| 871 | // |
| 872 | //======================================================================= |
| 873 | //function : Tolerance |
| 874 | //purpose : returns the tolerance used in the approximation. |
| 875 | //======================================================================= |
| 876 | // |
| 877 | Standard_Real AppDef_Variational::Tolerance() const |
| 878 | { |
| 879 | return myTolerance; |
| 880 | } |
| 881 | // |
| 882 | //======================================================================= |
| 883 | //function : NbIterations |
| 884 | //purpose : returns the number of iterations used in the approximation. |
| 885 | //======================================================================= |
| 886 | // |
| 887 | Standard_Integer AppDef_Variational::NbIterations() const |
| 888 | { |
| 889 | return myNbIterations; |
| 890 | } |
| 891 | // |
| 892 | //======================================================================= |
| 893 | //function : Dump |
| 894 | //purpose : Prints on the stream o information on the current state |
| 895 | // of the object. |
| 896 | //======================================================================= |
| 897 | // |
| 898 | void AppDef_Variational::Dump(Standard_OStream& o) const |
| 899 | { |
| 900 | o << " \nVariational Smoothing " << endl; |
| 901 | o << " Number of multipoints " << myNbPoints << endl; |
| 902 | o << " Number of 2d par multipoint " << myNbP2d << endl; |
| 903 | o << " Nombre of 3d par multipoint " << myNbP3d << endl; |
| 904 | o << " Number of PassagePoint " << myNbPassPoints << endl; |
| 905 | o << " Number of TangencyPoints " << myNbTangPoints << endl; |
| 906 | o << " Number of CurvaturePoints " << myNbCurvPoints << endl; |
| 907 | o << " \nTolerance " << o.setf(ios::scientific) << setprecision(3) << setw(9) << myTolerance; |
| 908 | if ( WithMinMax()) { o << " as Max Error." << endl;} |
| 909 | else { o << " as size Error." << endl;} |
| 910 | o << "CriteriumWeights : " << myPercent[0] << " , " |
| 911 | << myPercent[1] << " , " << myPercent[2] << endl; |
| 912 | |
| 913 | if (myIsDone ) { |
| 914 | o << " MaxError " << setprecision(3) << setw(9) << myMaxError << endl; |
| 915 | o << " Index of MaxError " << myMaxErrorIndex << endl; |
| 916 | o << " Average Error " << setprecision(3) << setw(9) << myAverageError << endl; |
| 917 | o << " Quadratic Error " << setprecision(3) << setw(9) << myCriterium[0] << endl; |
| 918 | o << " Tension Criterium " << setprecision(3) << setw(9) << myCriterium[1] << endl; |
| 919 | o << " Flexion Criterium " << setprecision(3) << setw(9) << myCriterium[2] << endl; |
| 920 | o << " Jerk Criterium " << setprecision(3) << setw(9) << myCriterium[3] << endl; |
| 921 | o << " NbSegments " << myKnots->Length()-1 << endl; |
| 922 | } |
| 923 | else |
| 924 | { if (myIsOverConstr) o << "The probleme is overconstraint " << endl; |
| 925 | else o << " Erreur dans l''approximation" << endl; |
| 926 | } |
| 927 | } |
| 928 | // |
| 929 | //======================================================================= |
| 930 | //function : SetConstraints |
| 931 | //purpose : Define the constraints to approximate |
| 932 | // If this value is incompatible with the others fields |
| 933 | // this method modify nothing and returns false |
| 934 | //======================================================================= |
| 935 | // |
| 936 | Standard_Boolean AppDef_Variational::SetConstraints(const Handle(AppParCurves_HArray1OfConstraintCouple)& aConstraint) |
| 937 | |
| 938 | { |
| 939 | myConstraints=aConstraint; |
| 940 | Init(); |
| 941 | if (myIsOverConstr ) return Standard_False; |
| 942 | else return Standard_True; |
| 943 | } |
| 944 | // |
| 945 | //======================================================================= |
| 946 | //function : SetParameters |
| 947 | //purpose : Defines the parameters used by the approximations. |
| 948 | //======================================================================= |
| 949 | // |
| 950 | void AppDef_Variational::SetParameters(const Handle(TColStd_HArray1OfReal)& param) |
| 951 | { |
| 952 | myParameters->ChangeArray1() = param->Array1(); |
| 953 | } |
| 954 | // |
| 955 | //======================================================================= |
| 956 | //function : SetKnots |
| 957 | //purpose : Defines the knots used by the approximations |
| 958 | // -- If this value is incompatible with the others fields |
| 959 | // -- this method modify nothing and returns false |
| 960 | //======================================================================= |
| 961 | // |
| 962 | Standard_Boolean AppDef_Variational::SetKnots(const Handle(TColStd_HArray1OfReal)& knots) |
| 963 | { |
| 964 | myKnots->ChangeArray1() = knots->Array1(); |
| 965 | return Standard_True; |
| 966 | } |
| 967 | // |
| 968 | //======================================================================= |
| 969 | //function : SetMaxDegree |
| 970 | //purpose : Define the Maximum Degree used in the approximation |
| 971 | // If this value is incompatible with the others fields |
| 972 | // this method modify nothing and returns false |
| 973 | //======================================================================= |
| 974 | // |
| 975 | Standard_Boolean AppDef_Variational::SetMaxDegree(const Standard_Integer Degree) |
| 976 | { |
| 977 | if (((Degree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 ) |
| 978 | return Standard_False; |
| 979 | else |
| 980 | { |
| 981 | myMaxDegree=Degree; |
| 982 | |
| 983 | InitSmoothCriterion(); |
| 984 | |
| 985 | return Standard_True; |
| 986 | } |
| 987 | |
| 988 | } |
| 989 | // |
| 990 | //======================================================================= |
| 991 | //function : SetMaxSegment |
| 992 | //purpose : Define the maximum number of segments used in the approximation |
| 993 | // If this value is incompatible with the others fields |
| 994 | // this method modify nothing and returns false |
| 995 | //======================================================================= |
| 996 | // |
| 997 | Standard_Boolean AppDef_Variational::SetMaxSegment(const Standard_Integer NbSegment) |
| 998 | { |
| 999 | if ( myWithCutting == Standard_True && |
| 1000 | ((myMaxDegree-myNivCont)*NbSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 ) |
| 1001 | return Standard_False; |
| 1002 | else |
| 1003 | { |
| 1004 | myMaxSegment=NbSegment; |
| 1005 | return Standard_True; |
| 1006 | } |
| 1007 | } |
| 1008 | // |
| 1009 | //======================================================================= |
| 1010 | //function : SetContinuity |
| 1011 | //purpose : Define the Continuity used in the approximation |
| 1012 | // If this value is incompatible with the others fields |
| 1013 | // this method modify nothing and returns false |
| 1014 | //======================================================================= |
| 1015 | // |
| 1016 | Standard_Boolean AppDef_Variational::SetContinuity(const GeomAbs_Shape C) |
| 1017 | { |
| 1018 | Standard_Integer NivCont=0; |
| 1019 | switch (C) { |
| 1020 | case GeomAbs_C0: |
| 1021 | NivCont=0; |
| 1022 | break ; |
| 1023 | case GeomAbs_C1: |
| 1024 | NivCont=1; |
| 1025 | break ; |
| 1026 | case GeomAbs_C2: |
| 1027 | NivCont=2; |
| 1028 | break ; |
| 1029 | default: |
| 1030 | Standard_ConstructionError::Raise(); |
| 1031 | } |
| 1032 | if (((myMaxDegree-NivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 ) |
| 1033 | return Standard_False; |
| 1034 | else |
| 1035 | { |
| 1036 | myContinuity= C; |
| 1037 | myNivCont=NivCont; |
| 1038 | |
| 1039 | InitSmoothCriterion(); |
| 1040 | return Standard_True; |
| 1041 | } |
| 1042 | } |
| 1043 | // |
| 1044 | //======================================================================= |
| 1045 | //function : SetWithMinMax |
| 1046 | //purpose : Define if the approximation search to minimize the |
| 1047 | // maximum Error or not. |
| 1048 | //======================================================================= |
| 1049 | // |
| 1050 | void AppDef_Variational::SetWithMinMax(const Standard_Boolean MinMax) |
| 1051 | { |
| 1052 | myWithMinMax=MinMax; |
| 1053 | |
| 1054 | InitSmoothCriterion(); |
| 1055 | } |
| 1056 | // |
| 1057 | //======================================================================= |
| 1058 | //function : SetWithCutting |
| 1059 | //purpose : Define if the approximation can insert new Knots or not. |
| 1060 | // If this value is incompatible with the others fields |
| 1061 | // this method modify nothing and returns false |
| 1062 | //======================================================================= |
| 1063 | // |
| 1064 | Standard_Boolean AppDef_Variational::SetWithCutting(const Standard_Boolean Cutting) |
| 1065 | { |
| 1066 | if (Cutting == Standard_False) |
| 1067 | { |
| 1068 | if (((myMaxDegree-myNivCont)*myKnots->Length()-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 ) |
| 1069 | return Standard_False; |
| 1070 | |
| 1071 | else |
| 1072 | { |
| 1073 | myWithCutting=Cutting; |
| 1074 | InitSmoothCriterion(); |
| 1075 | return Standard_True; |
| 1076 | } |
| 1077 | } |
| 1078 | else |
| 1079 | { |
| 1080 | if (((myMaxDegree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 ) |
| 1081 | return Standard_False; |
| 1082 | |
| 1083 | else |
| 1084 | { |
| 1085 | myWithCutting=Cutting; |
| 1086 | InitSmoothCriterion(); |
| 1087 | return Standard_True; |
| 1088 | } |
| 1089 | } |
| 1090 | } |
| 1091 | // |
| 1092 | //======================================================================= |
| 1093 | //function : SetCriteriumWeight |
| 1094 | //purpose : define the Weights (as percent) associed to the criterium used in |
| 1095 | // the optimization. |
| 1096 | //======================================================================= |
| 1097 | // |
| 1098 | void AppDef_Variational::SetCriteriumWeight(const Standard_Real Percent1, const Standard_Real Percent2, const Standard_Real Percent3) |
| 1099 | { |
| 1100 | if (Percent1 < 0 || Percent2 < 0 || Percent3 < 0 ) Standard_DomainError::Raise(); |
| 1101 | Standard_Real Total = Percent1 + Percent2 + Percent3; |
| 1102 | myPercent[0] = Percent1/Total; |
| 1103 | myPercent[1] = Percent2/Total; |
| 1104 | myPercent[2] = Percent3/Total; |
| 1105 | |
| 1106 | InitSmoothCriterion(); |
| 1107 | } |
| 1108 | // |
| 1109 | //======================================================================= |
| 1110 | //function : SetCriteriumWeight |
| 1111 | //purpose : define the Weight (as percent) associed to the |
| 1112 | // criterium Order used in the optimization : Others |
| 1113 | // weights are updated. |
| 1114 | //======================================================================= |
| 1115 | // |
| 1116 | void AppDef_Variational::SetCriteriumWeight(const Standard_Integer Order, const Standard_Real Percent) |
| 1117 | { |
| 1118 | if ( Percent < 0 ) Standard_DomainError::Raise(); |
| 1119 | if ( Order < 1 || Order > 3 ) Standard_ConstructionError::Raise(); |
| 1120 | myPercent[Order-1] = Percent; |
| 1121 | Standard_Real Total = myPercent[0] + myPercent[1] + myPercent[2]; |
| 1122 | myPercent[0] = myPercent[0] / Total; |
| 1123 | myPercent[1] = myPercent[1] / Total; |
| 1124 | myPercent[2] = myPercent[2] / Total; |
| 1125 | |
| 1126 | InitSmoothCriterion(); |
| 1127 | |
| 1128 | } |
| 1129 | // |
| 1130 | //======================================================================= |
| 1131 | //function : SetTolerance |
| 1132 | //purpose : define the tolerance used in the approximation. |
| 1133 | //======================================================================= |
| 1134 | // |
| 1135 | void AppDef_Variational::SetTolerance(const Standard_Real Tol) |
| 1136 | { |
| 1137 | myTolerance=Tol; |
| 1138 | InitSmoothCriterion(); |
| 1139 | } |
| 1140 | // |
| 1141 | //======================================================================= |
| 1142 | //function : SetNbIterations |
| 1143 | //purpose : define the number of iterations used in the approximation. |
| 1144 | //======================================================================= |
| 1145 | // |
| 1146 | void AppDef_Variational::SetNbIterations(const Standard_Integer Iter) |
| 1147 | { |
| 1148 | myNbIterations=Iter; |
| 1149 | } |
| 1150 | |
| 1151 | //====================== Private Methodes =============================// |
| 1152 | //======================================================================= |
| 1153 | //function : TheMotor |
| 1154 | //purpose : Smoothing's motor like STRIM routine "MOTLIS" |
| 1155 | //======================================================================= |
| 1156 | void AppDef_Variational::TheMotor( |
| 1157 | Handle(AppDef_SmoothCriterion)& J, |
| 1158 | // const Standard_Real WQuadratic, |
| 1159 | const Standard_Real , |
| 1160 | const Standard_Real WQuality, |
| 1161 | Handle(FEmTool_Curve)& TheCurve, |
| 1162 | TColStd_Array1OfReal& Ecarts) |
| 1163 | { |
| 1164 | // ... |
| 1165 | |
| 1166 | const Standard_Real BigValue = 1.e37, SmallValue = 1.e-6, SmallestValue = 1.e-9; |
| 1167 | |
| 1168 | Handle(TColStd_HArray1OfReal) CurrentTi, NewTi, OldTi; |
| 1169 | Handle(TColStd_HArray2OfInteger) Dependence; |
| 1170 | Standard_Boolean lestim, lconst, ToOptim, iscut; |
| 1171 | Standard_Boolean isnear = Standard_False, again = Standard_True; |
| 1172 | Standard_Integer NbEst, ICDANA, NumPnt, Iter; |
| 1173 | Standard_Integer MaxNbEst =5; |
| 1174 | Standard_Real VOCRI[3] = {BigValue, BigValue, BigValue}, EROLD = BigValue, |
| 1175 | VALCRI[3], ERRMAX = BigValue, ERRMOY, ERRQUA; |
| 1176 | Standard_Real CBLONG, LNOLD; |
| 1177 | Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1; |
| 1178 | Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
| 1179 | Handle(FEmTool_Curve) CCurrent, COld, CNew; |
| 1180 | Standard_Real EpsLength = SmallValue; |
| 1181 | Standard_Real EpsDeg; |
| 1182 | |
| 1183 | Standard_Real e1, e2, e3; |
| 1184 | Standard_Real J1min, J2min, J3min; |
| 1185 | Standard_Integer iprog; |
| 1186 | |
| 1187 | // (0) Init |
| 1188 | |
| 1189 | J->GetEstimation(e1, e2, e3); |
| 1190 | J1min = 1.e-8; J2min = J3min = (e1 + 1.e-8) * 1.e-6; |
| 1191 | |
| 1192 | if(e1 < J1min) e1 = J1min;// Like in |
| 1193 | if(e2 < J2min) e2 = J2min;// MOTLIS |
| 1194 | if(e3 < J3min) e3 = J3min; |
| 1195 | |
| 1196 | |
| 1197 | J->SetEstimation(e1, e2, e3); |
| 1198 | |
| 1199 | CCurrent = TheCurve; |
| 1200 | CurrentTi = new TColStd_HArray1OfReal(1, myParameters->Length()); |
| 1201 | CurrentTi->ChangeArray1() = myParameters->Array1(); |
| 1202 | OldTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length()); |
| 1203 | OldTi->ChangeArray1() = CurrentTi->Array1(); |
| 1204 | COld = CCurrent; |
| 1205 | LNOLD = CBLONG = J->EstLength(); |
| 1206 | Dependence = J->DependenceTable(); |
| 1207 | |
| 1208 | J->SetCurve(CCurrent); |
| 1209 | FEmTool_Assembly * TheAssembly = |
| 1210 | new FEmTool_Assembly (Dependence->Array2(), J->AssemblyTable()); |
| 1211 | |
| 1212 | //============ Optimization ============================ |
| 1213 | // Standard_Integer inagain = 0; |
| 1214 | while (again) { |
| 1215 | |
| 1216 | // (1) Loop Optimization / Estimation |
| 1217 | lestim = Standard_True; |
| 1218 | lconst = Standard_True; |
| 1219 | NbEst = 0; |
| 1220 | |
| 1221 | J->SetCurve(CCurrent); |
| 1222 | |
| 1223 | while(lestim) { |
| 1224 | |
| 1225 | // (1.1) Curve's Optimization. |
| 1226 | EpsLength = SmallValue * CBLONG / NbrPnt; |
| 1227 | CNew = new (FEmTool_Curve) (CCurrent->Dimension(), |
| 1228 | CCurrent->NbElements(), CCurrent->Base(), EpsLength); |
| 1229 | CNew->Knots() = CCurrent->Knots(); |
| 1230 | |
| 1231 | J->SetParameters(CurrentTi); |
| 1232 | EpsDeg = Min(WQuality * .1, CBLONG * .001); |
| 1233 | |
| 1234 | Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1()); |
| 1235 | |
| 1236 | lconst = Standard_False; |
| 1237 | |
| 1238 | // (1.2) calcul des criteres de qualites et amelioration |
| 1239 | // des estimation. |
| 1240 | ICDANA = J->QualityValues(J1min, J2min, J3min, |
| 1241 | VALCRI[0], VALCRI[1], VALCRI[2]); |
| 1242 | |
| 1243 | if(ICDANA > 0) lconst = Standard_True; |
| 1244 | |
| 1245 | J->ErrorValues(ERRMAX, ERRQUA, ERRMOY); |
| 1246 | |
| 1247 | isnear = ((Sqrt(ERRQUA / NbrPnt) < 2*WQuality) && |
| 1248 | (myNbIterations > 1)); |
| 1249 | |
| 1250 | // (1.3) Optimisation des ti par proj orthogonale |
| 1251 | // et calcul de l'erreur aux points. |
| 1252 | |
| 1253 | if (isnear) { |
| 1254 | NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length()); |
| 1255 | Project(CNew, CurrentTi->Array1(), |
| 1256 | NewTi->ChangeArray1(), |
| 1257 | Ecarts, NumPnt, |
| 1258 | ERRMAX, ERRQUA, ERRMOY, 2); |
| 1259 | } |
| 1260 | else NewTi = CurrentTi; |
| 1261 | |
| 1262 | |
| 1263 | // (1.4) Progression's test |
| 1264 | iprog = 0; |
| 1265 | if ((EROLD > WQuality) && (ERRMAX < 0.95*EROLD)) iprog++; |
| 1266 | if ((EROLD > WQuality) && (ERRMAX < 0.8*EROLD)) iprog++; |
| 1267 | if ((EROLD > WQuality) && (ERRMAX < WQuality)) iprog++; |
| 1268 | if ((EROLD > WQuality) && (ERRMAX < 0.99*EROLD) |
| 1269 | && (ERRMAX < 1.1*WQuality)) iprog++; |
| 1270 | if ( VALCRI[0] < 0.975 * VOCRI[0]) iprog++; |
| 1271 | if ( VALCRI[0] < 0.9 * VOCRI[0]) iprog++; |
| 1272 | if ( VALCRI[1] < 0.95 * VOCRI[1]) iprog++; |
| 1273 | if ( VALCRI[1] < 0.8 * VOCRI[1]) iprog++; |
| 1274 | if ( VALCRI[2] < 0.95 * VOCRI[2]) iprog++; |
| 1275 | if ( VALCRI[2] < 0.8 * VOCRI[2]) iprog++; |
| 1276 | if ((VOCRI[1]>SmallestValue)&&(VOCRI[2]>SmallestValue)) { |
| 1277 | if ((VALCRI[1]/VOCRI[1] + 2*VALCRI[2]/VOCRI[2]) < 2.8) iprog++; |
| 1278 | } |
| 1279 | |
| 1280 | if (iprog < 2 && NbEst == 0 ) { |
| 1281 | // (1.5) Invalidation of new knots. |
| 1282 | VALCRI[0] = VOCRI[0]; |
| 1283 | VALCRI[1] = VOCRI[1]; |
| 1284 | VALCRI[2] = VOCRI[2]; |
| 1285 | ERRMAX = EROLD; |
| 1286 | CBLONG = LNOLD; |
| 1287 | CCurrent = COld; |
| 1288 | CurrentTi = OldTi; |
| 1289 | |
| 1290 | goto L8000; // exit |
| 1291 | } |
| 1292 | |
| 1293 | VOCRI[0] = VALCRI[0]; |
| 1294 | VOCRI[1] = VALCRI[1]; |
| 1295 | VOCRI[2] = VALCRI[2]; |
| 1296 | LNOLD = CBLONG; |
| 1297 | EROLD = ERRMAX; |
| 1298 | |
| 1299 | CCurrent = CNew; |
| 1300 | CurrentTi = NewTi; |
| 1301 | |
| 1302 | // (1.6) Test if the Estimations seems OK, else repeat |
| 1303 | NbEst++; |
| 1304 | lestim = ( (NbEst<MaxNbEst) && (ICDANA == 2)&& (iprog > 0) ); |
| 1305 | |
| 1306 | if (lestim && isnear) { |
| 1307 | // (1.7) Optimization of ti by ACR. |
| 1308 | |
| 1309 | std::stable_sort (CurrentTi->begin(), CurrentTi->end()); |
| 1310 | |
| 1311 | Standard_Integer Decima = 4; |
| 1312 | |
| 1313 | CCurrent->Length(0., 1., CBLONG); |
| 1314 | J->EstLength() = CBLONG; |
| 1315 | |
| 1316 | ACR(CCurrent, CurrentTi->ChangeArray1(), Decima); |
| 1317 | lconst = Standard_True; |
| 1318 | |
| 1319 | } |
| 1320 | } |
| 1321 | |
| 1322 | |
| 1323 | // (2) loop of parametric / geometric optimization |
| 1324 | |
| 1325 | Iter = 1; |
| 1326 | ToOptim = ((Iter < myNbIterations) && (isnear)); |
| 1327 | |
| 1328 | while(ToOptim) { |
| 1329 | Iter++; |
| 1330 | // (2.1) Save curent results |
| 1331 | VOCRI[0] = VALCRI[0]; |
| 1332 | VOCRI[1] = VALCRI[1]; |
| 1333 | VOCRI[2] = VALCRI[2]; |
| 1334 | EROLD = ERRMAX; |
| 1335 | LNOLD = CBLONG; |
| 1336 | COld = CCurrent; |
| 1337 | OldTi->ChangeArray1() = CurrentTi->Array1(); |
| 1338 | |
| 1339 | // (2.2) Optimization des ti by ACR. |
| 1340 | |
| 1341 | std::stable_sort (CurrentTi->begin(), CurrentTi->end()); |
| 1342 | |
| 1343 | Standard_Integer Decima = 4; |
| 1344 | |
| 1345 | CCurrent->Length(0., 1., CBLONG); |
| 1346 | J->EstLength() = CBLONG; |
| 1347 | |
| 1348 | ACR(CCurrent, CurrentTi->ChangeArray1(), Decima); |
| 1349 | lconst = Standard_True; |
| 1350 | |
| 1351 | // (2.3) Optimisation des courbes |
| 1352 | EpsLength = SmallValue * CBLONG / NbrPnt; |
| 1353 | |
| 1354 | CNew = new (FEmTool_Curve) (CCurrent->Dimension(), |
| 1355 | CCurrent->NbElements(), CCurrent->Base(), EpsLength); |
| 1356 | CNew->Knots() = CCurrent->Knots(); |
| 1357 | |
| 1358 | J->SetParameters(CurrentTi); |
| 1359 | |
| 1360 | EpsDeg = Min(WQuality * .1, CBLONG * .001); |
| 1361 | Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1()); |
| 1362 | |
| 1363 | CCurrent = CNew; |
| 1364 | |
| 1365 | // (2.4) calcul des criteres de qualites et amelioration |
| 1366 | // des estimation. |
| 1367 | |
| 1368 | ICDANA = J->QualityValues(J1min, J2min, J3min, VALCRI[0], VALCRI[1], VALCRI[2]); |
| 1369 | if(ICDANA > 0) lconst = Standard_True; |
| 1370 | |
| 1371 | J->GetEstimation(e1, e2, e3); |
| 1372 | // (2.5) Optimisation des ti par proj orthogonale |
| 1373 | |
| 1374 | NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length()); |
| 1375 | Project(CCurrent, CurrentTi->Array1(), |
| 1376 | NewTi->ChangeArray1(), |
| 1377 | Ecarts, NumPnt, |
| 1378 | ERRMAX, ERRQUA, ERRMOY, 2); |
| 1379 | |
| 1380 | // (2.6) Test de non regression |
| 1381 | |
| 1382 | Standard_Integer iregre = 0; |
| 1383 | if (NbrConstraint < NbrPnt) { |
| 1384 | if ( (ERRMAX > WQuality) && (ERRMAX > 1.05*EROLD)) iregre++; |
| 1385 | if ( (ERRMAX > WQuality) && (ERRMAX > 2*EROLD)) iregre++; |
| 1386 | if ( (EROLD > WQuality) && (ERRMAX <= 0.5*EROLD)) iregre--; |
| 1387 | } |
| 1388 | Standard_Real E1, E2, E3; |
| 1389 | J->GetEstimation(E1, E2, E3); |
| 1390 | if ( (VALCRI[0] > E1) && (VALCRI[0] > 1.1*VOCRI[0])) iregre++; |
| 1391 | if ( (VALCRI[1] > E2) && (VALCRI[1] > 1.1*VOCRI[1])) iregre++; |
| 1392 | if ( (VALCRI[2] > E3) && (VALCRI[2] > 1.1*VOCRI[2])) iregre++; |
| 1393 | |
| 1394 | |
| 1395 | if (iregre >= 2) { |
| 1396 | // if (iregre >= 1) { |
| 1397 | // (2.7) on restaure l'iteration precedente |
| 1398 | VALCRI[0] = VOCRI[0]; |
| 1399 | VALCRI[1] = VOCRI[1]; |
| 1400 | VALCRI[2] = VOCRI[2]; |
| 1401 | ERRMAX = EROLD; |
| 1402 | CBLONG = LNOLD; |
| 1403 | CCurrent = COld; |
| 1404 | CurrentTi->ChangeArray1() = OldTi->Array1(); |
| 1405 | ToOptim = Standard_False; |
| 1406 | } |
| 1407 | else { // Iteration is Ok. |
| 1408 | CCurrent = CNew; |
| 1409 | CurrentTi = NewTi; |
| 1410 | } |
| 1411 | if (Iter >= myNbIterations) ToOptim = Standard_False; |
| 1412 | } |
| 1413 | |
| 1414 | // (3) Decoupe eventuelle |
| 1415 | |
| 1416 | if ( (CCurrent->NbElements() < myMaxSegment) && myWithCutting ) { |
| 1417 | |
| 1418 | // (3.1) Sauvgarde de l'etat precedent |
| 1419 | VOCRI[0] = VALCRI[0]; |
| 1420 | VOCRI[1] = VALCRI[1]; |
| 1421 | VOCRI[2] = VALCRI[2]; |
| 1422 | EROLD = ERRMAX; |
| 1423 | COld = CCurrent; |
| 1424 | OldTi->ChangeArray1() = CurrentTi->Array1(); |
| 1425 | |
| 1426 | // (3.2) On arrange les ti : Trie + recadrage sur (0,1) |
| 1427 | // ---> On trie, afin d'assurer l'ordre par la suite. |
| 1428 | |
| 1429 | std::stable_sort (CurrentTi->begin(), CurrentTi->end()); |
| 1430 | |
| 1431 | if ((CurrentTi->Value(1)!= 0.) || |
| 1432 | (CurrentTi->Value(NbrPnt)!= 1.)) { |
| 1433 | Standard_Real t, DelatT = |
| 1434 | 1.0 /(CurrentTi->Value(NbrPnt)-CurrentTi->Value(1)); |
| 1435 | for (Standard_Integer ii=2; ii<NbrPnt; ii++) { |
| 1436 | t = (CurrentTi->Value(ii)-CurrentTi->Value(1))*DelatT; |
| 1437 | CurrentTi->SetValue(ii, t); |
| 1438 | } |
| 1439 | CurrentTi->SetValue(1, 0.); |
| 1440 | CurrentTi->SetValue(NbrPnt, 1.); |
| 1441 | } |
| 1442 | |
| 1443 | // (3.3) Insert new Knots |
| 1444 | |
| 1445 | SplitCurve(CCurrent, CurrentTi->Array1(), EpsLength, CNew, iscut); |
| 1446 | if (!iscut) again = Standard_False; |
| 1447 | else { |
| 1448 | CCurrent = CNew; |
| 1449 | // New Knots => New Assembly. |
| 1450 | J->SetCurve(CNew); |
| 1451 | delete TheAssembly; |
| 1452 | TheAssembly = new FEmTool_Assembly (Dependence->Array2(), |
| 1453 | J->AssemblyTable()); |
| 1454 | } |
| 1455 | } |
| 1456 | else { again = Standard_False;} |
| 1457 | } |
| 1458 | |
| 1459 | // ================ Great loop end =================== |
| 1460 | |
| 1461 | L8000: |
| 1462 | |
| 1463 | // (4) Compute the best Error. |
| 1464 | NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length()); |
| 1465 | Project(CCurrent, CurrentTi->Array1(), |
| 1466 | NewTi->ChangeArray1(), |
| 1467 | Ecarts, NumPnt, |
| 1468 | ERRMAX, ERRQUA, ERRMOY, 10); |
| 1469 | |
| 1470 | // (5) field's update |
| 1471 | |
| 1472 | TheCurve = CCurrent; |
| 1473 | J->EstLength() = CBLONG; |
| 1474 | myParameters->ChangeArray1() = NewTi->Array1(); |
| 1475 | myCriterium[0] = ERRQUA; |
| 1476 | myCriterium[1] = Sqrt(VALCRI[0]); |
| 1477 | myCriterium[2] = Sqrt(VALCRI[1]); |
| 1478 | myCriterium[3] = Sqrt(VALCRI[2]); |
| 1479 | myMaxError = ERRMAX; |
| 1480 | myMaxErrorIndex = NumPnt; |
| 1481 | if(NbrPnt > NbrConstraint) |
| 1482 | myAverageError = ERRMOY / (NbrPnt - NbrConstraint); |
| 1483 | else |
| 1484 | myAverageError = ERRMOY / NbrConstraint; |
| 1485 | |
| 1486 | delete TheAssembly; |
| 1487 | } |
| 1488 | |
| 1489 | //======================================================================= |
| 1490 | //function : Optimization |
| 1491 | //purpose : (like FORTRAN subroutine MINIMI) |
| 1492 | //======================================================================= |
| 1493 | void AppDef_Variational::Optimization(Handle(AppDef_SmoothCriterion)& J, |
| 1494 | FEmTool_Assembly& A, |
| 1495 | const Standard_Boolean ToAssemble, |
| 1496 | const Standard_Real EpsDeg, |
| 1497 | Handle(FEmTool_Curve)& Curve, |
| 1498 | const TColStd_Array1OfReal& Parameters) const |
| 1499 | { |
| 1500 | Standard_Integer MxDeg = Curve->Base()->WorkDegree(), |
| 1501 | NbElm = Curve->NbElements(), |
| 1502 | NbDim = Curve->Dimension(); |
| 1503 | |
| 1504 | Handle(FEmTool_HAssemblyTable) AssTable; |
| 1505 | |
| 1506 | math_Matrix H(0, MxDeg, 0, MxDeg); |
| 1507 | math_Vector G(0, MxDeg), Sol(1, A.NbGlobVar()); |
| 1508 | |
| 1509 | Standard_Integer el, dim; |
| 1510 | |
| 1511 | A.GetAssemblyTable(AssTable); |
| 1512 | Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
| 1513 | |
| 1514 | Standard_Real CBLONG = J->EstLength(); |
| 1515 | |
| 1516 | // Updating Assembly |
| 1517 | if (ToAssemble) |
| 1518 | A.NullifyMatrix(); |
| 1519 | A.NullifyVector(); |
| 1520 | |
| 1521 | |
| 1522 | for (el = 1; el <= NbElm; el++) { |
| 1523 | if (ToAssemble) { |
| 1524 | J->Hessian(el, 1, 1, H); |
| 1525 | for(dim = 1; dim <= NbDim; dim++) |
| 1526 | A.AddMatrix(el, dim, dim, H); |
| 1527 | } |
| 1528 | |
| 1529 | for(dim = 1; dim <= NbDim; dim++) { |
| 1530 | J->Gradient(el, dim, G); |
| 1531 | A.AddVector(el, dim, G); |
| 1532 | } |
| 1533 | } |
| 1534 | |
| 1535 | |
| 1536 | // Solution of system |
| 1537 | if (ToAssemble) { |
| 1538 | if(NbConstr != 0) { //Treatment of constraints |
| 1539 | AssemblingConstraints(Curve, Parameters, CBLONG, A); |
| 1540 | } |
| 1541 | A.Solve(); |
| 1542 | } |
| 1543 | A.Solution(Sol); |
| 1544 | |
| 1545 | |
| 1546 | // Updating J |
| 1547 | J->SetCurve(Curve); |
| 1548 | J->InputVector(Sol, AssTable); |
| 1549 | |
| 1550 | // Updating Curve and reduction of degree |
| 1551 | |
| 1552 | Standard_Integer Newdeg; |
| 1553 | Standard_Real MaxError; |
| 1554 | |
| 1555 | if(NbConstr == 0) { |
| 1556 | for(el = 1; el <= NbElm; el++) { |
| 1557 | Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError); |
| 1558 | } |
| 1559 | } |
| 1560 | else { |
| 1561 | |
| 1562 | TColStd_Array1OfReal& TabInt = Curve->Knots(); |
| 1563 | Standard_Integer Icnt = 1, p0 = Parameters.Lower() - myFirstPoint, point; |
| 1564 | for(el = 1; el <= NbElm; el++) { |
| 1565 | while((Icnt < NbConstr) && |
| 1566 | (Parameters(p0 + myTypConstraints->Value(2 * Icnt - 1)) <= TabInt(el))) Icnt++; |
| 1567 | point = p0 + myTypConstraints->Value(2 * Icnt - 1); |
| 1568 | if(Parameters(point) <= TabInt(el) || Parameters(point) >= TabInt(el + 1)) |
| 1569 | Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError); |
| 1570 | else |
| 1571 | if(Curve->Degree(el) < MxDeg) Curve->SetDegree(el, MxDeg); |
| 1572 | } |
| 1573 | } |
| 1574 | } |
| 1575 | |
| 1576 | void AppDef_Variational::Project(const Handle(FEmTool_Curve)& C, |
| 1577 | const TColStd_Array1OfReal& Ti, |
| 1578 | TColStd_Array1OfReal& ProjTi, |
| 1579 | TColStd_Array1OfReal& Distance, |
| 1580 | Standard_Integer& NumPoints, |
| 1581 | Standard_Real& MaxErr, |
| 1582 | Standard_Real& QuaErr, |
| 1583 | Standard_Real& AveErr, |
| 1584 | const Standard_Integer NbIterations) const |
| 1585 | |
| 1586 | { |
| 1587 | // Initialisation |
| 1588 | |
| 1589 | const Standard_Real Seuil = 1.e-9, Eps = 1.e-12; |
| 1590 | |
| 1591 | MaxErr = QuaErr = AveErr = 0.; |
| 1592 | |
| 1593 | Standard_Integer Ipnt, NItCv, Iter, i, i0 = -myDimension, d0 = Distance.Lower() - 1; |
| 1594 | |
| 1595 | Standard_Real TNew, Dist, T0, Dist0, F1, F2, Aux, DF, Ecart; |
| 1596 | |
| 1597 | Standard_Boolean EnCour; |
| 1598 | |
| 1599 | TColStd_Array1OfReal ValOfC(1, myDimension), FirstDerOfC(1, myDimension), |
| 1600 | SecndDerOfC(1, myDimension); |
| 1601 | |
| 1602 | for(Ipnt = 1; Ipnt <= ProjTi.Length(); Ipnt++) { |
| 1603 | |
| 1604 | i0 += myDimension; |
| 1605 | |
| 1606 | TNew = Ti(Ipnt); |
| 1607 | |
| 1608 | EnCour = Standard_True; |
| 1609 | NItCv = 0; |
| 1610 | Iter = 0; |
| 1611 | C->D0(TNew, ValOfC); |
| 1612 | |
| 1613 | Dist = 0; |
| 1614 | for(i = 1; i <= myDimension; i++) { |
| 1615 | Aux = ValOfC(i) - myTabPoints->Value(i0 + i); |
| 1616 | Dist += Aux * Aux; |
| 1617 | } |
| 1618 | Dist = Sqrt(Dist); |
| 1619 | |
| 1620 | // ------- Newton's method for solving (C'(t),C(t) - P) = 0 |
| 1621 | |
| 1622 | while( EnCour ) { |
| 1623 | |
| 1624 | Iter++; |
| 1625 | T0 = TNew; |
| 1626 | Dist0 = Dist; |
| 1627 | |
| 1628 | C->D2(TNew, SecndDerOfC); |
| 1629 | C->D1(TNew, FirstDerOfC); |
| 1630 | |
| 1631 | F1 = F2 = 0.; |
| 1632 | for(i = 1; i <= myDimension; i++) { |
| 1633 | Aux = ValOfC(i) - myTabPoints->Value(i0 + i); |
| 1634 | DF = FirstDerOfC(i); |
| 1635 | F1 += Aux*DF; // (C'(t),C(t) - P) |
| 1636 | F2 += DF*DF + Aux * SecndDerOfC(i); // ((C'(t),C(t) - P))' |
| 1637 | } |
| 1638 | |
| 1639 | if(Abs(F2) < Eps) |
| 1640 | EnCour = Standard_False; |
| 1641 | else { |
| 1642 | // Formula of Newton x(k+1) = x(k) - F(x(k))/F'(x(k)) |
| 1643 | TNew -= F1 / F2; |
| 1644 | if(TNew < 0.) TNew = 0.; |
| 1645 | if(TNew > 1.) TNew = 1.; |
| 1646 | |
| 1647 | |
| 1648 | // Analysis of result |
| 1649 | |
| 1650 | C->D0(TNew, ValOfC); |
| 1651 | |
| 1652 | Dist = 0; |
| 1653 | for(i = 1; i <= myDimension; i++) { |
| 1654 | Aux = ValOfC(i) - myTabPoints->Value(i0 + i); |
| 1655 | Dist += Aux * Aux; |
| 1656 | } |
| 1657 | Dist = Sqrt(Dist); |
| 1658 | |
| 1659 | Ecart = Dist0 - Dist; |
| 1660 | |
| 1661 | if(Ecart <= -Seuil) { |
| 1662 | // Pas d'amelioration on s'arrete |
| 1663 | EnCour = Standard_False; |
| 1664 | TNew = T0; |
| 1665 | Dist = Dist0; |
| 1666 | } |
| 1667 | else if(Ecart <= Seuil) |
| 1668 | // Convergence |
| 1669 | NItCv++; |
| 1670 | else |
| 1671 | NItCv = 0; |
| 1672 | |
| 1673 | if((NItCv >= 2) || (Iter >= NbIterations)) EnCour = Standard_False; |
| 1674 | |
| 1675 | } |
| 1676 | } |
| 1677 | |
| 1678 | |
| 1679 | ProjTi(Ipnt) = TNew; |
| 1680 | Distance(d0 + Ipnt) = Dist; |
| 1681 | if(Dist > MaxErr) { |
| 1682 | MaxErr = Dist; |
| 1683 | NumPoints = Ipnt; |
| 1684 | } |
| 1685 | QuaErr += Dist * Dist; |
| 1686 | AveErr += Dist; |
| 1687 | } |
| 1688 | |
| 1689 | NumPoints = NumPoints + myFirstPoint - 1;// Setting NumPoints to interval [myFirstPoint, myLastPoint] |
| 1690 | |
| 1691 | } |
| 1692 | |
| 1693 | void AppDef_Variational::ACR(Handle(FEmTool_Curve)& Curve, |
| 1694 | TColStd_Array1OfReal& Ti, |
| 1695 | const Standard_Integer Decima) const |
| 1696 | { |
| 1697 | |
| 1698 | const Standard_Real Eps = 1.e-8; |
| 1699 | |
| 1700 | TColStd_Array1OfReal& Knots = Curve->Knots(); |
| 1701 | Standard_Integer NbrPnt = Ti.Length(), TiFirst = Ti.Lower(), TiLast = Ti.Upper(), |
| 1702 | KFirst = Knots.Lower(), KLast = Knots.Upper(); |
| 1703 | |
| 1704 | Standard_Real CbLong, DeltaT, VTest, UNew, UOld, DU, TPara, TOld, DTInv, Ratio; |
| 1705 | Standard_Integer ipnt, ii, IElm, IOld, POld, PCnt, ICnt=0; |
| 1706 | Standard_Integer NbCntr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
| 1707 | |
| 1708 | // (1) Calcul de la longueur de courbe |
| 1709 | |
| 1710 | Curve->Length(Ti(TiFirst), Ti(TiLast), CbLong); |
| 1711 | |
| 1712 | // (2) Mise de l'acr dans Ti |
| 1713 | |
| 1714 | if(NbrPnt >= 2) { |
| 1715 | |
| 1716 | // (2.0) Initialisation |
| 1717 | DeltaT = (Ti(TiLast) - Ti(TiFirst)) / Decima; |
| 1718 | VTest = Ti(TiFirst) + DeltaT; |
| 1719 | |
| 1720 | if(NbCntr > 0) { |
| 1721 | PCnt = myTypConstraints->Value(1) - myFirstPoint + TiFirst; |
| 1722 | ICnt = 1; |
| 1723 | } |
| 1724 | else |
| 1725 | PCnt = TiLast + 1; |
| 1726 | |
| 1727 | UOld = 0.; |
| 1728 | |
| 1729 | TOld = Ti(TiFirst); |
| 1730 | POld = TiFirst; |
| 1731 | |
| 1732 | IElm = KFirst; |
| 1733 | IOld = IElm; |
| 1734 | |
| 1735 | Ti(TiFirst) = 0.; |
| 1736 | |
| 1737 | for(ipnt = TiFirst + 1; ipnt <= TiLast; ipnt++) { |
| 1738 | |
| 1739 | while((ICnt <= NbCntr) && (PCnt < ipnt)) { |
| 1740 | ICnt++; |
| 1741 | PCnt = myTypConstraints->Value(2*ICnt-1) - myFirstPoint + TiFirst; |
| 1742 | } |
| 1743 | |
| 1744 | TPara = Ti(ipnt); |
| 1745 | |
| 1746 | if(TPara >= VTest || PCnt == ipnt) { |
| 1747 | |
| 1748 | if ( Ti(TiLast) - TPara <= 1.e-2*DeltaT) { |
| 1749 | ipnt = TiLast; |
| 1750 | TPara = Ti(ipnt); |
| 1751 | } |
| 1752 | // (2.2), (2.3) Cacul la longueur de courbe |
| 1753 | Curve->Length(Ti(TiFirst), TPara, UNew); |
| 1754 | |
| 1755 | UNew /= CbLong; |
| 1756 | |
| 1757 | while(Knots(IElm + 1) < TPara && IElm < KLast - 1) IElm++; |
| 1758 | |
| 1759 | // (2.4) Mise a jours des parametres de decoupe |
| 1760 | DTInv = 1. / (TPara - TOld); |
| 1761 | DU = UNew - UOld; |
| 1762 | |
| 1763 | for(ii = IOld+1; ii <= IElm; ii++) { |
| 1764 | Ratio = (Knots(ii) - TOld) * DTInv; |
| 1765 | Knots(ii) = UOld + Ratio * DU; |
| 1766 | } |
| 1767 | |
| 1768 | // (2.5) Mise a jours des parametres de points. |
| 1769 | |
| 1770 | //Very strange loop, because it never works (POld+1 > ipnt-1) |
| 1771 | for(ii = POld+1; ii <= ipnt-1; ii++) { |
| 1772 | Ratio = ( Ti(ii) - TOld ) * DTInv; |
| 1773 | Ti(ii) = UOld + Ratio * DU; |
| 1774 | } |
| 1775 | |
| 1776 | Ti(ipnt) = UNew; |
| 1777 | |
| 1778 | UOld = UNew; |
| 1779 | IOld = IElm; |
| 1780 | TOld = TPara; |
| 1781 | POld = ipnt; |
| 1782 | |
| 1783 | } |
| 1784 | // --> Nouveau seuil parametrique pour le decimage |
| 1785 | |
| 1786 | if(TPara >= VTest) { |
| 1787 | // ii = RealToInt((TPara - VTest + Eps) / DeltaT); |
| 1788 | // VTest += (ii + 1) * DeltaT; |
| 1789 | VTest += Ceiling((TPara - VTest + Eps) / DeltaT) * DeltaT; |
| 1790 | if(VTest > 1. - Eps) VTest = 1.; |
| 1791 | } |
| 1792 | } |
| 1793 | } |
| 1794 | |
| 1795 | // --- On ajuste les valeurs extremes |
| 1796 | |
| 1797 | Ti(TiFirst) = 0.; |
| 1798 | Ti(TiLast) = 1.; |
| 1799 | ii = TiLast - 1; |
| 1800 | while ( Ti(ii) > Knots(KLast) ) { |
| 1801 | Ti(ii) = 1.; |
| 1802 | --ii; |
| 1803 | } |
| 1804 | Knots(KFirst) = 0.; |
| 1805 | Knots(KLast) = 1.; |
| 1806 | |
| 1807 | } |
| 1808 | |
| 1809 | //----------------------------------------------------------// |
| 1810 | // Standard_Integer NearIndex // |
| 1811 | // Purpose: searching nearest index of TabPar corresponding // |
| 1812 | // given value T (is similar to fortran routine MSRRE2) // |
| 1813 | //----------------------------------------------------------// |
| 1814 | |
| 1815 | static Standard_Integer NearIndex(const Standard_Real T, |
| 1816 | const TColStd_Array1OfReal& TabPar, |
| 1817 | const Standard_Real Eps, Standard_Integer& Flag) |
| 1818 | { |
| 1819 | Standard_Integer Loi = TabPar.Lower(), Upi = TabPar.Upper(); |
| 1820 | |
| 1821 | Flag = 0; |
| 1822 | |
| 1823 | if(T < TabPar(Loi)) { Flag = -1; return Loi;} |
| 1824 | if(T > TabPar(Upi)) { Flag = 1; return Upi;} |
| 1825 | |
| 1826 | Standard_Integer Ibeg = Loi, Ifin = Upi, Imidl; |
| 1827 | |
| 1828 | while(Ibeg + 1 != Ifin) { |
| 1829 | Imidl = (Ibeg + Ifin) / 2; |
| 1830 | if((T >= TabPar(Ibeg)) && (T <= TabPar(Imidl))) |
| 1831 | Ifin = Imidl; |
| 1832 | else |
| 1833 | Ibeg = Imidl; |
| 1834 | } |
| 1835 | |
| 1836 | if(Abs(T - TabPar(Ifin)) < Eps) return Ifin; |
| 1837 | |
| 1838 | return Ibeg; |
| 1839 | } |
| 1840 | |
| 1841 | |
| 1842 | //----------------------------------------------------------// |
| 1843 | // void GettingKnots // |
| 1844 | // Purpose: calculating values of new Knots for elements // |
| 1845 | // with degree that is equal Deg // |
| 1846 | //----------------------------------------------------------// |
| 1847 | |
| 1848 | static void GettingKnots(const TColStd_Array1OfReal& TabPar, |
| 1849 | const Handle(FEmTool_Curve)& InCurve, |
| 1850 | const Standard_Integer Deg, |
| 1851 | Standard_Integer& NbElm, |
| 1852 | TColStd_Array1OfReal& NewKnots) |
| 1853 | |
| 1854 | { |
| 1855 | |
| 1856 | const Standard_Real Eps = 1.e-12; |
| 1857 | |
| 1858 | TColStd_Array1OfReal& OldKnots = InCurve->Knots(); |
| 1859 | Standard_Integer NbMaxOld = InCurve->NbElements(); |
| 1860 | Standard_Integer NbMax = NewKnots.Upper(), Ipt, Ipt1, Ipt2; |
| 1861 | Standard_Integer el = 0, i1 = OldKnots.Lower(), i0 = i1 - 1, Flag; |
| 1862 | Standard_Real TPar; |
| 1863 | |
| 1864 | while((NbElm < NbMax) && (el < NbMaxOld)) { |
| 1865 | |
| 1866 | el++; i0++; i1++; // i0, i1 are indexes of left and right knots of element el |
| 1867 | |
| 1868 | if(InCurve->Degree(el) == Deg) { |
| 1869 | |
| 1870 | NbElm++; |
| 1871 | |
| 1872 | Ipt1 = NearIndex(OldKnots(i0), TabPar, Eps, Flag); |
| 1873 | if(Flag != 0) Ipt1 = TabPar.Lower(); |
| 1874 | Ipt2 = NearIndex(OldKnots(i1), TabPar, Eps, Flag); |
| 1875 | if(Flag != 0) Ipt2 = TabPar.Upper(); |
| 1876 | |
| 1877 | if(Ipt2 - Ipt1 >= 1) { |
| 1878 | |
| 1879 | Ipt = (Ipt1 + Ipt2) / 2; |
| 1880 | if(2 * Ipt == Ipt1 + Ipt2) |
| 1881 | TPar = 2. * TabPar(Ipt); |
| 1882 | else |
| 1883 | TPar = TabPar(Ipt) + TabPar(Ipt + 1); |
| 1884 | |
| 1885 | NewKnots(NbElm) = (OldKnots(i0) + OldKnots(i1) + TPar) / 4.; |
| 1886 | } |
| 1887 | else |
| 1888 | NewKnots(NbElm) = (OldKnots(i0) + OldKnots(i1)) / 2.; |
| 1889 | } |
| 1890 | } |
| 1891 | } |
| 1892 | |
| 1893 | void AppDef_Variational::SplitCurve(const Handle(FEmTool_Curve)& InCurve, |
| 1894 | const TColStd_Array1OfReal& Ti, |
| 1895 | const Standard_Real CurveTol, |
| 1896 | Handle(FEmTool_Curve)& OutCurve, |
| 1897 | Standard_Boolean& iscut) const |
| 1898 | { |
| 1899 | Standard_Integer NbElmOld = InCurve->NbElements(); |
| 1900 | |
| 1901 | if(NbElmOld >= myMaxSegment) {iscut = Standard_False; return;} |
| 1902 | #ifdef OCCT_DEBUG |
| 1903 | Standard_Integer MaxDegree = |
| 1904 | #endif |
| 1905 | InCurve->Base()->WorkDegree(); |
| 1906 | Standard_Integer NbElm = NbElmOld; |
| 1907 | TColStd_Array1OfReal NewKnots(NbElm + 1, myMaxSegment); |
| 1908 | #ifndef OCCT_DEBUG |
| 1909 | GettingKnots(Ti, InCurve, InCurve->Base()->WorkDegree(), NbElm, NewKnots); |
| 1910 | GettingKnots(Ti, InCurve, InCurve->Base()->WorkDegree() - 1, NbElm, NewKnots); |
| 1911 | #else |
| 1912 | GettingKnots(Ti, InCurve, MaxDegree, NbElm, NewKnots); |
| 1913 | GettingKnots(Ti, InCurve, MaxDegree - 1, NbElm, NewKnots); |
| 1914 | |
| 1915 | #endif |
| 1916 | if(NbElm > NbElmOld) { |
| 1917 | |
| 1918 | iscut = Standard_True; |
| 1919 | |
| 1920 | OutCurve = new FEmTool_Curve(InCurve->Dimension(), NbElm, InCurve->Base(), CurveTol); |
| 1921 | TColStd_Array1OfReal& OutKnots = OutCurve->Knots(); |
| 1922 | TColStd_Array1OfReal& InKnots = InCurve->Knots(); |
| 1923 | |
| 1924 | Standard_Integer i, i0 = OutKnots.Lower(); |
| 1925 | for(i = InKnots.Lower(); i <= InKnots.Upper(); i++) OutKnots(i) = InKnots(i); |
| 1926 | for(i = NbElmOld + 1; i <= NbElm; i++) OutKnots(i + i0) = NewKnots(i); |
| 1927 | |
| 1928 | std::sort (OutKnots.begin(), OutKnots.end()); |
| 1929 | } |
| 1930 | else |
| 1931 | iscut = Standard_False; |
| 1932 | |
| 1933 | } |
| 1934 | |
| 1935 | //======================================================================= |
| 1936 | //function : InitSmoothCriterion |
| 1937 | //purpose : Initializes the SmoothCriterion |
| 1938 | //======================================================================= |
| 1939 | void AppDef_Variational::InitSmoothCriterion() |
| 1940 | { |
| 1941 | |
| 1942 | const Standard_Real Eps2 = 1.e-6, Eps3 = 1.e-9; |
| 1943 | // const Standard_Real J1 = .01, J2 = .001, J3 = .001; |
| 1944 | |
| 1945 | |
| 1946 | |
| 1947 | Standard_Real Length; |
| 1948 | |
| 1949 | InitParameters(Length); |
| 1950 | |
| 1951 | mySmoothCriterion->SetParameters(myParameters); |
| 1952 | |
| 1953 | Standard_Real E1, E2, E3; |
| 1954 | |
| 1955 | InitCriterionEstimations(Length, E1, E2, E3); |
| 1956 | /* |
| 1957 | J1 = 1.e-8; J2 = J3 = (E1 + 1.e-8) * 1.e-6; |
| 1958 | |
| 1959 | if(E1 < J1) E1 = J1; |
| 1960 | if(E2 < J2) E2 = J2; |
| 1961 | if(E3 < J3) E3 = J3; |
| 1962 | */ |
| 1963 | mySmoothCriterion->EstLength() = Length; |
| 1964 | mySmoothCriterion->SetEstimation(E1, E2, E3); |
| 1965 | |
| 1966 | Standard_Real WQuadratic, WQuality; |
| 1967 | |
| 1968 | if(!myWithMinMax && myTolerance != 0.) |
| 1969 | WQuality = myTolerance; |
| 1970 | else if (myTolerance == 0.) |
| 1971 | WQuality = 1.; |
| 1972 | else |
| 1973 | WQuality = Max(myTolerance, Eps2 * Length); |
| 1974 | |
| 1975 | Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
| 1976 | WQuadratic = Sqrt((Standard_Real)(myNbPoints - NbConstr)) * WQuality; |
| 1977 | if(WQuadratic > Eps3) WQuadratic = 1./ WQuadratic; |
| 1978 | |
| 1979 | if(WQuadratic == 0.) WQuadratic = Max(Sqrt(E1), 1.); |
| 1980 | |
| 1981 | |
| 1982 | mySmoothCriterion->SetWeight(WQuadratic, WQuality, |
| 1983 | myPercent[0], myPercent[1], myPercent[2]); |
| 1984 | |
| 1985 | |
| 1986 | Handle(PLib_Base) TheBase = new PLib_HermitJacobi(myMaxDegree, myContinuity); |
| 1987 | Handle(FEmTool_Curve) TheCurve; |
| 1988 | Standard_Integer NbElem; |
| 1989 | Standard_Real CurvTol = Eps2 * Length / myNbPoints; |
| 1990 | |
| 1991 | // Decoupe de l'intervalle en fonction des contraintes |
| 1992 | if(myWithCutting == Standard_True && NbConstr != 0) { |
| 1993 | |
| 1994 | InitCutting(TheBase, CurvTol, TheCurve); |
| 1995 | |
| 1996 | } |
| 1997 | else { |
| 1998 | |
| 1999 | NbElem = 1; |
| 2000 | TheCurve = new FEmTool_Curve(myDimension, NbElem, TheBase, CurvTol); |
| 2001 | TheCurve->Knots().SetValue(TheCurve->Knots().Lower(), |
| 2002 | myParameters->Value(myFirstPoint)); |
| 2003 | TheCurve->Knots().SetValue(TheCurve->Knots().Upper(), |
| 2004 | myParameters->Value(myLastPoint)); |
| 2005 | |
| 2006 | } |
| 2007 | |
| 2008 | mySmoothCriterion->SetCurve(TheCurve); |
| 2009 | |
| 2010 | return; |
| 2011 | |
| 2012 | } |
| 2013 | |
| 2014 | // |
| 2015 | //======================================================================= |
| 2016 | //function : InitParameters |
| 2017 | //purpose : Calculation of initial estimation of the multicurve's length |
| 2018 | // and parameters for multipoints. |
| 2019 | //======================================================================= |
| 2020 | // |
| 2021 | void AppDef_Variational::InitParameters(Standard_Real& Length) |
| 2022 | { |
| 2023 | |
| 2024 | const Standard_Real Eps1 = Precision::Confusion() * .01; |
| 2025 | |
| 2026 | Standard_Real aux, dist; |
| 2027 | Standard_Integer i, i0, i1 = 0, ipoint; |
| 2028 | |
| 2029 | |
| 2030 | Length = 0.; |
| 2031 | myParameters->SetValue(myFirstPoint, Length); |
| 2032 | |
| 2033 | for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint; ipoint++) { |
| 2034 | i0 = i1; i1 += myDimension; |
| 2035 | dist = 0; |
| 2036 | for(i = 1; i <= myDimension; i++) { |
| 2037 | aux = myTabPoints->Value(i1 + i) - myTabPoints->Value(i0 + i); |
| 2038 | dist += aux * aux; |
| 2039 | } |
| 2040 | Length += Sqrt(dist); |
| 2041 | myParameters->SetValue(ipoint, Length); |
| 2042 | } |
| 2043 | |
| 2044 | |
| 2045 | if(Length <= Eps1) |
| 2046 | Standard_ConstructionError::Raise("AppDef_Variational::InitParameters"); |
| 2047 | |
| 2048 | |
| 2049 | for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint - 1; ipoint++) |
| 2050 | myParameters->SetValue(ipoint, myParameters->Value(ipoint) / Length); |
| 2051 | |
| 2052 | myParameters->SetValue(myLastPoint, 1.); |
| 2053 | |
| 2054 | // Avec peu de point il y a sans doute sous estimation ... |
| 2055 | if(myNbPoints < 10) Length *= (1. + 0.1 / (myNbPoints - 1)); |
| 2056 | } |
| 2057 | |
| 2058 | // |
| 2059 | //======================================================================= |
| 2060 | //function : InitCriterionEstimations |
| 2061 | //purpose : Calculation of initial estimation of the linear criteria |
| 2062 | // |
| 2063 | //======================================================================= |
| 2064 | // |
| 2065 | void AppDef_Variational::InitCriterionEstimations(const Standard_Real Length, |
| 2066 | Standard_Real& E1, |
| 2067 | Standard_Real& E2, |
| 2068 | Standard_Real& E3) const |
| 2069 | { |
| 2070 | E1 = Length * Length; |
| 2071 | |
| 2072 | const Standard_Real Eps1 = Precision::Confusion() * .01; |
| 2073 | |
| 2074 | math_Vector VTang1(1, myDimension), VTang2(1, myDimension), VTang3(1, myDimension), |
| 2075 | VScnd1(1, myDimension), VScnd2(1, myDimension), VScnd3(1, myDimension); |
| 2076 | |
| 2077 | // ========== Treatment of first point ================= |
| 2078 | |
| 2079 | Standard_Integer ipnt = myFirstPoint; |
| 2080 | |
| 2081 | EstTangent(ipnt, VTang1); |
| 2082 | ipnt++; |
| 2083 | EstTangent(ipnt, VTang2); |
| 2084 | ipnt++; |
| 2085 | EstTangent(ipnt, VTang3); |
| 2086 | |
| 2087 | ipnt = myFirstPoint; |
| 2088 | EstSecnd(ipnt, VTang1, VTang2, Length, VScnd1); |
| 2089 | ipnt++; |
| 2090 | EstSecnd(ipnt, VTang1, VTang3, Length, VScnd2); |
| 2091 | |
| 2092 | // Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 Begin |
| 2093 | // Standard_Real Delta = .5 * (myParameters->Value(ipnt) - myParameters->Value(--ipnt)); |
| 2094 | Standard_Integer anInd = ipnt; |
| 2095 | Standard_Real Delta = .5 * (myParameters->Value(anInd) - myParameters->Value(--ipnt)); |
| 2096 | // Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 End |
| 2097 | |
| 2098 | if(Delta <= Eps1) Delta = 1.; |
| 2099 | |
| 2100 | E2 = VScnd1.Norm2() * Delta; |
| 2101 | |
| 2102 | E3 = (Delta > Eps1) ? VScnd2.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.; |
| 2103 | // ========== Treatment of internal points ================= |
| 2104 | |
| 2105 | Standard_Integer CurrPoint = 2; |
| 2106 | |
| 2107 | for(ipnt = myFirstPoint + 1; ipnt < myLastPoint; ipnt++) { |
| 2108 | |
| 2109 | Delta = .5 * (myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1)); |
| 2110 | |
| 2111 | if(CurrPoint == 1) { |
| 2112 | if(ipnt + 1 != myLastPoint) { |
| 2113 | EstTangent(ipnt + 2, VTang3); |
| 2114 | EstSecnd(ipnt + 1, VTang1, VTang3, Length, VScnd2); |
| 2115 | } |
| 2116 | else |
| 2117 | EstSecnd(ipnt + 1, VTang1, VTang2, Length, VScnd2); |
| 2118 | |
| 2119 | E2 += VScnd1.Norm2() * Delta; |
| 2120 | E3 += (Delta > Eps1) ? VScnd2.Subtracted(VScnd3).Norm2() / (4. * Delta) : 0.; |
| 2121 | |
| 2122 | } |
| 2123 | else if(CurrPoint == 2) { |
| 2124 | if(ipnt + 1 != myLastPoint) { |
| 2125 | EstTangent(ipnt + 2, VTang1); |
| 2126 | EstSecnd(ipnt + 1, VTang2, VTang1, Length, VScnd3); |
| 2127 | } |
| 2128 | else |
| 2129 | EstSecnd(ipnt + 1, VTang2, VTang3, Length, VScnd3); |
| 2130 | |
| 2131 | E2 += VScnd2.Norm2() * Delta; |
| 2132 | E3 += (Delta > Eps1) ? VScnd3.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.; |
| 2133 | |
| 2134 | } |
| 2135 | else { |
| 2136 | if(ipnt + 1 != myLastPoint) { |
| 2137 | EstTangent(ipnt + 2, VTang2); |
| 2138 | EstSecnd(ipnt + 1, VTang3, VTang2, Length, VScnd1); |
| 2139 | } |
| 2140 | else |
| 2141 | EstSecnd(ipnt + 1, VTang3, VTang1, Length, VScnd1); |
| 2142 | |
| 2143 | E2 += VScnd3.Norm2() * Delta; |
| 2144 | E3 += (Delta > Eps1) ? VScnd1.Subtracted(VScnd2).Norm2() / (4. * Delta) : 0.; |
| 2145 | |
| 2146 | } |
| 2147 | |
| 2148 | CurrPoint++; if(CurrPoint == 4) CurrPoint = 1; |
| 2149 | } |
| 2150 | |
| 2151 | // ========== Treatment of last point ================= |
| 2152 | |
| 2153 | Delta = .5 * (myParameters->Value(myLastPoint) - myParameters->Value(myLastPoint - 1)); |
| 2154 | if(Delta <= Eps1) Delta = 1.; |
| 2155 | |
| 2156 | Standard_Real aux; |
| 2157 | |
| 2158 | if(CurrPoint == 1) { |
| 2159 | |
| 2160 | E2 += VScnd1.Norm2() * Delta; |
| 2161 | aux = VScnd1.Subtracted(VScnd3).Norm2(); |
| 2162 | E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux; |
| 2163 | |
| 2164 | } |
| 2165 | else if(CurrPoint == 2) { |
| 2166 | |
| 2167 | E2 += VScnd2.Norm2() * Delta; |
| 2168 | aux = VScnd2.Subtracted(VScnd1).Norm2(); |
| 2169 | E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux; |
| 2170 | |
| 2171 | } |
| 2172 | else { |
| 2173 | |
| 2174 | E2 += VScnd3.Norm2() * Delta; |
| 2175 | aux = VScnd3.Subtracted(VScnd2).Norm2(); |
| 2176 | E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux; |
| 2177 | |
| 2178 | } |
| 2179 | |
| 2180 | aux = Length * Length; |
| 2181 | |
| 2182 | E2 *= aux; E3 *= aux; |
| 2183 | |
| 2184 | |
| 2185 | } |
| 2186 | |
| 2187 | // |
| 2188 | //======================================================================= |
| 2189 | //function : EstTangent |
| 2190 | //purpose : Calculation of estimation of the Tangent |
| 2191 | // (see fortran routine MMLIPRI) |
| 2192 | //======================================================================= |
| 2193 | // |
| 2194 | |
| 2195 | void AppDef_Variational::EstTangent(const Standard_Integer ipnt, |
| 2196 | math_Vector& VTang) const |
| 2197 | |
| 2198 | { |
| 2199 | Standard_Integer i ; |
| 2200 | const Standard_Real Eps1 = Precision::Confusion() * .01; |
| 2201 | const Standard_Real EpsNorm = 1.e-9; |
| 2202 | |
| 2203 | Standard_Real Wpnt = 1.; |
| 2204 | |
| 2205 | |
| 2206 | if(ipnt == myFirstPoint) { |
| 2207 | // Estimation at first point |
| 2208 | if(myNbPoints < 3) |
| 2209 | Wpnt = 0.; |
| 2210 | else { |
| 2211 | |
| 2212 | Standard_Integer adr1 = 1, adr2 = adr1 + myDimension, |
| 2213 | adr3 = adr2 + myDimension; |
| 2214 | |
| 2215 | math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension); |
| 2216 | math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension); |
| 2217 | math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension); |
| 2218 | |
| 2219 | // Parabolic interpolation |
| 2220 | // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t, |
| 2221 | // first derivative for t=0 is A1 = ((d2-1)*P1 + P2 - d2*P3)/(d*(1-d)) |
| 2222 | // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*d |
| 2223 | Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm(); |
| 2224 | Standard_Real V2 = 0.; |
| 2225 | if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm(); |
| 2226 | if(V2 > Eps1) { |
| 2227 | Standard_Real d = V1 / (V1 + V2), d1; |
| 2228 | d1 = 1. / (d * (1 - d)); d *= d; |
| 2229 | VTang = ((d - 1.) * Pnt1 + Pnt2 - d * Pnt3) * d1; |
| 2230 | } |
| 2231 | else { |
| 2232 | // Simple 2-point estimation |
| 2233 | |
| 2234 | VTang = Pnt2 - Pnt1; |
| 2235 | } |
| 2236 | } |
| 2237 | } |
| 2238 | else if (ipnt == myLastPoint) { |
| 2239 | // Estimation at last point |
| 2240 | if(myNbPoints < 3) |
| 2241 | Wpnt = 0.; |
| 2242 | else { |
| 2243 | |
| 2244 | Standard_Integer adr1 = (myLastPoint - 3) * myDimension + 1, |
| 2245 | adr2 = adr1 + myDimension, |
| 2246 | adr3 = adr2 + myDimension; |
| 2247 | |
| 2248 | math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension); |
| 2249 | math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension); |
| 2250 | math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension); |
| 2251 | |
| 2252 | // Parabolic interpolation |
| 2253 | // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t, |
| 2254 | // first derivative for t=1 is 2*A2 + A1 = ((d2+1)*P1 - P2 - d2*P3)/(d*(1-d)) |
| 2255 | // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*(d-2) |
| 2256 | Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm(); |
| 2257 | Standard_Real V2 = 0.; |
| 2258 | if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm(); |
| 2259 | if(V2 > Eps1) { |
| 2260 | Standard_Real d = V1 / (V1 + V2), d1; |
| 2261 | d1 = 1. / (d * (1 - d)); d *= d - 2; |
| 2262 | VTang = ((d + 1.) * Pnt1 - Pnt2 - d * Pnt3) * d1; |
| 2263 | } |
| 2264 | else { |
| 2265 | // Simple 2-point estimation |
| 2266 | |
| 2267 | VTang = Pnt3 - Pnt2; |
| 2268 | } |
| 2269 | } |
| 2270 | } |
| 2271 | else { |
| 2272 | |
| 2273 | Standard_Integer adr1 = (ipnt - myFirstPoint - 1) * myDimension + 1, |
| 2274 | adr2 = adr1 + 2 * myDimension; |
| 2275 | |
| 2276 | math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension); |
| 2277 | math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension); |
| 2278 | |
| 2279 | VTang = Pnt2 - Pnt1; |
| 2280 | |
| 2281 | } |
| 2282 | |
| 2283 | Standard_Real Vnorm = VTang.Norm(); |
| 2284 | |
| 2285 | if(Vnorm <= EpsNorm) |
| 2286 | VTang.Init(0.); |
| 2287 | else |
| 2288 | VTang /= Vnorm; |
| 2289 | |
| 2290 | // Estimation with constraints |
| 2291 | |
| 2292 | Standard_Real Wcnt = 0.; |
| 2293 | Standard_Integer IdCnt = 1; |
| 2294 | |
| 2295 | // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint] |
| 2296 | |
| 2297 | Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
| 2298 | math_Vector VCnt(1, myDimension, 0.); |
| 2299 | |
| 2300 | if(NbConstr > 0) { |
| 2301 | |
| 2302 | while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt && |
| 2303 | IdCnt <= NbConstr) IdCnt++; |
| 2304 | if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) && |
| 2305 | (myTypConstraints->Value(2 * IdCnt) >= 1)) { |
| 2306 | Wcnt = 1.; |
| 2307 | Standard_Integer i0 = 2 * myDimension * (IdCnt - 1), k = 0; |
| 2308 | for( i = 1; i <= myNbP3d; i++) { |
| 2309 | for(Standard_Integer j = 1; j <= 3; j++) |
| 2310 | VCnt(++k) = myTabConstraints->Value(++i0); |
| 2311 | i0 += 3; |
| 2312 | } |
| 2313 | for(i = 1; i <= myNbP2d; i++) { |
| 2314 | for(Standard_Integer j = 1; j <= 2; j++) |
| 2315 | VCnt(++k) = myTabConstraints->Value(++i0); |
| 2316 | i0 += 2; |
| 2317 | } |
| 2318 | } |
| 2319 | } |
| 2320 | |
| 2321 | // Averaging of estimation |
| 2322 | |
| 2323 | Standard_Real Denom = Wpnt + Wcnt; |
| 2324 | if(Denom == 0.) Denom = 1.; |
| 2325 | else Denom = 1. / Denom; |
| 2326 | |
| 2327 | VTang = (Wpnt * VTang + Wcnt * VCnt) * Denom; |
| 2328 | |
| 2329 | Vnorm = VTang.Norm(); |
| 2330 | |
| 2331 | if(Vnorm <= EpsNorm) |
| 2332 | VTang.Init(0.); |
| 2333 | else |
| 2334 | VTang /= Vnorm; |
| 2335 | |
| 2336 | |
| 2337 | } |
| 2338 | |
| 2339 | // |
| 2340 | //======================================================================= |
| 2341 | //function : EstSecnd |
| 2342 | //purpose : Calculation of estimation of the second derivative |
| 2343 | // (see fortran routine MLIMSCN) |
| 2344 | //======================================================================= |
| 2345 | // |
| 2346 | void AppDef_Variational::EstSecnd(const Standard_Integer ipnt, |
| 2347 | const math_Vector& VTang1, |
| 2348 | const math_Vector& VTang2, |
| 2349 | const Standard_Real Length, |
| 2350 | math_Vector& VScnd) const |
| 2351 | { |
| 2352 | Standard_Integer i ; |
| 2353 | |
| 2354 | const Standard_Real Eps = 1.e-9; |
| 2355 | |
| 2356 | Standard_Real Wpnt = 1.; |
| 2357 | |
| 2358 | Standard_Real aux; |
| 2359 | |
| 2360 | if(ipnt == myFirstPoint) |
| 2361 | aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt); |
| 2362 | else if(ipnt == myLastPoint) |
| 2363 | aux = myParameters->Value(ipnt) - myParameters->Value(ipnt - 1); |
| 2364 | else |
| 2365 | aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1); |
| 2366 | |
| 2367 | if(aux <= Eps) |
| 2368 | aux = 1.; |
| 2369 | else |
| 2370 | aux = 1. / aux; |
| 2371 | |
| 2372 | VScnd = (VTang2 - VTang1) * aux; |
| 2373 | |
| 2374 | // Estimation with constraints |
| 2375 | |
| 2376 | Standard_Real Wcnt = 0.; |
| 2377 | Standard_Integer IdCnt = 1; |
| 2378 | |
| 2379 | // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint] |
| 2380 | |
| 2381 | Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
| 2382 | math_Vector VCnt(1, myDimension, 0.); |
| 2383 | |
| 2384 | if(NbConstr > 0) { |
| 2385 | |
| 2386 | while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt && |
| 2387 | IdCnt <= NbConstr) IdCnt++; |
| 2388 | |
| 2389 | if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) && |
| 2390 | (myTypConstraints->Value(2 * IdCnt) >= 2)) { |
| 2391 | Wcnt = 1.; |
| 2392 | Standard_Integer i0 = 2 * myDimension * (IdCnt - 1) + 3, k = 0; |
| 2393 | for( i = 1; i <= myNbP3d; i++) { |
| 2394 | for(Standard_Integer j = 1; j <= 3; j++) |
| 2395 | VCnt(++k) = myTabConstraints->Value(++i0); |
| 2396 | i0 += 3; |
| 2397 | } |
| 2398 | i0--; |
| 2399 | for(i = 1; i <= myNbP2d; i++) { |
| 2400 | for(Standard_Integer j = 1; j <= 2; j++) |
| 2401 | VCnt(++k) = myTabConstraints->Value(++i0); |
| 2402 | i0 += 2; |
| 2403 | } |
| 2404 | } |
| 2405 | } |
| 2406 | |
| 2407 | // Averaging of estimation |
| 2408 | |
| 2409 | Standard_Real Denom = Wpnt + Wcnt; |
| 2410 | if(Denom == 0.) Denom = 1.; |
| 2411 | else Denom = 1. / Denom; |
| 2412 | |
| 2413 | VScnd = (Wpnt * VScnd + (Wcnt * Length) * VCnt) * Denom; |
| 2414 | |
| 2415 | } |
| 2416 | |
| 2417 | |
| 2418 | // |
| 2419 | //======================================================================= |
| 2420 | //function : InitCutting |
| 2421 | //purpose : Realisation of curve's cutting a priory accordingly to |
| 2422 | // constraints (see fortran routine MLICUT) |
| 2423 | //======================================================================= |
| 2424 | // |
| 2425 | void AppDef_Variational::InitCutting(const Handle(PLib_Base)& aBase, |
| 2426 | const Standard_Real CurvTol, |
| 2427 | Handle(FEmTool_Curve)& aCurve) const |
| 2428 | { |
| 2429 | |
| 2430 | // Definition of number of elements |
| 2431 | Standard_Integer ORCMx = -1, NCont = 0, i, kk, NbElem; |
| 2432 | Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
| 2433 | |
| 2434 | for(i = 1; i <= NbConstr; i++) { |
| 2435 | kk = Abs(myTypConstraints->Value(2 * i)) + 1; |
| 2436 | ORCMx = Max(ORCMx, kk); |
| 2437 | NCont += kk; |
| 2438 | } |
| 2439 | |
| 2440 | if(ORCMx > myMaxDegree - myNivCont) |
| 2441 | Standard_ConstructionError::Raise("AppDef_Variational::InitCutting"); |
| 2442 | |
| 2443 | Standard_Integer NLibre = Max(myMaxDegree - myNivCont - (myMaxDegree + 1) / 4, |
| 2444 | myNivCont + 1); |
| 2445 | |
| 2446 | NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1; |
| 2447 | |
| 2448 | while((NbElem > myMaxSegment) && (NLibre < myMaxDegree - myNivCont)) { |
| 2449 | |
| 2450 | NLibre++; |
| 2451 | NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1; |
| 2452 | |
| 2453 | } |
| 2454 | |
| 2455 | |
| 2456 | if(NbElem > myMaxSegment) |
| 2457 | Standard_ConstructionError::Raise("AppDef_Variational::InitCutting"); |
| 2458 | |
| 2459 | |
| 2460 | aCurve = new FEmTool_Curve(myDimension, NbElem, aBase, CurvTol); |
| 2461 | |
| 2462 | Standard_Integer NCnt = (NCont - 1) / NbElem + 1; |
| 2463 | Standard_Integer NPlus = NbElem - (NCnt * NbElem - NCont); |
| 2464 | |
| 2465 | TColStd_Array1OfReal& Knot = aCurve->Knots(); |
| 2466 | |
| 2467 | Standard_Integer IDeb = 0, IFin = NbConstr + 1, |
| 2468 | NDeb = 0, NFin = 0, |
| 2469 | IndEl = Knot.Lower(), IUpper = Knot.Upper(), |
| 2470 | NbEl = 0; |
| 2471 | |
| 2472 | |
| 2473 | Knot(IndEl) = myParameters->Value(myFirstPoint); |
| 2474 | Knot(IUpper) = myParameters->Value(myLastPoint); |
| 2475 | |
| 2476 | while(NbElem - NbEl > 1) { |
| 2477 | |
| 2478 | IndEl++; NbEl++; |
| 2479 | if(NPlus == 0) NCnt--; |
| 2480 | |
| 2481 | while(NDeb < NCnt && IDeb < IFin) { |
| 2482 | IDeb++; |
| 2483 | NDeb += Abs(myTypConstraints->Value(2 * IDeb)) + 1; |
| 2484 | } |
| 2485 | |
| 2486 | if(NDeb == NCnt) { |
| 2487 | NDeb = 0; |
| 2488 | if(NPlus == 1 && |
| 2489 | myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) > Knot(IndEl - 1)) |
| 2490 | |
| 2491 | Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)); |
| 2492 | else |
| 2493 | Knot(IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) + |
| 2494 | myParameters->Value(myTypConstraints->Value(2 * IDeb + 1))) / 2; |
| 2495 | |
| 2496 | } |
| 2497 | else { |
| 2498 | NDeb -= NCnt; |
| 2499 | Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)); |
| 2500 | } |
| 2501 | |
| 2502 | NPlus--; |
| 2503 | if(NPlus == 0) NCnt--; |
| 2504 | |
| 2505 | if(NbElem - NbEl == 1) break; |
| 2506 | |
| 2507 | NbEl++; |
| 2508 | |
| 2509 | while(NFin < NCnt && IDeb < IFin) { |
| 2510 | IFin--; |
| 2511 | NFin += Abs(myTypConstraints->Value(2 * IFin)) + 1; |
| 2512 | } |
| 2513 | |
| 2514 | if(NFin == NCnt) { |
| 2515 | NFin = 0; |
| 2516 | Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) + |
| 2517 | myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2; |
| 2518 | } |
| 2519 | else { |
| 2520 | NFin -= NCnt; |
| 2521 | if(myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) < Knot(IUpper - IndEl + 1)) |
| 2522 | Knot(IUpper + 1 - IndEl) = myParameters->Value(myTypConstraints->Value(2 * IFin - 1)); |
| 2523 | else |
| 2524 | Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) + |
| 2525 | myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2; |
| 2526 | } |
| 2527 | } |
| 2528 | } |
| 2529 | |
| 2530 | //======================================================================= |
| 2531 | //function : Adjusting |
| 2532 | //purpose : Smoothing's adjusting like STRIM routine "MAJLIS" |
| 2533 | //======================================================================= |
| 2534 | void AppDef_Variational::Adjusting( |
| 2535 | Handle(AppDef_SmoothCriterion)& J, |
| 2536 | Standard_Real& WQuadratic, |
| 2537 | Standard_Real& WQuality, |
| 2538 | Handle(FEmTool_Curve)& TheCurve, |
| 2539 | TColStd_Array1OfReal& Ecarts) |
| 2540 | { |
| 2541 | |
| 2542 | // cout << "=========== Adjusting =============" << endl; |
| 2543 | |
| 2544 | /* Initialized data */ |
| 2545 | |
| 2546 | const Standard_Integer mxiter = 2; |
| 2547 | const Standard_Real eps1 = 1e-6; |
| 2548 | Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1; |
| 2549 | Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
| 2550 | Standard_Real CurvTol = eps1 * J->EstLength() / NbrPnt; |
| 2551 | |
| 2552 | |
| 2553 | /* Local variables */ |
| 2554 | Standard_Integer iter, ipnt; |
| 2555 | Standard_Real ecart, emold, erold, tpara; |
| 2556 | Standard_Real vocri[4], j1cibl, vtest, vseuil; |
| 2557 | Standard_Integer i, numint, flag; |
| 2558 | TColStd_Array1OfReal tbpoid(myFirstPoint, myLastPoint); |
| 2559 | Standard_Boolean loptim, lrejet; |
| 2560 | Handle(AppDef_SmoothCriterion) JNew; |
| 2561 | Handle(FEmTool_Curve) CNew; |
| 2562 | Standard_Real E1, E2, E3; |
| 2563 | |
| 2564 | |
| 2565 | /* (0.b) Initialisations */ |
| 2566 | |
| 2567 | loptim = Standard_True; |
| 2568 | iter = 0; |
| 2569 | tbpoid.Init(1.); |
| 2570 | |
| 2571 | |
| 2572 | /* ============ boucle sur le moteur de lissage ============== */ |
| 2573 | |
| 2574 | vtest = WQuality * .9; |
| 2575 | j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint)); |
| 2576 | |
| 2577 | while(loptim) { |
| 2578 | |
| 2579 | ++iter; |
| 2580 | |
| 2581 | /* (1) Sauvegarde de l'etat precedents */ |
| 2582 | |
| 2583 | vocri[0] = myCriterium[0]; |
| 2584 | vocri[1] = myCriterium[1]; |
| 2585 | vocri[2] = myCriterium[2]; |
| 2586 | vocri[3] = myCriterium[3]; |
| 2587 | erold = myMaxError; |
| 2588 | emold = myAverageError; |
| 2589 | |
| 2590 | /* (2) Augmentation du poids des moindre carre */ |
| 2591 | |
| 2592 | if (j1cibl > vtest) { |
| 2593 | WQuadratic = j1cibl / vtest * WQuadratic; |
| 2594 | } |
| 2595 | |
| 2596 | /* (3) Augmentation du poid associe aux points a problemes */ |
| 2597 | |
| 2598 | vseuil = WQuality * .88; |
| 2599 | |
| 2600 | for (ipnt = myFirstPoint; ipnt <= myLastPoint; ++ipnt) { |
| 2601 | if (Ecarts(ipnt) > vtest) { |
| 2602 | ecart = (Ecarts(ipnt) - vseuil) / WQuality; |
| 2603 | tbpoid(ipnt) = (ecart * 3 + 1.) * tbpoid(ipnt); |
| 2604 | } |
| 2605 | } |
| 2606 | |
| 2607 | /* (4) Decoupe force */ |
| 2608 | |
| 2609 | if (TheCurve->NbElements() < myMaxSegment && myWithCutting) { |
| 2610 | |
| 2611 | numint = NearIndex(myParameters->Value(myMaxErrorIndex), TheCurve->Knots(), 0, flag); |
| 2612 | |
| 2613 | tpara = (TheCurve->Knots()(numint) + TheCurve->Knots()(numint + 1) + |
| 2614 | myParameters->Value(myMaxErrorIndex) * 2) / 4; |
| 2615 | |
| 2616 | CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements() + 1, |
| 2617 | TheCurve->Base(), CurvTol); |
| 2618 | |
| 2619 | for(i = 1; i <= numint; i++) CNew->Knots()(i) = TheCurve->Knots()(i); |
| 2620 | for(i = numint + 1; i <= TheCurve->Knots().Length(); i++) |
| 2621 | CNew->Knots()(i + 1) = TheCurve->Knots()(i); |
| 2622 | |
| 2623 | CNew->Knots()(numint + 1) = tpara; |
| 2624 | |
| 2625 | } else { |
| 2626 | |
| 2627 | CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements(), TheCurve->Base(), CurvTol); |
| 2628 | |
| 2629 | CNew->Knots() = TheCurve->Knots(); |
| 2630 | } |
| 2631 | |
| 2632 | |
| 2633 | JNew = new AppDef_LinearCriteria(mySSP, myFirstPoint, myLastPoint); |
| 2634 | |
| 2635 | JNew->EstLength() = J->EstLength(); |
| 2636 | |
| 2637 | J->GetEstimation(E1, E2, E3); |
| 2638 | |
| 2639 | JNew->SetEstimation(E1, E2, E3); |
| 2640 | |
| 2641 | JNew->SetParameters(myParameters); |
| 2642 | |
| 2643 | JNew->SetWeight(WQuadratic, WQuality, myPercent[0], myPercent[1], myPercent[2]); |
| 2644 | |
| 2645 | JNew->SetWeight(tbpoid); |
| 2646 | |
| 2647 | JNew->SetCurve(CNew); |
| 2648 | |
| 2649 | /* (5) Relissage */ |
| 2650 | |
| 2651 | TheMotor(JNew, WQuadratic, WQuality, CNew, Ecarts); |
| 2652 | |
| 2653 | /* (6) Tests de rejet */ |
| 2654 | |
| 2655 | j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint)); |
| 2656 | vseuil = Sqrt(vocri[1]) + (erold - myMaxError) * 4; |
| 2657 | |
| 2658 | lrejet = ((myMaxError > WQuality) && (myMaxError > erold * 1.01)) || (Sqrt(myCriterium[1]) > vseuil * 1.05); |
| 2659 | |
| 2660 | if (lrejet) { |
| 2661 | myCriterium[0] = vocri[0]; |
| 2662 | myCriterium[1] = vocri[1]; |
| 2663 | myCriterium[2] = vocri[2]; |
| 2664 | myCriterium[3] = vocri[3]; |
| 2665 | myMaxError = erold; |
| 2666 | myAverageError = emold; |
| 2667 | |
| 2668 | loptim = Standard_False; |
| 2669 | } |
| 2670 | else { |
| 2671 | J = JNew; |
| 2672 | TheCurve = CNew; |
| 2673 | J->SetCurve(TheCurve); |
| 2674 | } |
| 2675 | |
| 2676 | /* (7) Test de convergence */ |
| 2677 | |
| 2678 | if (((iter >= mxiter) && (myMaxSegment == CNew->NbElements())) || myMaxError < WQuality) { |
| 2679 | loptim = Standard_False; |
| 2680 | } |
| 2681 | } |
| 2682 | } |
| 2683 | |
| 2684 | static Standard_Boolean NotParallel(gp_Vec& T, gp_Vec& V) |
| 2685 | { |
| 2686 | V = T; |
| 2687 | V.SetX(V.X() + 1.); |
| 2688 | if (V.CrossMagnitude(T) > 1.e-12) |
| 2689 | return Standard_True; |
| 2690 | V.SetY(V.Y() + 1.); |
| 2691 | if (V.CrossMagnitude(T) > 1.e-12) |
| 2692 | return Standard_True; |
| 2693 | V.SetZ(V.Z() + 1.); |
| 2694 | if (V.CrossMagnitude(T) > 1.e-12) |
| 2695 | return Standard_True; |
| 2696 | return Standard_False; |
| 2697 | } |
| 2698 | |
| 2699 | void AppDef_Variational::AssemblingConstraints(const Handle(FEmTool_Curve)& Curve, |
| 2700 | const TColStd_Array1OfReal& Parameters, |
| 2701 | const Standard_Real CBLONG, |
| 2702 | FEmTool_Assembly& A ) const |
| 2703 | { |
| 2704 | |
| 2705 | Standard_Integer MxDeg = Curve->Base()->WorkDegree(), |
| 2706 | NbElm = Curve->NbElements(), |
| 2707 | NbDim = Curve->Dimension(); |
| 2708 | |
| 2709 | TColStd_Array1OfReal G0(0, MxDeg), G1(0, MxDeg), G2(0, MxDeg); |
| 2710 | math_Vector V0((Standard_Real*)&G0(0), 0, MxDeg), |
| 2711 | V1((Standard_Real*)&G1(0), 0, MxDeg), |
| 2712 | V2((Standard_Real*)&G2(0), 0, MxDeg); |
| 2713 | |
| 2714 | Standard_Integer IndexOfConstraint, Ng3d, Ng2d, NBeg2d, NPass, NgPC1, |
| 2715 | NTang3d, NTang2d, |
| 2716 | Point, TypOfConstr, |
| 2717 | p0 = Parameters.Lower() - myFirstPoint, |
| 2718 | curel = 1, el, i, ipnt, ityp, j, k, pnt, curdim, |
| 2719 | jt, Ntheta = 6 * myNbP3d + 2 * myNbP2d; |
| 2720 | Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
| 2721 | |
| 2722 | // Ng3d = 3 * NbConstr + 2 * myNbTangPoints + 5 * myNbCurvPoints; |
| 2723 | // Ng2d = 2 * NbConstr + 1 * myNbTangPoints + 3 * myNbCurvPoints; |
| 2724 | Ng3d = 3 * NbConstr + 3 * myNbTangPoints + 5 * myNbCurvPoints; |
| 2725 | Ng2d = 2 * NbConstr + 2 * myNbTangPoints + 3 * myNbCurvPoints; |
| 2726 | NBeg2d = Ng3d * myNbP3d; |
| 2727 | // NgPC1 = NbConstr + myNbCurvPoints; |
| 2728 | NgPC1 = NbConstr + myNbTangPoints + myNbCurvPoints; |
| 2729 | NPass = 0; |
| 2730 | NTang3d = 3 * NgPC1; |
| 2731 | NTang2d = 2 * NgPC1; |
| 2732 | |
| 2733 | TColStd_Array1OfReal& Intervals = Curve->Knots(); |
| 2734 | |
| 2735 | Standard_Real t, R1, R2; |
| 2736 | |
| 2737 | Handle(PLib_Base) myBase = Curve->Base(); |
| 2738 | Handle(PLib_HermitJacobi) myHermitJacobi = Handle(PLib_HermitJacobi)::DownCast (myBase); |
| 2739 | Standard_Integer Order = myHermitJacobi->NivConstr() + 1; |
| 2740 | |
| 2741 | Standard_Real UFirst, ULast, coeff, c0, mfact, mfact1; |
| 2742 | |
| 2743 | A.NullifyConstraint(); |
| 2744 | |
| 2745 | ipnt = -1; |
| 2746 | ityp = 0; |
| 2747 | for(i = 1; i <= NbConstr; i++) { |
| 2748 | |
| 2749 | ipnt += 2; ityp += 2; |
| 2750 | |
| 2751 | Point = myTypConstraints->Value(ipnt); |
| 2752 | TypOfConstr = myTypConstraints->Value(ityp); |
| 2753 | |
| 2754 | t = Parameters(p0 + Point); |
| 2755 | |
| 2756 | for(el = curel; el <= NbElm; ) { |
| 2757 | if( t <= Intervals(++el) ) { |
| 2758 | curel = el - 1; |
| 2759 | break; |
| 2760 | } |
| 2761 | } |
| 2762 | |
| 2763 | |
| 2764 | UFirst = Intervals(curel); ULast = Intervals(curel + 1); |
| 2765 | coeff = (ULast - UFirst)/2.; c0 = (ULast + UFirst)/2.; |
| 2766 | |
| 2767 | t = (t - c0) / coeff; |
| 2768 | |
| 2769 | if(TypOfConstr == 0) { |
| 2770 | myBase->D0(t, G0); |
| 2771 | for(k = 1; k < Order; k++) { |
| 2772 | mfact = Pow(coeff, k); |
| 2773 | G0(k) *= mfact; |
| 2774 | G0(k + Order) *= mfact; |
| 2775 | } |
| 2776 | } |
| 2777 | else if(TypOfConstr == 1) { |
| 2778 | myBase->D1(t, G0, G1); |
| 2779 | for(k = 1; k < Order; k++) { |
| 2780 | mfact = Pow(coeff, k); |
| 2781 | G0(k) *= mfact; |
| 2782 | G0(k + Order) *= mfact; |
| 2783 | G1(k) *= mfact; |
| 2784 | G1(k + Order) *= mfact; |
| 2785 | } |
| 2786 | mfact = 1./coeff; |
| 2787 | for(k = 0; k <= MxDeg; k++) { |
| 2788 | G1(k) *= mfact; |
| 2789 | } |
| 2790 | } |
| 2791 | else { |
| 2792 | myBase->D2(t, G0, G1, G2); |
| 2793 | for(k = 1; k < Order; k++) { |
| 2794 | mfact = Pow(coeff, k); |
| 2795 | G0(k) *= mfact; |
| 2796 | G0(k + Order) *= mfact; |
| 2797 | G1(k) *= mfact; |
| 2798 | G1(k + Order) *= mfact; |
| 2799 | G2(k) *= mfact; |
| 2800 | G2(k + Order) *= mfact; |
| 2801 | } |
| 2802 | mfact = 1. / coeff; |
| 2803 | mfact1 = mfact / coeff; |
| 2804 | for(k = 0; k <= MxDeg; k++) { |
| 2805 | G1(k) *= mfact; |
| 2806 | G2(k) *= mfact1; |
| 2807 | } |
| 2808 | } |
| 2809 | |
| 2810 | NPass++; |
| 2811 | |
| 2812 | j = NbDim * (Point - myFirstPoint); |
| 2813 | Standard_Integer n0 = NPass; |
| 2814 | curdim = 0; |
| 2815 | for(pnt = 1; pnt <= myNbP3d; pnt++) { |
| 2816 | IndexOfConstraint = n0; |
| 2817 | for(k = 1; k <= 3; k++) { |
| 2818 | curdim++; |
| 2819 | A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k)); |
| 2820 | IndexOfConstraint += NgPC1; |
| 2821 | } |
| 2822 | j +=3; |
| 2823 | n0 += Ng3d; |
| 2824 | } |
| 2825 | |
| 2826 | n0 = NPass + NBeg2d; |
| 2827 | for(pnt = 1; pnt <= myNbP2d; pnt++) { |
| 2828 | IndexOfConstraint = n0; |
| 2829 | for(k = 1; k <= 2; k++) { |
| 2830 | curdim++; |
| 2831 | A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k)); |
| 2832 | IndexOfConstraint += NgPC1; |
| 2833 | } |
| 2834 | j +=2; |
| 2835 | n0 += Ng2d; |
| 2836 | } |
| 2837 | |
| 2838 | /* if(TypOfConstr == 1) { |
| 2839 | |
| 2840 | IndexOfConstraint = NTang3d + 1; |
| 2841 | jt = Ntheta * (i - 1); |
| 2842 | for(pnt = 1; pnt <= myNbP3d; pnt++) { |
| 2843 | for(k = 1; k <= 3; k++) { |
| 2844 | A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.); |
| 2845 | A.AddConstraint(IndexOfConstraint + 1, curel, k, myTtheta->Value(jt + 3 + k) * V1, 0.); |
| 2846 | } |
| 2847 | IndexOfConstraint += Ng3d; |
| 2848 | jt += 6; |
| 2849 | } |
| 2850 | |
| 2851 | IndexOfConstraint = NBeg2d + NTang2d + 1; |
| 2852 | for(pnt = 1; pnt <= myNbP2d; pnt++) { |
| 2853 | for(k = 1; k <= 2; k++) { |
| 2854 | A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.); |
| 2855 | } |
| 2856 | IndexOfConstraint += Ng2d; |
| 2857 | jt += 2; |
| 2858 | } |
| 2859 | |
| 2860 | NTang3d += 2; |
| 2861 | NTang2d += 1; |
| 2862 | } */ |
| 2863 | if(TypOfConstr == 1) { |
| 2864 | |
| 2865 | NPass++; |
| 2866 | n0 = NPass; |
| 2867 | j = 2 * NbDim * (i - 1); |
| 2868 | curdim = 0; |
| 2869 | for(pnt = 1; pnt <= myNbP3d; pnt++) { |
| 2870 | IndexOfConstraint = n0; |
| 2871 | for(k = 1; k <= 3; k++) { |
| 2872 | curdim++; |
| 2873 | A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k)); |
| 2874 | IndexOfConstraint += NgPC1; |
| 2875 | } |
| 2876 | n0 += Ng3d; |
| 2877 | j += 6; |
| 2878 | } |
| 2879 | |
| 2880 | n0 = NPass + NBeg2d; |
| 2881 | for(pnt = 1; pnt <= myNbP2d; pnt++) { |
| 2882 | IndexOfConstraint = n0; |
| 2883 | for(k = 1; k <= 2; k++) { |
| 2884 | curdim++; |
| 2885 | A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k)); |
| 2886 | IndexOfConstraint += NgPC1; |
| 2887 | } |
| 2888 | n0 += Ng2d; |
| 2889 | j += 4; |
| 2890 | } |
| 2891 | } |
| 2892 | if(TypOfConstr == 2) { |
| 2893 | |
| 2894 | NPass++; |
| 2895 | n0 = NPass; |
| 2896 | j = 2 * NbDim * (i - 1); |
| 2897 | curdim = 0; |
| 2898 | for(pnt = 1; pnt <= myNbP3d; pnt++) { |
| 2899 | IndexOfConstraint = n0; |
| 2900 | for(k = 1; k <= 3; k++) { |
| 2901 | curdim++; |
| 2902 | A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k)); |
| 2903 | IndexOfConstraint += NgPC1; |
| 2904 | } |
| 2905 | n0 += Ng3d; |
| 2906 | j += 6; |
| 2907 | } |
| 2908 | |
| 2909 | n0 = NPass + NBeg2d; |
| 2910 | for(pnt = 1; pnt <= myNbP2d; pnt++) { |
| 2911 | IndexOfConstraint = n0; |
| 2912 | for(k = 1; k <= 2; k++) { |
| 2913 | curdim++; |
| 2914 | A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k)); |
| 2915 | IndexOfConstraint += NgPC1; |
| 2916 | } |
| 2917 | n0 += Ng2d; |
| 2918 | j += 4; |
| 2919 | } |
| 2920 | |
| 2921 | j = 2 * NbDim * (i - 1) + 3; |
| 2922 | jt = Ntheta * (i - 1); |
| 2923 | IndexOfConstraint = NTang3d + 1; |
| 2924 | curdim = 0; |
| 2925 | for(pnt = 1; pnt <= myNbP3d; pnt++) { |
| 2926 | R1 = 0.; R2 = 0.; |
| 2927 | for(k = 1; k <= 3; k++) { |
| 2928 | R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k); |
| 2929 | R2 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + 3 + k); |
| 2930 | } |
| 2931 | R1 *= CBLONG * CBLONG; |
| 2932 | R2 *= CBLONG * CBLONG; |
| 2933 | for(k = 1; k <= 3; k++) { |
| 2934 | curdim++; |
| 2935 | if(k > 1) R1 = R2 = 0.; |
| 2936 | A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1); |
| 2937 | A.AddConstraint(IndexOfConstraint + 1, curel, curdim, myTfthet->Value(jt + 3 + k) * V2, R2); |
| 2938 | } |
| 2939 | IndexOfConstraint += Ng3d; |
| 2940 | j += 6; |
| 2941 | jt += 6; |
| 2942 | } |
| 2943 | |
| 2944 | j--; |
| 2945 | IndexOfConstraint = NBeg2d + NTang2d + 1; |
| 2946 | for(pnt = 1; pnt <= myNbP2d; pnt++) { |
| 2947 | R1 = 0.; |
| 2948 | for(k = 1; k <= 2; k++) { |
| 2949 | R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k); |
| 2950 | } |
| 2951 | R1 *= CBLONG * CBLONG; |
| 2952 | for(k = 1; k <= 2; k++) { |
| 2953 | curdim++; |
| 2954 | if(k > 1) R1 = 0.; |
| 2955 | A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1); |
| 2956 | } |
| 2957 | IndexOfConstraint += Ng2d; |
| 2958 | j += 4; |
| 2959 | jt += 2; |
| 2960 | } |
| 2961 | |
| 2962 | NTang3d += 2; |
| 2963 | NTang2d += 1; |
| 2964 | } |
| 2965 | |
| 2966 | } |
| 2967 | |
| 2968 | } |
| 2969 | |
| 2970 | Standard_Boolean AppDef_Variational::InitTthetaF(const Standard_Integer ndimen, |
| 2971 | const AppParCurves_Constraint typcon, |
| 2972 | const Standard_Integer begin, |
| 2973 | const Standard_Integer jndex) |
| 2974 | { |
| 2975 | if ((ndimen < 2)||(ndimen >3)) |
| 2976 | return Standard_False; |
| 2977 | gp_Vec T, V; |
| 2978 | gp_Vec theta1, theta2; |
| 2979 | gp_Vec F; |
| 2980 | Standard_Real XX, XY, YY, XZ, YZ, ZZ; |
| 2981 | |
| 2982 | if ((typcon == AppParCurves_TangencyPoint)||(typcon == AppParCurves_CurvaturePoint)) |
| 2983 | { |
| 2984 | T.SetX(myTabConstraints->Value(jndex)); |
| 2985 | T.SetY(myTabConstraints->Value(jndex + 1)); |
| 2986 | if (ndimen == 3) |
| 2987 | T.SetZ(myTabConstraints->Value(jndex + 2)); |
| 2988 | else |
| 2989 | T.SetZ(0.); |
| 2990 | if (ndimen == 2) |
| 2991 | { |
| 2992 | V.SetX(0.); |
| 2993 | V.SetY(0.); |
| 2994 | V.SetZ(1.); |
| 2995 | } |
| 2996 | if (ndimen == 3) |
| 2997 | if (!NotParallel(T, V)) |
| 2998 | return Standard_False; |
| 2999 | theta1 = V ^ T; |
| 3000 | theta1.Normalize(); |
| 3001 | myTtheta->SetValue(begin, theta1.X()); |
| 3002 | myTtheta->SetValue(begin + 1, theta1.Y()); |
| 3003 | if (ndimen == 3) |
| 3004 | { |
| 3005 | theta2 = T ^ theta1; |
| 3006 | theta2.Normalize(); |
| 3007 | myTtheta->SetValue(begin + 2, theta1.Z()); |
| 3008 | myTtheta->SetValue(begin + 3, theta2.X()); |
| 3009 | myTtheta->SetValue(begin + 4, theta2.Y()); |
| 3010 | myTtheta->SetValue(begin + 5, theta2.Z()); |
| 3011 | } |
| 3012 | |
| 3013 | // Calculation of myTfthet |
| 3014 | if (typcon == AppParCurves_CurvaturePoint) |
| 3015 | { |
| 3016 | XX = Pow(T.X(), 2); |
| 3017 | XY = T.X() * T.Y(); |
| 3018 | YY = Pow(T.Y(), 2); |
| 3019 | if (ndimen == 2) |
| 3020 | { |
| 3021 | F.SetX(YY * theta1.X() - XY * theta1.Y()); |
| 3022 | F.SetY(XX * theta1.Y() - XY * theta1.X()); |
| 3023 | myTfthet->SetValue(begin, F.X()); |
| 3024 | myTfthet->SetValue(begin + 1, F.Y()); |
| 3025 | } |
| 3026 | if (ndimen == 3) |
| 3027 | { |
| 3028 | XZ = T.X() * T.Z(); |
| 3029 | YZ = T.Y() * T.Z(); |
| 3030 | ZZ = Pow(T.Z(), 2); |
| 3031 | |
| 3032 | F.SetX((ZZ + YY) * theta1.X() - XY * theta1.Y() - XZ * theta1.Z()); |
| 3033 | F.SetY((XX + ZZ) * theta1.Y() - XY * theta1.X() - YZ * theta1.Z()); |
| 3034 | F.SetZ((XX + YY) * theta1.Z() - XZ * theta1.X() - YZ * theta1.Y()); |
| 3035 | myTfthet->SetValue(begin, F.X()); |
| 3036 | myTfthet->SetValue(begin + 1, F.Y()); |
| 3037 | myTfthet->SetValue(begin + 2, F.Z()); |
| 3038 | F.SetX((ZZ + YY) * theta2.X() - XY * theta2.Y() - XZ * theta2.Z()); |
| 3039 | F.SetY((XX + ZZ) * theta2.Y() - XY * theta2.X() - YZ * theta2.Z()); |
| 3040 | F.SetZ((XX + YY) * theta2.Z() - XZ * theta2.X() - YZ * theta2.Y()); |
| 3041 | myTfthet->SetValue(begin + 3, F.X()); |
| 3042 | myTfthet->SetValue(begin + 4, F.Y()); |
| 3043 | myTfthet->SetValue(begin + 5, F.Z()); |
| 3044 | } |
| 3045 | } |
| 3046 | } |
| 3047 | return Standard_True; |
| 3048 | } |
| 3049 | |