1 // Copyright (C) 2007-2008 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.
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
22 // File : SMESH_Pattern.hxx
23 // Created : Mon Aug 2 10:30:00 2004
24 // Author : Edward AGAPOV (eap)
26 #include "SMESH_Block.hxx"
28 #include <BRepAdaptor_Curve.hxx>
29 #include <BRepAdaptor_Curve2d.hxx>
30 #include <BRepAdaptor_Surface.hxx>
31 #include <BRepTools.hxx>
32 #include <BRepTools_WireExplorer.hxx>
33 #include <BRep_Builder.hxx>
34 #include <BRep_Tool.hxx>
35 #include <Bnd_Box.hxx>
36 #include <Extrema_ExtPC.hxx>
37 #include <Extrema_POnCurv.hxx>
38 #include <Geom2d_Curve.hxx>
39 #include <ShapeAnalysis.hxx>
40 #include <TopExp_Explorer.hxx>
41 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
42 #include <TopTools_ListIteratorOfListOfShape.hxx>
43 #include <TopTools_ListOfShape.hxx>
45 #include <TopoDS_Compound.hxx>
46 #include <TopoDS_Face.hxx>
47 #include <TopoDS_Iterator.hxx>
48 #include <TopoDS_Wire.hxx>
49 #include <gp_Trsf.hxx>
51 #include <math_FunctionSetRoot.hxx>
52 #include <math_Matrix.hxx>
53 #include <math_Vector.hxx>
55 #include "SMDS_MeshNode.hxx"
56 #include "SMDS_MeshVolume.hxx"
57 #include "SMDS_VolumeTool.hxx"
58 #include "utilities.h"
64 //#define DEBUG_PARAM_COMPUTE
66 //================================================================================
68 * \brief Set edge data
69 * \param edgeID - block subshape ID
70 * \param curve - edge geometry
71 * \param isForward - is curve orientation coincides with edge orientation in the block
73 //================================================================================
75 void SMESH_Block::TEdge::Set( const int edgeID, Adaptor3d_Curve* curve, const bool isForward )
77 myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID );
78 if ( myC3d ) delete myC3d;
80 myFirst = curve->FirstParameter();
81 myLast = curve->LastParameter();
83 std::swap( myFirst, myLast );
86 //================================================================================
88 * \brief Set coordinates of nodes at edge ends to work with mesh block
89 * \param edgeID - block subshape ID
90 * \param node1 - coordinates of node with lower ID
91 * \param node2 - coordinates of node with upper ID
93 //================================================================================
95 void SMESH_Block::TEdge::Set( const int edgeID, const gp_XYZ& node1, const gp_XYZ& node2 )
97 myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID );
101 if ( myC3d ) delete myC3d;
105 //=======================================================================
106 //function : SMESH_Block::TEdge::GetU
108 //=======================================================================
110 double SMESH_Block::TEdge::GetU( const gp_XYZ& theParams ) const
112 double u = theParams.Coord( myCoordInd );
113 if ( !myC3d ) // if mesh block
115 return ( 1 - u ) * myFirst + u * myLast;
118 //=======================================================================
119 //function : SMESH_Block::TEdge::Point
121 //=======================================================================
123 gp_XYZ SMESH_Block::TEdge::Point( const gp_XYZ& theParams ) const
125 double u = GetU( theParams );
126 if ( myC3d ) return myC3d->Value( u ).XYZ();
128 return myNodes[0] * ( 1 - u ) + myNodes[1] * u;
131 //================================================================================
135 //================================================================================
137 SMESH_Block::TEdge::~TEdge()
139 if ( myC3d ) delete myC3d;
142 //================================================================================
144 * \brief Set face data
145 * \param faceID - block subshape ID
146 * \param S - face surface geometry
147 * \param c2d - 4 pcurves in the order as returned by GetFaceEdgesIDs(faceID)
148 * \param isForward - orientation of pcurves comparing with block edge direction
150 //================================================================================
152 void SMESH_Block::TFace::Set( const int faceID,
153 Adaptor3d_Surface* S,
154 Adaptor2d_Curve2d* c2D[4],
155 const bool isForward[4] )
157 if ( myS ) delete myS;
160 vector< int > edgeIdVec;
161 GetFaceEdgesIDs( faceID, edgeIdVec );
162 for ( int iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges
164 myCoordInd[ iE ] = GetCoordIndOnEdge( edgeIdVec[ iE ] );
165 if ( myC2d[ iE ]) delete myC2d[ iE ];
166 myC2d[ iE ] = c2D[ iE ];
167 myFirst[ iE ] = myC2d[ iE ]->FirstParameter();
168 myLast [ iE ] = myC2d[ iE ]->LastParameter();
169 if ( !isForward[ iE ])
170 std::swap( myFirst[ iE ], myLast[ iE ] );
173 myCorner[ 0 ] = myC2d[ 0 ]->Value( myFirst[0] ).XY();
174 myCorner[ 1 ] = myC2d[ 0 ]->Value( myLast[0] ).XY();
175 myCorner[ 2 ] = myC2d[ 1 ]->Value( myLast[1] ).XY();
176 myCorner[ 3 ] = myC2d[ 1 ]->Value( myFirst[1] ).XY();
179 //================================================================================
181 * \brief Set face data to work with mesh block
182 * \param faceID - block subshape ID
183 * \param edgeU0 - filled data of edge u0 = GetFaceEdgesIDs(faceID)[ 0 ]
184 * \param edgeU1 - filled data of edge u1 = GetFaceEdgesIDs(faceID)[ 1 ]
186 //================================================================================
188 void SMESH_Block::TFace::Set( const int faceID, const TEdge& edgeU0, const TEdge& edgeU1 )
190 vector< int > edgeIdVec;
191 GetFaceEdgesIDs( faceID, edgeIdVec );
192 myNodes[ 0 ] = edgeU0.NodeXYZ( 1 );
193 myNodes[ 1 ] = edgeU0.NodeXYZ( 0 );
194 myNodes[ 2 ] = edgeU1.NodeXYZ( 0 );
195 myNodes[ 3 ] = edgeU1.NodeXYZ( 1 );
196 myCoordInd[ 0 ] = GetCoordIndOnEdge( edgeIdVec[ 0 ] );
197 myCoordInd[ 1 ] = GetCoordIndOnEdge( edgeIdVec[ 1 ] );
198 myCoordInd[ 2 ] = GetCoordIndOnEdge( edgeIdVec[ 2 ] );
199 myCoordInd[ 3 ] = GetCoordIndOnEdge( edgeIdVec[ 3 ] );
200 if ( myS ) delete myS;
204 //================================================================================
208 //================================================================================
210 SMESH_Block::TFace::~TFace()
212 if ( myS ) delete myS;
213 for ( int i = 0 ; i < 4; ++i )
214 if ( myC2d[ i ]) delete myC2d[ i ];
217 //=======================================================================
218 //function : SMESH_Block::TFace::GetCoefs
219 //purpose : return coefficients for addition of [0-3]-th edge and vertex
220 //=======================================================================
222 void SMESH_Block::TFace::GetCoefs(int iE,
223 const gp_XYZ& theParams,
225 double& Vcoef ) const
227 double dU = theParams.Coord( GetUInd() );
228 double dV = theParams.Coord( GetVInd() );
231 Ecoef = ( 1 - dV ); // u0
232 Vcoef = ( 1 - dU ) * ( 1 - dV ); break; // 00
235 Vcoef = dU * ( 1 - dV ); break; // 10
237 Ecoef = ( 1 - dU ); // 0v
238 Vcoef = dU * dV ; break; // 11
241 Vcoef = ( 1 - dU ) * dV ; break; // 01
246 //=======================================================================
247 //function : SMESH_Block::TFace::GetUV
249 //=======================================================================
251 gp_XY SMESH_Block::TFace::GetUV( const gp_XYZ& theParams ) const
254 for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
256 double Ecoef = 0, Vcoef = 0;
257 GetCoefs( iE, theParams, Ecoef, Vcoef );
259 double u = theParams.Coord( myCoordInd[ iE ] );
260 u = ( 1 - u ) * myFirst[ iE ] + u * myLast[ iE ];
261 uv += Ecoef * myC2d[ iE ]->Value( u ).XY();
263 uv -= Vcoef * myCorner[ iE ];
268 //=======================================================================
269 //function : SMESH_Block::TFace::Point
271 //=======================================================================
273 gp_XYZ SMESH_Block::TFace::Point( const gp_XYZ& theParams ) const
276 if ( !myS ) // if mesh block
278 for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
280 double Ecoef = 0, Vcoef = 0;
281 GetCoefs( iE, theParams, Ecoef, Vcoef );
283 double u = theParams.Coord( myCoordInd[ iE ] );
286 case 1: i1 = 3; i2 = 2; break;
287 case 2: i1 = 1; i2 = 2; break;
288 case 3: i1 = 0; i2 = 3; break;
290 p += Ecoef * ( myNodes[ i1 ] * ( 1 - u ) + myNodes[ i2 ] * u );
292 p -= Vcoef * myNodes[ iE ];
298 gp_XY uv = GetUV( theParams );
299 p = myS->Value( uv.X(), uv.Y() ).XYZ();
304 //=======================================================================
305 //function : GetShapeCoef
307 //=======================================================================
309 double* SMESH_Block::GetShapeCoef (const int theShapeID)
311 static double shapeCoef[][3] = {
312 // V000, V100, V010, V110
313 { -1,-1,-1 }, { 1,-1,-1 }, { -1, 1,-1 }, { 1, 1,-1 },
314 // V001, V101, V011, V111,
315 { -1,-1, 1 }, { 1,-1, 1 }, { -1, 1, 1 }, { 1, 1, 1 },
316 // Ex00, Ex10, Ex01, Ex11,
317 { 0,-1,-1 }, { 0, 1,-1 }, { 0,-1, 1 }, { 0, 1, 1 },
318 // E0y0, E1y0, E0y1, E1y1,
319 { -1, 0,-1 }, { 1, 0,-1 }, { -1, 0, 1 }, { 1, 0, 1 },
320 // E00z, E10z, E01z, E11z,
321 { -1,-1, 0 }, { 1,-1, 0 }, { -1, 1, 0 }, { 1, 1, 0 },
322 // Fxy0, Fxy1, Fx0z, Fx1z, F0yz, F1yz,
323 { 0, 0,-1 }, { 0, 0, 1 }, { 0,-1, 0 }, { 0, 1, 0 }, { -1, 0, 0 }, { 1, 0, 0 },
327 if ( theShapeID < ID_V000 || theShapeID > ID_F1yz )
328 return shapeCoef[ ID_Shell - 1 ];
330 return shapeCoef[ theShapeID - 1 ];
333 //=======================================================================
334 //function : ShellPoint
335 //purpose : return coordinates of a point in shell
336 //=======================================================================
338 bool SMESH_Block::ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const
340 thePoint.SetCoord( 0., 0., 0. );
341 for ( int shapeID = ID_V000; shapeID < ID_Shell; shapeID++ )
344 double* coefs = GetShapeCoef( shapeID );
346 for ( int iCoef = 0; iCoef < 3; iCoef++ ) {
347 if ( coefs[ iCoef ] != 0 ) {
348 if ( coefs[ iCoef ] < 0 )
349 k *= ( 1. - theParams.Coord( iCoef + 1 ));
351 k *= theParams.Coord( iCoef + 1 );
354 // add point on a shape
355 if ( fabs( k ) > DBL_MIN )
358 if ( shapeID < ID_Ex00 ) // vertex
359 VertexPoint( shapeID, Ps );
360 else if ( shapeID < ID_Fxy0 ) { // edge
361 EdgePoint( shapeID, theParams, Ps );
364 FacePoint( shapeID, theParams, Ps );
372 //=======================================================================
373 //function : ShellPoint
374 //purpose : computes coordinates of a point in shell by points on sub-shapes;
375 // thePointOnShape[ subShapeID ] must be a point on a subShape
376 //=======================================================================
378 bool SMESH_Block::ShellPoint(const gp_XYZ& theParams,
379 const vector<gp_XYZ>& thePointOnShape,
382 if ( thePointOnShape.size() < ID_F1yz )
385 double x = theParams.X(), y = theParams.Y(), z = theParams.Z();
386 double x1 = 1. - x, y1 = 1. - y, z1 = 1. - z;
387 const vector<gp_XYZ>& p = thePointOnShape;
390 x1 * p[ID_F0yz] + x * p[ID_F1yz] +
391 y1 * p[ID_Fx0z] + y * p[ID_Fx1z] +
392 z1 * p[ID_Fxy0] + z * p[ID_Fxy1] +
393 x1 * (y1 * (z1 * p[ID_V000] + z * p[ID_V001]) +
394 y * (z1 * p[ID_V010] + z * p[ID_V011])) +
395 x * (y1 * (z1 * p[ID_V100] + z * p[ID_V101]) +
396 y * (z1 * p[ID_V110] + z * p[ID_V111]));
398 x1 * (y1 * p[ID_E00z] + y * p[ID_E01z]) +
399 x * (y1 * p[ID_E10z] + y * p[ID_E11z]) +
400 y1 * (z1 * p[ID_Ex00] + z * p[ID_Ex01]) +
401 y * (z1 * p[ID_Ex10] + z * p[ID_Ex11]) +
402 z1 * (x1 * p[ID_E0y0] + x * p[ID_E1y0]) +
403 z * (x1 * p[ID_E0y1] + x * p[ID_E1y1]);
408 //=======================================================================
409 //function : Constructor
411 //=======================================================================
413 SMESH_Block::SMESH_Block():
416 myTolerance(-1.) // to be re-initialized
421 //=======================================================================
422 //function : NbVariables
424 //=======================================================================
426 Standard_Integer SMESH_Block::NbVariables() const
431 //=======================================================================
432 //function : NbEquations
434 //=======================================================================
436 Standard_Integer SMESH_Block::NbEquations() const
441 //=======================================================================
444 //=======================================================================
446 Standard_Boolean SMESH_Block::Value(const math_Vector& theXYZ, math_Vector& theFxyz)
448 gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
449 if ( params.IsEqual( myParam, DBL_MIN )) { // same param
450 theFxyz( 1 ) = funcValue( myValues[ SQUARE_DIST ]);
453 ShellPoint( params, P );
454 gp_Vec dP( P - myPoint );
455 theFxyz(1) = funcValue( dP.SquareMagnitude() );
460 //=======================================================================
461 //function : Derivatives
463 //=======================================================================
465 Standard_Boolean SMESH_Block::Derivatives(const math_Vector& XYZ,math_Matrix& Df)
468 return Values(XYZ,F,Df);
471 //=======================================================================
472 //function : GetStateNumber
474 //=======================================================================
476 Standard_Integer SMESH_Block::GetStateNumber ()
478 return 0; //myValues[0] < 1e-1;
481 //=======================================================================
484 //=======================================================================
486 Standard_Boolean SMESH_Block::Values(const math_Vector& theXYZ,
487 math_Vector& theFxyz,
490 gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
491 if ( params.IsEqual( myParam, DBL_MIN )) { // same param
492 theFxyz( 1 ) = funcValue( myValues[ SQUARE_DIST ] );
493 theDf( 1, DRV_1 ) = myValues[ DRV_1 ];
494 theDf( 1, DRV_2 ) = myValues[ DRV_2 ];
495 theDf( 1, DRV_3 ) = myValues[ DRV_3 ];
498 #ifdef DEBUG_PARAM_COMPUTE
499 MESSAGE ( "PARAM GUESS: " << params.X() << " "<< params.Y() << " "<< params.X() );
500 myNbIterations++; // how many times call ShellPoint()
502 ShellPoint( params, P );
504 gp_Vec dP( myPoint, P );
505 double sqDist = dP.SquareMagnitude();
506 theFxyz(1) = funcValue( sqDist );
508 if ( sqDist < myTolerance * myTolerance ) { // a solution found
510 myValues[ SQUARE_DIST ] = sqDist;
511 theFxyz(1) = theDf( 1,1 ) = theDf( 1,2 ) = theDf( 1,3 ) = 0;
515 if ( sqDist < myValues[ SQUARE_DIST ] ) // a better guess
517 // 3 partial derivatives
518 gp_Vec drv[ 3 ]; // where we move with a small step in each direction
519 for ( int iP = 1; iP <= 3; iP++ ) {
520 if ( iP == myFaceIndex ) {
521 drv[ iP - 1 ] = gp_Vec(0,0,0);
525 bool onEdge = ( theXYZ( iP ) + 0.001 > 1. );
527 params.SetCoord( iP, theXYZ( iP ) - 0.001 );
529 params.SetCoord( iP, theXYZ( iP ) + 0.001 );
530 ShellPoint( params, Pi );
531 params.SetCoord( iP, theXYZ( iP ) ); // restore params
532 gp_Vec dPi ( P, Pi );
533 if ( onEdge ) dPi *= -1.;
534 double mag = dPi.Magnitude();
539 for ( int iP = 0; iP < 3; iP++ ) {
541 theDf( 1, iP + 1 ) = dP * drv[iP];
543 // Distance from P to plane passing through myPoint and defined
544 // by the 2 other derivative directions:
545 // like IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P)
546 // where L is (P -> myPoint), P is defined by the 2 other derivative direction
547 int iPrev = ( iP ? iP - 1 : 2 );
548 int iNext = ( iP == 2 ? 0 : iP + 1 );
549 gp_Vec plnNorm = drv[ iPrev ].Crossed( drv [ iNext ] );
550 double Direc = plnNorm * drv[ iP ];
551 if ( Abs(Direc) <= DBL_MIN )
552 theDf( 1, iP + 1 ) = dP * drv[ iP ];
554 double Dis = plnNorm * P - plnNorm * myPoint;
555 theDf( 1, iP + 1 ) = Dis/Direc;
559 #ifdef DEBUG_PARAM_COMPUTE
560 MESSAGE ( "F = " << theFxyz(1) << " DRV: " << theDf(1,1) << " " << theDf(1,2) << " " << theDf(1,3) );
561 myNbIterations +=3; // how many times call ShellPoint()
564 // store better values
566 myValues[SQUARE_DIST]= sqDist;
567 myValues[DRV_1] = theDf(1,DRV_1);
568 myValues[DRV_2] = theDf(1,DRV_2);
569 myValues[DRV_3] = theDf(1,DRV_3);
575 //============================================================================
576 //function : computeParameters
577 //purpose : compute point parameters in the block using math_FunctionSetRoot
578 //============================================================================
580 bool SMESH_Block::computeParameters(const gp_Pnt& thePoint,
582 const gp_XYZ& theParamsHint)
584 myPoint = thePoint.XYZ();
586 myParam.SetCoord( -1,-1,-1 );
587 myValues[ SQUARE_DIST ] = 1e100;
589 math_Vector low ( 1, 3, 0.0 );
590 math_Vector up ( 1, 3, 1.0 );
591 math_Vector tol ( 1, 3, 1e-4 );
592 math_Vector start( 1, 3, 0.0 );
593 start( 1 ) = theParamsHint.X();
594 start( 2 ) = theParamsHint.Y();
595 start( 3 ) = theParamsHint.Z();
597 math_FunctionSetRoot paramSearch( *this, tol );
599 mySquareFunc = 0; // large approaching steps
600 //if ( hasHint ) mySquareFunc = 1; // small approaching steps
602 double loopTol = 10 * myTolerance;
604 while ( distance() > loopTol && nbLoops <= 3 )
606 paramSearch.Perform ( *static_cast<math_FunctionSetWithDerivatives*>(this),
608 start( 1 ) = myParam.X();
609 start( 2 ) = myParam.Y();
610 start( 3 ) = myParam.Z();
611 mySquareFunc = !mySquareFunc;
614 #ifdef DEBUG_PARAM_COMPUTE
615 mySumDist += distance();
616 MESSAGE ( " ------ SOLUTION: ( "<< myParam.X() <<" "<< myParam.Y() <<" "<< myParam.Z() <<" )"<<endl
617 << " ------ DIST : " << distance() << "\t Tol=" << myTolerance << "\t Nb LOOPS=" << nbLoops << endl
618 << " ------ NB IT: " << myNbIterations << ", SUM DIST: " << mySumDist );
623 if ( myFaceIndex > 0 )
624 theParams.SetCoord( myFaceIndex, myFaceParam );
629 //=======================================================================
630 //function : ComputeParameters
631 //purpose : compute point parameters in the block
632 //=======================================================================
634 bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint,
636 const int theShapeID,
637 const gp_XYZ& theParamsHint)
639 if ( VertexParameters( theShapeID, theParams ))
642 if ( IsEdgeID( theShapeID )) {
643 TEdge& e = myEdge[ theShapeID - ID_FirstE ];
644 Adaptor3d_Curve* curve = e.GetCurve();
645 Extrema_ExtPC anExtPC( thePoint, *curve, curve->FirstParameter(), curve->LastParameter() );
646 int i, nb = anExtPC.IsDone() ? anExtPC.NbExt() : 0;
647 for ( i = 1; i <= nb; i++ ) {
648 if ( anExtPC.IsMin( i ))
649 return EdgeParameters( theShapeID, anExtPC.Point( i ).Parameter(), theParams );
654 const bool isOnFace = IsFaceID( theShapeID );
655 double * coef = GetShapeCoef( theShapeID );
657 // Find the first guess paremeters
659 gp_XYZ start(0, 0, 0);
661 bool hasHint = ( 0 <= theParamsHint.X() && theParamsHint.X() <= 1 &&
662 0 <= theParamsHint.Y() && theParamsHint.Y() <= 1 &&
663 0 <= theParamsHint.Y() && theParamsHint.Y() <= 1 );
664 if ( !hasHint && !myGridComputed )
666 // define the first guess by thePoint projection on lines
667 // connecting vertices
668 bool needGrid = false;
669 gp_XYZ par000( 0, 0, 0 ), par111( 1, 1, 1 );
670 double zero = DBL_MIN * DBL_MIN;
671 for ( int iEdge = 0, iParam = 1; iParam <= 3 && !needGrid; iParam++ )
673 if ( isOnFace && coef[ iParam - 1 ] != 0 ) {
678 for ( int iE = 0; iE < 4; iE++, iEdge++ ) { // loop on 4 parallel edges
679 gp_Pnt p0 = myEdge[ iEdge ].Point( par000 );
680 gp_Pnt p1 = myEdge[ iEdge ].Point( par111 );
681 gp_Vec v01( p0, p1 ), v0P( p0, thePoint );
682 double len2 = v01.SquareMagnitude();
685 par = v0P.Dot( v01 ) / len2;
686 if ( par < 0 || par > 1 ) { // projection falls out of line ends => needGrid
693 start.SetCoord( iParam, sumParam / 4.);
696 // compute nodes of 3 x 3 x 3 grid
699 for ( double x = 0.25; x < 0.9; x += 0.25 )
700 for ( double y = 0.25; y < 0.9; y += 0.25 )
701 for ( double z = 0.25; z < 0.9; z += 0.25 ) {
702 TxyzPair & prmPtn = my3x3x3GridNodes[ iNode++ ];
703 prmPtn.first.SetCoord( x, y, z );
704 ShellPoint( prmPtn.first, prmPtn.second );
705 box.Add( gp_Pnt( prmPtn.second ));
707 myGridComputed = true;
708 myTolerance = sqrt( box.SquareExtent() ) * 1e-5;
714 start = theParamsHint;
716 else if ( myGridComputed )
718 double minDist = DBL_MAX;
719 gp_XYZ* bestParam = 0;
720 for ( int iNode = 0; iNode < 27; iNode++ ) {
721 TxyzPair & prmPtn = my3x3x3GridNodes[ iNode ];
722 double dist = ( thePoint.XYZ() - prmPtn.second ).SquareModulus();
723 if ( dist < minDist ) {
725 bestParam = & prmPtn.first;
734 // put a point on the face
735 for ( int iCoord = 0; iCoord < 3; iCoord++ )
736 if ( coef[ iCoord ] ) {
737 myFaceIndex = iCoord + 1;
738 myFaceParam = ( coef[ iCoord ] < 0.5 ) ? 0.0 : 1.0;
739 start.SetCoord( myFaceIndex, myFaceParam );
743 #ifdef DEBUG_PARAM_COMPUTE
744 MESSAGE ( " #### POINT " <<thePoint.X()<<" "<<thePoint.Y()<<" "<<thePoint.Z()<<" ####" );
747 if ( myTolerance < 0 ) myTolerance = 1e-6;
749 const double parDelta = 1e-4;
750 const double sqTolerance = myTolerance * myTolerance;
752 gp_XYZ solution = start, params = start;
753 double sqDistance = 1e100;
754 int nbLoops = 0, nbGetWorst = 0;
756 while ( nbLoops <= 100 )
759 ShellPoint( params, P );
761 gp_Vec dP( thePoint, P );
762 double sqDist = dP.SquareMagnitude();
764 if ( sqDist > sqDistance ) { // solution get worse
765 if ( ++nbGetWorst > 2 )
766 return computeParameters( thePoint, theParams, solution );
768 #ifdef DEBUG_PARAM_COMPUTE
769 MESSAGE ( "PARAMS: ( " << params.X() <<" "<< params.Y() <<" "<< params.Z() <<" )" );
770 MESSAGE ( "DIST: " << sqrt( sqDist ) );
773 if ( sqDist < sqDistance ) { // get better
777 if ( sqDistance < sqTolerance ) // a solution found
781 // look for a next better solution
782 for ( int iP = 1; iP <= 3; iP++ ) {
783 if ( iP == myFaceIndex )
785 // see where we move with a small (=parDelta) step in this direction
786 gp_XYZ nearParams = params;
787 bool onEdge = ( params.Coord( iP ) + parDelta > 1. );
789 nearParams.SetCoord( iP, params.Coord( iP ) - parDelta );
791 nearParams.SetCoord( iP, params.Coord( iP ) + parDelta );
792 ShellPoint( nearParams, Pi );
793 gp_Vec dPi ( P, Pi );
794 if ( onEdge ) dPi *= -1.;
795 // modify a parameter
796 double mag = dPi.Magnitude();
799 gp_Vec dir = dPi / mag; // dir we move modifying the parameter
800 double dist = dir * dP; // where we should get to
801 double dPar = dist / mag * parDelta; // predict parameter change
802 double curPar = params.Coord( iP );
803 double par = curPar - dPar; // new parameter value
804 while ( par > 1 || par < 0 ) {
808 params.SetCoord( iP, par );
813 #ifdef DEBUG_PARAM_COMPUTE
814 myNbIterations += nbLoops*4; // how many times ShellPoint called
815 mySumDist += sqrt( sqDistance );
816 MESSAGE ( " ------ SOLUTION: ( "<<solution.X()<<" "<<solution.Y()<<" "<<solution.Z()<<" )"<< std::endl
817 << " ------ DIST : " << sqrt( sqDistance ) << "\t Tol=" << myTolerance << "\t Nb LOOPS=" << nbLoops << std::endl
818 << " ------ NB IT: " << myNbIterations << ", SUM DIST: " << mySumDist );
821 theParams = solution;
823 if ( myFaceIndex > 0 )
824 theParams.SetCoord( myFaceIndex, myFaceParam );
829 //=======================================================================
830 //function : VertexParameters
831 //purpose : return parameters of a vertex given by TShapeID
832 //=======================================================================
834 bool SMESH_Block::VertexParameters(const int theVertexID, gp_XYZ& theParams)
836 switch ( theVertexID ) {
837 case ID_V000: theParams.SetCoord(0., 0., 0.); return true;
838 case ID_V100: theParams.SetCoord(1., 0., 0.); return true;
839 case ID_V110: theParams.SetCoord(1., 1., 0.); return true;
840 case ID_V010: theParams.SetCoord(0., 1., 0.); return true;
846 //=======================================================================
847 //function : EdgeParameters
848 //purpose : return parameters of a point given by theU on edge
849 //=======================================================================
851 bool SMESH_Block::EdgeParameters(const int theEdgeID, const double theU, gp_XYZ& theParams)
853 if ( IsEdgeID( theEdgeID )) {
854 vector< int > vertexVec;
855 GetEdgeVertexIDs( theEdgeID, vertexVec );
856 VertexParameters( vertexVec[0], theParams );
857 TEdge& e = myEdge[ theEdgeID - ID_Ex00 ];
858 double param = ( theU - e.EndParam(0) ) / ( e.EndParam(1) - e.EndParam(0) );
859 theParams.SetCoord( e.CoordInd(), param );
865 //=======================================================================
866 //function : DumpShapeID
867 //purpose : debug an id of a block sub-shape
868 //=======================================================================
870 #define CASEDUMP(id,strm) case id: strm << #id; break;
872 ostream& SMESH_Block::DumpShapeID (const int id, ostream& stream)
875 CASEDUMP( ID_V000, stream );
876 CASEDUMP( ID_V100, stream );
877 CASEDUMP( ID_V010, stream );
878 CASEDUMP( ID_V110, stream );
879 CASEDUMP( ID_V001, stream );
880 CASEDUMP( ID_V101, stream );
881 CASEDUMP( ID_V011, stream );
882 CASEDUMP( ID_V111, stream );
883 CASEDUMP( ID_Ex00, stream );
884 CASEDUMP( ID_Ex10, stream );
885 CASEDUMP( ID_Ex01, stream );
886 CASEDUMP( ID_Ex11, stream );
887 CASEDUMP( ID_E0y0, stream );
888 CASEDUMP( ID_E1y0, stream );
889 CASEDUMP( ID_E0y1, stream );
890 CASEDUMP( ID_E1y1, stream );
891 CASEDUMP( ID_E00z, stream );
892 CASEDUMP( ID_E10z, stream );
893 CASEDUMP( ID_E01z, stream );
894 CASEDUMP( ID_E11z, stream );
895 CASEDUMP( ID_Fxy0, stream );
896 CASEDUMP( ID_Fxy1, stream );
897 CASEDUMP( ID_Fx0z, stream );
898 CASEDUMP( ID_Fx1z, stream );
899 CASEDUMP( ID_F0yz, stream );
900 CASEDUMP( ID_F1yz, stream );
901 CASEDUMP( ID_Shell, stream );
902 default: stream << "ID_INVALID";
907 //=======================================================================
908 //function : GetShapeIDByParams
909 //purpose : define an id of the block sub-shape by normlized point coord
910 //=======================================================================
912 int SMESH_Block::GetShapeIDByParams ( const gp_XYZ& theCoord )
914 // id ( 0 - 26 ) computation:
916 // vertex ( 0 - 7 ) : id = 1*x + 2*y + 4*z
918 // edge || X ( 8 - 11 ) : id = 8 + 1*y + 2*z
919 // edge || Y ( 12 - 15 ): id = 1*x + 12 + 2*z
920 // edge || Z ( 16 - 19 ): id = 1*x + 2*y + 16
922 // face || XY ( 20 - 21 ): id = 8 + 12 + 1*z - 0
923 // face || XZ ( 22 - 23 ): id = 8 + 1*y + 16 - 2
924 // face || YZ ( 24 - 25 ): id = 1*x + 12 + 16 - 4
926 static int iAddBnd[] = { 1, 2, 4 };
927 static int iAddNotBnd[] = { 8, 12, 16 };
928 static int iFaceSubst[] = { 0, 2, 4 };
932 for ( int iCoord = 0; iCoord < 3; iCoord++ )
934 double val = theCoord.Coord( iCoord + 1 );
937 else if ( val == 1.0 )
938 id += iAddBnd[ iOnBoundary++ ];
940 id += iAddNotBnd[ iCoord ];
942 if ( iOnBoundary == 1 ) // face
943 id -= iFaceSubst[ (id - 20) / 4 ];
944 else if ( iOnBoundary == 0 ) // shell
947 if ( id > 26 || id < 0 ) {
948 MESSAGE( "GetShapeIDByParams() = " << id
949 <<" "<< theCoord.X() <<" "<< theCoord.Y() <<" "<< theCoord.Z() );
952 return id + 1; // shape ids start at 1
955 //=======================================================================
956 //function : GetOrderedEdges
957 //purpose : return nb wires and a list of oredered edges
958 //=======================================================================
960 int SMESH_Block::GetOrderedEdges (const TopoDS_Face& theFace,
961 TopoDS_Vertex theFirstVertex,
962 list< TopoDS_Edge >& theEdges,
963 list< int > & theNbVertexInWires)
965 // put wires in a list, so that an outer wire comes first
966 list<TopoDS_Wire> aWireList;
967 //TopoDS_Wire anOuterWire = BRepTools::OuterWire( theFace ); ### issue 0020184
968 TopoDS_Wire anOuterWire = ShapeAnalysis::OuterWire( theFace );
969 //aWireList.push_back( anOuterWire ); ### issue 0020184
970 for ( TopoDS_Iterator wIt (theFace); wIt.More(); wIt.Next() )
971 if ( !anOuterWire.IsSame( wIt.Value() ))
972 aWireList.push_back( TopoDS::Wire( wIt.Value() ));
974 aWireList.push_front( TopoDS::Wire( wIt.Value() ));// ### issue 0020184
976 // loop on edges of wires
977 theNbVertexInWires.clear();
978 list<TopoDS_Wire>::iterator wlIt = aWireList.begin();
979 for ( ; wlIt != aWireList.end(); wlIt++ )
982 BRepTools_WireExplorer wExp( *wlIt, theFace );
983 for ( iE = 0; wExp.More(); wExp.Next(), iE++ )
985 TopoDS_Edge edge = wExp.Current();
986 edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
987 theEdges.push_back( edge );
989 theNbVertexInWires.push_back( iE );
991 if ( wlIt == aWireList.begin() && theEdges.size() > 1 ) { // the outer wire
992 // orient closed edges
993 list< TopoDS_Edge >::iterator eIt, eIt2;
994 for ( eIt = theEdges.begin(); eIt != theEdges.end(); eIt++ )
996 TopoDS_Edge& edge = *eIt;
997 if ( TopExp::FirstVertex( edge ).IsSame( TopExp::LastVertex( edge ) ))
1000 bool isNext = ( eIt2 == theEdges.begin() );
1001 TopoDS_Edge edge2 = isNext ? *(++eIt2) : *(--eIt2);
1003 Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( edge, theFace, f1,l1 );
1004 Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( edge2, theFace, f2,l2 );
1005 gp_Pnt2d pf = c1->Value( edge.Orientation() == TopAbs_FORWARD ? f1 : l1 );
1006 gp_Pnt2d pl = c1->Value( edge.Orientation() == TopAbs_FORWARD ? l1 : f1 );
1007 bool isFirst = ( edge2.Orientation() == TopAbs_FORWARD ? isNext : !isNext );
1008 gp_Pnt2d p2 = c2->Value( isFirst ? f2 : l2 );
1009 isFirst = ( p2.SquareDistance( pf ) < p2.SquareDistance( pl ));
1010 if ( isNext ? isFirst : !isFirst )
1012 // to make a seam go first
1013 if ( theFirstVertex.IsNull() )
1014 theFirstVertex = TopExp::FirstVertex( edge, true );
1017 // rotate theEdges until it begins from theFirstVertex
1018 if ( ! theFirstVertex.IsNull() ) {
1019 TopoDS_Vertex vv[2];
1020 TopExp::Vertices( theEdges.front(), vv[0], vv[1], true );
1021 // on closed face, make seam edge the first in the list
1022 while ( !vv[0].IsSame( theFirstVertex ) || vv[0].IsSame( vv[1] ))
1024 theEdges.splice(theEdges.end(), theEdges,
1025 theEdges.begin(), ++theEdges.begin());
1026 TopExp::Vertices( theEdges.front(), vv[0], vv[1], true );
1027 if ( iE++ > theNbVertexInWires.back() ) {
1029 gp_Pnt p = BRep_Tool::Pnt( theFirstVertex );
1030 MESSAGE ( " : Warning : vertex "<< theFirstVertex.TShape().operator->()
1031 << " ( " << p.X() << " " << p.Y() << " " << p.Z() << " )"
1032 << " not found in outer wire of face "<< theFace.TShape().operator->()
1033 << " with vertices: " );
1034 wExp.Init( *wlIt, theFace );
1035 for ( int i = 0; wExp.More(); wExp.Next(), i++ )
1037 TopoDS_Edge edge = wExp.Current();
1038 edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
1039 TopoDS_Vertex v = TopExp::FirstVertex( edge, true );
1040 gp_Pnt p = BRep_Tool::Pnt( v );
1041 MESSAGE_ADD ( i << " " << v.TShape().operator->() << " "
1042 << p.X() << " " << p.Y() << " " << p.Z() << " " << std::endl );
1045 break; // break infinite loop
1052 return aWireList.size();
1054 //================================================================================
1056 * \brief Call it after geometry initialisation
1058 //================================================================================
1060 void SMESH_Block::init()
1064 myGridComputed = false;
1067 //=======================================================================
1068 //function : LoadMeshBlock
1069 //purpose : prepare to work with theVolume
1070 //=======================================================================
1072 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
1074 bool SMESH_Block::LoadMeshBlock(const SMDS_MeshVolume* theVolume,
1075 const int theNode000Index,
1076 const int theNode001Index,
1077 vector<const SMDS_MeshNode*>& theOrderedNodes)
1079 MESSAGE(" ::LoadMeshBlock()");
1082 SMDS_VolumeTool vTool;
1083 if (!vTool.Set( theVolume ) || vTool.NbNodes() != 8 ||
1084 !vTool.IsLinked( theNode000Index, theNode001Index )) {
1085 MESSAGE(" Bad arguments ");
1088 vTool.SetExternalNormal();
1089 // In terms of indices used for access to nodes and faces in SMDS_VolumeTool:
1090 int V000, V100, V010, V110, V001, V101, V011, V111; // 8 vertices
1091 int Fxy0, Fxy1; // bottom and top faces
1092 // vertices of faces
1093 vector<int> vFxy0, vFxy1;
1095 V000 = theNode000Index;
1096 V001 = theNode001Index;
1098 // get faces sharing V000 and V001
1099 list<int> fV000, fV001;
1101 for ( iF = 0; iF < vTool.NbFaces(); ++iF ) {
1102 const int* nid = vTool.GetFaceNodesIndices( iF );
1103 for ( iN = 0; iN < 4; ++iN )
1104 if ( nid[ iN ] == V000 ) {
1105 fV000.push_back( iF );
1106 } else if ( nid[ iN ] == V001 ) {
1107 fV001.push_back( iF );
1111 // find the bottom (Fxy0), the top (Fxy1) faces
1112 list<int>::iterator fIt1, fIt2, Fxy0Pos;
1113 for ( fIt1 = fV000.begin(); fIt1 != fV000.end(); fIt1++) {
1114 fIt2 = std::find( fV001.begin(), fV001.end(), *fIt1 );
1115 if ( fIt2 != fV001.end() ) { // *fIt1 is in the both lists
1116 fV001.erase( fIt2 ); // erase Fx0z or F0yz from fV001
1117 } else { // *fIt1 is in fV000 only
1118 Fxy0Pos = fIt1; // points to Fxy0
1122 Fxy1 = fV001.front();
1123 const SMDS_MeshNode** nn = vTool.GetNodes();
1125 // find bottom veritices, their order is that a face normal is external
1127 const int* nid = vTool.GetFaceNodesIndices( Fxy0 );
1128 for ( i = 0; i < 4; ++i )
1129 if ( nid[ i ] == V000 )
1131 for ( iN = 0; iN < 4; ++iN, ++i ) {
1132 if ( i == 4 ) i = 0;
1133 vFxy0[ iN ] = nid[ i ];
1135 // find top veritices, their order is that a face normal is external
1137 nid = vTool.GetFaceNodesIndices( Fxy1 );
1138 for ( i = 0; i < 4; ++i )
1139 if ( nid[ i ] == V001 )
1141 for ( iN = 0; iN < 4; ++iN, ++i ) {
1142 if ( i == 4 ) i = 0;
1143 vFxy1[ iN ] = nid[ i ];
1145 // find indices of the rest veritices
1153 // set points coordinates
1154 myPnt[ ID_V000 - 1 ] = gpXYZ( nn[ V000 ] );
1155 myPnt[ ID_V100 - 1 ] = gpXYZ( nn[ V100 ] );
1156 myPnt[ ID_V010 - 1 ] = gpXYZ( nn[ V010 ] );
1157 myPnt[ ID_V110 - 1 ] = gpXYZ( nn[ V110 ] );
1158 myPnt[ ID_V001 - 1 ] = gpXYZ( nn[ V001 ] );
1159 myPnt[ ID_V101 - 1 ] = gpXYZ( nn[ V101 ] );
1160 myPnt[ ID_V011 - 1 ] = gpXYZ( nn[ V011 ] );
1161 myPnt[ ID_V111 - 1 ] = gpXYZ( nn[ V111 ] );
1163 // fill theOrderedNodes
1164 theOrderedNodes.resize( 8 );
1165 theOrderedNodes[ 0 ] = nn[ V000 ];
1166 theOrderedNodes[ 1 ] = nn[ V100 ];
1167 theOrderedNodes[ 2 ] = nn[ V010 ];
1168 theOrderedNodes[ 3 ] = nn[ V110 ];
1169 theOrderedNodes[ 4 ] = nn[ V001 ];
1170 theOrderedNodes[ 5 ] = nn[ V101 ];
1171 theOrderedNodes[ 6 ] = nn[ V011 ];
1172 theOrderedNodes[ 7 ] = nn[ V111 ];
1175 vector< int > vertexVec;
1176 for ( iE = 0; iE < NbEdges(); ++iE ) {
1177 GetEdgeVertexIDs(( iE + ID_FirstE ), vertexVec );
1178 myEdge[ iE ].Set(( iE + ID_FirstE ),
1179 myPnt[ vertexVec[0] - 1 ],
1180 myPnt[ vertexVec[1] - 1 ]);
1183 // fill faces' corners
1184 for ( iF = ID_Fxy0; iF < ID_Shell; ++iF )
1186 TFace& tFace = myFace[ iF - ID_FirstF ];
1187 vector< int > edgeIdVec(4, -1);
1188 GetFaceEdgesIDs( iF, edgeIdVec );
1189 tFace.Set( iF, myEdge[ edgeIdVec [ 0 ] - ID_Ex00], myEdge[ edgeIdVec [ 1 ] - ID_Ex00]);
1195 //=======================================================================
1196 //function : LoadBlockShapes
1197 //purpose : Initialize block geometry with theShell,
1198 // add sub-shapes of theBlock to theShapeIDMap so that they get
1199 // IDs acoording to enum TShapeID
1200 //=======================================================================
1202 bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell,
1203 const TopoDS_Vertex& theVertex000,
1204 const TopoDS_Vertex& theVertex001,
1205 TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
1207 MESSAGE(" ::LoadBlockShapes()");
1208 return ( FindBlockShapes( theShell, theVertex000, theVertex001, theShapeIDMap ) &&
1209 LoadBlockShapes( theShapeIDMap ));
1212 //=======================================================================
1213 //function : LoadBlockShapes
1214 //purpose : add sub-shapes of theBlock to theShapeIDMap so that they get
1215 // IDs acoording to enum TShapeID
1216 //=======================================================================
1218 bool SMESH_Block::FindBlockShapes(const TopoDS_Shell& theShell,
1219 const TopoDS_Vertex& theVertex000,
1220 const TopoDS_Vertex& theVertex001,
1221 TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
1223 MESSAGE(" ::FindBlockShapes()");
1226 TopoDS_Shape V000, V100, V010, V110, V001, V101, V011, V111;
1228 TopoDS_Shape Ex00, Ex10, Ex01, Ex11;
1229 TopoDS_Shape E0y0, E1y0, E0y1, E1y1;
1230 TopoDS_Shape E00z, E10z, E01z, E11z;
1232 TopoDS_Shape Fxy0, Fx0z, F0yz, Fxy1, Fx1z, F1yz;
1234 // nb of faces bound to a vertex in TopTools_IndexedDataMapOfShapeListOfShape
1235 // filled by TopExp::MapShapesAndAncestors()
1236 const int NB_FACES_BY_VERTEX = 6;
1238 TopTools_IndexedDataMapOfShapeListOfShape vfMap;
1239 TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_FACE, vfMap );
1240 if ( vfMap.Extent() != 8 ) {
1241 MESSAGE(" Wrong nb of vertices in the block: " << vfMap.Extent() );
1245 V000 = theVertex000;
1246 V001 = theVertex001;
1248 if ( V000.IsNull() ) {
1249 // find vertex 000 - the one with smallest coordinates
1250 double minVal = DBL_MAX, minX, val;
1251 for ( int i = 1; i <= 8; i++ ) {
1252 const TopoDS_Vertex& v = TopoDS::Vertex( vfMap.FindKey( i ));
1253 gp_Pnt P = BRep_Tool::Pnt( v );
1254 val = P.X() + P.Y() + P.Z();
1255 if ( val < minVal || ( val == minVal && P.X() < minX )) {
1261 // find vertex 001 - the one on the most vertical edge passing through V000
1262 TopTools_IndexedDataMapOfShapeListOfShape veMap;
1263 TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_EDGE, veMap );
1264 gp_Vec dir001 = gp::DZ();
1265 gp_Pnt p000 = BRep_Tool::Pnt( TopoDS::Vertex( V000 ));
1266 double maxVal = -DBL_MAX;
1267 TopTools_ListIteratorOfListOfShape eIt ( veMap.FindFromKey( V000 ));
1268 for ( ; eIt.More(); eIt.Next() ) {
1269 const TopoDS_Edge& e = TopoDS::Edge( eIt.Value() );
1270 TopoDS_Vertex v = TopExp::FirstVertex( e );
1271 if ( v.IsSame( V000 ))
1272 v = TopExp::LastVertex( e );
1273 val = dir001 * gp_Vec( p000, BRep_Tool::Pnt( v )).Normalized();
1274 if ( val > maxVal ) {
1281 // find the bottom (Fxy0), Fx0z and F0yz faces
1283 const TopTools_ListOfShape& f000List = vfMap.FindFromKey( V000 );
1284 const TopTools_ListOfShape& f001List = vfMap.FindFromKey( V001 );
1285 if (f000List.Extent() != NB_FACES_BY_VERTEX ||
1286 f001List.Extent() != NB_FACES_BY_VERTEX ) {
1287 MESSAGE(" LoadBlockShapes() " << f000List.Extent() << " " << f001List.Extent());
1290 TopTools_ListIteratorOfListOfShape f001It, f000It ( f000List );
1291 int i, j, iFound1, iFound2;
1292 for ( j = 0; f000It.More(); f000It.Next(), j++ )
1294 if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1295 const TopoDS_Shape& F = f000It.Value();
1296 for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1297 if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1298 if ( F.IsSame( f001It.Value() ))
1301 if ( f001It.More() ) // Fx0z or F0yz found
1302 if ( Fx0z.IsNull() ) {
1309 else // F is the bottom face
1312 if ( Fxy0.IsNull() || Fx0z.IsNull() || F0yz.IsNull() ) {
1313 MESSAGE( Fxy0.IsNull() <<" "<< Fx0z.IsNull() <<" "<< F0yz.IsNull() );
1317 // choose the top face (Fxy1)
1318 for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1319 if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1320 if ( i != iFound1 && i != iFound2 )
1323 Fxy1 = f001It.Value();
1324 if ( Fxy1.IsNull() ) {
1325 MESSAGE(" LoadBlockShapes() error ");
1329 // find bottom edges and veritices
1330 list< TopoDS_Edge > eList;
1331 list< int > nbVertexInWires;
1332 GetOrderedEdges( TopoDS::Face( Fxy0 ), TopoDS::Vertex( V000 ), eList, nbVertexInWires );
1333 if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1334 MESSAGE(" LoadBlockShapes() error ");
1337 list< TopoDS_Edge >::iterator elIt = eList.begin();
1338 for ( i = 0; elIt != eList.end(); elIt++, i++ )
1340 case 0: E0y0 = *elIt; V010 = TopExp::LastVertex( *elIt, true ); break;
1341 case 1: Ex10 = *elIt; V110 = TopExp::LastVertex( *elIt, true ); break;
1342 case 2: E1y0 = *elIt; V100 = TopExp::LastVertex( *elIt, true ); break;
1343 case 3: Ex00 = *elIt; break;
1346 if ( i != 4 || E0y0.IsNull() || Ex10.IsNull() || E1y0.IsNull() || Ex00.IsNull() ) {
1347 MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1352 // find top edges and veritices
1354 GetOrderedEdges( TopoDS::Face( Fxy1 ), TopoDS::Vertex( V001 ), eList, nbVertexInWires );
1355 if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1356 MESSAGE(" LoadBlockShapes() error ");
1359 for ( i = 0, elIt = eList.begin(); elIt != eList.end(); elIt++, i++ )
1361 case 0: Ex01 = *elIt; V101 = TopExp::LastVertex( *elIt, true ); break;
1362 case 1: E1y1 = *elIt; V111 = TopExp::LastVertex( *elIt, true ); break;
1363 case 2: Ex11 = *elIt; V011 = TopExp::LastVertex( *elIt, true ); break;
1364 case 3: E0y1 = *elIt; break;
1367 if ( i != 4 || Ex01.IsNull() || E1y1.IsNull() || Ex11.IsNull() || E0y1.IsNull() ) {
1368 MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1372 // swap Fx0z and F0yz if necessary
1373 TopExp_Explorer exp( Fx0z, TopAbs_VERTEX );
1374 for ( ; exp.More(); exp.Next() ) // Fx0z shares V101 and V100
1375 if ( V101.IsSame( exp.Current() ) || V100.IsSame( exp.Current() ))
1376 break; // V101 or V100 found
1377 if ( !exp.More() ) { // not found
1378 std::swap( Fx0z, F0yz);
1381 // find Fx1z and F1yz faces
1382 const TopTools_ListOfShape& f111List = vfMap.FindFromKey( V111 );
1383 const TopTools_ListOfShape& f110List = vfMap.FindFromKey( V110 );
1384 if (f111List.Extent() != NB_FACES_BY_VERTEX ||
1385 f110List.Extent() != NB_FACES_BY_VERTEX ) {
1386 MESSAGE(" LoadBlockShapes() " << f111List.Extent() << " " << f110List.Extent());
1389 TopTools_ListIteratorOfListOfShape f111It, f110It ( f110List);
1390 for ( j = 0 ; f110It.More(); f110It.Next(), j++ ) {
1391 if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1392 const TopoDS_Shape& F = f110It.Value();
1393 for ( i = 0, f111It.Initialize( f111List ); f111It.More(); f111It.Next(), i++ ) {
1394 if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1395 if ( F.IsSame( f111It.Value() )) { // Fx1z or F1yz found
1396 if ( Fx1z.IsNull() )
1403 if ( Fx1z.IsNull() || F1yz.IsNull() ) {
1404 MESSAGE(" LoadBlockShapes() error ");
1408 // swap Fx1z and F1yz if necessary
1409 for ( exp.Init( Fx1z, TopAbs_VERTEX ); exp.More(); exp.Next() )
1410 if ( V010.IsSame( exp.Current() ) || V011.IsSame( exp.Current() ))
1412 if ( !exp.More() ) {
1413 std::swap( Fx1z, F1yz);
1416 // find vertical edges
1417 for ( exp.Init( Fx0z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1418 const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1419 const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1420 if ( vFirst.IsSame( V001 ))
1422 else if ( vFirst.IsSame( V100 ))
1425 if ( E00z.IsNull() || E10z.IsNull() ) {
1426 MESSAGE(" LoadBlockShapes() error ");
1429 for ( exp.Init( Fx1z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1430 const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1431 const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1432 if ( vFirst.IsSame( V111 ))
1434 else if ( vFirst.IsSame( V010 ))
1437 if ( E01z.IsNull() || E11z.IsNull() ) {
1438 MESSAGE(" LoadBlockShapes() error ");
1442 // load shapes in theShapeIDMap
1444 theShapeIDMap.Clear();
1446 theShapeIDMap.Add(V000.Oriented( TopAbs_FORWARD ));
1447 theShapeIDMap.Add(V100.Oriented( TopAbs_FORWARD ));
1448 theShapeIDMap.Add(V010.Oriented( TopAbs_FORWARD ));
1449 theShapeIDMap.Add(V110.Oriented( TopAbs_FORWARD ));
1450 theShapeIDMap.Add(V001.Oriented( TopAbs_FORWARD ));
1451 theShapeIDMap.Add(V101.Oriented( TopAbs_FORWARD ));
1452 theShapeIDMap.Add(V011.Oriented( TopAbs_FORWARD ));
1453 theShapeIDMap.Add(V111.Oriented( TopAbs_FORWARD ));
1455 theShapeIDMap.Add(Ex00);
1456 theShapeIDMap.Add(Ex10);
1457 theShapeIDMap.Add(Ex01);
1458 theShapeIDMap.Add(Ex11);
1460 theShapeIDMap.Add(E0y0);
1461 theShapeIDMap.Add(E1y0);
1462 theShapeIDMap.Add(E0y1);
1463 theShapeIDMap.Add(E1y1);
1465 theShapeIDMap.Add(E00z);
1466 theShapeIDMap.Add(E10z);
1467 theShapeIDMap.Add(E01z);
1468 theShapeIDMap.Add(E11z);
1470 theShapeIDMap.Add(Fxy0);
1471 theShapeIDMap.Add(Fxy1);
1472 theShapeIDMap.Add(Fx0z);
1473 theShapeIDMap.Add(Fx1z);
1474 theShapeIDMap.Add(F0yz);
1475 theShapeIDMap.Add(F1yz);
1477 theShapeIDMap.Add(theShell);
1482 //================================================================================
1484 * \brief Initialize block geometry with shapes from theShapeIDMap
1485 * \param theShapeIDMap - map of block subshapes
1486 * \retval bool - is a success
1488 //================================================================================
1490 bool SMESH_Block::LoadBlockShapes(const TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
1494 // store shapes geometry
1495 for ( int shapeID = 1; shapeID < theShapeIDMap.Extent(); shapeID++ )
1497 const TopoDS_Shape& S = theShapeIDMap( shapeID );
1498 switch ( S.ShapeType() )
1500 case TopAbs_VERTEX: {
1502 if ( !IsVertexID( ID_V111 )) return false;
1503 myPnt[ shapeID - ID_V000 ] = BRep_Tool::Pnt( TopoDS::Vertex( S )).XYZ();
1508 if ( !IsEdgeID( shapeID )) return false;
1509 const TopoDS_Edge& edge = TopoDS::Edge( S );
1510 TEdge& tEdge = myEdge[ shapeID - ID_FirstE ];
1512 new BRepAdaptor_Curve( edge ),
1513 IsForwardEdge( edge, theShapeIDMap ));
1518 if ( !LoadFace( TopoDS::Face( S ), shapeID, theShapeIDMap ))
1524 } // loop on shapes in theShapeIDMap
1529 //================================================================================
1531 * \brief Load face geometry
1532 * \param theFace - face
1533 * \param theFaceID - face in-block ID
1534 * \param theShapeIDMap - map of block subshapes
1535 * \retval bool - is a success
1537 * It is enough to compute params or coordinates on the face.
1538 * Face subshapes must be loaded into theShapeIDMap before
1540 //================================================================================
1542 bool SMESH_Block::LoadFace(const TopoDS_Face& theFace,
1543 const int theFaceID,
1544 const TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
1546 if ( !IsFaceID( theFaceID ) ) return false;
1548 Adaptor2d_Curve2d* c2d[4];
1550 vector< int > edgeIdVec;
1551 GetFaceEdgesIDs( theFaceID, edgeIdVec );
1552 for ( int iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges
1554 if ( edgeIdVec[ iE ] > theShapeIDMap.Extent() )
1556 const TopoDS_Edge& edge = TopoDS::Edge( theShapeIDMap( edgeIdVec[ iE ]));
1557 c2d[ iE ] = new BRepAdaptor_Curve2d( edge, theFace );
1558 isForward[ iE ] = IsForwardEdge( edge, theShapeIDMap );
1560 TFace& tFace = myFace[ theFaceID - ID_FirstF ];
1561 tFace.Set( theFaceID, new BRepAdaptor_Surface( theFace ), c2d, isForward );
1565 //================================================================================
1567 * \brief/ Insert theShape into theShapeIDMap with theShapeID
1568 * \param theShape - shape to insert
1569 * \param theShapeID - shape in-block ID
1570 * \param theShapeIDMap - map of block subshapes
1572 //================================================================================
1574 bool SMESH_Block::Insert(const TopoDS_Shape& theShape,
1575 const int theShapeID,
1576 TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
1578 if ( !theShape.IsNull() && theShapeID > 0 )
1580 if ( theShapeIDMap.Contains( theShape ))
1581 return ( theShapeIDMap.FindIndex( theShape ) == theShapeID );
1583 if ( theShapeID <= theShapeIDMap.Extent() ) {
1584 theShapeIDMap.Substitute( theShapeID, theShape );
1587 while ( theShapeIDMap.Extent() < theShapeID - 1 ) {
1588 TopoDS_Compound comp;
1589 BRep_Builder().MakeCompound( comp );
1590 theShapeIDMap.Add( comp );
1592 theShapeIDMap.Add( theShape );
1599 //=======================================================================
1600 //function : GetFaceEdgesIDs
1601 //purpose : return edges IDs in the order u0, u1, 0v, 1v
1602 // u0 means "|| u, v == 0"
1603 //=======================================================================
1605 void SMESH_Block::GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec )
1607 edgeVec.resize( 4 );
1610 edgeVec[ 0 ] = ID_Ex00;
1611 edgeVec[ 1 ] = ID_Ex10;
1612 edgeVec[ 2 ] = ID_E0y0;
1613 edgeVec[ 3 ] = ID_E1y0;
1616 edgeVec[ 0 ] = ID_Ex01;
1617 edgeVec[ 1 ] = ID_Ex11;
1618 edgeVec[ 2 ] = ID_E0y1;
1619 edgeVec[ 3 ] = ID_E1y1;
1622 edgeVec[ 0 ] = ID_Ex00;
1623 edgeVec[ 1 ] = ID_Ex01;
1624 edgeVec[ 2 ] = ID_E00z;
1625 edgeVec[ 3 ] = ID_E10z;
1628 edgeVec[ 0 ] = ID_Ex10;
1629 edgeVec[ 1 ] = ID_Ex11;
1630 edgeVec[ 2 ] = ID_E01z;
1631 edgeVec[ 3 ] = ID_E11z;
1634 edgeVec[ 0 ] = ID_E0y0;
1635 edgeVec[ 1 ] = ID_E0y1;
1636 edgeVec[ 2 ] = ID_E00z;
1637 edgeVec[ 3 ] = ID_E01z;
1640 edgeVec[ 0 ] = ID_E1y0;
1641 edgeVec[ 1 ] = ID_E1y1;
1642 edgeVec[ 2 ] = ID_E10z;
1643 edgeVec[ 3 ] = ID_E11z;
1646 MESSAGE(" GetFaceEdgesIDs(), wrong face ID: " << faceID );
1650 //=======================================================================
1651 //function : GetEdgeVertexIDs
1652 //purpose : return vertex IDs of an edge
1653 //=======================================================================
1655 void SMESH_Block::GetEdgeVertexIDs (const int edgeID, vector< int >& vertexVec )
1657 vertexVec.resize( 2 );
1661 vertexVec[ 0 ] = ID_V000;
1662 vertexVec[ 1 ] = ID_V100;
1665 vertexVec[ 0 ] = ID_V010;
1666 vertexVec[ 1 ] = ID_V110;
1669 vertexVec[ 0 ] = ID_V001;
1670 vertexVec[ 1 ] = ID_V101;
1673 vertexVec[ 0 ] = ID_V011;
1674 vertexVec[ 1 ] = ID_V111;
1678 vertexVec[ 0 ] = ID_V000;
1679 vertexVec[ 1 ] = ID_V010;
1682 vertexVec[ 0 ] = ID_V100;
1683 vertexVec[ 1 ] = ID_V110;
1686 vertexVec[ 0 ] = ID_V001;
1687 vertexVec[ 1 ] = ID_V011;
1690 vertexVec[ 0 ] = ID_V101;
1691 vertexVec[ 1 ] = ID_V111;
1695 vertexVec[ 0 ] = ID_V000;
1696 vertexVec[ 1 ] = ID_V001;
1699 vertexVec[ 0 ] = ID_V100;
1700 vertexVec[ 1 ] = ID_V101;
1703 vertexVec[ 0 ] = ID_V010;
1704 vertexVec[ 1 ] = ID_V011;
1707 vertexVec[ 0 ] = ID_V110;
1708 vertexVec[ 1 ] = ID_V111;
1711 vertexVec.resize(0);
1712 MESSAGE(" GetEdgeVertexIDs(), wrong edge ID: " << edgeID );