1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 // File : SMESH_Pattern.hxx
24 // Created : Mon Aug 2 10:30:00 2004
25 // Author : Edward AGAPOV (eap)
27 #include "SMESH_Block.hxx"
29 #include "SMDS_MeshNode.hxx"
30 #include "SMDS_MeshVolume.hxx"
31 #include "SMDS_VolumeTool.hxx"
32 #include "SMESH_MeshAlgos.hxx"
34 #include <BRepAdaptor_Curve.hxx>
35 #include <BRepAdaptor_Curve2d.hxx>
36 #include <BRepAdaptor_Surface.hxx>
37 #include <BRepTools.hxx>
38 #include <BRepTools_WireExplorer.hxx>
39 #include <BRep_Builder.hxx>
40 #include <BRep_Tool.hxx>
41 #include <Bnd_B2d.hxx>
42 #include <Bnd_Box.hxx>
43 #include <Extrema_ExtPC.hxx>
44 #include <Extrema_ExtPS.hxx>
45 #include <Extrema_POnCurv.hxx>
46 #include <Extrema_POnSurf.hxx>
47 #include <Geom2d_Curve.hxx>
48 #include <ShapeAnalysis.hxx>
49 #include <TopExp_Explorer.hxx>
50 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
51 #include <TopTools_ListIteratorOfListOfShape.hxx>
52 #include <TopTools_ListOfShape.hxx>
54 #include <TopoDS_Compound.hxx>
55 #include <TopoDS_Face.hxx>
56 #include <TopoDS_Iterator.hxx>
57 #include <TopoDS_Wire.hxx>
58 #include <gp_Trsf.hxx>
60 #include <math_FunctionSetRoot.hxx>
61 #include <math_Matrix.hxx>
62 #include <math_Vector.hxx>
64 #include <utilities.h>
71 //#define DEBUG_PARAM_COMPUTE
73 //================================================================================
75 * \brief Set edge data
76 * \param edgeID - block sub-shape ID
77 * \param curve - edge geometry
78 * \param isForward - is curve orientation coincides with edge orientation in the block
80 //================================================================================
82 void SMESH_Block::TEdge::Set( const int edgeID, Adaptor3d_Curve* curve, const bool isForward )
84 myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID );
85 if ( myC3d ) delete myC3d;
87 myFirst = curve->FirstParameter();
88 myLast = curve->LastParameter();
90 std::swap( myFirst, myLast );
93 //================================================================================
95 * \brief Set coordinates of nodes at edge ends to work with mesh block
96 * \param edgeID - block sub-shape ID
97 * \param node1 - coordinates of node with lower ID
98 * \param node2 - coordinates of node with upper ID
100 //================================================================================
102 void SMESH_Block::TEdge::Set( const int edgeID, const gp_XYZ& node1, const gp_XYZ& node2 )
104 myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID );
105 myNodes[ 0 ] = node1;
106 myNodes[ 1 ] = node2;
108 if ( myC3d ) delete myC3d;
112 //=======================================================================
113 //function : SMESH_Block::TEdge::GetU
115 //=======================================================================
117 double SMESH_Block::TEdge::GetU( const gp_XYZ& theParams ) const
119 double u = theParams.Coord( myCoordInd );
120 if ( !myC3d ) // if mesh block
122 return ( 1 - u ) * myFirst + u * myLast;
125 //=======================================================================
126 //function : SMESH_Block::TEdge::Point
128 //=======================================================================
130 gp_XYZ SMESH_Block::TEdge::Point( const gp_XYZ& theParams ) const
132 double u = GetU( theParams );
133 if ( myC3d ) return myC3d->Value( u ).XYZ();
135 return myNodes[0] * ( 1 - u ) + myNodes[1] * u;
138 //================================================================================
142 //================================================================================
144 SMESH_Block::TEdge::~TEdge()
146 if ( myC3d ) delete myC3d;
149 //================================================================================
151 * \brief Set face data
152 * \param faceID - block sub-shape ID
153 * \param S - face surface geometry
154 * \param c2d - 4 pcurves in the order as returned by GetFaceEdgesIDs(faceID)
155 * \param isForward - orientation of pcurves comparing with block edge direction
157 //================================================================================
159 void SMESH_Block::TFace::Set( const int faceID,
160 Adaptor3d_Surface* S,
161 Adaptor2d_Curve2d* c2D[4],
162 const bool isForward[4] )
164 if ( myS ) delete myS;
167 vector< int > edgeIdVec;
168 GetFaceEdgesIDs( faceID, edgeIdVec );
169 for ( size_t iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges
171 myCoordInd[ iE ] = GetCoordIndOnEdge( edgeIdVec[ iE ] );
172 if ( myC2d[ iE ]) delete myC2d[ iE ];
173 myC2d[ iE ] = c2D[ iE ];
174 myFirst[ iE ] = myC2d[ iE ]->FirstParameter();
175 myLast [ iE ] = myC2d[ iE ]->LastParameter();
176 if ( !isForward[ iE ])
177 std::swap( myFirst[ iE ], myLast[ iE ] );
180 myCorner[ 0 ] = myC2d[ 0 ]->Value( myFirst[0] ).XY();
181 myCorner[ 1 ] = myC2d[ 0 ]->Value( myLast[0] ).XY();
182 myCorner[ 2 ] = myC2d[ 1 ]->Value( myLast[1] ).XY();
183 myCorner[ 3 ] = myC2d[ 1 ]->Value( myFirst[1] ).XY();
186 //================================================================================
188 * \brief Set face data to work with mesh block
189 * \param faceID - block sub-shape ID
190 * \param edgeU0 - filled data of edge u0 = GetFaceEdgesIDs(faceID)[ 0 ]
191 * \param edgeU1 - filled data of edge u1 = GetFaceEdgesIDs(faceID)[ 1 ]
193 //================================================================================
195 void SMESH_Block::TFace::Set( const int faceID, const TEdge& edgeU0, const TEdge& edgeU1 )
197 vector< int > edgeIdVec;
198 GetFaceEdgesIDs( faceID, edgeIdVec );
199 myNodes[ 0 ] = edgeU0.NodeXYZ( 1 );
200 myNodes[ 1 ] = edgeU0.NodeXYZ( 0 );
201 myNodes[ 2 ] = edgeU1.NodeXYZ( 0 );
202 myNodes[ 3 ] = edgeU1.NodeXYZ( 1 );
203 myCoordInd[ 0 ] = GetCoordIndOnEdge( edgeIdVec[ 0 ] );
204 myCoordInd[ 1 ] = GetCoordIndOnEdge( edgeIdVec[ 1 ] );
205 myCoordInd[ 2 ] = GetCoordIndOnEdge( edgeIdVec[ 2 ] );
206 myCoordInd[ 3 ] = GetCoordIndOnEdge( edgeIdVec[ 3 ] );
207 if ( myS ) delete myS;
211 //================================================================================
215 //================================================================================
217 SMESH_Block::TFace::~TFace()
219 if ( myS ) delete myS;
220 for ( int i = 0 ; i < 4; ++i )
221 if ( myC2d[ i ]) delete myC2d[ i ];
224 //=======================================================================
225 //function : SMESH_Block::TFace::GetCoefs
226 //purpose : return coefficients for addition of [0-3]-th edge and vertex
227 //=======================================================================
229 void SMESH_Block::TFace::GetCoefs(int iE,
230 const gp_XYZ& theParams,
232 double& Vcoef ) const
234 double dU = theParams.Coord( GetUInd() );
235 double dV = theParams.Coord( GetVInd() );
238 Ecoef = ( 1 - dV ); // u0
239 Vcoef = ( 1 - dU ) * ( 1 - dV ); break; // 00
242 Vcoef = dU * ( 1 - dV ); break; // 10
244 Ecoef = ( 1 - dU ); // 0v
245 Vcoef = dU * dV ; break; // 11
248 Vcoef = ( 1 - dU ) * dV ; break; // 01
253 //=======================================================================
254 //function : SMESH_Block::TFace::GetUV
256 //=======================================================================
258 gp_XY SMESH_Block::TFace::GetUV( const gp_XYZ& theParams ) const
261 for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
263 double Ecoef = 0, Vcoef = 0;
264 GetCoefs( iE, theParams, Ecoef, Vcoef );
266 double u = theParams.Coord( myCoordInd[ iE ] );
267 u = ( 1 - u ) * myFirst[ iE ] + u * myLast[ iE ];
268 uv += Ecoef * myC2d[ iE ]->Value( u ).XY();
270 uv -= Vcoef * myCorner[ iE ];
275 //=======================================================================
276 //function : SMESH_Block::TFace::Point
278 //=======================================================================
280 gp_XYZ SMESH_Block::TFace::Point( const gp_XYZ& theParams ) const
283 if ( !myS ) // if mesh block
285 for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
287 double Ecoef = 0, Vcoef = 0;
288 GetCoefs( iE, theParams, Ecoef, Vcoef );
290 double u = theParams.Coord( myCoordInd[ iE ] );
293 case 1: i1 = 3; i2 = 2; break;
294 case 2: i1 = 1; i2 = 2; break;
295 case 3: i1 = 0; i2 = 3; break;
297 p += Ecoef * ( myNodes[ i1 ] * ( 1 - u ) + myNodes[ i2 ] * u );
299 p -= Vcoef * myNodes[ iE ];
305 gp_XY uv = GetUV( theParams );
306 p = myS->Value( uv.X(), uv.Y() ).XYZ();
315 bool isPntInTria( const gp_XY& p, const gp_XY& t0, const gp_XY& t1, const gp_XY& t2 )
318 SMESH_MeshAlgos::GetBarycentricCoords( p, t0, t1, t2, bc0, bc1 );
319 return ( bc0 >= 0. && bc1 >= 0. && bc0 + bc1 <= 1. );
323 bool isPntInQuad( const gp_XY& p,
324 const gp_XY& q0, const gp_XY& q1, const gp_XY& q2, const gp_XY& q3 )
326 const int in1 = isPntInTria( p, q0, q1, q2 );
327 const int in2 = isPntInTria( p, q0, q2, q3 );
328 return in1 + in2 == 1;
332 //=======================================================================
333 //function : IsUVInQuad
334 //purpose : Checks if UV is in a quardilateral defined by 4 nornalized points
335 //=======================================================================
337 bool SMESH_Block::TFace::IsUVInQuad( const gp_XY& uv,
338 const gp_XYZ& param0, const gp_XYZ& param1,
339 const gp_XYZ& param2, const gp_XYZ& param3 ) const
341 gp_XY q0 = GetUV( param0 );
342 gp_XY q1 = GetUV( param1 );
343 gp_XY q2 = GetUV( param2 );
344 gp_XY q3 = GetUV( param3 );
345 return isPntInQuad( uv, q0,q1,q2,q3);
348 //=======================================================================
349 //function : GetUVRange
350 //purpose : returns UV range of the face
351 //=======================================================================
353 gp_XY SMESH_Block::TFace::GetUVRange() const
355 if ( !myS ) return gp_XY(1.,1.);
358 for ( int iE = 0; iE < 4; ++iE )
360 //TColStd_Array1OfReal T(1,
362 return bb.CornerMax() - bb.CornerMin();
365 //=======================================================================
366 //function : GetShapeCoef
368 //=======================================================================
370 double* SMESH_Block::GetShapeCoef (const int theShapeID)
372 static double shapeCoef[][3] = {
373 // V000, V100, V010, V110
374 { -1,-1,-1 }, { 1,-1,-1 }, { -1, 1,-1 }, { 1, 1,-1 },
375 // V001, V101, V011, V111,
376 { -1,-1, 1 }, { 1,-1, 1 }, { -1, 1, 1 }, { 1, 1, 1 },
377 // Ex00, Ex10, Ex01, Ex11,
378 { 0,-1,-1 }, { 0, 1,-1 }, { 0,-1, 1 }, { 0, 1, 1 },
379 // E0y0, E1y0, E0y1, E1y1,
380 { -1, 0,-1 }, { 1, 0,-1 }, { -1, 0, 1 }, { 1, 0, 1 },
381 // E00z, E10z, E01z, E11z,
382 { -1,-1, 0 }, { 1,-1, 0 }, { -1, 1, 0 }, { 1, 1, 0 },
383 // Fxy0, Fxy1, Fx0z, Fx1z, F0yz, F1yz,
384 { 0, 0,-1 }, { 0, 0, 1 }, { 0,-1, 0 }, { 0, 1, 0 }, { -1, 0, 0 }, { 1, 0, 0 },
388 if ( theShapeID < ID_V000 || theShapeID > ID_F1yz )
389 return shapeCoef[ ID_Shell - 1 ];
391 return shapeCoef[ theShapeID - 1 ];
394 //=======================================================================
395 //function : ShellPoint
396 //purpose : return coordinates of a point in shell
397 //=======================================================================
399 bool SMESH_Block::ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const
401 thePoint.SetCoord( 0., 0., 0. );
402 for ( int shapeID = ID_V000; shapeID < ID_Shell; shapeID++ )
405 double* coefs = GetShapeCoef( shapeID );
407 for ( int iCoef = 0; iCoef < 3; iCoef++ ) {
408 if ( coefs[ iCoef ] != 0 ) {
409 if ( coefs[ iCoef ] < 0 )
410 k *= ( 1. - theParams.Coord( iCoef + 1 ));
412 k *= theParams.Coord( iCoef + 1 );
415 // add point on a shape
416 if ( fabs( k ) > DBL_MIN )
419 if ( shapeID < ID_Ex00 ) // vertex
420 VertexPoint( shapeID, Ps );
421 else if ( shapeID < ID_Fxy0 ) { // edge
422 EdgePoint( shapeID, theParams, Ps );
425 FacePoint( shapeID, theParams, Ps );
433 //=======================================================================
434 //function : ShellPoint
435 //purpose : computes coordinates of a point in shell by points on sub-shapes;
436 // thePointOnShape[ subShapeID ] must be a point on a sub-shape
437 //=======================================================================
439 bool SMESH_Block::ShellPoint(const gp_XYZ& theParams,
440 const vector<gp_XYZ>& thePointOnShape,
443 if ( thePointOnShape.size() < ID_F1yz )
446 const double x = theParams.X(), y = theParams.Y(), z = theParams.Z();
447 const double x1 = 1. - x, y1 = 1. - y, z1 = 1. - z;
448 const vector<gp_XYZ>& p = thePointOnShape;
451 x1 * p[ID_F0yz] + x * p[ID_F1yz] +
452 y1 * p[ID_Fx0z] + y * p[ID_Fx1z] +
453 z1 * p[ID_Fxy0] + z * p[ID_Fxy1] +
454 x1 * (y1 * (z1 * p[ID_V000] + z * p[ID_V001]) +
455 y * (z1 * p[ID_V010] + z * p[ID_V011])) +
456 x * (y1 * (z1 * p[ID_V100] + z * p[ID_V101]) +
457 y * (z1 * p[ID_V110] + z * p[ID_V111]));
459 x1 * (y1 * p[ID_E00z] + y * p[ID_E01z]) +
460 x * (y1 * p[ID_E10z] + y * p[ID_E11z]) +
461 y1 * (z1 * p[ID_Ex00] + z * p[ID_Ex01]) +
462 y * (z1 * p[ID_Ex10] + z * p[ID_Ex11]) +
463 z1 * (x1 * p[ID_E0y0] + x * p[ID_E1y0]) +
464 z * (x1 * p[ID_E0y1] + x * p[ID_E1y1]);
469 //=======================================================================
470 //function : Constructor
472 //=======================================================================
474 SMESH_Block::SMESH_Block():
477 myTolerance(-1.) // to be re-initialized
482 //=======================================================================
483 //function : NbVariables
485 //=======================================================================
487 Standard_Integer SMESH_Block::NbVariables() const
492 //=======================================================================
493 //function : NbEquations
495 //=======================================================================
497 Standard_Integer SMESH_Block::NbEquations() const
502 //=======================================================================
505 //=======================================================================
507 Standard_Boolean SMESH_Block::Value(const math_Vector& theXYZ, math_Vector& theFxyz)
509 gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
510 if ( params.IsEqual( myParam, DBL_MIN )) { // same param
511 theFxyz( 1 ) = funcValue( myValues[ SQUARE_DIST ]);
514 ShellPoint( params, P );
515 gp_Vec dP( P - myPoint );
516 theFxyz(1) = funcValue( dP.SquareMagnitude() );
521 //=======================================================================
522 //function : Derivatives
524 //=======================================================================
526 Standard_Boolean SMESH_Block::Derivatives(const math_Vector& XYZ,math_Matrix& Df)
529 return Values(XYZ,F,Df);
532 //=======================================================================
533 //function : GetStateNumber
535 //=======================================================================
537 Standard_Integer SMESH_Block::GetStateNumber ()
539 return 0; //myValues[0] < 1e-1;
542 //=======================================================================
545 //=======================================================================
547 Standard_Boolean SMESH_Block::Values(const math_Vector& theXYZ,
548 math_Vector& theFxyz,
551 gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
552 if ( params.IsEqual( myParam, DBL_MIN )) { // same param
553 theFxyz( 1 ) = funcValue( myValues[ SQUARE_DIST ] );
554 theDf( 1, DRV_1 ) = myValues[ DRV_1 ];
555 theDf( 1, DRV_2 ) = myValues[ DRV_2 ];
556 theDf( 1, DRV_3 ) = myValues[ DRV_3 ];
559 #ifdef DEBUG_PARAM_COMPUTE
560 MESSAGE ( "PARAM GUESS: " << params.X() << " "<< params.Y() << " "<< params.X() );
561 myNbIterations++; // how many times call ShellPoint()
563 ShellPoint( params, P );
565 gp_Vec dP( myPoint, P );
566 double sqDist = dP.SquareMagnitude();
567 theFxyz(1) = funcValue( sqDist );
569 if ( sqDist < myTolerance * myTolerance ) { // a solution found
571 myValues[ SQUARE_DIST ] = sqDist;
572 theFxyz(1) = theDf( 1,1 ) = theDf( 1,2 ) = theDf( 1,3 ) = 0;
576 if ( sqDist < myValues[ SQUARE_DIST ] ) // a better guess
578 // 3 partial derivatives
579 gp_Vec drv[ 3 ]; // where we move with a small step in each direction
580 for ( int iP = 1; iP <= 3; iP++ ) {
581 if ( iP == myFaceIndex ) {
582 drv[ iP - 1 ] = gp_Vec(0,0,0);
586 bool onEdge = ( theXYZ( iP ) + 0.001 > 1. );
588 params.SetCoord( iP, theXYZ( iP ) - 0.001 );
590 params.SetCoord( iP, theXYZ( iP ) + 0.001 );
591 ShellPoint( params, Pi );
592 params.SetCoord( iP, theXYZ( iP ) ); // restore params
593 gp_Vec dPi ( P, Pi );
594 if ( onEdge ) dPi *= -1.;
595 double mag = dPi.Magnitude();
599 // drv[ iP - 1 ] = dPi / 0.001;
601 for ( int iP = 0; iP < 3; iP++ ) {
603 theDf( 1, iP + 1 ) = dP * drv[iP];
605 // Distance from P to plane passing through myPoint and defined
606 // by the 2 other derivative directions:
607 // like IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P)
608 // where L is (P -> myPoint), P is defined by the 2 other derivative direction
609 int iPrev = ( iP ? iP - 1 : 2 );
610 int iNext = ( iP == 2 ? 0 : iP + 1 );
611 gp_Vec plnNorm = drv[ iPrev ].Crossed( drv [ iNext ] );
612 double Direc = plnNorm * drv[ iP ];
613 if ( Abs(Direc) <= DBL_MIN )
614 theDf( 1, iP + 1 ) = dP * drv[ iP ];
616 double Dis = plnNorm * P - plnNorm * myPoint;
617 theDf( 1, iP + 1 ) = Dis/Direc;
621 #ifdef DEBUG_PARAM_COMPUTE
622 MESSAGE ( "F = " << theFxyz(1) << " DRV: " << theDf(1,1) << " " << theDf(1,2) << " " << theDf(1,3) );
623 myNbIterations +=3; // how many times call ShellPoint()
626 // store better values
628 myValues[SQUARE_DIST]= sqDist;
629 myValues[DRV_1] = theDf(1,DRV_1);
630 myValues[DRV_2] = theDf(1,DRV_2);
631 myValues[DRV_3] = theDf(1,DRV_3);
637 //============================================================================
638 //function : computeParameters
639 //purpose : compute point parameters in the block using math_FunctionSetRoot
640 //============================================================================
642 bool SMESH_Block::computeParameters(const gp_Pnt& thePoint,
644 const gp_XYZ& theParamsHint,
647 myPoint = thePoint.XYZ();
649 myParam.SetCoord( -1,-1,-1 );
650 myValues[ SQUARE_DIST ] = 1e100;
652 math_Vector low ( 1, 3, 0.0 );
653 math_Vector up ( 1, 3, 1.0 );
654 math_Vector tol ( 1, 3, 1e-4 );
655 math_Vector start( 1, 3, 0.0 );
656 start( 1 ) = theParamsHint.X();
657 start( 2 ) = theParamsHint.Y();
658 start( 3 ) = theParamsHint.Z();
660 math_FunctionSetRoot paramSearch( *this, tol );
662 mySquareFunc = 0; // large approaching steps
663 //if ( hasHint ) mySquareFunc = 1; // small approaching steps
665 double loopTol = 10 * myTolerance;
667 while ( distance() > loopTol && nbLoops <= 3 )
669 paramSearch.Perform ( *static_cast<math_FunctionSetWithDerivatives*>(this),
671 start( 1 ) = myParam.X();
672 start( 2 ) = myParam.Y();
673 start( 3 ) = myParam.Z();
674 mySquareFunc = !mySquareFunc;
677 #ifdef DEBUG_PARAM_COMPUTE
678 mySumDist += distance();
679 MESSAGE ( " ------ SOLUTION: ( "<< myParam.X() <<" "<< myParam.Y() <<" "<< myParam.Z() <<" )"<<endl
680 << " ------ DIST : " << distance() << "\t Tol=" << myTolerance << "\t Nb LOOPS=" << nbLoops << endl
681 << " ------ NB IT: " << myNbIterations << ", SUM DIST: " << mySumDist );
686 if ( myFaceIndex > 0 )
688 theParams.SetCoord( myFaceIndex, myFaceParam );
689 if ( distance() > loopTol )
690 refineParametersOnFace( thePoint, theParams, theShapeID );
695 //=======================================================================
696 //function : ComputeParameters
697 //purpose : compute point parameters in the block
698 //=======================================================================
700 bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint,
702 const int theShapeID,
703 const gp_XYZ& theParamsHint)
705 if ( VertexParameters( theShapeID, theParams ))
708 if ( IsEdgeID( theShapeID )) {
709 TEdge& e = myEdge[ theShapeID - ID_FirstE ];
710 Adaptor3d_Curve* curve = e.GetCurve();
711 Extrema_ExtPC anExtPC( thePoint, *curve, curve->FirstParameter(), curve->LastParameter() );
712 int i, nb = anExtPC.IsDone() ? anExtPC.NbExt() : 0;
713 for ( i = 1; i <= nb; i++ ) {
714 if ( anExtPC.IsMin( i ))
715 return EdgeParameters( theShapeID, anExtPC.Point( i ).Parameter(), theParams );
720 const bool isOnFace = IsFaceID( theShapeID );
721 double * coef = GetShapeCoef( theShapeID );
723 // Find the first guess paremeters
725 gp_XYZ start(0, 0, 0);
727 bool hasHint = ( 0 <= theParamsHint.X() && theParamsHint.X() <= 1 &&
728 0 <= theParamsHint.Y() && theParamsHint.Y() <= 1 &&
729 0 <= theParamsHint.Z() && theParamsHint.Z() <= 1 );
730 if ( !hasHint && !myGridComputed )
732 // define the first guess by thePoint projection on lines
733 // connecting vertices
734 bool needGrid = false;
735 gp_XYZ par000( 0, 0, 0 ), par111( 1, 1, 1 );
736 double zero = DBL_MIN * DBL_MIN;
737 for ( int iEdge = 0, iParam = 1; iParam <= 3 && !needGrid; iParam++ )
739 if ( isOnFace && coef[ iParam - 1 ] != 0 ) {
744 for ( int iE = 0; iE < 4; iE++, iEdge++ ) { // loop on 4 parallel edges
745 gp_Pnt p0 = myEdge[ iEdge ].Point( par000 );
746 gp_Pnt p1 = myEdge[ iEdge ].Point( par111 );
747 gp_Vec v01( p0, p1 ), v0P( p0, thePoint );
748 double len2 = v01.SquareMagnitude();
751 par = v0P.Dot( v01 ) / len2;
752 if ( par < 0 || par > 1 ) { // projection falls out of line ends => needGrid
759 start.SetCoord( iParam, sumParam / 4.);
762 // compute nodes of 10 x 10 x 10 grid
765 for ( double x = 0.05; x < 1.; x += 0.1 )
766 for ( double y = 0.05; y < 1.; y += 0.1 )
767 for ( double z = 0.05; z < 1.; z += 0.1 ) {
768 TxyzPair & prmPtn = my3x3x3GridNodes[ iNode++ ];
769 prmPtn.first.SetCoord( x, y, z );
770 ShellPoint( prmPtn.first, prmPtn.second );
771 box.Add( gp_Pnt( prmPtn.second ));
773 myGridComputed = true;
774 if ( myTolerance < 0 )
775 myTolerance = sqrt( box.SquareExtent() ) * 1e-5;
781 start = theParamsHint;
783 if ( myGridComputed )
785 double minDist = DBL_MAX;
789 if ( ShellPoint( start, p ))
790 minDist = thePoint.SquareDistance( p );
792 gp_XYZ* bestParam = 0;
793 for ( int iNode = 0; iNode < 1000; iNode++ ) {
794 TxyzPair & prmPtn = my3x3x3GridNodes[ iNode ];
795 double dist = ( thePoint.XYZ() - prmPtn.second ).SquareModulus();
796 if ( dist < minDist ) {
798 bestParam = & prmPtn.first;
808 // put a point on the face
809 for ( int iCoord = 0; iCoord < 3; iCoord++ )
810 if ( coef[ iCoord ] ) {
811 myFaceIndex = iCoord + 1;
812 myFaceParam = ( coef[ iCoord ] < 0.5 ) ? 0.0 : 1.0;
813 start.SetCoord( myFaceIndex, myFaceParam );
817 #ifdef DEBUG_PARAM_COMPUTE
818 MESSAGE ( " #### POINT " <<thePoint.X()<<" "<<thePoint.Y()<<" "<<thePoint.Z()<<" ####" );
821 if ( myTolerance < 0 ) myTolerance = 1e-6;
823 const double parDelta = 1e-4;
824 const double sqTolerance = myTolerance * myTolerance;
826 gp_XYZ solution = start, params = start;
827 double sqDistance = 1e100;
828 int nbLoops = 0, nbGetWorst = 0;
830 while ( nbLoops <= 100 )
833 ShellPoint( params, P );
835 gp_Vec dP( thePoint, P );
836 double sqDist = dP.SquareMagnitude();
838 if ( sqDist > sqDistance ) { // solution get worse
839 if ( ++nbGetWorst > 2 )
840 return computeParameters( thePoint, theParams, solution, theShapeID );
842 #ifdef DEBUG_PARAM_COMPUTE
843 MESSAGE ( "PARAMS: ( " << params.X() <<" "<< params.Y() <<" "<< params.Z() <<" )" );
844 MESSAGE ( "DIST: " << sqrt( sqDist ) );
847 if ( sqDist < sqDistance ) { // get better
851 if ( sqDistance < sqTolerance ) // a solution found
855 // look for a next better solution
856 for ( int iP = 1; iP <= 3; iP++ ) {
857 if ( iP == myFaceIndex )
859 // see where we move with a small (=parDelta) step in this direction
860 gp_XYZ nearParams = params;
861 bool onEdge = ( params.Coord( iP ) + parDelta > 1. );
863 nearParams.SetCoord( iP, params.Coord( iP ) - parDelta );
865 nearParams.SetCoord( iP, params.Coord( iP ) + parDelta );
866 ShellPoint( nearParams, Pi );
867 gp_Vec dPi ( P, Pi );
868 if ( onEdge ) dPi *= -1.;
869 // modify a parameter
870 double mag = dPi.Magnitude();
873 gp_Vec dir = dPi / mag; // dir we move modifying the parameter
874 double dist = dir * dP; // where we should get to
875 double dPar = dist / mag * parDelta; // predict parameter change
876 double curPar = params.Coord( iP );
877 double par = curPar - dPar; // new parameter value
878 while ( par > 1 || par < 0 ) {
882 params.SetCoord( iP, par );
887 #ifdef DEBUG_PARAM_COMPUTE
888 myNbIterations += nbLoops*4; // how many times ShellPoint called
889 mySumDist += sqrt( sqDistance );
890 MESSAGE ( " ------ SOLUTION: ( "<<solution.X()<<" "<<solution.Y()<<" "<<solution.Z()<<" )"<< std::endl
891 << " ------ DIST : " << sqrt( sqDistance ) << "\t Tol=" << myTolerance << "\t Nb LOOPS=" << nbLoops << std::endl
892 << " ------ NB IT: " << myNbIterations << ", SUM DIST: " << mySumDist );
895 if ( myFaceIndex > 0 )
896 theParams.SetCoord( myFaceIndex, myFaceParam );
898 const double reachedDist = sqrt( sqDistance );
899 // if ( reachedDist > 1000 * myTolerance &&
900 // computeParameters( thePoint, theParams, solution ) &&
901 // reachedDist > distance() )
904 if ( reachedDist > 10 * myTolerance &&
905 computeParameters( thePoint, theParams, solution, theShapeID ) &&
906 reachedDist > distance() )
909 theParams = solution;
910 myValues[ SQUARE_DIST ] = sqDistance;
912 if ( reachedDist > 10 * myTolerance && myFaceIndex > 0 )
913 refineParametersOnFace( thePoint, theParams, theShapeID );
918 //================================================================================
920 * \brief Find more precise solution
921 * \param [in] thePoint - the point
922 * \param [in,out] theParams - solution to precise
923 * \param [in] theFaceID - FACE ID
925 //================================================================================
927 void SMESH_Block::refineParametersOnFace( const gp_Pnt& thePoint,
931 // find UV of thePoint on the FACE
935 const TFace& tface = myFace[ theFaceID - ID_FirstF ];
936 if ( !tface.Surface() ) return;
938 Extrema_ExtPS extPS( thePoint, *tface.Surface(),
939 tface.Surface()->UResolution( myTolerance ),
940 tface.Surface()->VResolution( myTolerance ));
941 if ( !extPS.IsDone() || extPS.NbExt() < 1 )
944 double minDist = 1e100;
945 for ( int i = 1; i <= extPS.NbExt(); ++i )
946 if ( extPS.SquareDistance( i ) < minDist )
948 minDist = extPS.SquareDistance( i );
949 extPS.Point( i ).Parameter( U,V );
951 if ( minDist > 100 * myTolerance * myTolerance )
955 if ( findUVByHalfDivision( thePoint, uv, tface, theParams))
958 int nbGetWorstLimit = 20;
959 if ( findUVAround( thePoint, uv, tface, theParams, nbGetWorstLimit ))
962 double dist2, prevSolDist = distance();
963 gp_XYZ sol = theParams;
964 for ( double delta = 1./10; delta > 0.001; delta /= 2.5, nbGetWorstLimit *= 2 )
966 for ( double y = delta; y < 1.; y += delta )
968 sol.SetCoord( tface.GetVInd(), y );
969 for ( double x = delta; x < 1.; x += delta )
971 sol.SetCoord( tface.GetUInd(), x );
972 dist2 = thePoint.SquareDistance( tface.Point( sol ));
973 if ( dist2 < prevSolDist * prevSolDist )
975 if ( findUVAround( thePoint, uv, tface, theParams, nbGetWorstLimit ))
977 if ( distance() < 1000 * myTolerance )
979 prevSolDist = distance();
986 //================================================================================
988 * \brief Finds parameters corresponding to a given UV of a given face using half-division
989 * \param [in] theUV - the UV to locate
990 * \param [in] tface - the face
991 * \param [in,out] theParams - the starting parameters to improve
992 * \return bool - \c true if found solution is within myTolerance
994 //================================================================================
996 bool SMESH_Block::findUVByHalfDivision( const gp_Pnt& thePoint,
998 const SMESH_Block::TFace& tface,
1001 int nbGetUV = 0; // just for statistics
1003 // find a range of parameters including the UV
1005 double xMin, xMax, yMin, yMax;
1006 //#define _DEBUG_REFINE_
1007 #ifdef _DEBUG_REFINE_
1008 cout << "SMESH_Block::refineParametersOnFace(): dividing Starts at dist " << distance()<< endl;
1010 double dx = 0.1, xSol = theParams.Coord( tface.GetUInd() );
1011 double dy = 0.1, ySol = theParams.Coord( tface.GetVInd() );
1012 gp_XYZ xXYZ( 0,0,0 ); xXYZ.SetCoord( tface.GetUInd(), 1 );
1013 gp_XYZ yXYZ( 0,0,0 ); yXYZ.SetCoord( tface.GetVInd(), 1 );
1014 gp_XYZ xy0,xy1,xy2,xy3;
1015 bool isInQuad = false;
1018 xMin = Max( 0., xSol - 0.5*dx ); xMax = Min( 1.0, xSol + 0.5*dx );
1019 yMin = Max( 0., ySol - 0.5*dy ); yMax = Min( 1.0, ySol + 0.5*dy );
1020 xy0.SetLinearForm( xMin, xXYZ, yMin, yXYZ );
1021 xy1.SetLinearForm( xMax, xXYZ, yMin, yXYZ );
1022 xy2.SetLinearForm( xMax, xXYZ, yMax, yXYZ );
1023 xy3.SetLinearForm( xMin, xXYZ, yMax, yXYZ );
1024 isInQuad = tface.IsUVInQuad( theUV, xy0,xy1,xy2,xy3 );
1030 xSol = 0.5 * (xMax + xMin) ;
1031 ySol = 0.5 * (yMax + yMin) ;
1032 if ( xMin == 0. && yMin == 0. && xMax == 1. && yMax == 1. ) // avoid infinite loop
1034 #ifdef _DEBUG_REFINE_
1035 cout << "SMESH_Block::refineParametersOnFace(): tface.IsUVInQuad() fails" << endl;
1036 cout << " nbGetUV = " << nbGetUV << endl;
1043 // refine solution using half-division technic
1045 gp_XYZ sol = theParams;
1047 const double paramTol = 0.001;
1048 while ( dx > paramTol || dy > paramTol )
1051 bool xDivided = ( dx > paramTol );
1054 double xMid = 0.5 * ( xMin + xMax );
1055 gp_XYZ parMid1 = xMid * xXYZ + yMin * yXYZ;
1056 gp_XYZ parMid2 = xMid * xXYZ + yMax * yXYZ;
1058 if ( tface.IsUVInQuad( theUV, xy0,parMid1,parMid2,xy3 ))
1061 xy1 = parMid1; xy2 = parMid2;
1063 else if ( tface.IsUVInQuad( theUV, parMid1,xy1,xy2,parMid2 ))
1067 xy0 = parMid1; xy3 = parMid2;
1077 bool yDivided = ( dy > paramTol );
1080 double yMid = 0.5 * ( yMin + yMax );
1081 gp_XYZ parMid2 = xMax * xXYZ + yMid * yXYZ;
1082 gp_XYZ parMid3 = xMin * xXYZ + yMid * yXYZ;
1084 if ( tface.IsUVInQuad( theUV, xy0,xy1,parMid2,parMid3 ))
1087 xy2 = parMid2; xy3 = parMid3;
1089 else if ( tface.IsUVInQuad( theUV, parMid3,parMid2,xy2,xy3 ))
1093 xy0 = parMid3; xy1 = parMid2;
1102 if ( !xDivided && !yDivided )
1104 #ifdef _DEBUG_REFINE_
1105 cout << "SMESH_Block::refineParametersOnFace(): nothing divided" << endl;
1106 cout << " nbGetUV = " << nbGetUV << endl;
1111 // evaluate reached distance to thePoint
1112 sol.SetCoord( tface.GetUInd(), 0.5 * ( xMin + xMax ));
1113 sol.SetCoord( tface.GetVInd(), 0.5 * ( yMin + yMax ));
1114 if ( saveBetterSolution( sol, theParams, thePoint.SquareDistance( tface.Point( sol ))))
1116 #ifdef _DEBUG_REFINE_
1117 cout << "SMESH_Block::refineParametersOnFace(): dividing succeeded" << endl;
1118 cout << " nbGetUV = " << nbGetUV << endl;
1123 #ifdef _DEBUG_REFINE_
1124 cout << "SMESH_Block::refineParametersOnFace(): dividing Ends at dist " << distance()<< endl;
1125 cout << " nbGetUV = " << nbGetUV << endl;
1131 //================================================================================
1133 * \brief Finds parameters corresponding to a given UV of a given face by searching
1134 * around the starting solution
1135 * \param [in] theUV - the UV to locate
1136 * \param [in] tface - the face
1137 * \param [in,out] theParams - the starting parameters to improve
1138 * \param [in] nbGetWorstLimit - nb of steps from the starting solution w/o improvement
1139 * to stop searching in this direction
1140 * \return bool - \c true if found solution is within myTolerance
1142 //================================================================================
1144 bool SMESH_Block::findUVAround( const gp_Pnt& thePoint,
1146 const SMESH_Block::TFace& tface,
1148 int nbGetWorstLimit )
1150 #ifdef _DEBUG_REFINE_
1151 cout << "SMESH_Block::refineParametersOnFace(): walk around Starts at dist " << distance()<< endl;
1152 cout << " nbGetUV = " << (nbGetUV=0) << endl;
1154 const double paramTol = 0.001;
1155 const double dx = 0.01, dy = 0.01;
1156 double xMin = theParams.Coord( tface.GetUInd() ), xMax;
1157 double yMin = theParams.Coord( tface.GetVInd() ), yMax;
1159 if ( xMin + dx < 1. )
1162 xMax = 1, xMin = 1 - dx;
1163 gp_XYZ sol = theParams;
1164 sol.SetCoord( tface.GetUInd(), xMax );
1165 sol.SetCoord( tface.GetVInd(), yMax );
1167 if ( saveBetterSolution( sol, theParams, thePoint.SquareDistance( tface.Point( sol ))))
1170 int xMaxNbGetWorst = 0, xMinNbGetWorst = 0, yMaxNbGetWorst = 0, yMinNbGetWorst = 0;
1171 double xMaxBestDist = 1e100, xMinBestDist = 1e100, yMaxBestDist = 1e100, yMinBestDist = 1e100;
1172 double x, y, bestDist, dist;
1173 while ( xMax - xMin < 1 || yMax - yMin < 1 )
1179 for ( x = Max(0.,xMin); x <= xMax+paramTol; x += dx )
1181 y = Max( 0., yMin - dy );
1182 sol.SetCoord( tface.GetUInd(), x );
1183 sol.SetCoord( tface.GetVInd(), y );
1185 dist = thePoint.SquareDistance( tface.Point( sol ));
1186 bestDist = Min( dist, bestDist );
1187 if ( saveBetterSolution( sol, theParams, dist ))
1189 sol.SetCoord( tface.GetUInd(), Min( 1., x + 0.5*dx ));
1190 sol.SetCoord( tface.GetVInd(), y + 0.5*dy );
1192 dist = thePoint.SquareDistance( tface.Point( sol ));
1193 bestDist = Min( dist, bestDist );
1194 if ( saveBetterSolution( sol, theParams, dist ))
1197 yMin = Max(0., yMin-dy );
1198 yMinNbGetWorst += ( yMinBestDist < bestDist );
1199 yMinBestDist = Min( yMinBestDist, bestDist );
1200 if ( yMinNbGetWorst > nbGetWorstLimit )
1206 for ( x = Max(0.,xMin); x <= xMax+paramTol; x += dx )
1208 y = Min( 1., yMax + dy );
1209 sol.SetCoord( tface.GetUInd(), x );
1210 sol.SetCoord( tface.GetVInd(), y );
1212 dist = thePoint.SquareDistance( tface.Point( sol ));
1213 bestDist = Min( dist, bestDist );
1214 if ( saveBetterSolution( sol, theParams, dist ))
1216 sol.SetCoord( tface.GetUInd(), Min( 1., x + 0.5*dx ));
1217 sol.SetCoord( tface.GetVInd(), y - 0.5*dy );
1219 dist = thePoint.SquareDistance( tface.Point( sol ));
1220 bestDist = Min( dist, bestDist );
1221 if ( saveBetterSolution( sol, theParams, dist ))
1224 yMax = Min(1., yMax+dy );
1225 yMaxNbGetWorst += ( yMaxBestDist < bestDist );
1226 yMaxBestDist = Min( yMaxBestDist, bestDist );
1227 if ( yMaxNbGetWorst > nbGetWorstLimit )
1234 for ( y = Max(0.,yMin); y <= yMax+paramTol; y += dy )
1236 x = Max( 0., xMin - dx );
1237 sol.SetCoord( tface.GetUInd(), x );
1238 sol.SetCoord( tface.GetVInd(), y );
1240 dist = thePoint.SquareDistance( tface.Point( sol ));
1241 bestDist = Min( dist, bestDist );
1242 if ( saveBetterSolution( sol, theParams, dist ))
1244 sol.SetCoord( tface.GetUInd(), x + 0.5*dx );
1245 sol.SetCoord( tface.GetVInd(), Min( 1., y + 0.5*dy ));
1247 dist = thePoint.SquareDistance( tface.Point( sol ));
1248 bestDist = Min( dist, bestDist );
1249 if ( saveBetterSolution( sol, theParams, dist ))
1252 xMin = Max(0., xMin-dx );
1253 xMinNbGetWorst += ( xMinBestDist < bestDist );
1254 xMinBestDist = Min( xMinBestDist, bestDist );
1255 if ( xMinNbGetWorst > nbGetWorstLimit )
1261 for ( y = Max(0.,yMin); y <= yMax+paramTol; y += dy )
1263 x = Min( 1., xMax + dx );
1264 sol.SetCoord( tface.GetUInd(), x );
1265 sol.SetCoord( tface.GetVInd(), y );
1267 dist = thePoint.SquareDistance( tface.Point( sol ));
1268 bestDist = Min( dist, bestDist );
1269 if ( saveBetterSolution( sol, theParams, dist ))
1271 sol.SetCoord( tface.GetUInd(), x - 0.5*dx);
1272 sol.SetCoord( tface.GetVInd(), Min( 1., y + 0.5*dy ));
1274 dist = thePoint.SquareDistance( tface.Point( sol ));
1275 bestDist = Min( dist, bestDist );
1276 if ( saveBetterSolution( sol, theParams, dist ))
1279 xMax = Min(1., xMax+dx );
1280 xMaxNbGetWorst += ( xMaxBestDist < bestDist );
1281 xMaxBestDist = Min( xMaxBestDist, bestDist );
1282 if ( xMaxNbGetWorst > nbGetWorstLimit )
1286 #ifdef _DEBUG_REFINE_
1287 cout << "SMESH_Block::refineParametersOnFace(): walk around failed at dist " << distance()<< endl;
1288 //cout << " nbGetUV = " << nbGetUV << endl;
1294 //================================================================================
1296 * \brief Store a solution if it's better than a previous one
1297 * \param [in] theNewParams - a new solution
1298 * \param [out] theParams - the parameters to store solution in
1299 * \param [in] sqDistance - a square distance reached at theNewParams
1300 * \return bool - true if the reached distance is within the tolerance
1302 //================================================================================
1304 bool SMESH_Block::saveBetterSolution( const gp_XYZ& theNewParams,
1308 if ( myValues[ SQUARE_DIST ] > sqDistance )
1310 myValues[ SQUARE_DIST ] = sqDistance;
1311 theParams = theNewParams;
1312 if ( distance() <= myTolerance )
1318 //=======================================================================
1319 //function : SetTolerance
1320 //purpose : set tolerance for ComputeParameters()
1321 //=======================================================================
1323 void SMESH_Block::SetTolerance(const double tol)
1329 //=======================================================================
1330 //function : IsToleranceReached
1331 //purpose : return true if solution found by ComputeParameters() is within the tolerance
1332 //=======================================================================
1334 bool SMESH_Block::IsToleranceReached() const
1336 return distance() < myTolerance;
1339 //=======================================================================
1340 //function : VertexParameters
1341 //purpose : return parameters of a vertex given by TShapeID
1342 //=======================================================================
1344 bool SMESH_Block::VertexParameters(const int theVertexID, gp_XYZ& theParams)
1346 switch ( theVertexID ) {
1347 case ID_V000: theParams.SetCoord(0., 0., 0.); return true;
1348 case ID_V100: theParams.SetCoord(1., 0., 0.); return true;
1349 case ID_V110: theParams.SetCoord(1., 1., 0.); return true;
1350 case ID_V010: theParams.SetCoord(0., 1., 0.); return true;
1356 //=======================================================================
1357 //function : EdgeParameters
1358 //purpose : return parameters of a point given by theU on edge
1359 //=======================================================================
1361 bool SMESH_Block::EdgeParameters(const int theEdgeID, const double theU, gp_XYZ& theParams)
1363 if ( IsEdgeID( theEdgeID )) {
1364 vector< int > vertexVec;
1365 GetEdgeVertexIDs( theEdgeID, vertexVec );
1366 VertexParameters( vertexVec[0], theParams );
1367 TEdge& e = myEdge[ theEdgeID - ID_Ex00 ];
1368 double param = ( theU - e.EndParam(0) ) / ( e.EndParam(1) - e.EndParam(0) );
1369 theParams.SetCoord( e.CoordInd(), param );
1375 //=======================================================================
1376 //function : DumpShapeID
1377 //purpose : debug an id of a block sub-shape
1378 //=======================================================================
1380 #define CASEDUMP(id,strm) case id: strm << #id; break;
1382 ostream& SMESH_Block::DumpShapeID (const int id, ostream& stream)
1385 CASEDUMP( ID_V000, stream );
1386 CASEDUMP( ID_V100, stream );
1387 CASEDUMP( ID_V010, stream );
1388 CASEDUMP( ID_V110, stream );
1389 CASEDUMP( ID_V001, stream );
1390 CASEDUMP( ID_V101, stream );
1391 CASEDUMP( ID_V011, stream );
1392 CASEDUMP( ID_V111, stream );
1393 CASEDUMP( ID_Ex00, stream );
1394 CASEDUMP( ID_Ex10, stream );
1395 CASEDUMP( ID_Ex01, stream );
1396 CASEDUMP( ID_Ex11, stream );
1397 CASEDUMP( ID_E0y0, stream );
1398 CASEDUMP( ID_E1y0, stream );
1399 CASEDUMP( ID_E0y1, stream );
1400 CASEDUMP( ID_E1y1, stream );
1401 CASEDUMP( ID_E00z, stream );
1402 CASEDUMP( ID_E10z, stream );
1403 CASEDUMP( ID_E01z, stream );
1404 CASEDUMP( ID_E11z, stream );
1405 CASEDUMP( ID_Fxy0, stream );
1406 CASEDUMP( ID_Fxy1, stream );
1407 CASEDUMP( ID_Fx0z, stream );
1408 CASEDUMP( ID_Fx1z, stream );
1409 CASEDUMP( ID_F0yz, stream );
1410 CASEDUMP( ID_F1yz, stream );
1411 CASEDUMP( ID_Shell, stream );
1412 default: stream << "ID_INVALID";
1417 //=======================================================================
1418 //function : GetShapeIDByParams
1419 //purpose : define an id of the block sub-shape by normlized point coord
1420 //=======================================================================
1422 int SMESH_Block::GetShapeIDByParams ( const gp_XYZ& theCoord )
1424 // id ( 0 - 26 ) computation:
1426 // vertex ( 0 - 7 ) : id = 1*x + 2*y + 4*z
1428 // edge || X ( 8 - 11 ) : id = 8 + 1*y + 2*z
1429 // edge || Y ( 12 - 15 ): id = 1*x + 12 + 2*z
1430 // edge || Z ( 16 - 19 ): id = 1*x + 2*y + 16
1432 // face || XY ( 20 - 21 ): id = 8 + 12 + 1*z - 0
1433 // face || XZ ( 22 - 23 ): id = 8 + 1*y + 16 - 2
1434 // face || YZ ( 24 - 25 ): id = 1*x + 12 + 16 - 4
1436 static int iAddBnd[] = { 1, 2, 4 };
1437 static int iAddNotBnd[] = { 8, 12, 16 };
1438 static int iFaceSubst[] = { 0, 2, 4 };
1441 int iOnBoundary = 0;
1442 for ( int iCoord = 0; iCoord < 3; iCoord++ )
1444 double val = theCoord.Coord( iCoord + 1 );
1447 else if ( val == 1.0 )
1448 id += iAddBnd[ iOnBoundary++ ];
1450 id += iAddNotBnd[ iCoord ];
1452 if ( iOnBoundary == 1 ) // face
1453 id -= iFaceSubst[ (id - 20) / 4 ];
1454 else if ( iOnBoundary == 0 ) // shell
1457 if ( id > 26 || id < 0 ) {
1458 MESSAGE( "GetShapeIDByParams() = " << id
1459 <<" "<< theCoord.X() <<" "<< theCoord.Y() <<" "<< theCoord.Z() );
1462 return id + 1; // shape ids start at 1
1465 //================================================================================
1467 * \brief Return number of wires and a list of oredered edges.
1468 * \param theFace - the face to process
1469 * \param theEdges - all ordered edges of theFace (outer edges go first).
1470 * \param theNbEdgesInWires - nb of edges (== nb of vertices in closed wire) in each wire
1471 * \param theFirstVertex - the vertex of the outer wire to set first in the returned
1472 * list ( theFirstVertex may be NULL )
1473 * \param theShapeAnalysisAlgo - if true, ShapeAnalysis::OuterWire() is used to find
1474 * the outer wire else BRepTools::OuterWire() is used.
1475 * \retval int - nb of wires
1477 * Always try to set a seam edge first.
1478 * BRepTools::OuterWire() fails e.g. in the case of issue 0020184,
1479 * ShapeAnalysis::OuterWire() fails in the case of issue 0020452
1481 //================================================================================
1483 int SMESH_Block::GetOrderedEdges (const TopoDS_Face& theFace,
1484 list< TopoDS_Edge >& theEdges,
1485 list< int > & theNbEdgesInWires,
1486 TopoDS_Vertex theFirstVertex,
1487 const bool theShapeAnalysisAlgo)
1489 // put wires in a list, so that an outer wire comes first
1490 list<TopoDS_Wire> aWireList;
1491 TopoDS_Wire anOuterWire =
1492 theShapeAnalysisAlgo ? ShapeAnalysis::OuterWire( theFace ) : BRepTools::OuterWire( theFace );
1493 for ( TopoDS_Iterator wIt (theFace); wIt.More(); wIt.Next() )
1494 if ( wIt.Value().ShapeType() == TopAbs_WIRE ) // it can be internal vertex!
1496 if ( !anOuterWire.IsSame( wIt.Value() ))
1497 aWireList.push_back( TopoDS::Wire( wIt.Value() ));
1499 aWireList.push_front( TopoDS::Wire( wIt.Value() ));
1502 // loop on edges of wires
1503 theNbEdgesInWires.clear();
1504 list<TopoDS_Wire>::iterator wlIt = aWireList.begin();
1505 for ( ; wlIt != aWireList.end(); wlIt++ )
1508 BRepTools_WireExplorer wExp( *wlIt, theFace );
1509 for ( iE = 0; wExp.More(); wExp.Next(), iE++ )
1511 TopoDS_Edge edge = wExp.Current();
1512 // commented for issue 0020557, other related ones: 0020526, PAL19080
1513 // edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
1514 theEdges.push_back( edge );
1516 if ( iE == 0 ) // wExp returns nothing if e.g. the wire contains one internal edge
1518 for ( TopoDS_Iterator e( *wlIt ); e.More(); e.Next(), ++iE )
1519 theEdges.push_back( TopoDS::Edge( e.Value() ));
1521 theNbEdgesInWires.push_back( iE );
1523 if ( wlIt == aWireList.begin() && theEdges.size() > 1 ) { // the outer wire
1524 // orient closed edges
1525 list< TopoDS_Edge >::iterator eIt, eIt2;
1526 for ( eIt = theEdges.begin(); eIt != theEdges.end(); eIt++ )
1528 TopoDS_Edge& edge = *eIt;
1529 if ( TopExp::FirstVertex( edge ).IsSame( TopExp::LastVertex( edge ) ))
1532 bool isNext = ( eIt2 == theEdges.begin() );
1533 TopoDS_Edge edge2 = isNext ? *(++eIt2) : *(--eIt2);
1535 Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( edge, theFace, f1,l1 );
1536 Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( edge2, theFace, f2,l2 );
1537 gp_Pnt2d pf = c1->Value( edge.Orientation() == TopAbs_FORWARD ? f1 : l1 );
1538 gp_Pnt2d pl = c1->Value( edge.Orientation() == TopAbs_FORWARD ? l1 : f1 );
1539 bool isFirst = ( edge2.Orientation() == TopAbs_FORWARD ? isNext : !isNext );
1540 gp_Pnt2d p2 = c2->Value( isFirst ? f2 : l2 );
1541 isFirst = ( p2.SquareDistance( pf ) < p2.SquareDistance( pl ));
1542 if ( isNext ? isFirst : !isFirst )
1544 // to make a seam go first
1545 if ( theFirstVertex.IsNull() )
1546 theFirstVertex = TopExp::FirstVertex( edge, true );
1549 // rotate theEdges until it begins from theFirstVertex
1550 if ( ! theFirstVertex.IsNull() ) {
1551 TopoDS_Vertex vv[2];
1552 TopExp::Vertices( theEdges.front(), vv[0], vv[1], true );
1553 // on closed face, make seam edge the first in the list
1554 while ( !vv[0].IsSame( theFirstVertex ) || vv[0].IsSame( vv[1] ))
1556 theEdges.splice(theEdges.end(), theEdges,
1557 theEdges.begin(), ++theEdges.begin());
1558 TopExp::Vertices( theEdges.front(), vv[0], vv[1], true );
1559 if ( iE++ > theNbEdgesInWires.back() ) {
1561 gp_Pnt p = BRep_Tool::Pnt( theFirstVertex );
1562 MESSAGE ( " : Warning : vertex "<< theFirstVertex.TShape().operator->()
1563 << " ( " << p.X() << " " << p.Y() << " " << p.Z() << " )"
1564 << " not found in outer wire of face "<< theFace.TShape().operator->()
1565 << " with vertices: " );
1566 wExp.Init( *wlIt, theFace );
1567 for ( int i = 0; wExp.More(); wExp.Next(), i++ )
1569 TopoDS_Edge edge = wExp.Current();
1570 edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
1571 TopoDS_Vertex v = TopExp::FirstVertex( edge, true );
1572 gp_Pnt p = BRep_Tool::Pnt( v );
1573 MESSAGE_ADD ( i << " " << v.TShape().operator->() << " "
1574 << p.X() << " " << p.Y() << " " << p.Z() << " " << std::endl );
1577 break; // break infinite loop
1584 return aWireList.size();
1586 //================================================================================
1588 * \brief Call it after geometry initialisation
1590 //================================================================================
1592 void SMESH_Block::init()
1596 myGridComputed = false;
1599 //=======================================================================
1600 //function : LoadMeshBlock
1601 //purpose : prepare to work with theVolume
1602 //=======================================================================
1604 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
1606 bool SMESH_Block::LoadMeshBlock(const SMDS_MeshVolume* theVolume,
1607 const int theNode000Index,
1608 const int theNode001Index,
1609 vector<const SMDS_MeshNode*>& theOrderedNodes)
1613 SMDS_VolumeTool vTool;
1614 if (!vTool.Set( theVolume ) || vTool.NbNodes() != 8 ||
1615 !vTool.IsLinked( theNode000Index, theNode001Index )) {
1616 MESSAGE(" Bad arguments ");
1619 vTool.SetExternalNormal();
1620 // In terms of indices used for access to nodes and faces in SMDS_VolumeTool:
1621 int V000, V100, V010, V110, V001, V101, V011, V111; // 8 vertices
1622 int Fxy0, Fxy1; // bottom and top faces
1623 // vertices of faces
1624 vector<int> vFxy0, vFxy1;
1626 V000 = theNode000Index;
1627 V001 = theNode001Index;
1629 // get faces sharing V000 and V001
1630 list<int> fV000, fV001;
1632 for ( iF = 0; iF < vTool.NbFaces(); ++iF ) {
1633 const int* nid = vTool.GetFaceNodesIndices( iF );
1634 for ( iN = 0; iN < 4; ++iN )
1635 if ( nid[ iN ] == V000 ) {
1636 fV000.push_back( iF );
1637 } else if ( nid[ iN ] == V001 ) {
1638 fV001.push_back( iF );
1642 // find the bottom (Fxy0), the top (Fxy1) faces
1643 list<int>::iterator fIt1, fIt2, Fxy0Pos;
1644 for ( fIt1 = fV000.begin(); fIt1 != fV000.end(); fIt1++) {
1645 fIt2 = std::find( fV001.begin(), fV001.end(), *fIt1 );
1646 if ( fIt2 != fV001.end() ) { // *fIt1 is in the both lists
1647 fV001.erase( fIt2 ); // erase Fx0z or F0yz from fV001
1648 } else { // *fIt1 is in fV000 only
1649 Fxy0Pos = fIt1; // points to Fxy0
1653 Fxy1 = fV001.front();
1654 const SMDS_MeshNode** nn = vTool.GetNodes();
1656 // find bottom veritices, their order is that a face normal is external
1658 const int* nid = vTool.GetFaceNodesIndices( Fxy0 );
1659 for ( i = 0; i < 4; ++i )
1660 if ( nid[ i ] == V000 )
1662 for ( iN = 0; iN < 4; ++iN, ++i ) {
1663 if ( i == 4 ) i = 0;
1664 vFxy0[ iN ] = nid[ i ];
1666 // find top veritices, their order is that a face normal is external
1668 nid = vTool.GetFaceNodesIndices( Fxy1 );
1669 for ( i = 0; i < 4; ++i )
1670 if ( nid[ i ] == V001 )
1672 for ( iN = 0; iN < 4; ++iN, ++i ) {
1673 if ( i == 4 ) i = 0;
1674 vFxy1[ iN ] = nid[ i ];
1676 // find indices of the rest veritices
1684 // set points coordinates
1685 myPnt[ ID_V000 - 1 ] = gpXYZ( nn[ V000 ] );
1686 myPnt[ ID_V100 - 1 ] = gpXYZ( nn[ V100 ] );
1687 myPnt[ ID_V010 - 1 ] = gpXYZ( nn[ V010 ] );
1688 myPnt[ ID_V110 - 1 ] = gpXYZ( nn[ V110 ] );
1689 myPnt[ ID_V001 - 1 ] = gpXYZ( nn[ V001 ] );
1690 myPnt[ ID_V101 - 1 ] = gpXYZ( nn[ V101 ] );
1691 myPnt[ ID_V011 - 1 ] = gpXYZ( nn[ V011 ] );
1692 myPnt[ ID_V111 - 1 ] = gpXYZ( nn[ V111 ] );
1694 // fill theOrderedNodes
1695 theOrderedNodes.resize( 8 );
1696 theOrderedNodes[ 0 ] = nn[ V000 ];
1697 theOrderedNodes[ 1 ] = nn[ V100 ];
1698 theOrderedNodes[ 2 ] = nn[ V010 ];
1699 theOrderedNodes[ 3 ] = nn[ V110 ];
1700 theOrderedNodes[ 4 ] = nn[ V001 ];
1701 theOrderedNodes[ 5 ] = nn[ V101 ];
1702 theOrderedNodes[ 6 ] = nn[ V011 ];
1703 theOrderedNodes[ 7 ] = nn[ V111 ];
1706 vector< int > vertexVec;
1707 for ( iE = 0; iE < NbEdges(); ++iE ) {
1708 GetEdgeVertexIDs(( iE + ID_FirstE ), vertexVec );
1709 myEdge[ iE ].Set(( iE + ID_FirstE ),
1710 myPnt[ vertexVec[0] - 1 ],
1711 myPnt[ vertexVec[1] - 1 ]);
1714 // fill faces' corners
1715 for ( iF = ID_Fxy0; iF < ID_Shell; ++iF )
1717 TFace& tFace = myFace[ iF - ID_FirstF ];
1718 vector< int > edgeIdVec(4, -1);
1719 GetFaceEdgesIDs( iF, edgeIdVec );
1720 tFace.Set( iF, myEdge[ edgeIdVec [ 0 ] - ID_Ex00], myEdge[ edgeIdVec [ 1 ] - ID_Ex00]);
1726 //=======================================================================
1727 //function : LoadBlockShapes
1728 //purpose : Initialize block geometry with theShell,
1729 // add sub-shapes of theBlock to theShapeIDMap so that they get
1730 // IDs acoording to enum TShapeID
1731 //=======================================================================
1733 bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell,
1734 const TopoDS_Vertex& theVertex000,
1735 const TopoDS_Vertex& theVertex001,
1736 TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
1738 return ( FindBlockShapes( theShell, theVertex000, theVertex001, theShapeIDMap ) &&
1739 LoadBlockShapes( theShapeIDMap ));
1742 //=======================================================================
1743 //function : LoadBlockShapes
1744 //purpose : add sub-shapes of theBlock to theShapeIDMap so that they get
1745 // IDs acoording to enum TShapeID
1746 //=======================================================================
1748 bool SMESH_Block::FindBlockShapes(const TopoDS_Shell& theShell,
1749 const TopoDS_Vertex& theVertex000,
1750 const TopoDS_Vertex& theVertex001,
1751 TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
1754 TopoDS_Shape V000, V100, V010, V110, V001, V101, V011, V111;
1756 TopoDS_Shape Ex00, Ex10, Ex01, Ex11;
1757 TopoDS_Shape E0y0, E1y0, E0y1, E1y1;
1758 TopoDS_Shape E00z, E10z, E01z, E11z;
1760 TopoDS_Shape Fxy0, Fx0z, F0yz, Fxy1, Fx1z, F1yz;
1762 // nb of faces bound to a vertex in TopTools_IndexedDataMapOfShapeListOfShape
1763 // filled by TopExp::MapShapesAndAncestors()
1764 const int NB_FACES_BY_VERTEX = 6;
1766 TopTools_IndexedDataMapOfShapeListOfShape vfMap;
1767 TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_FACE, vfMap );
1768 if ( vfMap.Extent() != 8 ) {
1769 MESSAGE(" Wrong nb of vertices in the block: " << vfMap.Extent() );
1773 V000 = theVertex000;
1774 V001 = theVertex001;
1776 if ( V000.IsNull() ) {
1777 // find vertex 000 - the one with smallest coordinates
1778 double minVal = DBL_MAX, minX = DBL_MAX, val;
1779 for ( int i = 1; i <= 8; i++ ) {
1780 const TopoDS_Vertex& v = TopoDS::Vertex( vfMap.FindKey( i ));
1781 gp_Pnt P = BRep_Tool::Pnt( v );
1782 val = P.X() + P.Y() + P.Z();
1783 if ( val < minVal || ( val == minVal && P.X() < minX )) {
1789 // find vertex 001 - the one on the most vertical edge passing through V000
1790 TopTools_IndexedDataMapOfShapeListOfShape veMap;
1791 TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_EDGE, veMap );
1792 gp_Vec dir001 = gp::DZ();
1793 gp_Pnt p000 = BRep_Tool::Pnt( TopoDS::Vertex( V000 ));
1794 double maxVal = -DBL_MAX;
1795 TopTools_ListIteratorOfListOfShape eIt ( veMap.FindFromKey( V000 ));
1796 for ( ; eIt.More(); eIt.Next() ) {
1797 const TopoDS_Edge& e = TopoDS::Edge( eIt.Value() );
1798 TopoDS_Vertex v = TopExp::FirstVertex( e );
1799 if ( v.IsSame( V000 )) {
1800 v = TopExp::LastVertex( e );
1801 if ( v.IsSame( V000 ))
1804 val = dir001 * gp_Vec( p000, BRep_Tool::Pnt( v )).Normalized();
1805 if ( val > maxVal ) {
1812 // find the bottom (Fxy0), Fx0z and F0yz faces
1814 const TopTools_ListOfShape& f000List = vfMap.FindFromKey( V000 );
1815 const TopTools_ListOfShape& f001List = vfMap.FindFromKey( V001 );
1816 if (f000List.Extent() != NB_FACES_BY_VERTEX ||
1817 f001List.Extent() != NB_FACES_BY_VERTEX ) {
1818 MESSAGE(" LoadBlockShapes() " << f000List.Extent() << " " << f001List.Extent());
1821 TopTools_ListIteratorOfListOfShape f001It, f000It ( f000List );
1822 int i, j, iFound1=0, iFound2=0;
1823 for ( j = 0; f000It.More(); f000It.Next(), j++ )
1825 if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1826 const TopoDS_Shape& F = f000It.Value();
1827 for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1828 if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1829 if ( F.IsSame( f001It.Value() ))
1832 if ( f001It.More() ) // Fx0z or F0yz found
1833 if ( Fx0z.IsNull() ) {
1840 else // F is the bottom face
1843 if ( Fxy0.IsNull() || Fx0z.IsNull() || F0yz.IsNull() ) {
1844 MESSAGE( Fxy0.IsNull() <<" "<< Fx0z.IsNull() <<" "<< F0yz.IsNull() );
1848 // choose the top face (Fxy1)
1849 for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1850 if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1851 if ( i != iFound1 && i != iFound2 )
1854 Fxy1 = f001It.Value();
1855 if ( Fxy1.IsNull() ) {
1856 MESSAGE(" LoadBlockShapes() error ");
1860 // find bottom edges and veritices
1861 list< TopoDS_Edge > eList;
1862 list< int > nbVertexInWires;
1863 GetOrderedEdges( TopoDS::Face( Fxy0 ), eList, nbVertexInWires, TopoDS::Vertex( V000 ) );
1864 if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1865 MESSAGE(" LoadBlockShapes() error ");
1868 list< TopoDS_Edge >::iterator elIt = eList.begin();
1869 for ( i = 0; elIt != eList.end(); elIt++, i++ )
1871 case 0: E0y0 = *elIt; V010 = TopExp::LastVertex( *elIt, true ); break;
1872 case 1: Ex10 = *elIt; V110 = TopExp::LastVertex( *elIt, true ); break;
1873 case 2: E1y0 = *elIt; V100 = TopExp::LastVertex( *elIt, true ); break;
1874 case 3: Ex00 = *elIt; break;
1877 if ( i != 4 || E0y0.IsNull() || Ex10.IsNull() || E1y0.IsNull() || Ex00.IsNull() ) {
1878 MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1883 // find top edges and veritices
1885 GetOrderedEdges( TopoDS::Face( Fxy1 ), eList, nbVertexInWires, TopoDS::Vertex( V001 ) );
1886 if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1887 MESSAGE(" LoadBlockShapes() error ");
1890 for ( i = 0, elIt = eList.begin(); elIt != eList.end(); elIt++, i++ )
1892 case 0: Ex01 = *elIt; V101 = TopExp::LastVertex( *elIt, true ); break;
1893 case 1: E1y1 = *elIt; V111 = TopExp::LastVertex( *elIt, true ); break;
1894 case 2: Ex11 = *elIt; V011 = TopExp::LastVertex( *elIt, true ); break;
1895 case 3: E0y1 = *elIt; break;
1898 if ( i != 4 || Ex01.IsNull() || E1y1.IsNull() || Ex11.IsNull() || E0y1.IsNull() ) {
1899 MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1903 // swap Fx0z and F0yz if necessary
1904 TopExp_Explorer exp( Fx0z, TopAbs_VERTEX );
1905 for ( ; exp.More(); exp.Next() ) // Fx0z shares V101 and V100
1906 if ( V101.IsSame( exp.Current() ) || V100.IsSame( exp.Current() ))
1907 break; // V101 or V100 found
1908 if ( !exp.More() ) { // not found
1909 std::swap( Fx0z, F0yz);
1912 // find Fx1z and F1yz faces
1913 const TopTools_ListOfShape& f111List = vfMap.FindFromKey( V111 );
1914 const TopTools_ListOfShape& f110List = vfMap.FindFromKey( V110 );
1915 if (f111List.Extent() != NB_FACES_BY_VERTEX ||
1916 f110List.Extent() != NB_FACES_BY_VERTEX ) {
1917 MESSAGE(" LoadBlockShapes() " << f111List.Extent() << " " << f110List.Extent());
1920 TopTools_ListIteratorOfListOfShape f111It, f110It ( f110List);
1921 for ( j = 0 ; f110It.More(); f110It.Next(), j++ ) {
1922 if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1923 const TopoDS_Shape& F = f110It.Value();
1924 for ( i = 0, f111It.Initialize( f111List ); f111It.More(); f111It.Next(), i++ ) {
1925 if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1926 if ( F.IsSame( f111It.Value() )) { // Fx1z or F1yz found
1927 if ( Fx1z.IsNull() )
1934 if ( Fx1z.IsNull() || F1yz.IsNull() ) {
1935 MESSAGE(" LoadBlockShapes() error ");
1939 // swap Fx1z and F1yz if necessary
1940 for ( exp.Init( Fx1z, TopAbs_VERTEX ); exp.More(); exp.Next() )
1941 if ( V010.IsSame( exp.Current() ) || V011.IsSame( exp.Current() ))
1943 if ( !exp.More() ) {
1944 std::swap( Fx1z, F1yz);
1947 // find vertical edges
1948 for ( exp.Init( Fx0z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1949 const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1950 const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1951 if ( vFirst.IsSame( V001 ))
1953 else if ( vFirst.IsSame( V100 ))
1956 if ( E00z.IsNull() || E10z.IsNull() ) {
1957 MESSAGE(" LoadBlockShapes() error ");
1960 for ( exp.Init( Fx1z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1961 const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1962 const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1963 if ( vFirst.IsSame( V111 ))
1965 else if ( vFirst.IsSame( V010 ))
1968 if ( E01z.IsNull() || E11z.IsNull() ) {
1969 MESSAGE(" LoadBlockShapes() error ");
1973 // load shapes in theShapeIDMap
1975 theShapeIDMap.Clear();
1977 theShapeIDMap.Add(V000.Oriented( TopAbs_FORWARD ));
1978 theShapeIDMap.Add(V100.Oriented( TopAbs_FORWARD ));
1979 theShapeIDMap.Add(V010.Oriented( TopAbs_FORWARD ));
1980 theShapeIDMap.Add(V110.Oriented( TopAbs_FORWARD ));
1981 theShapeIDMap.Add(V001.Oriented( TopAbs_FORWARD ));
1982 theShapeIDMap.Add(V101.Oriented( TopAbs_FORWARD ));
1983 theShapeIDMap.Add(V011.Oriented( TopAbs_FORWARD ));
1984 theShapeIDMap.Add(V111.Oriented( TopAbs_FORWARD ));
1986 theShapeIDMap.Add(Ex00);
1987 theShapeIDMap.Add(Ex10);
1988 theShapeIDMap.Add(Ex01);
1989 theShapeIDMap.Add(Ex11);
1991 theShapeIDMap.Add(E0y0);
1992 theShapeIDMap.Add(E1y0);
1993 theShapeIDMap.Add(E0y1);
1994 theShapeIDMap.Add(E1y1);
1996 theShapeIDMap.Add(E00z);
1997 theShapeIDMap.Add(E10z);
1998 theShapeIDMap.Add(E01z);
1999 theShapeIDMap.Add(E11z);
2001 theShapeIDMap.Add(Fxy0);
2002 theShapeIDMap.Add(Fxy1);
2003 theShapeIDMap.Add(Fx0z);
2004 theShapeIDMap.Add(Fx1z);
2005 theShapeIDMap.Add(F0yz);
2006 theShapeIDMap.Add(F1yz);
2008 theShapeIDMap.Add(theShell);
2013 //================================================================================
2015 * \brief Initialize block geometry with shapes from theShapeIDMap
2016 * \param theShapeIDMap - map of block sub-shapes
2017 * \retval bool - is a success
2019 //================================================================================
2021 bool SMESH_Block::LoadBlockShapes(const TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
2025 // store shapes geometry
2026 for ( int shapeID = 1; shapeID < theShapeIDMap.Extent(); shapeID++ )
2028 const TopoDS_Shape& S = theShapeIDMap( shapeID );
2029 switch ( S.ShapeType() )
2031 case TopAbs_VERTEX: {
2033 if ( !IsVertexID( ID_V111 )) return false;
2034 myPnt[ shapeID - ID_V000 ] = BRep_Tool::Pnt( TopoDS::Vertex( S )).XYZ();
2039 if ( !IsEdgeID( shapeID )) return false;
2040 const TopoDS_Edge& edge = TopoDS::Edge( S );
2041 TEdge& tEdge = myEdge[ shapeID - ID_FirstE ];
2043 new BRepAdaptor_Curve( edge ),
2044 IsForwardEdge( edge, theShapeIDMap ));
2049 if ( !LoadFace( TopoDS::Face( S ), shapeID, theShapeIDMap ))
2055 } // loop on shapes in theShapeIDMap
2060 //================================================================================
2062 * \brief Load face geometry
2063 * \param theFace - face
2064 * \param theFaceID - face in-block ID
2065 * \param theShapeIDMap - map of block sub-shapes
2066 * \retval bool - is a success
2068 * It is enough to compute params or coordinates on the face.
2069 * Face sub-shapes must be loaded into theShapeIDMap before
2071 //================================================================================
2073 bool SMESH_Block::LoadFace(const TopoDS_Face& theFace,
2074 const int theFaceID,
2075 const TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
2077 if ( !IsFaceID( theFaceID ) ) return false;
2079 Adaptor2d_Curve2d* c2d[4];
2081 vector< int > edgeIdVec;
2082 GetFaceEdgesIDs( theFaceID, edgeIdVec );
2083 for ( size_t iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges
2085 if ( edgeIdVec[ iE ] > theShapeIDMap.Extent() )
2087 const TopoDS_Edge& edge = TopoDS::Edge( theShapeIDMap( edgeIdVec[ iE ]));
2088 c2d[ iE ] = new BRepAdaptor_Curve2d( edge, theFace );
2089 isForward[ iE ] = IsForwardEdge( edge, theShapeIDMap );
2091 TFace& tFace = myFace[ theFaceID - ID_FirstF ];
2092 tFace.Set( theFaceID, new BRepAdaptor_Surface( theFace ), c2d, isForward );
2096 //================================================================================
2098 * \brief/ Insert theShape into theShapeIDMap with theShapeID
2099 * \param theShape - shape to insert
2100 * \param theShapeID - shape in-block ID
2101 * \param theShapeIDMap - map of block sub-shapes
2103 //================================================================================
2105 bool SMESH_Block::Insert(const TopoDS_Shape& theShape,
2106 const int theShapeID,
2107 TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
2109 if ( !theShape.IsNull() && theShapeID > 0 )
2111 if ( theShapeIDMap.Contains( theShape ))
2112 return ( theShapeIDMap.FindIndex( theShape ) == theShapeID );
2114 if ( theShapeID <= theShapeIDMap.Extent() ) {
2115 theShapeIDMap.Substitute( theShapeID, theShape );
2118 while ( theShapeIDMap.Extent() < theShapeID - 1 ) {
2119 TopoDS_Compound comp;
2120 BRep_Builder().MakeCompound( comp );
2121 theShapeIDMap.Add( comp );
2123 theShapeIDMap.Add( theShape );
2130 //=======================================================================
2131 //function : GetFaceEdgesIDs
2132 //purpose : return edges IDs in the order u0, u1, 0v, 1v
2133 // u0 means "|| u, v == 0"
2134 //=======================================================================
2136 void SMESH_Block::GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec )
2138 edgeVec.resize( 4 );
2141 edgeVec[ 0 ] = ID_Ex00;
2142 edgeVec[ 1 ] = ID_Ex10;
2143 edgeVec[ 2 ] = ID_E0y0;
2144 edgeVec[ 3 ] = ID_E1y0;
2147 edgeVec[ 0 ] = ID_Ex01;
2148 edgeVec[ 1 ] = ID_Ex11;
2149 edgeVec[ 2 ] = ID_E0y1;
2150 edgeVec[ 3 ] = ID_E1y1;
2153 edgeVec[ 0 ] = ID_Ex00;
2154 edgeVec[ 1 ] = ID_Ex01;
2155 edgeVec[ 2 ] = ID_E00z;
2156 edgeVec[ 3 ] = ID_E10z;
2159 edgeVec[ 0 ] = ID_Ex10;
2160 edgeVec[ 1 ] = ID_Ex11;
2161 edgeVec[ 2 ] = ID_E01z;
2162 edgeVec[ 3 ] = ID_E11z;
2165 edgeVec[ 0 ] = ID_E0y0;
2166 edgeVec[ 1 ] = ID_E0y1;
2167 edgeVec[ 2 ] = ID_E00z;
2168 edgeVec[ 3 ] = ID_E01z;
2171 edgeVec[ 0 ] = ID_E1y0;
2172 edgeVec[ 1 ] = ID_E1y1;
2173 edgeVec[ 2 ] = ID_E10z;
2174 edgeVec[ 3 ] = ID_E11z;
2177 MESSAGE(" GetFaceEdgesIDs(), wrong face ID: " << faceID );
2181 //=======================================================================
2182 //function : GetEdgeVertexIDs
2183 //purpose : return vertex IDs of an edge
2184 //=======================================================================
2186 void SMESH_Block::GetEdgeVertexIDs (const int edgeID, vector< int >& vertexVec )
2188 vertexVec.resize( 2 );
2192 vertexVec[ 0 ] = ID_V000;
2193 vertexVec[ 1 ] = ID_V100;
2196 vertexVec[ 0 ] = ID_V010;
2197 vertexVec[ 1 ] = ID_V110;
2200 vertexVec[ 0 ] = ID_V001;
2201 vertexVec[ 1 ] = ID_V101;
2204 vertexVec[ 0 ] = ID_V011;
2205 vertexVec[ 1 ] = ID_V111;
2209 vertexVec[ 0 ] = ID_V000;
2210 vertexVec[ 1 ] = ID_V010;
2213 vertexVec[ 0 ] = ID_V100;
2214 vertexVec[ 1 ] = ID_V110;
2217 vertexVec[ 0 ] = ID_V001;
2218 vertexVec[ 1 ] = ID_V011;
2221 vertexVec[ 0 ] = ID_V101;
2222 vertexVec[ 1 ] = ID_V111;
2226 vertexVec[ 0 ] = ID_V000;
2227 vertexVec[ 1 ] = ID_V001;
2230 vertexVec[ 0 ] = ID_V100;
2231 vertexVec[ 1 ] = ID_V101;
2234 vertexVec[ 0 ] = ID_V010;
2235 vertexVec[ 1 ] = ID_V011;
2238 vertexVec[ 0 ] = ID_V110;
2239 vertexVec[ 1 ] = ID_V111;
2242 vertexVec.resize(0);
2243 MESSAGE(" GetEdgeVertexIDs(), wrong edge ID: " << edgeID );