+
+ //================================================================================
+ /*!
+ * \brief Return coordinate system for z-th layer of nodes
+ */
+ //================================================================================
+
+ gp_Ax2 getLayerCoordSys(const int z,
+ const vector< const TNodeColumn* >& columns,
+ int& xColumn)
+ {
+ // gravity center of a layer
+ gp_XYZ O(0,0,0);
+ int vertexCol = -1;
+ for ( int i = 0; i < columns.size(); ++i )
+ {
+ O += gpXYZ( (*columns[ i ])[ z ]);
+ if ( vertexCol < 0 &&
+ columns[ i ]->front()->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
+ vertexCol = i;
+ }
+ O /= columns.size();
+
+ // Z axis
+ gp_Vec Z(0,0,0);
+ int iPrev = columns.size()-1;
+ for ( int i = 0; i < columns.size(); ++i )
+ {
+ gp_Vec v1( O, gpXYZ( (*columns[ iPrev ])[ z ]));
+ gp_Vec v2( O, gpXYZ( (*columns[ i ] )[ z ]));
+ Z += v1 ^ v2;
+ iPrev = i;
+ }
+
+ if ( vertexCol >= 0 )
+ {
+ O = gpXYZ( (*columns[ vertexCol ])[ z ]);
+ }
+ if ( xColumn < 0 || xColumn >= columns.size() )
+ {
+ // select a column for X dir
+ double maxDist = 0;
+ for ( int i = 0; i < columns.size(); ++i )
+ {
+ double dist = ( O - gpXYZ((*columns[ i ])[ z ])).SquareModulus();
+ if ( dist > maxDist )
+ {
+ xColumn = i;
+ maxDist = dist;
+ }
+ }
+ }
+
+ // X axis
+ gp_Vec X( O, gpXYZ( (*columns[ xColumn ])[ z ]));
+
+ return gp_Ax2( O, Z, X);
+ }
+
+ //================================================================================
+ /*!
+ * \brief Removes submeshes meshed with regular grid from given list
+ * \retval int - nb of removed submeshes
+ */
+ //================================================================================
+
+ int removeQuasiQuads(list< SMESH_subMesh* >& notQuadSubMesh)
+ {
+ int oldNbSM = notQuadSubMesh.size();
+ SMESHDS_Mesh* mesh = notQuadSubMesh.front()->GetFather()->GetMeshDS();
+ list< SMESH_subMesh* >::iterator smIt = notQuadSubMesh.begin();
+#define __NEXT_SM { ++smIt; continue; }
+ while ( smIt != notQuadSubMesh.end() )
+ {
+ SMESH_subMesh* faceSm = *smIt;
+ SMESHDS_SubMesh* faceSmDS = faceSm->GetSubMeshDS();
+ int nbQuads = faceSmDS->NbElements();
+ if ( nbQuads == 0 ) __NEXT_SM;
+
+ // get oredered edges
+ list< TopoDS_Edge > orderedEdges;
+ list< int > nbEdgesInWires;
+ TopoDS_Vertex V000;
+ int nbWires = SMESH_Block::GetOrderedEdges( TopoDS::Face( faceSm->GetSubShape() ),
+ V000, orderedEdges, nbEdgesInWires );
+ if ( nbWires != 1 || nbEdgesInWires.front() <= 4 )
+ __NEXT_SM;
+
+ // get nb of segements on edges
+ list<int> nbSegOnEdge;
+ list< TopoDS_Edge >::iterator edge = orderedEdges.begin();
+ for ( ; edge != orderedEdges.end(); ++edge )
+ {
+ if ( SMESHDS_SubMesh* edgeSmDS = mesh->MeshElements( *edge ))
+ nbSegOnEdge.push_back( edgeSmDS->NbElements() );
+ else
+ nbSegOnEdge.push_back(0);
+ }
+
+ // unite nbSegOnEdge of continues edges
+ int nbEdges = nbEdgesInWires.front();
+ list<int>::iterator nbSegIt = nbSegOnEdge.begin();
+ for ( edge = orderedEdges.begin(); edge != orderedEdges.end(); )
+ {
+ const TopoDS_Edge& e1 = *edge++;
+ const TopoDS_Edge& e2 = ( edge == orderedEdges.end() ? orderedEdges.front() : *edge );
+ if ( SMESH_Algo::IsContinuous( e1, e2 ))
+ {
+ // common vertex of continues edges must be shared by two 2D mesh elems of geom face
+ TopoDS_Vertex vCommon = TopExp::LastVertex( e1, true );
+ const SMDS_MeshNode* vNode = SMESH_Algo::VertexNode( vCommon, mesh );
+ int nbF = 0;
+ if ( vNode )
+ {
+ SMDS_ElemIteratorPtr fIt = vNode->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() )
+ nbF += faceSmDS->Contains( fIt->next() );
+ }
+ list<int>::iterator nbSegIt1 = nbSegIt++;
+ if ( !vNode || nbF == 2 ) // !vNode - two edges can be meshed as one
+ {
+ // unite
+ if ( nbSegIt == nbSegOnEdge.end() ) nbSegIt = nbSegOnEdge.begin();
+ *nbSegIt += *nbSegIt1;
+ nbSegOnEdge.erase( nbSegIt1 );
+ --nbEdges;
+ }
+ }
+ else
+ {
+ ++nbSegIt;
+ }
+ }
+ vector<int> nbSegVec( nbSegOnEdge.begin(), nbSegOnEdge.end());
+ if ( nbSegVec.size() == 4 &&
+ nbSegVec[0] == nbSegVec[2] &&
+ nbSegVec[1] == nbSegVec[3] &&
+ nbSegVec[0] * nbSegVec[1] == nbQuads
+ )
+ smIt = notQuadSubMesh.erase( smIt );
+ else
+ __NEXT_SM;
+ }
+
+ return oldNbSM - notQuadSubMesh.size();
+ }