+ // columns for vertices
+ const SMDS_MeshNode* n0 = cols->begin()->second.front();
+ id = n0->getshapeId();
+ myShapeIndex2ColumnMap[ id ] = make_pair( cols, isForward );
+
+ const SMDS_MeshNode* n1 = cols->rbegin()->second.front();
+ id = n1->getshapeId();
+ myShapeIndex2ColumnMap[ id ] = make_pair( cols, !isForward );
+ }
+ }
+
+// #define SHOWYXZ(msg, xyz) { gp_Pnt p(xyz); cout << msg << " ("<< p.X() << "; " <<p.Y() << "; " <<p.Z() << ") " <<endl; }
+
+// double _u[]={ 0.1, 0.1, 0.9, 0.9 };
+// double _v[]={ 0.1, 0.9, 0.1, 0.9 };
+// for ( int z = 0; z < 2; ++z )
+// for ( int i = 0; i < 4; ++i )
+// {
+// //gp_XYZ testPar(0.25, 0.25, 0), testCoord;
+// int iFace = (z ? ID_TOP_FACE : ID_BOT_FACE);
+// gp_XYZ testPar(_u[i], _v[i], z), testCoord;
+// if ( !FacePoint( iFace, testPar, testCoord ))
+// RETURN_BAD_RESULT("TEST FacePoint() FAILED");
+// SHOWYXZ("IN TEST PARAM" , testPar);
+// SHOWYXZ("OUT TEST CORD" , testCoord);
+// if ( !ComputeParameters( testCoord, testPar , iFace))
+// RETURN_BAD_RESULT("TEST ComputeParameters() FAILED");
+// SHOWYXZ("OUT TEST PARAM" , testPar);
+// }
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Return pointer to column of nodes
+ * \param node - bottom node from which the returned column goes up
+ * \retval const TNodeColumn* - the found column
+ */
+//================================================================================
+
+const TNodeColumn* StdMeshers_PrismAsBlock::GetNodeColumn(const SMDS_MeshNode* node) const
+{
+ int sID = node->getshapeId();
+
+ map<int, pair< TParam2ColumnMap*, bool > >::const_iterator col_frw =
+ myShapeIndex2ColumnMap.find( sID );
+ if ( col_frw != myShapeIndex2ColumnMap.end() ) {
+ const TParam2ColumnMap* cols = col_frw->second.first;
+ TParam2ColumnIt u_col = cols->begin();
+ for ( ; u_col != cols->end(); ++u_col )
+ if ( u_col->second[ 0 ] == node )
+ return & u_col->second;
+ }
+ return 0;
+}
+
+//=======================================================================
+//function : GetLayersTransformation
+//purpose : Return transformations to get coordinates of nodes of each layer
+// by nodes of the bottom. Layer is a set of nodes at a certain step
+// from bottom to top.
+// Transformation to get top node from bottom ones is computed
+// only if the top FACE is not meshed.
+//=======================================================================
+
+bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector<gp_Trsf> & trsf,
+ const Prism_3D::TPrismTopo& prism) const
+{
+ const bool itTopMeshed = !SubMesh( ID_BOT_FACE )->IsEmpty();
+ const size_t zSize = VerticalSize();
+ if ( zSize < 3 && !itTopMeshed ) return true;
+ trsf.resize( zSize - 1 );
+
+ // Select some node columns by which we will define coordinate system of layers
+
+ vector< const TNodeColumn* > columns;
+ {
+ bool isReverse;
+ list< TopoDS_Edge >::const_iterator edgeIt = prism.myBottomEdges.begin();
+ for ( int iE = 0; iE < prism.myNbEdgesInWires.front(); ++iE, ++edgeIt )
+ {
+ if ( SMESH_Algo::isDegenerated( *edgeIt )) continue;
+ const TParam2ColumnMap* u2colMap =
+ GetParam2ColumnMap( MeshDS()->ShapeToIndex( *edgeIt ), isReverse );
+ if ( !u2colMap ) return false;
+ double f = u2colMap->begin()->first, l = u2colMap->rbegin()->first;
+ //isReverse = ( edgeIt->Orientation() == TopAbs_REVERSED );
+ //if ( isReverse ) swap ( f, l ); -- u2colMap takes orientation into account
+ const int nbCol = 5;
+ for ( int i = 0; i < nbCol; ++i )
+ {
+ double u = f + i/double(nbCol) * ( l - f );
+ const TNodeColumn* col = & getColumn( u2colMap, u )->second;
+ if ( columns.empty() || col != columns.back() )
+ columns.push_back( col );
+ }
+ }
+ }
+
+ // Find tolerance to check transformations
+
+ double tol2;
+ {
+ Bnd_B3d bndBox;
+ for ( size_t i = 0; i < columns.size(); ++i )
+ bndBox.Add( gpXYZ( columns[i]->front() ));
+ tol2 = bndBox.SquareExtent() * 1e-5;
+ }
+
+ // Compute transformations
+
+ int xCol = -1;
+ gp_Trsf fromCsZ, toCs0;
+ gp_Ax3 cs0 = getLayerCoordSys(0, columns, xCol );
+ //double dist0 = cs0.Location().Distance( gpXYZ( (*columns[0])[0]));
+ toCs0.SetTransformation( cs0 );
+ for ( size_t z = 1; z < zSize; ++z )
+ {
+ gp_Ax3 csZ = getLayerCoordSys(z, columns, xCol );
+ //double distZ = csZ.Location().Distance( gpXYZ( (*columns[0])[z]));
+ fromCsZ.SetTransformation( csZ );
+ fromCsZ.Invert();
+ gp_Trsf& t = trsf[ z-1 ];
+ t = fromCsZ * toCs0;
+ //t.SetScaleFactor( distZ/dist0 ); - it does not work properly, wrong base point
+
+ // check a transformation
+ for ( size_t i = 0; i < columns.size(); ++i )
+ {
+ gp_Pnt p0 = gpXYZ( (*columns[i])[0] );
+ gp_Pnt pz = gpXYZ( (*columns[i])[z] );
+ t.Transforms( p0.ChangeCoord() );
+ if ( p0.SquareDistance( pz ) > tol2 )
+ {
+ t = gp_Trsf();
+ return ( z == zSize - 1 ); // OK if fails only bottom->top trsf
+ }
+ }
+ }
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Check curve orientation of a bottom edge
+ * \param meshDS - mesh DS
+ * \param columnsMap - node columns map of side face
+ * \param bottomEdge - the bottom edge
+ * \param sideFaceID - side face in-block ID
+ * \retval bool - true if orientation coincide with in-block forward orientation
+ */
+//================================================================================
+
+bool StdMeshers_PrismAsBlock::IsForwardEdge(SMESHDS_Mesh* meshDS,
+ const TParam2ColumnMap& columnsMap,
+ const TopoDS_Edge & bottomEdge,
+ const int sideFaceID)
+{
+ bool isForward = false;
+ if ( SMESH_MesherHelper::IsClosedEdge( bottomEdge ))
+ {
+ isForward = ( bottomEdge.Orientation() == TopAbs_FORWARD );
+ }
+ else
+ {
+ const TNodeColumn& firstCol = columnsMap.begin()->second;
+ const SMDS_MeshNode* bottomNode = firstCol[0];
+ TopoDS_Shape firstVertex = SMESH_MesherHelper::GetSubShapeByNode( bottomNode, meshDS );
+ isForward = ( firstVertex.IsSame( TopExp::FirstVertex( bottomEdge, true )));
+ }
+ // on 2 of 4 sides first vertex is end
+ if ( sideFaceID == ID_Fx1z || sideFaceID == ID_F0yz )
+ isForward = !isForward;
+ return isForward;
+}
+
+//=======================================================================
+//function : faceGridToPythonDump
+//purpose : Prints a script creating a normal grid on the prism side
+//=======================================================================
+
+void StdMeshers_PrismAsBlock::faceGridToPythonDump(const SMESH_Block::TShapeID face,
+ const int nb)
+{
+#ifdef _DEBUG_
+ gp_XYZ pOnF[6] = { gp_XYZ(0,0,0), gp_XYZ(0,0,1),
+ gp_XYZ(0,0,0), gp_XYZ(0,1,0),
+ gp_XYZ(0,0,0), gp_XYZ(1,0,0) };
+ gp_XYZ p2;
+ cout << "mesh = smesh.Mesh( 'Face " << face << "')" << endl;
+ SMESH_Block::TFace& f = myFace[ face - ID_FirstF ];
+ gp_XYZ params = pOnF[ face - ID_FirstF ];
+ //const int nb = 10; // nb face rows
+ for ( int j = 0; j <= nb; ++j )
+ {
+ params.SetCoord( f.GetVInd(), double( j )/ nb );
+ for ( int i = 0; i <= nb; ++i )
+ {
+ params.SetCoord( f.GetUInd(), double( i )/ nb );
+ gp_XYZ p = f.Point( params );
+ gp_XY uv = f.GetUV( params );
+ cout << "mesh.AddNode( " << p.X() << ", " << p.Y() << ", " << p.Z() << " )"
+ << " # " << 1 + i + j * ( nb + 1 )
+ << " ( " << i << ", " << j << " ) "
+ << " UV( " << uv.X() << ", " << uv.Y() << " )" << endl;
+ ShellPoint( params, p2 );
+ double dist = ( p2 - p ).Modulus();
+ if ( dist > 1e-4 )
+ cout << "#### dist from ShellPoint " << dist
+ << " (" << p2.X() << ", " << p2.Y() << ", " << p2.Z() << " ) " << endl;
+ }
+ }
+ for ( int j = 0; j < nb; ++j )
+ for ( int i = 0; i < nb; ++i )
+ {
+ int n = 1 + i + j * ( nb + 1 );
+ cout << "mesh.AddFace([ "
+ << n << ", " << n+1 << ", "
+ << n+nb+2 << ", " << n+nb+1 << "]) " << endl;
+ }
+#else
+ (void)face; // unused in release mode
+ (void)nb; // unused in release mode
+#endif
+}
+
+//================================================================================
+/*!
+ * \brief Constructor
+ * \param faceID - in-block ID
+ * \param face - geom FACE
+ * \param baseEdge - EDGE proreply oriented in the bottom EDGE !!!
+ * \param columnsMap - map of node columns
+ * \param first - first normalized param
+ * \param last - last normalized param
+ */
+//================================================================================
+
+StdMeshers_PrismAsBlock::TSideFace::TSideFace(SMESH_Mesh& mesh,
+ const int faceID,
+ const Prism_3D::TQuadList& quadList,
+ const TopoDS_Edge& baseEdge,
+ TParam2ColumnMap* columnsMap,
+ const double first,
+ const double last):
+ myID( faceID ),
+ myParamToColumnMap( columnsMap ),
+ myHelper( mesh )
+{
+ myParams.resize( 1 );
+ myParams[ 0 ] = make_pair( first, last );
+ mySurface = PSurface( new BRepAdaptor_Surface( quadList.front()->face ));
+ myBaseEdge = baseEdge;
+ myIsForward = StdMeshers_PrismAsBlock::IsForwardEdge( myHelper.GetMeshDS(),
+ *myParamToColumnMap,
+ myBaseEdge, myID );
+ myHelper.SetSubShape( quadList.front()->face );
+
+ if ( quadList.size() > 1 ) // side is vertically composite
+ {
+ // fill myShapeID2Surf map to enable finding a right surface by any sub-shape ID
+
+ SMESHDS_Mesh* meshDS = myHelper.GetMeshDS();
+
+ TopTools_IndexedDataMapOfShapeListOfShape subToFaces;
+ Prism_3D::TQuadList::const_iterator quad = quadList.begin();
+ for ( ; quad != quadList.end(); ++quad )
+ {
+ const TopoDS_Face& face = (*quad)->face;
+ TopExp::MapShapesAndAncestors( face, TopAbs_VERTEX, TopAbs_FACE, subToFaces );
+ TopExp::MapShapesAndAncestors( face, TopAbs_EDGE, TopAbs_FACE, subToFaces );
+ myShapeID2Surf.insert( make_pair( meshDS->ShapeToIndex( face ),
+ PSurface( new BRepAdaptor_Surface( face ))));
+ }
+ for ( int i = 1; i <= subToFaces.Extent(); ++i )
+ {
+ const TopoDS_Shape& sub = subToFaces.FindKey( i );
+ TopTools_ListOfShape& faces = subToFaces( i );
+ int subID = meshDS->ShapeToIndex( sub );
+ int faceID = meshDS->ShapeToIndex( faces.First() );
+ myShapeID2Surf.insert ( make_pair( subID, myShapeID2Surf[ faceID ]));
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Constructor of a complex side face
+ */
+//================================================================================
+
+StdMeshers_PrismAsBlock::TSideFace::
+TSideFace(SMESH_Mesh& mesh,
+ const vector< TSideFace* >& components,
+ const vector< pair< double, double> > & params)
+ :myID( components[0] ? components[0]->myID : 0 ),
+ myParamToColumnMap( 0 ),
+ myParams( params ),
+ myIsForward( true ),
+ myComponents( components ),
+ myHelper( mesh )
+{
+ if ( myID == ID_Fx1z || myID == ID_F0yz )
+ {
+ // reverse components
+ std::reverse( myComponents.begin(), myComponents.end() );
+ std::reverse( myParams.begin(), myParams.end() );
+ for ( size_t i = 0; i < myParams.size(); ++i )
+ {
+ const double f = myParams[i].first;
+ const double l = myParams[i].second;
+ myParams[i] = make_pair( 1. - l, 1. - f );
+ }
+ }
+}
+//================================================================================
+/*!
+ * \brief Copy constructor
+ * \param other - other side
+ */
+//================================================================================
+
+StdMeshers_PrismAsBlock::TSideFace::TSideFace( const TSideFace& other ):
+ myID ( other.myID ),
+ myParamToColumnMap ( other.myParamToColumnMap ),
+ mySurface ( other.mySurface ),
+ myBaseEdge ( other.myBaseEdge ),
+ myShapeID2Surf ( other.myShapeID2Surf ),
+ myParams ( other.myParams ),
+ myIsForward ( other.myIsForward ),
+ myComponents ( other.myComponents.size() ),
+ myHelper ( *other.myHelper.GetMesh() )
+{
+ for ( size_t i = 0 ; i < myComponents.size(); ++i )
+ myComponents[ i ] = new TSideFace( *other.myComponents[ i ]);
+}
+
+//================================================================================
+/*!
+ * \brief Deletes myComponents
+ */
+//================================================================================
+
+StdMeshers_PrismAsBlock::TSideFace::~TSideFace()
+{
+ for ( size_t i = 0 ; i < myComponents.size(); ++i )
+ if ( myComponents[ i ] )
+ delete myComponents[ i ];
+}
+
+//================================================================================
+/*!
+ * \brief Return geometry of the vertical curve
+ * \param isMax - true means curve located closer to (1,1,1) block point
+ * \retval Adaptor3d_Curve* - curve adaptor
+ */
+//================================================================================
+
+Adaptor3d_Curve* StdMeshers_PrismAsBlock::TSideFace::VertiCurve(const bool isMax) const
+{
+ if ( !myComponents.empty() ) {
+ if ( isMax )
+ return myComponents.back()->VertiCurve(isMax);
+ else
+ return myComponents.front()->VertiCurve(isMax);
+ }
+ double f = myParams[0].first, l = myParams[0].second;
+ if ( !myIsForward ) std::swap( f, l );
+ return new TVerticalEdgeAdaptor( myParamToColumnMap, isMax ? l : f );
+}
+
+//================================================================================
+/*!
+ * \brief Return geometry of the top or bottom curve
+ * \param isTop -
+ * \retval Adaptor3d_Curve* -
+ */
+//================================================================================
+
+Adaptor3d_Curve* StdMeshers_PrismAsBlock::TSideFace::HorizCurve(const bool isTop) const
+{
+ return new THorizontalEdgeAdaptor( this, isTop );
+}
+
+//================================================================================
+/*!
+ * \brief Return pcurves
+ * \param pcurv - array of 4 pcurves
+ * \retval bool - is a success
+ */
+//================================================================================
+
+bool StdMeshers_PrismAsBlock::TSideFace::GetPCurves(Adaptor2d_Curve2d* pcurv[4]) const
+{
+ int iEdge[ 4 ] = { BOTTOM_EDGE, TOP_EDGE, V0_EDGE, V1_EDGE };
+
+ for ( int i = 0 ; i < 4 ; ++i ) {
+ Handle(Geom2d_Line) line;
+ switch ( iEdge[ i ] ) {
+ case TOP_EDGE:
+ line = new Geom2d_Line( gp_Pnt2d( 0, 1 ), gp::DX2d() ); break;
+ case BOTTOM_EDGE:
+ line = new Geom2d_Line( gp::Origin2d(), gp::DX2d() ); break;
+ case V0_EDGE:
+ line = new Geom2d_Line( gp::Origin2d(), gp::DY2d() ); break;
+ case V1_EDGE:
+ line = new Geom2d_Line( gp_Pnt2d( 1, 0 ), gp::DY2d() ); break;
+ }
+ pcurv[ i ] = new Geom2dAdaptor_Curve( line, 0, 1 );
+ }
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Returns geometry of pcurve on a horizontal face
+ * \param isTop - is top or bottom face
+ * \param horFace - a horizontal face
+ * \retval Adaptor2d_Curve2d* - curve adaptor
+ */
+//================================================================================
+
+Adaptor2d_Curve2d*
+StdMeshers_PrismAsBlock::TSideFace::HorizPCurve(const bool isTop,
+ const TopoDS_Face& horFace) const
+{
+ return new TPCurveOnHorFaceAdaptor( this, isTop, horFace );
+}
+
+//================================================================================
+/*!
+ * \brief Return a component corresponding to parameter
+ * \param U - parameter along a horizontal size
+ * \param localU - parameter along a horizontal size of a component
+ * \retval TSideFace* - found component
+ */
+//================================================================================