1 // Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License.
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 // File : SMESH_Pattern.hxx
21 // Created : Mon Aug 2 10:30:00 2004
22 // Author : Edward AGAPOV (eap)
24 #include "SMESH_Block.hxx"
26 #include <BRepAdaptor_Curve.hxx>
27 #include <BRepAdaptor_Curve2d.hxx>
28 #include <BRepAdaptor_Surface.hxx>
29 #include <BRepTools.hxx>
30 #include <BRepTools_WireExplorer.hxx>
31 #include <BRep_Builder.hxx>
32 #include <BRep_Tool.hxx>
33 #include <Extrema_ExtPC.hxx>
34 #include <Extrema_POnCurv.hxx>
35 #include <Geom2d_Curve.hxx>
36 #include <TopExp_Explorer.hxx>
37 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
38 #include <TopTools_ListIteratorOfListOfShape.hxx>
39 #include <TopTools_ListOfShape.hxx>
41 #include <TopoDS_Compound.hxx>
42 #include <TopoDS_Face.hxx>
43 #include <TopoDS_Iterator.hxx>
44 #include <TopoDS_Wire.hxx>
45 #include <gp_Trsf.hxx>
47 #include <math_FunctionSetRoot.hxx>
49 #include "SMDS_MeshNode.hxx"
50 #include "SMDS_MeshVolume.hxx"
51 #include "SMDS_VolumeTool.hxx"
52 #include "utilities.h"
60 //================================================================================
62 * \brief Set edge data
63 * \param edgeID - block subshape ID
64 * \param curve - edge geometry
65 * \param isForward - is curve orientation coincides with edge orientation in the block
67 //================================================================================
69 void SMESH_Block::TEdge::Set( const int edgeID, Adaptor3d_Curve* curve, const bool isForward )
71 myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID );
72 if ( myC3d ) delete myC3d;
74 myFirst = curve->FirstParameter();
75 myLast = curve->LastParameter();
77 std::swap( myFirst, myLast );
80 //================================================================================
82 * \brief Set coordinates of nodes at edge ends to work with mesh block
83 * \param edgeID - block subshape ID
84 * \param node1 - coordinates of node with lower ID
85 * \param node2 - coordinates of node with upper ID
87 //================================================================================
89 void SMESH_Block::TEdge::Set( const int edgeID, const gp_XYZ& node1, const gp_XYZ& node2 )
91 myCoordInd = SMESH_Block::GetCoordIndOnEdge( edgeID );
95 if ( myC3d ) delete myC3d;
99 //=======================================================================
100 //function : SMESH_Block::TEdge::GetU
102 //=======================================================================
104 double SMESH_Block::TEdge::GetU( const gp_XYZ& theParams ) const
106 double u = theParams.Coord( myCoordInd );
107 if ( !myC3d ) // if mesh block
109 return ( 1 - u ) * myFirst + u * myLast;
112 //=======================================================================
113 //function : SMESH_Block::TEdge::Point
115 //=======================================================================
117 gp_XYZ SMESH_Block::TEdge::Point( const gp_XYZ& theParams ) const
119 double u = GetU( theParams );
120 if ( myC3d ) return myC3d->Value( u ).XYZ();
122 return myNodes[0] * ( 1 - u ) + myNodes[1] * u;
125 //================================================================================
129 //================================================================================
131 SMESH_Block::TEdge::~TEdge()
133 if ( myC3d ) delete myC3d;
136 //================================================================================
138 * \brief Set face data
139 * \param faceID - block subshape ID
140 * \param S - face surface geometry
141 * \param c2d - 4 pcurves in the order as returned by GetFaceEdgesIDs(faceID)
142 * \param isForward - orientation of pcurves comparing with block edge direction
144 //================================================================================
146 void SMESH_Block::TFace::Set( const int faceID,
147 Adaptor3d_Surface* S,
148 Adaptor2d_Curve2d* c2D[4],
149 const bool isForward[4] )
151 if ( myS ) delete myS;
154 vector< int > edgeIdVec;
155 GetFaceEdgesIDs( faceID, edgeIdVec );
156 for ( int iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges
158 myCoordInd[ iE ] = GetCoordIndOnEdge( edgeIdVec[ iE ] );
159 if ( myC2d[ iE ]) delete myC2d[ iE ];
160 myC2d[ iE ] = c2D[ iE ];
161 myFirst[ iE ] = myC2d[ iE ]->FirstParameter();
162 myLast [ iE ] = myC2d[ iE ]->LastParameter();
163 if ( !isForward[ iE ])
164 std::swap( myFirst[ iE ], myLast[ iE ] );
167 myCorner[ 0 ] = myC2d[ 0 ]->Value( myFirst[0] ).XY();
168 myCorner[ 1 ] = myC2d[ 0 ]->Value( myLast[0] ).XY();
169 myCorner[ 2 ] = myC2d[ 1 ]->Value( myLast[1] ).XY();
170 myCorner[ 3 ] = myC2d[ 1 ]->Value( myFirst[1] ).XY();
173 //================================================================================
175 * \brief Set face data to work with mesh block
176 * \param faceID - block subshape ID
177 * \param edgeU0 - filled data of edge u0 = GetFaceEdgesIDs(faceID)[ 0 ]
178 * \param edgeU1 - filled data of edge u1 = GetFaceEdgesIDs(faceID)[ 1 ]
180 //================================================================================
182 void SMESH_Block::TFace::Set( const int faceID, const TEdge& edgeU0, const TEdge& edgeU1 )
184 vector< int > edgeIdVec;
185 GetFaceEdgesIDs( faceID, edgeIdVec );
186 myNodes[ 0 ] = edgeU0.NodeXYZ( 1 );
187 myNodes[ 1 ] = edgeU0.NodeXYZ( 0 );
188 myNodes[ 2 ] = edgeU1.NodeXYZ( 0 );
189 myNodes[ 3 ] = edgeU1.NodeXYZ( 1 );
190 myCoordInd[ 0 ] = GetCoordIndOnEdge( edgeIdVec[ 0 ] );
191 myCoordInd[ 1 ] = GetCoordIndOnEdge( edgeIdVec[ 1 ] );
192 myCoordInd[ 2 ] = GetCoordIndOnEdge( edgeIdVec[ 2 ] );
193 myCoordInd[ 3 ] = GetCoordIndOnEdge( edgeIdVec[ 3 ] );
194 if ( myS ) delete myS;
198 //================================================================================
202 //================================================================================
204 SMESH_Block::TFace::~TFace()
206 if ( myS ) delete myS;
207 for ( int i = 0 ; i < 4; ++i )
208 if ( myC2d[ i ]) delete myC2d[ i ];
211 //=======================================================================
212 //function : SMESH_Block::TFace::GetCoefs
213 //purpose : return coefficients for addition of [0-3]-th edge and vertex
214 //=======================================================================
216 void SMESH_Block::TFace::GetCoefs(int iE,
217 const gp_XYZ& theParams,
219 double& Vcoef ) const
221 double dU = theParams.Coord( GetUInd() );
222 double dV = theParams.Coord( GetVInd() );
225 Ecoef = ( 1 - dV ); // u0
226 Vcoef = ( 1 - dU ) * ( 1 - dV ); break; // 00
229 Vcoef = dU * ( 1 - dV ); break; // 10
231 Ecoef = ( 1 - dU ); // 0v
232 Vcoef = dU * dV ; break; // 11
235 Vcoef = ( 1 - dU ) * dV ; break; // 01
240 //=======================================================================
241 //function : SMESH_Block::TFace::GetUV
243 //=======================================================================
245 gp_XY SMESH_Block::TFace::GetUV( const gp_XYZ& theParams ) const
248 for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
250 double Ecoef = 0, Vcoef = 0;
251 GetCoefs( iE, theParams, Ecoef, Vcoef );
253 double u = theParams.Coord( myCoordInd[ iE ] );
254 u = ( 1 - u ) * myFirst[ iE ] + u * myLast[ iE ];
255 uv += Ecoef * myC2d[ iE ]->Value( u ).XY();
257 uv -= Vcoef * myCorner[ iE ];
262 //=======================================================================
263 //function : SMESH_Block::TFace::Point
265 //=======================================================================
267 gp_XYZ SMESH_Block::TFace::Point( const gp_XYZ& theParams ) const
270 if ( !myS ) // if mesh block
272 for ( int iE = 0; iE < 4; iE++ ) // loop on 4 edges
274 double Ecoef = 0, Vcoef = 0;
275 GetCoefs( iE, theParams, Ecoef, Vcoef );
277 double u = theParams.Coord( myCoordInd[ iE ] );
280 case 1: i1 = 3; i2 = 2; break;
281 case 2: i1 = 1; i2 = 2; break;
282 case 3: i1 = 0; i2 = 3; break;
284 p += Ecoef * ( myNodes[ i1 ] * ( 1 - u ) + myNodes[ i2 ] * u );
286 p -= Vcoef * myNodes[ iE ];
292 gp_XY uv = GetUV( theParams );
293 p = myS->Value( uv.X(), uv.Y() ).XYZ();
298 //=======================================================================
299 //function : GetShapeCoef
301 //=======================================================================
303 double* SMESH_Block::GetShapeCoef (const int theShapeID)
305 static double shapeCoef[][3] = {
306 // V000, V100, V010, V110
307 { -1,-1,-1 }, { 1,-1,-1 }, { -1, 1,-1 }, { 1, 1,-1 },
308 // V001, V101, V011, V111,
309 { -1,-1, 1 }, { 1,-1, 1 }, { -1, 1, 1 }, { 1, 1, 1 },
310 // Ex00, Ex10, Ex01, Ex11,
311 { 0,-1,-1 }, { 0, 1,-1 }, { 0,-1, 1 }, { 0, 1, 1 },
312 // E0y0, E1y0, E0y1, E1y1,
313 { -1, 0,-1 }, { 1, 0,-1 }, { -1, 0, 1 }, { 1, 0, 1 },
314 // E00z, E10z, E01z, E11z,
315 { -1,-1, 0 }, { 1,-1, 0 }, { -1, 1, 0 }, { 1, 1, 0 },
316 // Fxy0, Fxy1, Fx0z, Fx1z, F0yz, F1yz,
317 { 0, 0,-1 }, { 0, 0, 1 }, { 0,-1, 0 }, { 0, 1, 0 }, { -1, 0, 0 }, { 1, 0, 0 },
321 if ( theShapeID < ID_V000 || theShapeID > ID_F1yz )
322 return shapeCoef[ ID_Shell - 1 ];
324 return shapeCoef[ theShapeID - 1 ];
327 //=======================================================================
328 //function : ShellPoint
329 //purpose : return coordinates of a point in shell
330 //=======================================================================
332 bool SMESH_Block::ShellPoint( const gp_XYZ& theParams, gp_XYZ& thePoint ) const
334 thePoint.SetCoord( 0., 0., 0. );
335 for ( int shapeID = ID_V000; shapeID < ID_Shell; shapeID++ )
338 double* coefs = GetShapeCoef( shapeID );
340 for ( int iCoef = 0; iCoef < 3; iCoef++ ) {
341 if ( coefs[ iCoef ] != 0 ) {
342 if ( coefs[ iCoef ] < 0 )
343 k *= ( 1. - theParams.Coord( iCoef + 1 ));
345 k *= theParams.Coord( iCoef + 1 );
348 // add point on a shape
349 if ( fabs( k ) > DBL_MIN )
352 if ( shapeID < ID_Ex00 ) // vertex
353 VertexPoint( shapeID, Ps );
354 else if ( shapeID < ID_Fxy0 ) { // edge
355 EdgePoint( shapeID, theParams, Ps );
358 FacePoint( shapeID, theParams, Ps );
366 //=======================================================================
367 //function : ShellPoint
368 //purpose : computes coordinates of a point in shell by points on sub-shapes;
369 // thePointOnShape[ subShapeID ] must be a point on a subShape
370 //=======================================================================
372 bool SMESH_Block::ShellPoint(const gp_XYZ& theParams,
373 const vector<gp_XYZ>& thePointOnShape,
376 if ( thePointOnShape.size() < ID_F1yz )
379 double x = theParams.X(), y = theParams.Y(), z = theParams.Z();
380 double x1 = 1. - x, y1 = 1. - y, z1 = 1. - z;
381 const vector<gp_XYZ>& p = thePointOnShape;
384 x1 * p[ID_F0yz] + x * p[ID_F1yz]
385 + y1 * p[ID_Fx0z] + y * p[ID_Fx1z]
386 + z1 * p[ID_Fxy0] + z * p[ID_Fxy1]
387 + x1 * (y1 * (z1 * p[ID_V000] + z * p[ID_V001])
388 + y * (z1 * p[ID_V010] + z * p[ID_V011]))
389 + x * (y1 * (z1 * p[ID_V100] + z * p[ID_V101])
390 + y * (z1 * p[ID_V110] + z * p[ID_V111]));
392 x1 * (y1 * p[ID_E00z] + y * p[ID_E01z])
393 + x * (y1 * p[ID_E10z] + y * p[ID_E11z])
394 + y1 * (z1 * p[ID_Ex00] + z * p[ID_Ex01])
395 + y * (z1 * p[ID_Ex10] + z * p[ID_Ex11])
396 + z1 * (x1 * p[ID_E0y0] + x * p[ID_E1y0])
397 + z * (x1 * p[ID_E0y1] + x * p[ID_E1y1]);
402 //=======================================================================
403 //function : NbVariables
405 //=======================================================================
407 Standard_Integer SMESH_Block::NbVariables() const
412 //=======================================================================
413 //function : NbEquations
415 //=======================================================================
417 Standard_Integer SMESH_Block::NbEquations() const
422 //=======================================================================
425 //=======================================================================
427 Standard_Boolean SMESH_Block::Value(const math_Vector& theXYZ, math_Vector& theFxyz)
429 gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
430 if ( params.IsEqual( myParam, DBL_MIN )) { // same param
431 theFxyz( 1 ) = myValues[ 0 ];
434 ShellPoint( params, P );
435 gp_Vec dP( P - myPoint );
436 theFxyz(1) = SQRT_FUNC ? dP.SquareMagnitude() : dP.Magnitude();
441 //=======================================================================
442 //function : Derivatives
444 //=======================================================================
446 Standard_Boolean SMESH_Block::Derivatives(const math_Vector& XYZ,math_Matrix& Df)
448 MESSAGE( "SMESH_Block::Derivatives()");
450 return Values(XYZ,F,Df);
453 //=======================================================================
456 //=======================================================================
458 //#define DEBUG_PARAM_COMPUTE
460 Standard_Boolean SMESH_Block::Values(const math_Vector& theXYZ,
461 math_Vector& theFxyz,
464 // MESSAGE( endl<<"SMESH_Block::Values( "<<theXYZ(1)<<", "<<theXYZ(2)<<", "<<theXYZ(3)<<")");
466 gp_XYZ P, params( theXYZ(1), theXYZ(2), theXYZ(3) );
467 if ( params.IsEqual( myParam, DBL_MIN )) { // same param
468 theFxyz( 1 ) = myValues[ 0 ];
469 theDf( 1,1 ) = myValues[ 1 ];
470 theDf( 1,2 ) = myValues[ 2 ];
471 theDf( 1,3 ) = myValues[ 3 ];
474 #ifdef DEBUG_PARAM_COMPUTE
475 cout << "PARAM GUESS: " << params.X() << " "<< params.Y() << " "<< params.X() << endl;
477 ShellPoint( params, P );
478 //myNbIterations++; // how many time call ShellPoint()
480 gp_Vec dP( myPoint, P );
481 theFxyz(1) = SQRT_FUNC ? dP.SquareMagnitude() : dP.Magnitude();
482 if ( theFxyz(1) < 1e-6 ) {
491 if ( theFxyz(1) < myValues[0] ) // a better guess
493 // 3 partial derivatives
494 gp_Vec drv[ 3 ]; // where we move with a small step in each direction
495 for ( int iP = 1; iP <= 3; iP++ ) {
496 if ( iP - 1 == myFaceIndex ) {
497 drv[ iP - 1 ] = gp_Vec(0,0,0);
501 bool onEdge = ( theXYZ( iP ) + 0.001 > 1. );
503 params.SetCoord( iP, theXYZ( iP ) - 0.001 );
505 params.SetCoord( iP, theXYZ( iP ) + 0.001 );
506 ShellPoint( params, Pi );
507 params.SetCoord( iP, theXYZ( iP ) ); // restore params
508 gp_Vec dPi ( P, Pi );
509 if ( onEdge ) dPi *= -1.;
510 double mag = dPi.Magnitude();
515 for ( int iP = 0; iP < 3; iP++ ) {
516 theDf( 1, iP + 1 ) = dP * drv[iP];
517 // Distance from P to plane passing through myPoint and defined
518 // by the 2 other derivative directions:
519 // like IntAna_IntConicQuad::Perform (const gp_Lin& L, const gp_Pln& P)
520 // where L is (P -> myPoint), P is defined by the 2 other derivative direction
521 // int iPrev = ( iP ? iP - 1 : 2 );
522 // int iNext = ( iP == 2 ? 0 : iP + 1 );
523 // gp_Vec plnNorm = drv[ iPrev ].Crossed( drv [ iNext ] );
524 // double Direc = plnNorm * drv[ iP ];
525 // if ( Abs(Direc) <= DBL_MIN )
526 // theDf( 1, iP + 1 ) = dP * drv[ iP ];
528 // double Dis = plnNorm * P - plnNorm * myPoint;
529 // theDf( 1, iP + 1 ) = Dis/Direc;
532 #ifdef DEBUG_PARAM_COMPUTE
533 cout << "F = " << theFxyz(1) <<
534 " DRV: " << theDf(1,1) << " " << theDf(1,2) << " " << theDf(1,3) << endl;
536 //myNbIterations +=3; // how many time call ShellPoint()
538 // store better values
540 myValues[0]= theFxyz(1);
541 myValues[1]= theDf(1,1);
542 myValues[2]= theDf(1,2);
543 myValues[3]= theDf(1,3);
545 // SCRUTE( theFxyz(1) );
546 // SCRUTE( theDf( 1,1 ));
547 // SCRUTE( theDf( 1,2 ));
548 // SCRUTE( theDf( 1,3 ));
554 //=======================================================================
555 //function : ComputeParameters
556 //purpose : compute point parameters in the block
557 //=======================================================================
559 bool SMESH_Block::ComputeParameters(const gp_Pnt& thePoint,
561 const int theShapeID)
563 if ( VertexParameters( theShapeID, theParams ))
566 if ( IsEdgeID( theShapeID )) {
567 TEdge& e = myEdge[ theShapeID - ID_FirstE ];
568 Adaptor3d_Curve* curve = e.GetCurve();
569 Extrema_ExtPC anExtPC( thePoint, *curve, curve->FirstParameter(), curve->LastParameter() );
570 int i, nb = anExtPC.IsDone() ? anExtPC.NbExt() : 0;
571 for ( i = 1; i <= nb; i++ ) {
572 if ( anExtPC.IsMin( i ))
573 return EdgeParameters( theShapeID, anExtPC.Point( i ).Parameter(), theParams );
578 // MESSAGE( endl<<"SMESH_Block::ComputeParameters( "
579 // <<thePoint.X()<<", "<<thePoint.Y()<<", "<<thePoint.Z()<<")");
580 myPoint = thePoint.XYZ();
582 myParam.SetCoord( -1,-1,-1 );
585 const bool isOnFace = IsFaceID( theShapeID );
586 double * coef = GetShapeCoef( theShapeID );
589 math_Vector start( 1, 3, 0.0 );
590 if ( !myGridComputed )
592 // define the first guess by thePoint projection on lines
593 // connecting vertices
594 bool needGrid = false;
595 gp_XYZ par000( 0, 0, 0 ), par111( 1, 1, 1 );
596 double zero = DBL_MIN * DBL_MIN;
597 for ( int iEdge = 0, iParam = 1; iParam <= 3 && !needGrid; iParam++ )
599 if ( isOnFace && coef[ iParam - 1 ] != 0 ) {
603 for ( int iE = 0; iE < 4; iE++, iEdge++ ) { // loop on 4 parallel edges
604 gp_Pnt p0 = myEdge[ iEdge ].Point( par000 );
605 gp_Pnt p1 = myEdge[ iEdge ].Point( par111 );
606 gp_Vec v01( p0, p1 ), v0P( p0, thePoint );
607 double len2 = v01.SquareMagnitude();
610 par = v0P.Dot( v01 ) / len2;
611 if ( par < 0 || par > 1 ) { // projection falls out of line ends => needGrid
616 start( iParam ) += par;
618 start( iParam ) /= 4.;
621 // compute nodes of 3 x 3 x 3 grid
623 for ( double x = 0.25; x < 0.9; x += 0.25 )
624 for ( double y = 0.25; y < 0.9; y += 0.25 )
625 for ( double z = 0.25; z < 0.9; z += 0.25 ) {
626 TxyzPair & prmPtn = my3x3x3GridNodes[ iNode++ ];
627 prmPtn.first.SetCoord( x, y, z );
628 ShellPoint( prmPtn.first, prmPtn.second );
630 myGridComputed = true;
633 if ( myGridComputed ) {
634 double minDist = DBL_MAX;
635 gp_XYZ* bestParam = 0;
636 for ( int iNode = 0; iNode < 27; iNode++ ) {
637 TxyzPair & prmPtn = my3x3x3GridNodes[ iNode ];
638 double dist = ( thePoint.XYZ() - prmPtn.second ).SquareModulus();
639 if ( dist < minDist ) {
641 bestParam = & prmPtn.first;
644 start( 1 ) = bestParam->X();
645 start( 2 ) = bestParam->Y();
646 start( 3 ) = bestParam->Z();
651 // put a point on the face
652 for ( int iCoord = 0; iCoord < 3; iCoord++ )
653 if ( coef[ iCoord ] ) {
654 myFaceIndex = iCoord;
655 myFaceParam = ( coef[ myFaceIndex ] < 0.5 ) ? 0.0 : 1.0;
656 start( iCoord + 1 ) = myFaceParam;
659 math_Vector low ( 1, 3, 0.0 );
660 math_Vector up ( 1, 3, 1.0 );
661 math_Vector tol ( 1, 3, 1e-4 );
662 math_FunctionSetRoot paramSearch( *this, tol );
664 #ifdef DEBUG_PARAM_COMPUTE
665 cout << " #### POINT " <<thePoint.X()<<" "<<thePoint.Y()<<" "<<thePoint.Z()<<" ####"<< endl;
666 cout << " ** START ** " << start(1) << " " << start(2) << " " << start(3) << " " << endl;
669 while ( myValues[0] > 1e-1 && nbLoops++ < 10 ) {
670 paramSearch.Perform ( *static_cast<math_FunctionSetWithDerivatives*>(this),
672 if ( !paramSearch.IsDone() ) {
673 //MESSAGE( " !paramSearch.IsDone() " );
676 //MESSAGE( " NB ITERATIONS: " << paramSearch.NbIterations() );
678 start( 1 ) = myParam.X();
679 start( 2 ) = myParam.Y();
680 start( 3 ) = myParam.Z();
681 //MESSAGE( "Distance: " << ( SQRT_FUNC ? sqrt(myValues[0]) : myValues[0] ));
683 #ifdef DEBUG_PARAM_COMPUTE
684 cout << "-------SOLUTION-------: " << endl
685 << myParam.X() << " " << myParam.Y() << " " << myParam.Z() << endl
686 << " ------ DIST :" << myValues[0] << endl;
688 // MESSAGE( endl << myParam.X() << " " << myParam.Y() << " " << myParam.Z() << endl);
689 // mySumDist += myValues[0];
690 // MESSAGE( " TOTAL NB ITERATIONS: " << myNbIterations <<
691 // " DIST: " << ( SQRT_FUNC ? sqrt(mySumDist) : mySumDist ));
693 if ( myFaceIndex >= 0 )
694 myParam.SetCoord( myFaceIndex + 1, myFaceParam );
701 //=======================================================================
702 //function : VertexParameters
703 //purpose : return parameters of a vertex given by TShapeID
704 //=======================================================================
706 bool SMESH_Block::VertexParameters(const int theVertexID, gp_XYZ& theParams)
708 switch ( theVertexID ) {
709 case ID_V000: theParams.SetCoord(0., 0., 0.); return true;
710 case ID_V100: theParams.SetCoord(1., 0., 0.); return true;
711 case ID_V110: theParams.SetCoord(1., 1., 0.); return true;
712 case ID_V010: theParams.SetCoord(0., 1., 0.); return true;
718 //=======================================================================
719 //function : EdgeParameters
720 //purpose : return parameters of a point given by theU on edge
721 //=======================================================================
723 bool SMESH_Block::EdgeParameters(const int theEdgeID, const double theU, gp_XYZ& theParams)
725 if ( IsEdgeID( theEdgeID )) {
726 vector< int > vertexVec;
727 GetEdgeVertexIDs( theEdgeID, vertexVec );
728 VertexParameters( vertexVec[0], theParams );
729 TEdge& e = myEdge[ theEdgeID - ID_Ex00 ];
730 double param = ( theU - e.EndParam(0) ) / ( e.EndParam(1) - e.EndParam(0) );
731 theParams.SetCoord( e.CoordInd(), param );
737 //=======================================================================
738 //function : GetStateNumber
740 //=======================================================================
742 Standard_Integer SMESH_Block::GetStateNumber ()
744 // MESSAGE( endl<<"SMESH_Block::GetStateNumber( "<<myParam.X()<<", "<<
745 // myParam.Y()<<", "<<myParam.Z()<<") DISTANCE: " << myValues[0]);
746 return myValues[0] < 1e-1;
749 //=======================================================================
750 //function : DumpShapeID
751 //purpose : debug an id of a block sub-shape
752 //=======================================================================
754 #define CASEDUMP(id,strm) case id: strm << #id; break;
756 ostream& SMESH_Block::DumpShapeID (const int id, ostream& stream)
759 CASEDUMP( ID_V000, stream );
760 CASEDUMP( ID_V100, stream );
761 CASEDUMP( ID_V010, stream );
762 CASEDUMP( ID_V110, stream );
763 CASEDUMP( ID_V001, stream );
764 CASEDUMP( ID_V101, stream );
765 CASEDUMP( ID_V011, stream );
766 CASEDUMP( ID_V111, stream );
767 CASEDUMP( ID_Ex00, stream );
768 CASEDUMP( ID_Ex10, stream );
769 CASEDUMP( ID_Ex01, stream );
770 CASEDUMP( ID_Ex11, stream );
771 CASEDUMP( ID_E0y0, stream );
772 CASEDUMP( ID_E1y0, stream );
773 CASEDUMP( ID_E0y1, stream );
774 CASEDUMP( ID_E1y1, stream );
775 CASEDUMP( ID_E00z, stream );
776 CASEDUMP( ID_E10z, stream );
777 CASEDUMP( ID_E01z, stream );
778 CASEDUMP( ID_E11z, stream );
779 CASEDUMP( ID_Fxy0, stream );
780 CASEDUMP( ID_Fxy1, stream );
781 CASEDUMP( ID_Fx0z, stream );
782 CASEDUMP( ID_Fx1z, stream );
783 CASEDUMP( ID_F0yz, stream );
784 CASEDUMP( ID_F1yz, stream );
785 CASEDUMP( ID_Shell, stream );
786 default: stream << "ID_INVALID";
791 //=======================================================================
792 //function : GetShapeIDByParams
793 //purpose : define an id of the block sub-shape by normlized point coord
794 //=======================================================================
796 int SMESH_Block::GetShapeIDByParams ( const gp_XYZ& theCoord )
798 // id ( 0 - 26 ) computation:
800 // vertex ( 0 - 7 ) : id = 1*x + 2*y + 4*z
802 // edge || X ( 8 - 11 ) : id = 8 + 1*y + 2*z
803 // edge || Y ( 12 - 15 ): id = 1*x + 12 + 2*z
804 // edge || Z ( 16 - 19 ): id = 1*x + 2*y + 16
806 // face || XY ( 20 - 21 ): id = 8 + 12 + 1*z - 0
807 // face || XZ ( 22 - 23 ): id = 8 + 1*y + 16 - 2
808 // face || YZ ( 24 - 25 ): id = 1*x + 12 + 16 - 4
810 static int iAddBnd[] = { 1, 2, 4 };
811 static int iAddNotBnd[] = { 8, 12, 16 };
812 static int iFaceSubst[] = { 0, 2, 4 };
816 for ( int iCoord = 0; iCoord < 3; iCoord++ )
818 double val = theCoord.Coord( iCoord + 1 );
821 else if ( val == 1.0 )
822 id += iAddBnd[ iOnBoundary++ ];
824 id += iAddNotBnd[ iCoord ];
826 if ( iOnBoundary == 1 ) // face
827 id -= iFaceSubst[ (id - 20) / 4 ];
828 else if ( iOnBoundary == 0 ) // shell
831 if ( id > 26 || id < 0 ) {
832 MESSAGE( "GetShapeIDByParams() = " << id
833 <<" "<< theCoord.X() <<" "<< theCoord.Y() <<" "<< theCoord.Z() );
836 return id + 1; // shape ids start at 1
839 //=======================================================================
840 //function : GetOrderedEdges
841 //purpose : return nb wires and a list of oredered edges
842 //=======================================================================
844 int SMESH_Block::GetOrderedEdges (const TopoDS_Face& theFace,
845 TopoDS_Vertex theFirstVertex,
846 list< TopoDS_Edge >& theEdges,
847 list< int > & theNbVertexInWires)
849 // put wires in a list, so that an outer wire comes first
850 list<TopoDS_Wire> aWireList;
851 TopoDS_Wire anOuterWire = BRepTools::OuterWire( theFace );
852 aWireList.push_back( anOuterWire );
853 for ( TopoDS_Iterator wIt (theFace); wIt.More(); wIt.Next() )
854 if ( !anOuterWire.IsSame( wIt.Value() ))
855 aWireList.push_back( TopoDS::Wire( wIt.Value() ));
857 // loop on edges of wires
858 theNbVertexInWires.clear();
859 list<TopoDS_Wire>::iterator wlIt = aWireList.begin();
860 for ( ; wlIt != aWireList.end(); wlIt++ )
863 BRepTools_WireExplorer wExp( *wlIt, theFace );
864 for ( iE = 0; wExp.More(); wExp.Next(), iE++ )
866 TopoDS_Edge edge = wExp.Current();
867 edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
868 theEdges.push_back( edge );
870 theNbVertexInWires.push_back( iE );
872 if ( wlIt == aWireList.begin() && theEdges.size() > 1 ) { // the outer wire
873 // orient closed edges
874 list< TopoDS_Edge >::iterator eIt, eIt2;
875 for ( eIt = theEdges.begin(); eIt != theEdges.end(); eIt++ )
877 TopoDS_Edge& edge = *eIt;
878 if ( TopExp::FirstVertex( edge ).IsSame( TopExp::LastVertex( edge ) ))
881 bool isNext = ( eIt2 == theEdges.begin() );
882 TopoDS_Edge edge2 = isNext ? *(++eIt2) : *(--eIt2);
884 Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( edge, theFace, f1,l1 );
885 Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( edge2, theFace, f2,l2 );
886 gp_Pnt2d pf = c1->Value( edge.Orientation() == TopAbs_FORWARD ? f1 : l1 );
887 gp_Pnt2d pl = c1->Value( edge.Orientation() == TopAbs_FORWARD ? l1 : f1 );
888 bool isFirst = ( edge2.Orientation() == TopAbs_FORWARD ? isNext : !isNext );
889 gp_Pnt2d p2 = c2->Value( isFirst ? f2 : l2 );
890 isFirst = ( p2.SquareDistance( pf ) < p2.SquareDistance( pl ));
891 if ( isNext ? isFirst : !isFirst )
893 // to make a seam go first
894 if ( theFirstVertex.IsNull() )
895 theFirstVertex = TopExp::FirstVertex( edge, true );
898 // rotate theEdges until it begins from theFirstVertex
899 if ( ! theFirstVertex.IsNull() ) {
901 TopExp::Vertices( theEdges.front(), vv[0], vv[1], true );
902 // on closed face, make seam edge the first in the list
903 while ( !vv[0].IsSame( theFirstVertex ) || vv[0].IsSame( vv[1] ))
905 theEdges.splice(theEdges.end(), theEdges,
906 theEdges.begin(), ++theEdges.begin());
907 TopExp::Vertices( theEdges.front(), vv[0], vv[1], true );
908 if ( iE++ > theNbVertexInWires.back() ) {
910 gp_Pnt p = BRep_Tool::Pnt( theFirstVertex );
911 cout << " : Warning : vertex "<< theFirstVertex.TShape().operator->()
912 << " ( " << p.X() << " " << p.Y() << " " << p.Z() << " )"
913 << " not found in outer wire of face "<< theFace.TShape().operator->()
914 << " with vertices: " << endl;
915 wExp.Init( *wlIt, theFace );
916 for ( int i = 0; wExp.More(); wExp.Next(), i++ )
918 TopoDS_Edge edge = wExp.Current();
919 edge = TopoDS::Edge( edge.Oriented( wExp.Orientation() ));
920 TopoDS_Vertex v = TopExp::FirstVertex( edge, true );
921 gp_Pnt p = BRep_Tool::Pnt( v );
922 cout << i << " " << v.TShape().operator->() << " "
923 << p.X() << " " << p.Y() << " " << p.Z() << " " << endl;
926 break; // break infinite loop
933 return aWireList.size();
935 //================================================================================
937 * \brief Call it after geometry initialisation
939 //================================================================================
941 void SMESH_Block::init()
945 myGridComputed = false;
948 //=======================================================================
949 //function : LoadMeshBlock
950 //purpose : prepare to work with theVolume
951 //=======================================================================
953 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
955 bool SMESH_Block::LoadMeshBlock(const SMDS_MeshVolume* theVolume,
956 const int theNode000Index,
957 const int theNode001Index,
958 vector<const SMDS_MeshNode*>& theOrderedNodes)
960 MESSAGE(" ::LoadMeshBlock()");
963 SMDS_VolumeTool vTool;
964 if (!vTool.Set( theVolume ) || vTool.NbNodes() != 8 ||
965 !vTool.IsLinked( theNode000Index, theNode001Index )) {
966 MESSAGE(" Bad arguments ");
969 vTool.SetExternalNormal();
970 // In terms of indices used for access to nodes and faces in SMDS_VolumeTool:
971 int V000, V100, V010, V110, V001, V101, V011, V111; // 8 vertices
972 int Fxy0, Fxy1; // bottom and top faces
974 vector<int> vFxy0, vFxy1;
976 V000 = theNode000Index;
977 V001 = theNode001Index;
979 // get faces sharing V000 and V001
980 list<int> fV000, fV001;
982 for ( iF = 0; iF < vTool.NbFaces(); ++iF ) {
983 const int* nid = vTool.GetFaceNodesIndices( iF );
984 for ( iN = 0; iN < 4; ++iN )
985 if ( nid[ iN ] == V000 ) {
986 fV000.push_back( iF );
987 } else if ( nid[ iN ] == V001 ) {
988 fV001.push_back( iF );
992 // find the bottom (Fxy0), the top (Fxy1) faces
993 list<int>::iterator fIt1, fIt2, Fxy0Pos;
994 for ( fIt1 = fV000.begin(); fIt1 != fV000.end(); fIt1++) {
995 fIt2 = std::find( fV001.begin(), fV001.end(), *fIt1 );
996 if ( fIt2 != fV001.end() ) { // *fIt1 is in the both lists
997 fV001.erase( fIt2 ); // erase Fx0z or F0yz from fV001
998 } else { // *fIt1 is in fV000 only
999 Fxy0Pos = fIt1; // points to Fxy0
1003 Fxy1 = fV001.front();
1004 const SMDS_MeshNode** nn = vTool.GetNodes();
1006 // find bottom veritices, their order is that a face normal is external
1008 const int* nid = vTool.GetFaceNodesIndices( Fxy0 );
1009 for ( i = 0; i < 4; ++i )
1010 if ( nid[ i ] == V000 )
1012 for ( iN = 0; iN < 4; ++iN, ++i ) {
1013 if ( i == 4 ) i = 0;
1014 vFxy0[ iN ] = nid[ i ];
1016 // find top veritices, their order is that a face normal is external
1018 nid = vTool.GetFaceNodesIndices( Fxy1 );
1019 for ( i = 0; i < 4; ++i )
1020 if ( nid[ i ] == V001 )
1022 for ( iN = 0; iN < 4; ++iN, ++i ) {
1023 if ( i == 4 ) i = 0;
1024 vFxy1[ iN ] = nid[ i ];
1026 // find indices of the rest veritices
1034 // set points coordinates
1035 myPnt[ ID_V000 - 1 ] = gpXYZ( nn[ V000 ] );
1036 myPnt[ ID_V100 - 1 ] = gpXYZ( nn[ V100 ] );
1037 myPnt[ ID_V010 - 1 ] = gpXYZ( nn[ V010 ] );
1038 myPnt[ ID_V110 - 1 ] = gpXYZ( nn[ V110 ] );
1039 myPnt[ ID_V001 - 1 ] = gpXYZ( nn[ V001 ] );
1040 myPnt[ ID_V101 - 1 ] = gpXYZ( nn[ V101 ] );
1041 myPnt[ ID_V011 - 1 ] = gpXYZ( nn[ V011 ] );
1042 myPnt[ ID_V111 - 1 ] = gpXYZ( nn[ V111 ] );
1044 // fill theOrderedNodes
1045 theOrderedNodes.resize( 8 );
1046 theOrderedNodes[ 0 ] = nn[ V000 ];
1047 theOrderedNodes[ 1 ] = nn[ V100 ];
1048 theOrderedNodes[ 2 ] = nn[ V010 ];
1049 theOrderedNodes[ 3 ] = nn[ V110 ];
1050 theOrderedNodes[ 4 ] = nn[ V001 ];
1051 theOrderedNodes[ 5 ] = nn[ V101 ];
1052 theOrderedNodes[ 6 ] = nn[ V011 ];
1053 theOrderedNodes[ 7 ] = nn[ V111 ];
1056 vector< int > vertexVec;
1057 for ( iE = 0; iE < NbEdges(); ++iE ) {
1058 GetEdgeVertexIDs(( iE + ID_FirstE ), vertexVec );
1059 myEdge[ iE ].Set(( iE + ID_FirstE ),
1060 myPnt[ vertexVec[0] - 1 ],
1061 myPnt[ vertexVec[1] - 1 ]);
1064 // fill faces' corners
1065 for ( iF = ID_Fxy0; iF < ID_Shell; ++iF )
1067 TFace& tFace = myFace[ iF - ID_FirstF ];
1068 vector< int > edgeIdVec(4, -1);
1069 GetFaceEdgesIDs( iF, edgeIdVec );
1070 tFace.Set( iF, myEdge[ edgeIdVec [ 0 ]], myEdge[ edgeIdVec [ 1 ]]);
1076 //=======================================================================
1077 //function : LoadBlockShapes
1078 //purpose : Initialize block geometry with theShell,
1079 // add sub-shapes of theBlock to theShapeIDMap so that they get
1080 // IDs acoording to enum TShapeID
1081 //=======================================================================
1083 bool SMESH_Block::LoadBlockShapes(const TopoDS_Shell& theShell,
1084 const TopoDS_Vertex& theVertex000,
1085 const TopoDS_Vertex& theVertex001,
1086 TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
1088 MESSAGE(" ::LoadBlockShapes()");
1089 return ( FindBlockShapes( theShell, theVertex000, theVertex001, theShapeIDMap ) &&
1090 LoadBlockShapes( theShapeIDMap ));
1093 //=======================================================================
1094 //function : LoadBlockShapes
1095 //purpose : add sub-shapes of theBlock to theShapeIDMap so that they get
1096 // IDs acoording to enum TShapeID
1097 //=======================================================================
1099 bool SMESH_Block::FindBlockShapes(const TopoDS_Shell& theShell,
1100 const TopoDS_Vertex& theVertex000,
1101 const TopoDS_Vertex& theVertex001,
1102 TopTools_IndexedMapOfOrientedShape& theShapeIDMap )
1104 MESSAGE(" ::FindBlockShapes()");
1107 TopoDS_Shape V000, V100, V010, V110, V001, V101, V011, V111;
1109 TopoDS_Shape Ex00, Ex10, Ex01, Ex11;
1110 TopoDS_Shape E0y0, E1y0, E0y1, E1y1;
1111 TopoDS_Shape E00z, E10z, E01z, E11z;
1113 TopoDS_Shape Fxy0, Fx0z, F0yz, Fxy1, Fx1z, F1yz;
1115 // nb of faces bound to a vertex in TopTools_IndexedDataMapOfShapeListOfShape
1116 // filled by TopExp::MapShapesAndAncestors()
1117 const int NB_FACES_BY_VERTEX = 6;
1119 TopTools_IndexedDataMapOfShapeListOfShape vfMap;
1120 TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_FACE, vfMap );
1121 if ( vfMap.Extent() != 8 ) {
1122 MESSAGE(" Wrong nb of vertices in the block: " << vfMap.Extent() );
1126 V000 = theVertex000;
1127 V001 = theVertex001;
1129 if ( V000.IsNull() ) {
1130 // find vertex 000 - the one with smallest coordinates
1131 double minVal = DBL_MAX, minX, val;
1132 for ( int i = 1; i <= 8; i++ ) {
1133 const TopoDS_Vertex& v = TopoDS::Vertex( vfMap.FindKey( i ));
1134 gp_Pnt P = BRep_Tool::Pnt( v );
1135 val = P.X() + P.Y() + P.Z();
1136 if ( val < minVal || ( val == minVal && P.X() < minX )) {
1142 // find vertex 001 - the one on the most vertical edge passing through V000
1143 TopTools_IndexedDataMapOfShapeListOfShape veMap;
1144 TopExp::MapShapesAndAncestors( theShell, TopAbs_VERTEX, TopAbs_EDGE, veMap );
1145 gp_Vec dir001 = gp::DZ();
1146 gp_Pnt p000 = BRep_Tool::Pnt( TopoDS::Vertex( V000 ));
1147 double maxVal = -DBL_MAX;
1148 TopTools_ListIteratorOfListOfShape eIt ( veMap.FindFromKey( V000 ));
1149 for ( ; eIt.More(); eIt.Next() ) {
1150 const TopoDS_Edge& e = TopoDS::Edge( eIt.Value() );
1151 TopoDS_Vertex v = TopExp::FirstVertex( e );
1152 if ( v.IsSame( V000 ))
1153 v = TopExp::LastVertex( e );
1154 val = dir001 * gp_Vec( p000, BRep_Tool::Pnt( v )).Normalized();
1155 if ( val > maxVal ) {
1162 // find the bottom (Fxy0), Fx0z and F0yz faces
1164 const TopTools_ListOfShape& f000List = vfMap.FindFromKey( V000 );
1165 const TopTools_ListOfShape& f001List = vfMap.FindFromKey( V001 );
1166 if (f000List.Extent() != NB_FACES_BY_VERTEX ||
1167 f001List.Extent() != NB_FACES_BY_VERTEX ) {
1168 MESSAGE(" LoadBlockShapes() " << f000List.Extent() << " " << f001List.Extent());
1171 TopTools_ListIteratorOfListOfShape f001It, f000It ( f000List );
1172 int i, j, iFound1, iFound2;
1173 for ( j = 0; f000It.More(); f000It.Next(), j++ )
1175 if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1176 const TopoDS_Shape& F = f000It.Value();
1177 for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1178 if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1179 if ( F.IsSame( f001It.Value() ))
1182 if ( f001It.More() ) // Fx0z or F0yz found
1183 if ( Fx0z.IsNull() ) {
1190 else // F is the bottom face
1193 if ( Fxy0.IsNull() || Fx0z.IsNull() || F0yz.IsNull() ) {
1194 MESSAGE( Fxy0.IsNull() <<" "<< Fx0z.IsNull() <<" "<< F0yz.IsNull() );
1198 // choose the top face (Fxy1)
1199 for ( i = 0, f001It.Initialize( f001List ); f001It.More(); f001It.Next(), i++ ) {
1200 if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1201 if ( i != iFound1 && i != iFound2 )
1204 Fxy1 = f001It.Value();
1205 if ( Fxy1.IsNull() ) {
1206 MESSAGE(" LoadBlockShapes() error ");
1210 // find bottom edges and veritices
1211 list< TopoDS_Edge > eList;
1212 list< int > nbVertexInWires;
1213 GetOrderedEdges( TopoDS::Face( Fxy0 ), TopoDS::Vertex( V000 ), eList, nbVertexInWires );
1214 if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1215 MESSAGE(" LoadBlockShapes() error ");
1218 list< TopoDS_Edge >::iterator elIt = eList.begin();
1219 for ( i = 0; elIt != eList.end(); elIt++, i++ )
1221 case 0: E0y0 = *elIt; V010 = TopExp::LastVertex( *elIt, true ); break;
1222 case 1: Ex10 = *elIt; V110 = TopExp::LastVertex( *elIt, true ); break;
1223 case 2: E1y0 = *elIt; V100 = TopExp::LastVertex( *elIt, true ); break;
1224 case 3: Ex00 = *elIt; break;
1227 if ( i != 4 || E0y0.IsNull() || Ex10.IsNull() || E1y0.IsNull() || Ex00.IsNull() ) {
1228 MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1233 // find top edges and veritices
1235 GetOrderedEdges( TopoDS::Face( Fxy1 ), TopoDS::Vertex( V001 ), eList, nbVertexInWires );
1236 if ( nbVertexInWires.size() != 1 || nbVertexInWires.front() != 4 ) {
1237 MESSAGE(" LoadBlockShapes() error ");
1240 for ( i = 0, elIt = eList.begin(); elIt != eList.end(); elIt++, i++ )
1242 case 0: Ex01 = *elIt; V101 = TopExp::LastVertex( *elIt, true ); break;
1243 case 1: E1y1 = *elIt; V111 = TopExp::LastVertex( *elIt, true ); break;
1244 case 2: Ex11 = *elIt; V011 = TopExp::LastVertex( *elIt, true ); break;
1245 case 3: E0y1 = *elIt; break;
1248 if ( i != 4 || Ex01.IsNull() || E1y1.IsNull() || Ex11.IsNull() || E0y1.IsNull() ) {
1249 MESSAGE(" LoadBlockShapes() error, eList.size()=" << eList.size());
1253 // swap Fx0z and F0yz if necessary
1254 TopExp_Explorer exp( Fx0z, TopAbs_VERTEX );
1255 for ( ; exp.More(); exp.Next() ) // Fx0z shares V101 and V100
1256 if ( V101.IsSame( exp.Current() ) || V100.IsSame( exp.Current() ))
1257 break; // V101 or V100 found
1258 if ( !exp.More() ) { // not found
1259 std::swap( Fx0z, F0yz);
1262 // find Fx1z and F1yz faces
1263 const TopTools_ListOfShape& f111List = vfMap.FindFromKey( V111 );
1264 const TopTools_ListOfShape& f110List = vfMap.FindFromKey( V110 );
1265 if (f111List.Extent() != NB_FACES_BY_VERTEX ||
1266 f110List.Extent() != NB_FACES_BY_VERTEX ) {
1267 MESSAGE(" LoadBlockShapes() " << f111List.Extent() << " " << f110List.Extent());
1270 TopTools_ListIteratorOfListOfShape f111It, f110It ( f110List);
1271 for ( j = 0 ; f110It.More(); f110It.Next(), j++ ) {
1272 if ( NB_FACES_BY_VERTEX == 6 && j % 2 ) continue; // each face encounters twice
1273 const TopoDS_Shape& F = f110It.Value();
1274 for ( i = 0, f111It.Initialize( f111List ); f111It.More(); f111It.Next(), i++ ) {
1275 if ( NB_FACES_BY_VERTEX == 6 && i % 2 ) continue; // each face encounters twice
1276 if ( F.IsSame( f111It.Value() )) { // Fx1z or F1yz found
1277 if ( Fx1z.IsNull() )
1284 if ( Fx1z.IsNull() || F1yz.IsNull() ) {
1285 MESSAGE(" LoadBlockShapes() error ");
1289 // swap Fx1z and F1yz if necessary
1290 for ( exp.Init( Fx1z, TopAbs_VERTEX ); exp.More(); exp.Next() )
1291 if ( V010.IsSame( exp.Current() ) || V011.IsSame( exp.Current() ))
1293 if ( !exp.More() ) {
1294 std::swap( Fx1z, F1yz);
1297 // find vertical edges
1298 for ( exp.Init( Fx0z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1299 const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1300 const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1301 if ( vFirst.IsSame( V001 ))
1303 else if ( vFirst.IsSame( V100 ))
1306 if ( E00z.IsNull() || E10z.IsNull() ) {
1307 MESSAGE(" LoadBlockShapes() error ");
1310 for ( exp.Init( Fx1z, TopAbs_EDGE ); exp.More(); exp.Next() ) {
1311 const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
1312 const TopoDS_Shape& vFirst = TopExp::FirstVertex( edge, true );
1313 if ( vFirst.IsSame( V111 ))
1315 else if ( vFirst.IsSame( V010 ))
1318 if ( E01z.IsNull() || E11z.IsNull() ) {
1319 MESSAGE(" LoadBlockShapes() error ");
1323 // load shapes in theShapeIDMap
1325 theShapeIDMap.Clear();
1327 theShapeIDMap.Add(V000.Oriented( TopAbs_FORWARD ));
1328 theShapeIDMap.Add(V100.Oriented( TopAbs_FORWARD ));
1329 theShapeIDMap.Add(V010.Oriented( TopAbs_FORWARD ));
1330 theShapeIDMap.Add(V110.Oriented( TopAbs_FORWARD ));
1331 theShapeIDMap.Add(V001.Oriented( TopAbs_FORWARD ));
1332 theShapeIDMap.Add(V101.Oriented( TopAbs_FORWARD ));
1333 theShapeIDMap.Add(V011.Oriented( TopAbs_FORWARD ));
1334 theShapeIDMap.Add(V111.Oriented( TopAbs_FORWARD ));
1336 theShapeIDMap.Add(Ex00);
1337 theShapeIDMap.Add(Ex10);
1338 theShapeIDMap.Add(Ex01);
1339 theShapeIDMap.Add(Ex11);
1341 theShapeIDMap.Add(E0y0);
1342 theShapeIDMap.Add(E1y0);
1343 theShapeIDMap.Add(E0y1);
1344 theShapeIDMap.Add(E1y1);
1346 theShapeIDMap.Add(E00z);
1347 theShapeIDMap.Add(E10z);
1348 theShapeIDMap.Add(E01z);
1349 theShapeIDMap.Add(E11z);
1351 theShapeIDMap.Add(Fxy0);
1352 theShapeIDMap.Add(Fxy1);
1353 theShapeIDMap.Add(Fx0z);
1354 theShapeIDMap.Add(Fx1z);
1355 theShapeIDMap.Add(F0yz);
1356 theShapeIDMap.Add(F1yz);
1358 theShapeIDMap.Add(theShell);
1363 //================================================================================
1365 * \brief Initialize block geometry with shapes from theShapeIDMap
1366 * \param theShapeIDMap - map of block subshapes
1367 * \retval bool - is a success
1369 //================================================================================
1371 bool SMESH_Block::LoadBlockShapes(const TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
1375 // store shapes geometry
1376 for ( int shapeID = 1; shapeID < theShapeIDMap.Extent(); shapeID++ )
1378 const TopoDS_Shape& S = theShapeIDMap( shapeID );
1379 switch ( S.ShapeType() )
1381 case TopAbs_VERTEX: {
1383 if ( !IsVertexID( ID_V111 )) return false;
1384 myPnt[ shapeID - ID_V000 ] = BRep_Tool::Pnt( TopoDS::Vertex( S )).XYZ();
1389 if ( !IsEdgeID( shapeID )) return false;
1390 const TopoDS_Edge& edge = TopoDS::Edge( S );
1391 TEdge& tEdge = myEdge[ shapeID - ID_FirstE ];
1393 new BRepAdaptor_Curve( edge ),
1394 IsForwardEdge( edge, theShapeIDMap ));
1399 if ( !LoadFace( TopoDS::Face( S ), shapeID, theShapeIDMap ))
1405 } // loop on shapes in theShapeIDMap
1410 //================================================================================
1412 * \brief Load face geometry
1413 * \param theFace - face
1414 * \param theFaceID - face in-block ID
1415 * \param theShapeIDMap - map of block subshapes
1416 * \retval bool - is a success
1418 * It is enough to compute params or coordinates on the face.
1419 * Face subshapes must be loaded into theShapeIDMap before
1421 //================================================================================
1423 bool SMESH_Block::LoadFace(const TopoDS_Face& theFace,
1424 const int theFaceID,
1425 const TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
1427 if ( !IsFaceID( theFaceID ) ) return false;
1429 Adaptor2d_Curve2d* c2d[4];
1431 vector< int > edgeIdVec;
1432 GetFaceEdgesIDs( theFaceID, edgeIdVec );
1433 for ( int iE = 0; iE < edgeIdVec.size(); iE++ ) // loop on 4 edges
1435 if ( edgeIdVec[ iE ] > theShapeIDMap.Extent() )
1437 const TopoDS_Edge& edge = TopoDS::Edge( theShapeIDMap( edgeIdVec[ iE ]));
1438 c2d[ iE ] = new BRepAdaptor_Curve2d( edge, theFace );
1439 isForward[ iE ] = IsForwardEdge( edge, theShapeIDMap );
1441 TFace& tFace = myFace[ theFaceID - ID_FirstF ];
1442 tFace.Set( theFaceID, new BRepAdaptor_Surface( theFace ), c2d, isForward );
1446 //================================================================================
1448 * \brief/ Insert theShape into theShapeIDMap with theShapeID
1449 * \param theShape - shape to insert
1450 * \param theShapeID - shape in-block ID
1451 * \param theShapeIDMap - map of block subshapes
1453 //================================================================================
1455 bool SMESH_Block::Insert(const TopoDS_Shape& theShape,
1456 const int theShapeID,
1457 TopTools_IndexedMapOfOrientedShape& theShapeIDMap)
1459 if ( !theShape.IsNull() && theShapeID > 0 )
1461 if ( theShapeIDMap.Contains( theShape ))
1462 return ( theShapeIDMap.FindIndex( theShape ) == theShapeID );
1464 if ( theShapeID <= theShapeIDMap.Extent() ) {
1465 theShapeIDMap.Substitute( theShapeID, theShape );
1468 while ( theShapeIDMap.Extent() < theShapeID - 1 ) {
1469 TopoDS_Compound comp;
1470 BRep_Builder().MakeCompound( comp );
1471 theShapeIDMap.Add( comp );
1473 theShapeIDMap.Add( theShape );
1480 //=======================================================================
1481 //function : GetFaceEdgesIDs
1482 //purpose : return edges IDs in the order u0, u1, 0v, 1v
1483 // u0 means "|| u, v == 0"
1484 //=======================================================================
1486 void SMESH_Block::GetFaceEdgesIDs (const int faceID, vector< int >& edgeVec )
1488 edgeVec.resize( 4 );
1491 edgeVec[ 0 ] = ID_Ex00;
1492 edgeVec[ 1 ] = ID_Ex10;
1493 edgeVec[ 2 ] = ID_E0y0;
1494 edgeVec[ 3 ] = ID_E1y0;
1497 edgeVec[ 0 ] = ID_Ex01;
1498 edgeVec[ 1 ] = ID_Ex11;
1499 edgeVec[ 2 ] = ID_E0y1;
1500 edgeVec[ 3 ] = ID_E1y1;
1503 edgeVec[ 0 ] = ID_Ex00;
1504 edgeVec[ 1 ] = ID_Ex01;
1505 edgeVec[ 2 ] = ID_E00z;
1506 edgeVec[ 3 ] = ID_E10z;
1509 edgeVec[ 0 ] = ID_Ex10;
1510 edgeVec[ 1 ] = ID_Ex11;
1511 edgeVec[ 2 ] = ID_E01z;
1512 edgeVec[ 3 ] = ID_E11z;
1515 edgeVec[ 0 ] = ID_E0y0;
1516 edgeVec[ 1 ] = ID_E0y1;
1517 edgeVec[ 2 ] = ID_E00z;
1518 edgeVec[ 3 ] = ID_E01z;
1521 edgeVec[ 0 ] = ID_E1y0;
1522 edgeVec[ 1 ] = ID_E1y1;
1523 edgeVec[ 2 ] = ID_E10z;
1524 edgeVec[ 3 ] = ID_E11z;
1527 MESSAGE(" GetFaceEdgesIDs(), wrong face ID: " << faceID );
1531 //=======================================================================
1532 //function : GetEdgeVertexIDs
1533 //purpose : return vertex IDs of an edge
1534 //=======================================================================
1536 void SMESH_Block::GetEdgeVertexIDs (const int edgeID, vector< int >& vertexVec )
1538 vertexVec.resize( 2 );
1542 vertexVec[ 0 ] = ID_V000;
1543 vertexVec[ 1 ] = ID_V100;
1546 vertexVec[ 0 ] = ID_V010;
1547 vertexVec[ 1 ] = ID_V110;
1550 vertexVec[ 0 ] = ID_V001;
1551 vertexVec[ 1 ] = ID_V101;
1554 vertexVec[ 0 ] = ID_V011;
1555 vertexVec[ 1 ] = ID_V111;
1559 vertexVec[ 0 ] = ID_V000;
1560 vertexVec[ 1 ] = ID_V010;
1563 vertexVec[ 0 ] = ID_V100;
1564 vertexVec[ 1 ] = ID_V110;
1567 vertexVec[ 0 ] = ID_V001;
1568 vertexVec[ 1 ] = ID_V011;
1571 vertexVec[ 0 ] = ID_V101;
1572 vertexVec[ 1 ] = ID_V111;
1576 vertexVec[ 0 ] = ID_V000;
1577 vertexVec[ 1 ] = ID_V001;
1580 vertexVec[ 0 ] = ID_V100;
1581 vertexVec[ 1 ] = ID_V101;
1584 vertexVec[ 0 ] = ID_V010;
1585 vertexVec[ 1 ] = ID_V011;
1588 vertexVec[ 0 ] = ID_V110;
1589 vertexVec[ 1 ] = ID_V111;
1592 vertexVec.resize(0);
1593 MESSAGE(" GetEdgeVertexIDs(), wrong edge ID: " << edgeID );