+ // 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() && 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();