+//=======================================================================
+//function : allVerticalEdgesStraight
+//purpose : Defines if all "vertical" EDGEs are straight
+//=======================================================================
+
+bool StdMeshers_Prism_3D::allVerticalEdgesStraight( const Prism_3D::TPrismTopo& thePrism )
+{
+ for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
+ {
+ const Prism_3D::TQuadList& quads = thePrism.myWallQuads[i];
+ Prism_3D::TQuadList::const_iterator quadIt = quads.begin();
+ TopoDS_Edge prevQuadEdge;
+ for ( ; quadIt != quads.end(); ++quadIt )
+ {
+ StdMeshers_FaceSidePtr rightSide = (*quadIt)->side[ QUAD_RIGHT_SIDE ];
+
+ if ( !prevQuadEdge.IsNull() &&
+ !SMESH_Algo::IsContinuous( rightSide->Edge( 0 ), prevQuadEdge ))
+ return false;
+
+ for ( int iE = 0; iE < rightSide->NbEdges(); ++iE )
+ {
+ const TopoDS_Edge & rightE = rightSide->Edge( iE );
+ if ( !SMESH_Algo::IsStraight( rightE, /*degenResult=*/true ))
+ return false;
+
+ if ( iE > 0 &&
+ !SMESH_Algo::IsContinuous( rightSide->Edge( iE-1 ), rightE ))
+ return false;
+
+ prevQuadEdge = rightE;
+ }
+ }
+ }
+ return true;
+}
+
+//=======================================================================
+//function : project2dMesh
+//purpose : Project mesh faces from a source FACE of one prism (theSrcFace)
+// to a source FACE of another prism (theTgtFace)
+//=======================================================================
+
+bool StdMeshers_Prism_3D::project2dMesh(const TopoDS_Face& theSrcFace,
+ const TopoDS_Face& theTgtFace)
+{
+ TProjction2dAlgo* projector2D = TProjction2dAlgo::instance( this );
+ projector2D->myHyp.SetSourceFace( theSrcFace );
+ bool ok = projector2D->Compute( *myHelper->GetMesh(), theTgtFace );
+
+ SMESH_subMesh* tgtSM = myHelper->GetMesh()->GetSubMesh( theTgtFace );
+ if ( !ok && tgtSM->GetSubMeshDS() ) {
+ //tgtSM->ComputeStateEngine( SMESH_subMesh::CLEAN ); -- avoid propagation of events
+ SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
+ SMESHDS_SubMesh* tgtSMDS = tgtSM->GetSubMeshDS();
+ for ( SMDS_ElemIteratorPtr eIt = tgtSMDS->GetElements(); eIt->more(); )
+ meshDS->RemoveFreeElement( eIt->next(), tgtSMDS, /*fromGroups=*/false );
+ for ( SMDS_NodeIteratorPtr nIt = tgtSMDS->GetNodes(); nIt->more(); )
+ meshDS->RemoveFreeNode( nIt->next(), tgtSMDS, /*fromGroups=*/false );
+ }
+ tgtSM->ComputeStateEngine ( SMESH_subMesh::CHECK_COMPUTE_STATE );
+ tgtSM->ComputeSubMeshStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+
+ projector2D->SetEventListener( tgtSM );
+
+ return ok;
+}
+
+//================================================================================
+/*!
+ * \brief Set projection coordinates of a node to a face and it's sub-shapes
+ * \param faceID - the face given by in-block ID
+ * \param params - node normalized parameters
+ * \retval bool - is a success
+ */
+//================================================================================
+
+bool StdMeshers_Prism_3D::setFaceAndEdgesXYZ( const int faceID, const gp_XYZ& params, int z )
+{
+ // find base and top edges of the face
+ enum { BASE = 0, TOP, LEFT, RIGHT };
+ vector< int > edgeVec; // 0-base, 1-top
+ SMESH_Block::GetFaceEdgesIDs( faceID, edgeVec );
+
+ myBlock.EdgePoint( edgeVec[ BASE ], params, myShapeXYZ[ edgeVec[ BASE ]]);
+ myBlock.EdgePoint( edgeVec[ TOP ], params, myShapeXYZ[ edgeVec[ TOP ]]);
+
+ SHOWYXZ("\nparams ", params);
+ SHOWYXZ("TOP is " <<edgeVec[ TOP ], myShapeXYZ[ edgeVec[ TOP]]);
+ SHOWYXZ("BASE is "<<edgeVec[ BASE], myShapeXYZ[ edgeVec[ BASE]]);
+
+ if ( faceID == SMESH_Block::ID_Fx0z || faceID == SMESH_Block::ID_Fx1z )
+ {
+ myBlock.EdgePoint( edgeVec[ LEFT ], params, myShapeXYZ[ edgeVec[ LEFT ]]);
+ myBlock.EdgePoint( edgeVec[ RIGHT ], params, myShapeXYZ[ edgeVec[ RIGHT ]]);
+
+ SHOWYXZ("VER "<<edgeVec[ LEFT], myShapeXYZ[ edgeVec[ LEFT]]);
+ SHOWYXZ("VER "<<edgeVec[ RIGHT], myShapeXYZ[ edgeVec[ RIGHT]]);
+ }
+ myBlock.FacePoint( faceID, params, myShapeXYZ[ faceID ]);
+ SHOWYXZ("FacePoint "<<faceID, myShapeXYZ[ faceID]);
+
+ return true;
+}
+
+//=======================================================================
+//function : toSM
+//purpose : If (!isOK), sets the error to a sub-mesh of a current SOLID
+//=======================================================================
+
+bool StdMeshers_Prism_3D::toSM( bool isOK )
+{
+ if ( mySetErrorToSM &&
+ !isOK &&
+ myHelper &&
+ !myHelper->GetSubShape().IsNull() &&
+ myHelper->GetSubShape().ShapeType() == TopAbs_SOLID)
+ {
+ SMESH_subMesh* sm = myHelper->GetMesh()->GetSubMesh( myHelper->GetSubShape() );
+ sm->GetComputeError() = this->GetComputeError();
+ // clear error in order not to return it twice
+ _error = COMPERR_OK;
+ _comment.clear();
+ }
+ return isOK;
+}
+
+//=======================================================================
+//function : shapeID
+//purpose : Return index of a shape
+//=======================================================================
+
+int StdMeshers_Prism_3D::shapeID( const TopoDS_Shape& S )
+{
+ if ( S.IsNull() ) return 0;
+ if ( !myHelper ) return -3;
+ return myHelper->GetMeshDS()->ShapeToIndex( S );
+}
+
+namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
+{
+ struct EdgeWithNeighbors
+ {
+ TopoDS_Edge _edge;
+ int _iBase; /* index in a WIRE with non-base EDGEs excluded */
+ int _iL, _iR; /* used to connect edges in a base FACE */
+ bool _isBase; /* is used in a base FACE */
+ EdgeWithNeighbors(const TopoDS_Edge& E, int iE, int nbE, int shift, bool isBase ):
+ _edge( E ), _iBase( iE + shift ),
+ _iL( SMESH_MesherHelper::WrapIndex( iE-1, Max( 1, nbE )) + shift ),
+ _iR( SMESH_MesherHelper::WrapIndex( iE+1, Max( 1, nbE )) + shift ),
+ _isBase( isBase )
+ {
+ }
+ EdgeWithNeighbors() {}
+ bool IsInternal() const { return !_edge.IsNull() && _edge.Orientation() == TopAbs_INTERNAL; }
+ };
+ // PrismSide contains all FACEs linking a bottom EDGE with a top one.
+ struct PrismSide
+ {
+ TopoDS_Face _face; // a currently treated upper FACE
+ TopTools_IndexedMapOfShape *_faces; // all FACEs (pointer because of a private copy constructor)
+ TopoDS_Edge _topEdge; // a current top EDGE
+ vector< EdgeWithNeighbors >*_edges; // all EDGEs of _face
+ int _iBotEdge; // index of _topEdge within _edges
+ vector< bool > _isCheckedEdge; // mark EDGEs whose two owner FACEs found
+ int _nbCheckedEdges; // nb of EDGEs whose location is defined
+ PrismSide *_leftSide; // neighbor sides
+ PrismSide *_rightSide;
+ bool _isInternal; // whether this side raises from an INTERNAL EDGE
+ void SetExcluded() { _leftSide = _rightSide = NULL; }
+ bool IsExcluded() const { return !_leftSide; }
+ const TopoDS_Edge& Edge( int i ) const
+ {
+ return (*_edges)[ i ]._edge;
+ }
+ int FindEdge( const TopoDS_Edge& E ) const
+ {
+ for ( size_t i = 0; i < _edges->size(); ++i )
+ if ( E.IsSame( Edge( i ))) return i;
+ return -1;
+ }
+ bool IsSideFace( const TopoDS_Shape& face, const bool checkNeighbors ) const
+ {
+ if ( _faces->Contains( face )) // avoid returning true for a prism top FACE
+ return ( !_face.IsNull() || !( face.IsSame( _faces->FindKey( _faces->Extent() ))));
+
+ if ( checkNeighbors )
+ return (( _leftSide && _leftSide->IsSideFace ( face, false )) ||
+ ( _rightSide && _rightSide->IsSideFace( face, false )));
+
+ return false;
+ }
+ };
+ //--------------------------------------------------------------------------------
+ /*!
+ * \brief Return another faces sharing an edge
+ */
+ const TopoDS_Face & getAnotherFace( const TopoDS_Face& face,
+ const TopoDS_Edge& edge,
+ TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge)
+ {
+ TopTools_ListIteratorOfListOfShape faceIt( facesOfEdge.FindFromKey( edge ));
+ for ( ; faceIt.More(); faceIt.Next() )
+ if ( !face.IsSame( faceIt.Value() ))
+ return TopoDS::Face( faceIt.Value() );
+ return face;
+ }
+
+ //--------------------------------------------------------------------------------
+ /*!
+ * \brief Return ordered edges of a face
+ */
+ bool getEdges( const TopoDS_Face& face,
+ vector< EdgeWithNeighbors > & edges,
+ TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge,
+ const bool noHolesAllowed)
+ {
+ TopoDS_Face f = face;
+ if ( f.Orientation() != TopAbs_FORWARD &&
+ f.Orientation() != TopAbs_REVERSED )
+ f.Orientation( TopAbs_FORWARD );
+ list< TopoDS_Edge > ee;
+ list< int > nbEdgesInWires;
+ int nbW = SMESH_Block::GetOrderedEdges( f, ee, nbEdgesInWires );
+ if ( nbW > 1 && noHolesAllowed )
+ return false;
+
+ int iE, nbTot = 0, nbBase, iBase;
+ list< TopoDS_Edge >::iterator e = ee.begin();
+ list< int >::iterator nbE = nbEdgesInWires.begin();
+ for ( ; nbE != nbEdgesInWires.end(); ++nbE )
+ for ( iE = 0; iE < *nbE; ++e, ++iE )
+ if ( SMESH_Algo::isDegenerated( *e )) // degenerated EDGE is never used
+ {
+ e = --ee.erase( e );
+ --(*nbE);
+ --iE;
+ }
+
+ vector<int> isBase;
+ edges.clear();
+ e = ee.begin();
+ for ( nbE = nbEdgesInWires.begin(); nbE != nbEdgesInWires.end(); ++nbE )
+ {
+ nbBase = 0;
+ isBase.resize( *nbE );
+ list< TopoDS_Edge >::iterator eIt = e;
+ for ( iE = 0; iE < *nbE; ++eIt, ++iE )
+ {
+ isBase[ iE ] = ( getAnotherFace( face, *eIt, facesOfEdge ) != face );
+ nbBase += isBase[ iE ];
+ }
+ for ( iBase = 0, iE = 0; iE < *nbE; ++e, ++iE )
+ {
+ edges.push_back( EdgeWithNeighbors( *e, iBase, nbBase, nbTot, isBase[ iE ] ));
+ iBase += isBase[ iE ];
+ }
+ nbTot += nbBase;
+ }
+ if ( nbTot == 0 )
+ return false;
+
+ // IPAL53099. Set correct neighbors to INTERNAL EDGEs, which can be connected to
+ // EDGEs of the outer WIRE but this fact can't be detected by their order.
+ if ( nbW > 1 )
+ {
+ int iFirst = 0, iLast;
+ for ( nbE = nbEdgesInWires.begin(); nbE != nbEdgesInWires.end(); ++nbE )
+ {
+ iLast = iFirst + *nbE - 1;
+ TopoDS_Vertex vv[2] = { SMESH_MesherHelper::IthVertex( 0, edges[ iFirst ]._edge ),
+ SMESH_MesherHelper::IthVertex( 1, edges[ iLast ]._edge ) };
+ bool isConnectOk = ( vv[0].IsSame( vv[1] ));
+ if ( !isConnectOk )
+ {
+ edges[ iFirst ]._iL = edges[ iFirst ]._iBase; // connect to self
+ edges[ iLast ]._iR = edges[ iLast ]._iBase;
+
+ // look for an EDGE of the outer WIREs connected to vv
+ TopoDS_Vertex v0, v1;
+ for ( iE = 0; iE < iFirst; ++iE )
+ {
+ v0 = SMESH_MesherHelper::IthVertex( 0, edges[ iE ]._edge );
+ v1 = SMESH_MesherHelper::IthVertex( 1, edges[ iE ]._edge );
+ if ( vv[0].IsSame( v0 ) || vv[0].IsSame( v1 ))
+ edges[ iFirst ]._iL = edges[ iE ]._iBase;
+ if ( vv[1].IsSame( v0 ) || vv[1].IsSame( v1 ))
+ edges[ iLast ]._iR = edges[ iE ]._iBase;
+ }
+ }
+ iFirst += *nbE;
+ }
+ }
+ return edges.size();
+ }
+
+ //--------------------------------------------------------------------------------
+ /*!
+ * \brief Return number of faces sharing given edges
+ */
+ // int nbAdjacentFaces( const std::vector< EdgeWithNeighbors >& edges,
+ // const TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge )
+ // {
+ // TopTools_MapOfShape adjFaces;
+
+ // for ( size_t i = 0; i < edges.size(); ++i )
+ // {
+ // TopTools_ListIteratorOfListOfShape faceIt( facesOfEdge.FindFromKey( edges[i]._edge ));
+ // for ( ; faceIt.More(); faceIt.Next() )
+ // adjFaces.Add( faceIt.Value() );
+ // }
+ // return adjFaces.Extent();
+ // }
+}
+
+//================================================================================
+/*!
+ * \brief Return true if the algorithm can mesh this shape
+ * \param [in] aShape - shape to check
+ * \param [in] toCheckAll - if true, this check returns OK if all shapes are OK,
+ * else, returns OK if at least one shape is OK
+ */
+//================================================================================
+
+bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckAll)
+{
+ TopExp_Explorer sExp( shape, TopAbs_SOLID );
+ if ( !sExp.More() )
+ return false;
+
+ for ( ; sExp.More(); sExp.Next() )
+ {
+ // check nb shells
+ TopoDS_Shape shell;
+ TopExp_Explorer shExp( sExp.Current(), TopAbs_SHELL );
+ while ( shExp.More() ) {
+ shell = shExp.Current();
+ shExp.Next();
+ if ( shExp.More() && BRep_Tool::IsClosed( shExp.Current() ))
+ shell.Nullify();
+ }
+ if ( shell.IsNull() ) {
+ if ( toCheckAll ) return false;
+ continue;
+ }
+ // get all faces
+ TopTools_IndexedMapOfShape allFaces;
+ TopExp::MapShapes( sExp.Current(), TopAbs_FACE, allFaces );
+ if ( allFaces.Extent() < 3 ) {
+ if ( toCheckAll ) return false;
+ continue;
+ }
+ // is a box?
+ if ( allFaces.Extent() == 6 )
+ {
+ TopTools_IndexedMapOfOrientedShape map;
+ bool isBox = SMESH_Block::FindBlockShapes( TopoDS::Shell( shell ),
+ TopoDS_Vertex(), TopoDS_Vertex(), map );
+ if ( isBox ) {
+ if ( !toCheckAll ) return true;
+ continue;
+ }
+ }
+#ifdef _DEBUG_
+ TopTools_IndexedMapOfShape allShapes; // usage: allShapes.FindIndex( s )
+ TopExp::MapShapes( shape, allShapes );
+#endif
+
+ TopTools_IndexedDataMapOfShapeListOfShape facesOfEdge;
+ TopTools_ListIteratorOfListOfShape faceIt;
+ TopExp::MapShapesAndAncestors( sExp.Current(), TopAbs_EDGE, TopAbs_FACE , facesOfEdge );
+ if ( facesOfEdge.IsEmpty() ) {
+ if ( toCheckAll ) return false;
+ continue;
+ }
+
+ typedef vector< EdgeWithNeighbors > TEdgeWithNeighborsVec;
+ vector< TEdgeWithNeighborsVec > faceEdgesVec( allFaces.Extent() + 1 );
+ const size_t nbEdgesMax = facesOfEdge.Extent() * 2; // there can be seam EDGEs
+ TopTools_IndexedMapOfShape* facesOfSide = new TopTools_IndexedMapOfShape[ nbEdgesMax ];
+ SMESHUtils::ArrayDeleter<TopTools_IndexedMapOfShape> delFacesOfSide( facesOfSide );
+
+ // try to use each face as a bottom one
+ bool prismDetected = false;
+ vector< PrismSide > sides;
+ for ( int iF = 1; iF < allFaces.Extent() && !prismDetected; ++iF )
+ {
+ const TopoDS_Face& botF = TopoDS::Face( allFaces( iF ));
+
+ TEdgeWithNeighborsVec& botEdges = faceEdgesVec[ iF ];
+ if ( botEdges.empty() )
+ if ( !getEdges( botF, botEdges, facesOfEdge, /*noHoles=*/false ))
+ break;
+
+ int nbBase = 0;
+ for ( size_t iS = 0; iS < botEdges.size(); ++iS )
+ nbBase += botEdges[ iS ]._isBase;
+
+ if ( allFaces.Extent()-1 <= nbBase )
+ continue; // all faces are adjacent to botF - no top FACE
+
+ // init data of side FACEs
+ sides.clear();
+ sides.resize( nbBase );
+ size_t iS = 0;
+ for ( size_t iE = 0; iE < botEdges.size(); ++iE )
+ {
+ if ( !botEdges[ iE ]._isBase )
+ continue;
+ sides[ iS ]._topEdge = botEdges[ iE ]._edge;
+ sides[ iS ]._face = botF;
+ sides[ iS ]._leftSide = & sides[ botEdges[ iE ]._iR ];
+ sides[ iS ]._rightSide = & sides[ botEdges[ iE ]._iL ];
+ sides[ iS ]._isInternal = botEdges[ iE ].IsInternal();
+ sides[ iS ]._faces = & facesOfSide[ iS ];
+ sides[ iS ]._faces->Clear();
+ ++iS;
+ }
+
+ bool isOK = true; // ok for a current botF
+ bool isAdvanced = true; // is new data found in a current loop
+ int nbFoundSideFaces = 0;
+ for ( int iLoop = 0; isOK && isAdvanced; ++iLoop )
+ {
+ isAdvanced = false;
+ for ( size_t iS = 0; iS < sides.size() && isOK; ++iS )
+ {
+ PrismSide& side = sides[ iS ];
+ if ( side._face.IsNull() )
+ continue; // probably the prism top face is the last of side._faces
+
+ if ( side._topEdge.IsNull() )
+ {
+ // find vertical EDGEs --- EGDEs shared with neighbor side FACEs
+ for ( int is2nd = 0; is2nd < 2 && isOK; ++is2nd ) // 2 adjacent neighbors
+ {
+ int di = is2nd ? 1 : -1;
+ const PrismSide* adjSide = is2nd ? side._rightSide : side._leftSide;
+ for ( size_t i = 1; i < side._edges->size(); ++i )
+ {
+ int iE = SMESH_MesherHelper::WrapIndex( i*di + side._iBotEdge, side._edges->size());
+ if ( side._isCheckedEdge[ iE ] ) continue;
+ const TopoDS_Edge& vertE = side.Edge( iE );
+ const TopoDS_Shape& neighborF = getAnotherFace( side._face, vertE, facesOfEdge );
+ bool isEdgeShared = (( adjSide->IsSideFace( neighborF, side._isInternal )) ||
+ ( adjSide == &side && neighborF.IsSame( side._face )) );
+ if ( isEdgeShared ) // vertE is shared with adjSide
+ {
+ isAdvanced = true;
+ side._isCheckedEdge[ iE ] = true;
+ side._nbCheckedEdges++;
+ int nbNotCheckedE = side._edges->size() - side._nbCheckedEdges;
+ if ( nbNotCheckedE == 1 )
+ break;
+ }
+ else
+ {
+ if ( i == 1 && iLoop == 0 ) isOK = false;
+ break;
+ }
+ }
+ }
+ // find a top EDGE
+ int nbNotCheckedE = side._edges->size() - side._nbCheckedEdges;
+ if ( nbNotCheckedE == 1 )
+ {
+ vector<bool>::iterator ii = std::find( side._isCheckedEdge.begin(),
+ side._isCheckedEdge.end(), false );
+ if ( ii != side._isCheckedEdge.end() )
+ {
+ size_t iE = std::distance( side._isCheckedEdge.begin(), ii );
+ side._topEdge = side.Edge( iE );
+ }
+ }
+ isOK = ( nbNotCheckedE >= 1 );
+ }
+ else //if ( !side._topEdge.IsNull() )
+ {
+ // get a next face of a side
+ const TopoDS_Shape& f = getAnotherFace( side._face, side._topEdge, facesOfEdge );
+ side._faces->Add( f );
+ bool stop = false;
+ if ( f.IsSame( side._face ) || // _topEdge is a seam
+ SMESH_MesherHelper::Count( f, TopAbs_WIRE, false ) != 1 )
+ {
+ stop = true;
+ }
+ else if ( side._leftSide != & side && // not closed side face
+ side._leftSide->_faces->Contains( f ))
+ {
+ stop = true; // probably f is the prism top face
+ side._leftSide->_face.Nullify();
+ side._leftSide->_topEdge.Nullify();
+ }
+ else if ( side._rightSide != & side &&
+ side._rightSide->_faces->Contains( f ))
+ {
+ stop = true; // probably f is the prism top face
+ side._rightSide->_face.Nullify();
+ side._rightSide->_topEdge.Nullify();
+ }
+ if ( stop )
+ {
+ side._face.Nullify();
+ side._topEdge.Nullify();
+ continue;
+ }
+ side._face = TopoDS::Face( f );
+ int faceID = allFaces.FindIndex( side._face );
+ side._edges = & faceEdgesVec[ faceID ];
+ if ( side._edges->empty() )
+ if ( !getEdges( side._face, * side._edges, facesOfEdge, /*noHoles=*/true ))
+ break;
+ const int nbE = side._edges->size();
+ if ( nbE >= 4 )
+ {
+ isAdvanced = true;
+ ++nbFoundSideFaces;
+ side._iBotEdge = side.FindEdge( side._topEdge );
+ side._isCheckedEdge.clear();
+ side._isCheckedEdge.resize( nbE, false );
+ side._isCheckedEdge[ side._iBotEdge ] = true;
+ side._nbCheckedEdges = 1; // bottom EDGE is known
+ }
+ else // probably a triangular top face found
+ {
+ side._face.Nullify();
+ }
+ side._topEdge.Nullify();
+ isOK = ( !side._edges->empty() || side._faces->Extent() > 1 );
+
+ } //if ( !side._topEdge.IsNull() )
+
+ } // loop on prism sides
+
+ if ( nbFoundSideFaces > allFaces.Extent() )
+ {
+ isOK = false;
+ }
+ if ( iLoop > allFaces.Extent() * 10 )
+ {
+ isOK = false;
+#ifdef _DEBUG_
+ cerr << "BUG: infinite loop in StdMeshers_Prism_3D::IsApplicable()" << endl;
+#endif
+ }
+ } // while isAdvanced
+
+ if ( isOK && sides[0]._faces->Extent() > 1 )
+ {
+ const int nbFaces = sides[0]._faces->Extent();
+ if ( botEdges.size() == 1 ) // cylinder
+ {
+ prismDetected = ( nbFaces == allFaces.Extent()-1 );
+ }
+ else
+ {
+ const TopoDS_Shape& topFace = sides[0]._faces->FindKey( nbFaces );
+ size_t iS;
+ for ( iS = 1; iS < sides.size(); ++iS )
+ if ( ! sides[ iS ]._faces->Contains( topFace ))
+ break;
+ prismDetected = ( iS == sides.size() );
+ }
+ }
+ } // loop on allFaces
+
+ if ( !prismDetected && toCheckAll ) return false;
+ if ( prismDetected && !toCheckAll ) return true;
+
+ } // loop on solids
+
+ return toCheckAll;
+}
+
+namespace Prism_3D
+{
+ //================================================================================
+ /*!
+ * \brief Return true if this node and other one belong to one face
+ */
+ //================================================================================
+
+ bool Prism_3D::TNode::IsNeighbor( const Prism_3D::TNode& other ) const
+ {
+ if ( !other.myNode || !myNode ) return false;
+
+ SMDS_ElemIteratorPtr fIt = other.myNode->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() )
+ if ( fIt->next()->GetNodeIndex( myNode ) >= 0 )
+ return true;
+ return false;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Prism initialization
+ */
+ //================================================================================
+
+ void TPrismTopo::Clear()
+ {
+ myShape3D.Nullify();
+ myTop.Nullify();
+ myBottom.Nullify();
+ myWallQuads.clear();
+ myBottomEdges.clear();
+ myNbEdgesInWires.clear();
+ myWallQuads.clear();
+ }
+
+ //================================================================================
+ /*!
+ * \brief Set upside-down
+ */
+ //================================================================================
+
+ void TPrismTopo::SetUpsideDown()
+ {
+ std::swap( myBottom, myTop );
+ myBottomEdges.clear();
+ std::reverse( myBottomEdges.begin(), myBottomEdges.end() );
+ for ( size_t i = 0; i < myWallQuads.size(); ++i )
+ {
+ myWallQuads[i].reverse();
+ TQuadList::iterator q = myWallQuads[i].begin();
+ for ( ; q != myWallQuads[i].end(); ++q )
+ {
+ (*q)->shift( 2, /*keepUnitOri=*/true );
+ }
+ myBottomEdges.push_back( myWallQuads[i].front()->side[ QUAD_BOTTOM_SIDE ].grid->Edge(0) );
+ }
+ }
+
+} // namespace Prism_3D
+
+//================================================================================
+/*!
+ * \brief Constructor. Initialization is needed
+ */
+//================================================================================
+
+StdMeshers_PrismAsBlock::StdMeshers_PrismAsBlock()
+{
+ mySide = 0;
+}
+
+StdMeshers_PrismAsBlock::~StdMeshers_PrismAsBlock()
+{
+ Clear();
+}
+void StdMeshers_PrismAsBlock::Clear()
+{
+ myHelper = 0;
+ myShapeIDMap.Clear();
+ myError.reset();
+
+ if ( mySide ) {
+ delete mySide; mySide = 0;
+ }
+ myParam2ColumnMaps.clear();
+ myShapeIndex2ColumnMap.clear();
+}
+
+//=======================================================================
+//function : initPrism
+//purpose : Analyse shape geometry and mesh.
+// If there are triangles on one of faces, it becomes 'bottom'.
+// thePrism.myBottom can be already set up.
+//=======================================================================
+
+bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism,
+ const TopoDS_Shape& theShape3D,
+ const bool selectBottom)
+{
+ myHelper->SetSubShape( theShape3D );
+
+ SMESH_subMesh* mainSubMesh = myHelper->GetMesh()->GetSubMeshContaining( theShape3D );
+ if ( !mainSubMesh ) return toSM( error(COMPERR_BAD_INPUT_MESH,"Null submesh of shape3D"));
+
+ // detect not-quad FACE sub-meshes of the 3D SHAPE
+ list< SMESH_subMesh* > notQuadGeomSubMesh;
+ list< SMESH_subMesh* > notQuadElemSubMesh;
+ list< SMESH_subMesh* > meshedSubMesh;
+ int nbFaces = 0;
+ //
+ SMESH_subMeshIteratorPtr smIt = mainSubMesh->getDependsOnIterator(false,true);
+ while ( smIt->more() )
+ {
+ SMESH_subMesh* sm = smIt->next();
+ const TopoDS_Shape& face = sm->GetSubShape();
+ if ( face.ShapeType() > TopAbs_FACE ) break;
+ else if ( face.ShapeType() < TopAbs_FACE ) continue;
+ nbFaces++;
+
+ // is quadrangle FACE?
+ list< TopoDS_Edge > orderedEdges;
+ list< int > nbEdgesInWires;
+ int nbWires = SMESH_Block::GetOrderedEdges( TopoDS::Face( face ), orderedEdges,
+ nbEdgesInWires );
+ if ( nbWires != 1 || nbEdgesInWires.front() != 4 )
+ notQuadGeomSubMesh.push_back( sm );
+
+ // look for a not structured sub-mesh
+ if ( !sm->IsEmpty() )
+ {
+ meshedSubMesh.push_back( sm );
+ if ( !myHelper->IsSameElemGeometry( sm->GetSubMeshDS(), SMDSGeom_QUADRANGLE ) ||
+ !myHelper->IsStructured ( sm ))
+ notQuadElemSubMesh.push_back( sm );
+ }
+ }
+
+ int nbNotQuadMeshed = notQuadElemSubMesh.size();
+ int nbNotQuad = notQuadGeomSubMesh.size();
+ bool hasNotQuad = ( nbNotQuad || nbNotQuadMeshed );
+
+ // detect bad cases
+ if ( nbNotQuadMeshed > 2 )
+ {
+ return toSM( error(COMPERR_BAD_INPUT_MESH,
+ TCom("More than 2 faces with not quadrangle elements: ")
+ <<nbNotQuadMeshed));
+ }
+ if ( nbNotQuad > 2 || !thePrism.myBottom.IsNull() )
+ {
+ // Issue 0020843 - one of side FACEs is quasi-quadrilateral (not 4 EDGEs).
+ // Remove from notQuadGeomSubMesh faces meshed with regular grid
+ int nbQuasiQuads = removeQuasiQuads( notQuadGeomSubMesh, myHelper,
+ TQuadrangleAlgo::instance(this,myHelper) );
+ nbNotQuad -= nbQuasiQuads;
+ if ( nbNotQuad > 2 )
+ return toSM( error(COMPERR_BAD_SHAPE,
+ TCom("More than 2 not quadrilateral faces: ") <<nbNotQuad));
+ hasNotQuad = ( nbNotQuad || nbNotQuadMeshed );
+ }
+
+ // Analyse mesh and topology of FACEs: choose the bottom sub-mesh.
+ // If there are not quadrangle FACEs, they are top and bottom ones.
+ // Not quadrangle FACEs must be only on top and bottom.
+
+ SMESH_subMesh * botSM = 0;
+ SMESH_subMesh * topSM = 0;
+
+ if ( hasNotQuad ) // can choose a bottom FACE
+ {
+ if ( nbNotQuadMeshed > 0 ) botSM = notQuadElemSubMesh.front();
+ else botSM = notQuadGeomSubMesh.front();
+ if ( nbNotQuadMeshed > 1 ) topSM = notQuadElemSubMesh.back();
+ else if ( nbNotQuad > 1 ) topSM = notQuadGeomSubMesh.back();
+
+ if ( topSM == botSM ) {
+ if ( nbNotQuadMeshed > 1 ) topSM = notQuadElemSubMesh.front();
+ else topSM = notQuadGeomSubMesh.front();
+ }
+
+ // detect mesh triangles on wall FACEs
+ if ( nbNotQuad == 2 && nbNotQuadMeshed > 0 ) {
+ bool ok = false;
+ if ( nbNotQuadMeshed == 1 )
+ ok = ( find( notQuadGeomSubMesh.begin(),
+ notQuadGeomSubMesh.end(), botSM ) != notQuadGeomSubMesh.end() );
+ else
+ ok = ( notQuadGeomSubMesh == notQuadElemSubMesh );
+ if ( !ok )
+ return toSM( error(COMPERR_BAD_INPUT_MESH,
+ "Side face meshed with not quadrangle elements"));
+ }
+ }
+
+ thePrism.myNotQuadOnTop = ( nbNotQuadMeshed > 1 );
+
+ // use thePrism.myBottom
+ if ( !thePrism.myBottom.IsNull() )
+ {
+ if ( botSM ) { // <-- not quad geom or mesh on botSM
+ if ( ! botSM->GetSubShape().IsSame( thePrism.myBottom )) {
+ std::swap( botSM, topSM );
+ if ( !botSM || ! botSM->GetSubShape().IsSame( thePrism.myBottom )) {
+ if ( !selectBottom )
+ return toSM( error( COMPERR_BAD_INPUT_MESH,
+ "Incompatible non-structured sub-meshes"));
+ std::swap( botSM, topSM );
+ thePrism.myBottom = TopoDS::Face( botSM->GetSubShape() );
+ }
+ }
+ }
+ else if ( !selectBottom ) {
+ botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom );
+ }
+ }
+ if ( !botSM ) // find a proper bottom
+ {
+ bool savedSetErrorToSM = mySetErrorToSM;
+ mySetErrorToSM = false; // ingore errors in initPrism()
+
+ // search among meshed FACEs
+ list< SMESH_subMesh* >::iterator sm = meshedSubMesh.begin();
+ for ( ; !botSM && sm != meshedSubMesh.end(); ++sm )
+ {
+ thePrism.Clear();
+ botSM = *sm;
+ thePrism.myBottom = TopoDS::Face( botSM->GetSubShape() );
+ if ( !initPrism( thePrism, theShape3D, /*selectBottom=*/false ))
+ botSM = NULL;
+ }
+ // search among all FACEs
+ for ( TopExp_Explorer f( theShape3D, TopAbs_FACE ); !botSM && f.More(); f.Next() )
+ {
+ int minNbFaces = 2 + myHelper->Count( f.Current(), TopAbs_EDGE, false);
+ if ( nbFaces < minNbFaces) continue;
+ thePrism.Clear();
+ thePrism.myBottom = TopoDS::Face( f.Current() );
+ botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom );
+ if ( !initPrism( thePrism, theShape3D, /*selectBottom=*/false ))
+ botSM = NULL;
+ }
+ mySetErrorToSM = savedSetErrorToSM;
+ return botSM ? true : toSM( error( COMPERR_BAD_SHAPE ));
+ }
+
+ // find vertex 000 - the one with smallest coordinates (for easy DEBUG :-)
+ TopoDS_Vertex V000;
+ double minVal = DBL_MAX, minX = 0, val;
+ for ( TopExp_Explorer exp( botSM->GetSubShape(), TopAbs_VERTEX );
+ exp.More(); exp.Next() )
+ {
+ const TopoDS_Vertex& v = TopoDS::Vertex( exp.Current() );
+ gp_Pnt P = BRep_Tool::Pnt( v );
+ val = P.X() + P.Y() + P.Z();
+ if ( val < minVal || ( val == minVal && P.X() < minX )) {
+ V000 = v;
+ minVal = val;
+ minX = P.X();
+ }
+ }
+
+ thePrism.myShape3D = theShape3D;
+ if ( thePrism.myBottom.IsNull() )
+ thePrism.myBottom = TopoDS::Face( botSM->GetSubShape() );
+ thePrism.myBottom.Orientation( myHelper->GetSubShapeOri( theShape3D, thePrism.myBottom ));
+ thePrism.myTop. Orientation( myHelper->GetSubShapeOri( theShape3D, thePrism.myTop ));
+
+ // Get ordered bottom edges
+ TopoDS_Face reverseBottom = // to have order of top EDGEs as in the top FACE
+ TopoDS::Face( thePrism.myBottom.Reversed() );
+ SMESH_Block::GetOrderedEdges( reverseBottom,
+ thePrism.myBottomEdges,
+ thePrism.myNbEdgesInWires, V000 );
+
+ // Get Wall faces corresponding to the ordered bottom edges and the top FACE
+ if ( !getWallFaces( thePrism, nbFaces )) // it also sets thePrism.myTop
+ return false; //toSM( error(COMPERR_BAD_SHAPE, "Can't find side faces"));
+
+ if ( topSM )
+ {
+ if ( !thePrism.myTop.IsSame( topSM->GetSubShape() ))
+ return toSM( error
+ (notQuadGeomSubMesh.empty() ? COMPERR_BAD_INPUT_MESH : COMPERR_BAD_SHAPE,
+ "Non-quadrilateral faces are not opposite"));
+
+ // check that the found top and bottom FACEs are opposite
+ TopTools_IndexedMapOfShape topEdgesMap( thePrism.myBottomEdges.size() );
+ TopExp::MapShapes( thePrism.myTop, topEdgesMap );
+ list< TopoDS_Edge >::iterator edge = thePrism.myBottomEdges.begin();
+ for ( ; edge != thePrism.myBottomEdges.end(); ++edge )
+ if ( topEdgesMap.Contains( *edge ))
+ return toSM( error
+ (notQuadGeomSubMesh.empty() ? COMPERR_BAD_INPUT_MESH : COMPERR_BAD_SHAPE,
+ "Non-quadrilateral faces are not opposite"));
+ }
+
+ if ( thePrism.myBottomEdges.size() > thePrism.myWallQuads.size() )
+ {
+ // composite bottom sides => set thePrism upside-down
+ thePrism.SetUpsideDown();
+ }
+
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Initialization.
+ * \param helper - helper loaded with mesh and 3D shape
+ * \param thePrism - a prism data
+ * \retval bool - false if a mesh or a shape are KO
+ */
+//================================================================================
+
+bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper,
+ const Prism_3D::TPrismTopo& thePrism)
+{
+ myHelper = helper;
+ SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
+ SMESH_Mesh* mesh = myHelper->GetMesh();
+
+ if ( mySide ) {
+ delete mySide; mySide = 0;
+ }
+ vector< TSideFace* > sideFaces( NB_WALL_FACES, 0 );
+ vector< pair< double, double> > params( NB_WALL_FACES );
+ mySide = new TSideFace( *mesh, sideFaces, params );
+
+
+ SMESH_Block::init();
+ myShapeIDMap.Clear();
+ myShapeIndex2ColumnMap.clear();
+
+ int wallFaceIds[ NB_WALL_FACES ] = { // to walk around a block
+ SMESH_Block::ID_Fx0z, SMESH_Block::ID_F1yz,
+ SMESH_Block::ID_Fx1z, SMESH_Block::ID_F0yz
+ };
+
+ myError = SMESH_ComputeError::New();
+
+ myNotQuadOnTop = thePrism.myNotQuadOnTop;
+
+ // Find columns of wall nodes and calculate edges' lengths
+ // --------------------------------------------------------
+
+ myParam2ColumnMaps.clear();
+ myParam2ColumnMaps.resize( thePrism.myBottomEdges.size() ); // total nb edges
+
+ size_t iE, nbEdges = thePrism.myNbEdgesInWires.front(); // nb outer edges
+ vector< double > edgeLength( nbEdges );
+ multimap< double, int > len2edgeMap;
+
+ // for each EDGE: either split into several parts, or join with several next EDGEs
+ vector<int> nbSplitPerEdge( nbEdges, 0 );
+ vector<int> nbUnitePerEdge( nbEdges, 0 ); // -1 means "joined to a previous"
+
+ // consider continuous straight EDGEs as one side
+ const int nbSides = countNbSides( thePrism, nbUnitePerEdge, edgeLength );
+
+ list< TopoDS_Edge >::const_iterator edgeIt = thePrism.myBottomEdges.begin();
+ for ( iE = 0; iE < nbEdges; ++iE, ++edgeIt )
+ {
+ TParam2ColumnMap & faceColumns = myParam2ColumnMaps[ iE ];
+
+ Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[ iE ].begin();
+ for ( ; quad != thePrism.myWallQuads[ iE ].end(); ++quad )
+ {
+ const TopoDS_Edge& quadBot = (*quad)->side[ QUAD_BOTTOM_SIDE ].grid->Edge( 0 );
+ if ( !myHelper->LoadNodeColumns( faceColumns, (*quad)->face, quadBot, meshDS ))
+ return error(COMPERR_BAD_INPUT_MESH, TCom("Can't find regular quadrangle mesh ")
+ << "on a side face #" << MeshDS()->ShapeToIndex( (*quad)->face ));
+ }
+ SHOWYXZ("\np1 F " <<iE, gpXYZ(faceColumns.begin()->second.front() ));
+ SHOWYXZ("p2 F " <<iE, gpXYZ(faceColumns.rbegin()->second.front() ));
+ SHOWYXZ("V First "<<iE, BRep_Tool::Pnt( TopExp::FirstVertex(*edgeIt,true )));
+
+ if ( nbSides < NB_WALL_FACES ) // fill map used to split faces
+ len2edgeMap.insert( make_pair( edgeLength[ iE ], iE )); // sort edges by length
+ }
+ // Load columns of internal edges (forming holes)
+ // and fill map ShapeIndex to TParam2ColumnMap for them
+ for ( ; edgeIt != thePrism.myBottomEdges.end() ; ++edgeIt, ++iE )
+ {
+ TParam2ColumnMap & faceColumns = myParam2ColumnMaps[ iE ];
+
+ Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[ iE ].begin();
+ for ( ; quad != thePrism.myWallQuads[ iE ].end(); ++quad )
+ {
+ const TopoDS_Edge& quadBot = (*quad)->side[ QUAD_BOTTOM_SIDE ].grid->Edge( 0 );
+ if ( !myHelper->LoadNodeColumns( faceColumns, (*quad)->face, quadBot, meshDS ))
+ return error(COMPERR_BAD_INPUT_MESH, TCom("Can't find regular quadrangle mesh ")
+ << "on a side face #" << MeshDS()->ShapeToIndex( (*quad)->face ));
+
+ if ( !faceColumns.empty() && (int)faceColumns.begin()->second.size() != VerticalSize() )
+ return error(COMPERR_BAD_INPUT_MESH, "Different 'vertical' discretization");
+ }
+ // edge columns
+ int id = MeshDS()->ShapeToIndex( *edgeIt );
+ bool isForward = true; // meaningless for intenal wires
+ myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
+ // columns for vertices
+ // 1
+ const SMDS_MeshNode* n0 = faceColumns.begin()->second.front();
+ id = n0->getshapeId();
+ myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
+ // 2
+ const SMDS_MeshNode* n1 = faceColumns.rbegin()->second.front();
+ id = n1->getshapeId();
+ myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
+
+ // SHOWYXZ("\np1 F " <<iE, gpXYZ(faceColumns.begin()->second.front() ));
+ // SHOWYXZ("p2 F " <<iE, gpXYZ(faceColumns.rbegin()->second.front() ));
+ // SHOWYXZ("V First "<<iE, BRep_Tool::Pnt( TopExp::FirstVertex(*edgeIt,true )));
+ }
+
+ // Create 4 wall faces of a block
+ // -------------------------------
+
+ if ( nbSides <= NB_WALL_FACES ) // ************* Split faces if necessary
+ {
+ if ( nbSides != NB_WALL_FACES ) // define how to split
+ {
+ if ( len2edgeMap.size() != nbEdges )
+ RETURN_BAD_RESULT("Uniqueness of edge lengths not assured");
+
+ multimap< double, int >::reverse_iterator maxLen_i = len2edgeMap.rbegin();
+ multimap< double, int >::reverse_iterator midLen_i = ++len2edgeMap.rbegin();
+
+ double maxLen = maxLen_i->first;
+ double midLen = ( len2edgeMap.size() == 1 ) ? 0 : midLen_i->first;
+ switch ( nbEdges ) {
+ case 1: // 0-th edge is split into 4 parts
+ nbSplitPerEdge[ 0 ] = 4;
+ break;
+ case 2: // either the longest edge is split into 3 parts, or both edges into halves
+ if ( maxLen / 3 > midLen / 2 ) {
+ nbSplitPerEdge[ maxLen_i->second ] = 3;
+ }
+ else {
+ nbSplitPerEdge[ maxLen_i->second ] = 2;
+ nbSplitPerEdge[ midLen_i->second ] = 2;
+ }
+ break;
+ case 3:
+ if ( nbSides == 2 )
+ // split longest into 3 parts
+ nbSplitPerEdge[ maxLen_i->second ] = 3;
+ else
+ // split longest into halves
+ nbSplitPerEdge[ maxLen_i->second ] = 2;
+ }
+ }
+ }
+ else // **************************** Unite faces
+ {
+ int nbExraFaces = nbSides - 4; // nb of faces to fuse
+ for ( iE = 0; iE < nbEdges; ++iE )
+ {
+ if ( nbUnitePerEdge[ iE ] < 0 )
+ continue;
+ // look for already united faces
+ for ( size_t i = iE; i < iE + nbExraFaces; ++i )
+ {
+ if ( nbUnitePerEdge[ i ] > 0 ) // a side including nbUnitePerEdge[i]+1 edge
+ nbExraFaces += nbUnitePerEdge[ i ];
+ nbUnitePerEdge[ i ] = -1;
+ }
+ nbUnitePerEdge[ iE ] = nbExraFaces;
+ break;
+ }
+ }
+
+ // Create TSideFace's
+ int iSide = 0;
+ list< TopoDS_Edge >::const_iterator botE = thePrism.myBottomEdges.begin();
+ for ( iE = 0; iE < nbEdges; ++iE, ++botE )
+ {
+ TFaceQuadStructPtr quad = thePrism.myWallQuads[ iE ].front();
+ const int nbSplit = nbSplitPerEdge[ iE ];
+ const int nbExraFaces = nbUnitePerEdge[ iE ] + 1;
+ if ( nbSplit > 0 ) // split
+ {
+ vector< double > params;
+ splitParams( nbSplit, &myParam2ColumnMaps[ iE ], params );
+ const bool isForward =
+ StdMeshers_PrismAsBlock::IsForwardEdge( myHelper->GetMeshDS(),
+ myParam2ColumnMaps[iE],
+ *botE, SMESH_Block::ID_Fx0z );
+ for ( int i = 0; i < nbSplit; ++i ) {
+ double f = ( isForward ? params[ i ] : params[ nbSplit - i-1 ]);
+ double l = ( isForward ? params[ i+1 ] : params[ nbSplit - i ]);
+ TSideFace* comp = new TSideFace( *mesh, wallFaceIds[ iSide ],
+ thePrism.myWallQuads[ iE ], *botE,
+ &myParam2ColumnMaps[ iE ], f, l );
+ mySide->SetComponent( iSide++, comp );
+ }
+ }
+ else if ( nbExraFaces > 1 ) // unite
+ {
+ double u0 = 0, sumLen = 0;
+ for ( size_t i = iE; i < iE + nbExraFaces; ++i )
+ sumLen += edgeLength[ i ];
+
+ vector< TSideFace* > components( nbExraFaces );
+ vector< pair< double, double> > params( nbExraFaces );
+ bool endReached = false;
+ for ( int i = 0; i < nbExraFaces; ++i, ++botE, ++iE )
+ {
+ if ( iE == nbEdges )
+ {
+ endReached = true;
+ botE = thePrism.myBottomEdges.begin();
+ iE = 0;
+ }
+ components[ i ] = new TSideFace( *mesh, wallFaceIds[ iSide ],
+ thePrism.myWallQuads[ iE ], *botE,
+ &myParam2ColumnMaps[ iE ]);
+ double u1 = u0 + edgeLength[ iE ] / sumLen;
+ params[ i ] = make_pair( u0 , u1 );
+ u0 = u1;
+ }
+ TSideFace* comp = new TSideFace( *mesh, components, params );
+ mySide->SetComponent( iSide++, comp );
+ if ( endReached )
+ break;
+ --iE; // for increment in an external loop on iE
+ --botE;
+ }
+ else if ( nbExraFaces < 0 ) // skip already united face
+ {
+ }
+ else // use as is
+ {
+ TSideFace* comp = new TSideFace( *mesh, wallFaceIds[ iSide ],
+ thePrism.myWallQuads[ iE ], *botE,
+ &myParam2ColumnMaps[ iE ]);
+ mySide->SetComponent( iSide++, comp );
+ }
+ }
+
+
+ // Fill geometry fields of SMESH_Block
+ // ------------------------------------
+
+ vector< int > botEdgeIdVec;
+ SMESH_Block::GetFaceEdgesIDs( ID_BOT_FACE, botEdgeIdVec );
+
+ bool isForward[NB_WALL_FACES] = { true, true, true, true };
+ Adaptor2d_Curve2d* botPcurves[NB_WALL_FACES];
+ Adaptor2d_Curve2d* topPcurves[NB_WALL_FACES];
+
+ for ( int iF = 0; iF < NB_WALL_FACES; ++iF )
+ {
+ TSideFace * sideFace = mySide->GetComponent( iF );
+ if ( !sideFace )
+ RETURN_BAD_RESULT("NULL TSideFace");
+ int fID = sideFace->FaceID(); // in-block ID
+
+ // fill myShapeIDMap
+ if ( sideFace->InsertSubShapes( myShapeIDMap ) != 8 &&
+ !sideFace->IsComplex())
+ MESSAGE( ": Warning : InsertSubShapes() < 8 on side " << iF );
+
+ // side faces geometry
+ Adaptor2d_Curve2d* pcurves[NB_WALL_FACES];
+ if ( !sideFace->GetPCurves( pcurves ))
+ RETURN_BAD_RESULT("TSideFace::GetPCurves() failed");
+
+ SMESH_Block::TFace& tFace = myFace[ fID - ID_FirstF ];
+ tFace.Set( fID, sideFace->Surface(), pcurves, isForward );
+
+ SHOWYXZ( endl<<"F "<< iF << " id " << fID << " FRW " << sideFace->IsForward(), sideFace->Value(0,0));
+ // edges 3D geometry
+ vector< int > edgeIdVec;
+ SMESH_Block::GetFaceEdgesIDs( fID, edgeIdVec );
+ for ( int isMax = 0; isMax < 2; ++isMax ) {
+ {
+ int eID = edgeIdVec[ isMax ];
+ SMESH_Block::TEdge& tEdge = myEdge[ eID - ID_FirstE ];
+ tEdge.Set( eID, sideFace->HorizCurve(isMax), true);
+ SHOWYXZ(eID<<" HOR"<<isMax<<"(0)", sideFace->HorizCurve(isMax)->Value(0));
+ SHOWYXZ(eID<<" HOR"<<isMax<<"(1)", sideFace->HorizCurve(isMax)->Value(1));
+ }
+ {
+ int eID = edgeIdVec[ isMax+2 ];
+ SMESH_Block::TEdge& tEdge = myEdge[ eID - ID_FirstE ];
+ tEdge.Set( eID, sideFace->VertiCurve(isMax), true);
+ SHOWYXZ(eID<<" VER"<<isMax<<"(0)", sideFace->VertiCurve(isMax)->Value(0));
+ SHOWYXZ(eID<<" VER"<<isMax<<"(1)", sideFace->VertiCurve(isMax)->Value(1));
+
+ // corner points
+ vector< int > vertexIdVec;
+ SMESH_Block::GetEdgeVertexIDs( eID, vertexIdVec );
+ myPnt[ vertexIdVec[0] - ID_FirstV ] = tEdge.GetCurve()->Value(0).XYZ();
+ myPnt[ vertexIdVec[1] - ID_FirstV ] = tEdge.GetCurve()->Value(1).XYZ();
+ }
+ }
+ // pcurves on horizontal faces
+ for ( iE = 0; iE < NB_WALL_FACES; ++iE ) {
+ if ( edgeIdVec[ BOTTOM_EDGE ] == botEdgeIdVec[ iE ] ) {
+ botPcurves[ iE ] = sideFace->HorizPCurve( false, thePrism.myBottom );
+ topPcurves[ iE ] = sideFace->HorizPCurve( true, thePrism.myTop );
+ break;
+ }
+ }
+ //sideFace->dumpNodes( 4 ); // debug
+ }
+ // horizontal faces geometry
+ {
+ SMESH_Block::TFace& tFace = myFace[ ID_BOT_FACE - ID_FirstF ];
+ tFace.Set( ID_BOT_FACE, new BRepAdaptor_Surface( thePrism.myBottom ), botPcurves, isForward );
+ SMESH_Block::Insert( thePrism.myBottom, ID_BOT_FACE, myShapeIDMap );
+ }
+ {
+ SMESH_Block::TFace& tFace = myFace[ ID_TOP_FACE - ID_FirstF ];
+ tFace.Set( ID_TOP_FACE, new BRepAdaptor_Surface( thePrism.myTop ), topPcurves, isForward );
+ SMESH_Block::Insert( thePrism.myTop, ID_TOP_FACE, myShapeIDMap );
+ }
+ //faceGridToPythonDump( SMESH_Block::ID_Fxy0, 50 );
+ //faceGridToPythonDump( SMESH_Block::ID_Fxy1 );
+
+ // Fill map ShapeIndex to TParam2ColumnMap
+ // ----------------------------------------
+
+ list< TSideFace* > fList;
+ list< TSideFace* >::iterator fListIt;
+ fList.push_back( mySide );
+ for ( fListIt = fList.begin(); fListIt != fList.end(); ++fListIt)
+ {
+ int nb = (*fListIt)->NbComponents();
+ for ( int i = 0; i < nb; ++i ) {
+ if ( TSideFace* comp = (*fListIt)->GetComponent( i ))
+ fList.push_back( comp );
+ }
+ if ( TParam2ColumnMap* cols = (*fListIt)->GetColumns()) {
+ // columns for a base edge
+ int id = MeshDS()->ShapeToIndex( (*fListIt)->BaseEdge() );
+ bool isForward = (*fListIt)->IsForward();
+ myShapeIndex2ColumnMap[ id ] = make_pair( cols, isForward );
+
+ // 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 int 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 ( int 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 botton->top trsf
+ }
+ }
+ }
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Check curve orientation of a bootom edge
+ * \param meshDS - mesh DS
+ * \param columnsMap - node columns map of side face
+ * \param bottomEdge - the bootom edge
+ * \param sideFaceID - side face in-block ID
+ * \retval bool - true if orientation coinside 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;
+ }
+
+#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
+ */
+//================================================================================
+
+StdMeshers_PrismAsBlock::TSideFace*
+StdMeshers_PrismAsBlock::TSideFace::GetComponent(const double U,double & localU) const
+{
+ localU = U;
+ if ( myComponents.empty() )
+ return const_cast<TSideFace*>( this );
+
+ size_t i;
+ for ( i = 0; i < myComponents.size(); ++i )
+ if ( U < myParams[ i ].second )
+ break;
+ if ( i >= myComponents.size() )
+ i = myComponents.size() - 1;
+
+ double f = myParams[ i ].first, l = myParams[ i ].second;
+ localU = ( U - f ) / ( l - f );
+ return myComponents[ i ];
+}
+
+//================================================================================
+/*!
+ * \brief Find node columns for a parameter
+ * \param U - parameter along a horizontal edge
+ * \param col1 - the 1st found column
+ * \param col2 - the 2nd found column
+ * \retval r - normalized position of U between the found columns
+ */
+//================================================================================
+
+double StdMeshers_PrismAsBlock::TSideFace::GetColumns(const double U,
+ TParam2ColumnIt & col1,
+ TParam2ColumnIt & col2) const
+{
+ double u = U, r = 0;
+ if ( !myComponents.empty() ) {
+ TSideFace * comp = GetComponent(U,u);
+ return comp->GetColumns( u, col1, col2 );
+ }
+
+ if ( !myIsForward )
+ u = 1 - u;
+ double f = myParams[0].first, l = myParams[0].second;
+ u = f + u * ( l - f );
+
+ col1 = col2 = getColumn( myParamToColumnMap, u );
+ if ( ++col2 == myParamToColumnMap->end() ) {
+ --col2;
+ r = 0.5;
+ }
+ else {
+ double uf = col1->first;
+ double ul = col2->first;
+ r = ( u - uf ) / ( ul - uf );
+ }
+ return r;
+}
+
+//================================================================================
+/*!
+ * \brief Return all nodes at a given height together with their normalized parameters
+ * \param [in] Z - the height of interest
+ * \param [out] nodes - map of parameter to node
+ */
+//================================================================================
+
+void StdMeshers_PrismAsBlock::
+TSideFace::GetNodesAtZ(const int Z,
+ map<double, const SMDS_MeshNode* >& nodes ) const
+{
+ if ( !myComponents.empty() )
+ {
+ double u0 = 0.;
+ for ( size_t i = 0; i < myComponents.size(); ++i )
+ {
+ map<double, const SMDS_MeshNode* > nn;
+ myComponents[i]->GetNodesAtZ( Z, nn );
+ map<double, const SMDS_MeshNode* >::iterator u2n = nn.begin();
+ if ( !nodes.empty() && nodes.rbegin()->second == u2n->second )
+ ++u2n;
+ const double uRange = myParams[i].second - myParams[i].first;
+ for ( ; u2n != nn.end(); ++u2n )
+ nodes.insert( nodes.end(), make_pair( u0 + uRange * u2n->first, u2n->second ));
+ u0 += uRange;
+ }
+ }
+ else
+ {
+ double f = myParams[0].first, l = myParams[0].second;
+ if ( !myIsForward )
+ std::swap( f, l );
+ const double uRange = l - f;
+ if ( Abs( uRange ) < std::numeric_limits<double>::min() )
+ return;
+ TParam2ColumnIt u2col = getColumn( myParamToColumnMap, myParams[0].first + 1e-3 );
+ for ( ; u2col != myParamToColumnMap->end(); ++u2col )
+ if ( u2col->first > myParams[0].second + 1e-9 )
+ break;
+ else
+ nodes.insert( nodes.end(),
+ make_pair( ( u2col->first - f ) / uRange, u2col->second[ Z ] ));
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Return coordinates by normalized params
+ * \param U - horizontal param
+ * \param V - vertical param
+ * \retval gp_Pnt - result point
+ */
+//================================================================================
+
+gp_Pnt StdMeshers_PrismAsBlock::TSideFace::Value(const Standard_Real U,
+ const Standard_Real V) const
+{
+ if ( !myComponents.empty() ) {
+ double u;
+ TSideFace * comp = GetComponent(U,u);
+ return comp->Value( u, V );
+ }
+
+ TParam2ColumnIt u_col1, u_col2;
+ double vR, hR = GetColumns( U, u_col1, u_col2 );
+
+ const SMDS_MeshNode* nn[4];
+
+ // BEGIN issue 0020680: Bad cell created by Radial prism in center of torus
+ // Workaround for a wrongly located point returned by mySurface.Value() for
+ // UV located near boundary of BSpline surface.
+ // To bypass the problem, we take point from 3D curve of EDGE.
+ // It solves pb of the bloc_fiss_new.py
+ const double tol = 1e-3;
+ if ( V < tol || V+tol >= 1. )
+ {
+ nn[0] = V < tol ? u_col1->second.front() : u_col1->second.back();
+ nn[2] = V < tol ? u_col2->second.front() : u_col2->second.back();
+ TopoDS_Edge edge;
+ if ( V < tol )
+ {
+ edge = myBaseEdge;
+ }
+ else
+ {
+ TopoDS_Shape s = myHelper.GetSubShapeByNode( nn[0], myHelper.GetMeshDS() );
+ if ( s.ShapeType() != TopAbs_EDGE )
+ s = myHelper.GetSubShapeByNode( nn[2], myHelper.GetMeshDS() );
+ if ( !s.IsNull() && s.ShapeType() == TopAbs_EDGE )
+ edge = TopoDS::Edge( s );
+ }
+ if ( !edge.IsNull() )
+ {
+ double u1 = myHelper.GetNodeU( edge, nn[0], nn[2] );
+ double u3 = myHelper.GetNodeU( edge, nn[2], nn[0] );
+ double u = u1 * ( 1 - hR ) + u3 * hR;
+ TopLoc_Location loc; double f,l;
+ Handle(Geom_Curve) curve = BRep_Tool::Curve( edge,loc,f,l );
+ return curve->Value( u ).Transformed( loc );
+ }
+ }
+ // END issue 0020680: Bad cell created by Radial prism in center of torus
+
+ vR = getRAndNodes( & u_col1->second, V, nn[0], nn[1] );
+ vR = getRAndNodes( & u_col2->second, V, nn[2], nn[3] );
+
+ if ( !myShapeID2Surf.empty() ) // side is vertically composite
+ {
+ // find a FACE on which the 4 nodes lie
+ TSideFace* me = (TSideFace*) this;
+ int notFaceID1 = 0, notFaceID2 = 0;
+ for ( int i = 0; i < 4; ++i )
+ if ( nn[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) // node on FACE
+ {
+ me->mySurface = me->myShapeID2Surf[ nn[i]->getshapeId() ];
+ notFaceID2 = 0;
+ break;
+ }
+ else if ( notFaceID1 == 0 ) // node on EDGE or VERTEX
+ {
+ me->mySurface = me->myShapeID2Surf[ nn[i]->getshapeId() ];
+ notFaceID1 = nn[i]->getshapeId();
+ }
+ else if ( notFaceID1 != nn[i]->getshapeId() ) // node on other EDGE or VERTEX
+ {
+ if ( mySurface != me->myShapeID2Surf[ nn[i]->getshapeId() ])
+ notFaceID2 = nn[i]->getshapeId();
+ }
+ if ( notFaceID2 ) // no nodes of FACE and nodes are on different FACEs
+ {
+ SMESHDS_Mesh* meshDS = myHelper.GetMeshDS();
+ TopoDS_Shape face = myHelper.GetCommonAncestor( meshDS->IndexToShape( notFaceID1 ),
+ meshDS->IndexToShape( notFaceID2 ),
+ *myHelper.GetMesh(),
+ TopAbs_FACE );
+ if ( face.IsNull() )
+ throw SALOME_Exception("StdMeshers_PrismAsBlock::TSideFace::Value() face.IsNull()");
+ int faceID = meshDS->ShapeToIndex( face );
+ me->mySurface = me->myShapeID2Surf[ faceID ];
+ if ( !mySurface )
+ throw SALOME_Exception("StdMeshers_PrismAsBlock::TSideFace::Value() !mySurface");
+ }
+ }
+ ((TSideFace*) this)->myHelper.SetSubShape( mySurface->Face() );
+
+ gp_XY uv1 = myHelper.GetNodeUV( mySurface->Face(), nn[0], nn[2]);
+ gp_XY uv2 = myHelper.GetNodeUV( mySurface->Face(), nn[1], nn[3]);
+ gp_XY uv12 = uv1 * ( 1 - vR ) + uv2 * vR;
+
+ gp_XY uv3 = myHelper.GetNodeUV( mySurface->Face(), nn[2], nn[0]);
+ gp_XY uv4 = myHelper.GetNodeUV( mySurface->Face(), nn[3], nn[1]);
+ gp_XY uv34 = uv3 * ( 1 - vR ) + uv4 * vR;
+
+ gp_XY uv = uv12 * ( 1 - hR ) + uv34 * hR;
+
+ gp_Pnt p = mySurface->Value( uv.X(), uv.Y() );
+ return p;
+}
+
+
+//================================================================================
+/*!
+ * \brief Return boundary edge
+ * \param edge - edge index
+ * \retval TopoDS_Edge - found edge
+ */
+//================================================================================
+
+TopoDS_Edge StdMeshers_PrismAsBlock::TSideFace::GetEdge(const int iEdge) const
+{
+ if ( !myComponents.empty() ) {
+ switch ( iEdge ) {
+ case V0_EDGE : return myComponents.front()->GetEdge( iEdge );
+ case V1_EDGE : return myComponents.back() ->GetEdge( iEdge );
+ default: return TopoDS_Edge();
+ }
+ }
+ TopoDS_Shape edge;
+ const SMDS_MeshNode* node = 0;
+ SMESHDS_Mesh * meshDS = myHelper.GetMesh()->GetMeshDS();
+ TNodeColumn* column;
+
+ switch ( iEdge ) {
+ case TOP_EDGE:
+ case BOTTOM_EDGE:
+ column = & (( ++myParamToColumnMap->begin())->second );
+ node = ( iEdge == TOP_EDGE ) ? column->back() : column->front();
+ edge = myHelper.GetSubShapeByNode ( node, meshDS );
+ if ( edge.ShapeType() == TopAbs_VERTEX ) {
+ column = & ( myParamToColumnMap->begin()->second );
+ node = ( iEdge == TOP_EDGE ) ? column->back() : column->front();
+ }
+ break;
+ case V0_EDGE:
+ case V1_EDGE: {
+ bool back = ( iEdge == V1_EDGE );
+ if ( !myIsForward ) back = !back;
+ if ( back )
+ column = & ( myParamToColumnMap->rbegin()->second );
+ else
+ column = & ( myParamToColumnMap->begin()->second );
+ if ( column->size() > 0 )
+ edge = myHelper.GetSubShapeByNode( (*column)[ 1 ], meshDS );
+ if ( edge.IsNull() || edge.ShapeType() == TopAbs_VERTEX )
+ node = column->front();
+ break;
+ }
+ default:;
+ }
+ if ( !edge.IsNull() && edge.ShapeType() == TopAbs_EDGE )
+ return TopoDS::Edge( edge );
+
+ // find edge by 2 vertices
+ TopoDS_Shape V1 = edge;
+ TopoDS_Shape V2 = myHelper.GetSubShapeByNode( node, meshDS );
+ if ( !V2.IsNull() && V2.ShapeType() == TopAbs_VERTEX && !V2.IsSame( V1 ))
+ {
+ TopoDS_Shape ancestor = myHelper.GetCommonAncestor( V1, V2, *myHelper.GetMesh(), TopAbs_EDGE);
+ if ( !ancestor.IsNull() )
+ return TopoDS::Edge( ancestor );
+ }
+ return TopoDS_Edge();
+}
+
+//================================================================================
+/*!
+ * \brief Fill block sub-shapes
+ * \param shapeMap - map to fill in
+ * \retval int - nb inserted sub-shapes
+ */
+//================================================================================
+
+int StdMeshers_PrismAsBlock::TSideFace::InsertSubShapes(TBlockShapes& shapeMap) const
+{
+ int nbInserted = 0;
+
+ // Insert edges
+ vector< int > edgeIdVec;
+ SMESH_Block::GetFaceEdgesIDs( myID, edgeIdVec );
+
+ for ( int i = BOTTOM_EDGE; i <=V1_EDGE ; ++i ) {
+ TopoDS_Edge e = GetEdge( i );
+ if ( !e.IsNull() ) {
+ nbInserted += SMESH_Block::Insert( e, edgeIdVec[ i ], shapeMap);
+ }
+ }
+
+ // Insert corner vertices
+
+ TParam2ColumnIt col1, col2 ;
+ vector< int > vertIdVec;
+
+ // from V0 column
+ SMESH_Block::GetEdgeVertexIDs( edgeIdVec[ V0_EDGE ], vertIdVec);
+ GetColumns(0, col1, col2 );
+ const SMDS_MeshNode* node0 = col1->second.front();
+ const SMDS_MeshNode* node1 = col1->second.back();
+ TopoDS_Shape v0 = myHelper.GetSubShapeByNode( node0, myHelper.GetMeshDS());
+ TopoDS_Shape v1 = myHelper.GetSubShapeByNode( node1, myHelper.GetMeshDS());
+ if ( v0.ShapeType() == TopAbs_VERTEX ) {
+ nbInserted += SMESH_Block::Insert( v0, vertIdVec[ 0 ], shapeMap);
+ }
+ if ( v1.ShapeType() == TopAbs_VERTEX ) {
+ nbInserted += SMESH_Block::Insert( v1, vertIdVec[ 1 ], shapeMap);
+ }
+
+ // from V1 column
+ SMESH_Block::GetEdgeVertexIDs( edgeIdVec[ V1_EDGE ], vertIdVec);
+ GetColumns(1, col1, col2 );
+ node0 = col2->second.front();
+ node1 = col2->second.back();
+ v0 = myHelper.GetSubShapeByNode( node0, myHelper.GetMeshDS());
+ v1 = myHelper.GetSubShapeByNode( node1, myHelper.GetMeshDS());
+ if ( v0.ShapeType() == TopAbs_VERTEX ) {
+ nbInserted += SMESH_Block::Insert( v0, vertIdVec[ 0 ], shapeMap);
+ }
+ if ( v1.ShapeType() == TopAbs_VERTEX ) {
+ nbInserted += SMESH_Block::Insert( v1, vertIdVec[ 1 ], shapeMap);
+ }
+
+// TopoDS_Vertex V0, V1, Vcom;
+// TopExp::Vertices( myBaseEdge, V0, V1, true );
+// if ( !myIsForward ) std::swap( V0, V1 );
+
+// // bottom vertex IDs
+// SMESH_Block::GetEdgeVertexIDs( edgeIdVec[ _u0 ], vertIdVec);
+// SMESH_Block::Insert( V0, vertIdVec[ 0 ], shapeMap);
+// SMESH_Block::Insert( V1, vertIdVec[ 1 ], shapeMap);
+
+// TopoDS_Edge sideEdge = GetEdge( V0_EDGE );
+// if ( sideEdge.IsNull() || !TopExp::CommonVertex( botEdge, sideEdge, Vcom ))
+// return false;
+
+// // insert one side edge
+// int edgeID;
+// if ( Vcom.IsSame( V0 )) edgeID = edgeIdVec[ _v0 ];
+// else edgeID = edgeIdVec[ _v1 ];
+// SMESH_Block::Insert( sideEdge, edgeID, shapeMap);
+
+// // top vertex of the side edge
+// SMESH_Block::GetEdgeVertexIDs( edgeID, vertIdVec);
+// TopoDS_Vertex Vtop = TopExp::FirstVertex( sideEdge );
+// if ( Vcom.IsSame( Vtop ))
+// Vtop = TopExp::LastVertex( sideEdge );
+// SMESH_Block::Insert( Vtop, vertIdVec[ 1 ], shapeMap);