X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Prism_3D.cxx;h=151926415fb5de5b57d312af62301607b0c192f0;hp=f53624f82b5b823776f6e3e92c82778580a726d7;hb=3b5470266af557acb0db5e05c674e386361ce55c;hpb=3611527175ebcbc2093a22dc8be29bd8ed866c8c diff --git a/src/StdMeshers/StdMeshers_Prism_3D.cxx b/src/StdMeshers/StdMeshers_Prism_3D.cxx index f53624f82..151926415 100644 --- a/src/StdMeshers/StdMeshers_Prism_3D.cxx +++ b/src/StdMeshers/StdMeshers_Prism_3D.cxx @@ -1,24 +1,25 @@ -// Copyright (C) 2007-2008 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2011 CEA/DEN, EDF R&D, OPEN CASCADE // -// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, -// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS // -// This library is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 2.1 of the License. +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. // -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. // -// You should have received a copy of the GNU Lesser General Public -// License along with this library; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // -// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // + // SMESH SMESH : implementaion of SMESH idl descriptions // File : StdMeshers_Prism_3D.cxx // Module : SMESH @@ -37,14 +38,18 @@ #include "utilities.h" #include +#include #include #include +#include #include #include #include -#include #include +#include #include +#include +#include using namespace std; @@ -151,6 +156,152 @@ namespace { } params.push_back( parLast ); // 1. } + + //================================================================================ + /*! + * \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 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::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::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 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(); + } } //======================================================================= @@ -162,7 +313,7 @@ StdMeshers_Prism_3D::StdMeshers_Prism_3D(int hypId, int studyId, SMESH_Gen* gen) :SMESH_3D_Algo(hypId, studyId, gen) { _name = "Prism_3D"; - _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID); // 1 bit per shape type + _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID); // 1 bit per shape type myProjectTriangles = false; } @@ -257,75 +408,115 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh // Projections on the top and bottom faces are taken from nodes existing // on these faces; find correspondence between bottom and top nodes myBotToColumnMap.clear(); - if ( !assocOrProjBottom2Top() ) // it also fill myBotToColumnMap + if ( !assocOrProjBottom2Top() ) // it also fills myBotToColumnMap return false; // Create nodes inside the block - // loop on nodes inside the bottom face - TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin(); - for ( ; bot_column != myBotToColumnMap.end(); ++bot_column ) + // try to use transformation (issue 0020680) + vector trsf; + if ( myBlock.GetLayersTransformation(trsf)) { - const TNode& tBotNode = bot_column->first; // bottom TNode - if ( tBotNode.GetPositionType() != SMDS_TOP_FACE ) - continue; // node is not inside face - - // column nodes; middle part of the column are zero pointers - TNodeColumn& column = bot_column->second; - - // bottom node parameters and coords - myShapeXYZ[ ID_BOT_FACE ] = tBotNode.GetCoords(); - gp_XYZ botParams = tBotNode.GetParams(); - - // compute top node parameters - myShapeXYZ[ ID_TOP_FACE ] = gpXYZ( column.back() ); - gp_XYZ topParams = botParams; - topParams.SetZ( 1 ); - if ( column.size() > 2 ) { - gp_Pnt topCoords = myShapeXYZ[ ID_TOP_FACE ]; - if ( !myBlock.ComputeParameters( topCoords, topParams, ID_TOP_FACE, topParams )) - return error(TCom("Can't compute normalized parameters ") - << "for node " << column.back()->GetID() - << " on the face #"<< column.back()->GetPosition()->GetShapeId() ); - } + // loop on nodes inside the bottom face + TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin(); + for ( ; bot_column != myBotToColumnMap.end(); ++bot_column ) + { + const TNode& tBotNode = bot_column->first; // bottom TNode + if ( tBotNode.GetPositionType() != SMDS_TOP_FACE ) + continue; // node is not inside face + + // column nodes; middle part of the column are zero pointers + TNodeColumn& column = bot_column->second; + TNodeColumn::iterator columnNodes = column.begin(); + for ( int z = 0; columnNodes != column.end(); ++columnNodes, ++z) + { + const SMDS_MeshNode* & node = *columnNodes; + if ( node ) continue; // skip bottom or top node - // vertical loop - TNodeColumn::iterator columnNodes = column.begin(); - for ( int z = 0; columnNodes != column.end(); ++columnNodes, ++z) + gp_XYZ coords = tBotNode.GetCoords(); + trsf[z-1].Transforms( coords ); + node = meshDS->AddNode( coords.X(), coords.Y(), coords.Z() ); + meshDS->SetNodeInVolume( node, volumeID ); + } + } // loop on bottom nodes + } + else // use block approach + { + // loop on nodes inside the bottom face + TNode prevBNode; + TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin(); + for ( ; bot_column != myBotToColumnMap.end(); ++bot_column ) { - const SMDS_MeshNode* & node = *columnNodes; - if ( node ) continue; // skip bottom or top node - - // params of a node to create - double rz = (double) z / (double) ( column.size() - 1 ); - gp_XYZ params = botParams * ( 1 - rz ) + topParams * rz; - - // set coords on all faces and nodes - const int nbSideFaces = 4; - int sideFaceIDs[nbSideFaces] = { SMESH_Block::ID_Fx0z, - SMESH_Block::ID_Fx1z, - SMESH_Block::ID_F0yz, - SMESH_Block::ID_F1yz }; - for ( int iF = 0; iF < nbSideFaces; ++iF ) - if ( !setFaceAndEdgesXYZ( sideFaceIDs[ iF ], params, z )) - return false; - - // compute coords for a new node - gp_XYZ coords; - if ( !SMESH_Block::ShellPoint( params, myShapeXYZ, coords )) - return error("Can't compute coordinates by normalized parameters"); - - SHOWYXZ("TOPFacePoint ",myShapeXYZ[ ID_TOP_FACE]); - SHOWYXZ("BOT Node "<< tBotNode.myNode->GetID(),gpXYZ(tBotNode.myNode)); - SHOWYXZ("ShellPoint ",coords); - - // create a node - node = meshDS->AddNode( coords.X(), coords.Y(), coords.Z() ); - meshDS->SetNodeInVolume( node, volumeID ); - } - } // loop on bottom nodes + const TNode& tBotNode = bot_column->first; // bottom TNode + if ( tBotNode.GetPositionType() != SMDS_TOP_FACE ) + continue; // node is not inside face + // column nodes; middle part of the column are zero pointers + TNodeColumn& column = bot_column->second; + + // compute bottom node parameters + gp_XYZ paramHint(-1,-1,-1); + if ( prevBNode.IsNeighbor( tBotNode )) + paramHint = prevBNode.GetParams(); + if ( !myBlock.ComputeParameters( tBotNode.GetCoords(), tBotNode.ChangeParams(), + ID_BOT_FACE, paramHint )) + return error(TCom("Can't compute normalized parameters for node ") + << tBotNode.myNode->GetID() << " on the face #" + << myBlock.SubMesh( ID_BOT_FACE )->GetId() ); + prevBNode = tBotNode; + + myShapeXYZ[ ID_BOT_FACE ] = tBotNode.GetCoords(); + gp_XYZ botParams = tBotNode.GetParams(); + + // compute top node parameters + myShapeXYZ[ ID_TOP_FACE ] = gpXYZ( column.back() ); + gp_XYZ topParams = botParams; + topParams.SetZ( 1 ); + if ( column.size() > 2 ) { + gp_Pnt topCoords = myShapeXYZ[ ID_TOP_FACE ]; + if ( !myBlock.ComputeParameters( topCoords, topParams, ID_TOP_FACE, topParams )) + return error(TCom("Can't compute normalized parameters ") + << "for node " << column.back()->GetID() + << " on the face #"<< column.back()->getshapeId() ); + } + + // vertical loop + TNodeColumn::iterator columnNodes = column.begin(); + for ( int z = 0; columnNodes != column.end(); ++columnNodes, ++z) + { + const SMDS_MeshNode* & node = *columnNodes; + if ( node ) continue; // skip bottom or top node + + // params of a node to create + double rz = (double) z / (double) ( column.size() - 1 ); + gp_XYZ params = botParams * ( 1 - rz ) + topParams * rz; + + // set coords on all faces and nodes + const int nbSideFaces = 4; + int sideFaceIDs[nbSideFaces] = { SMESH_Block::ID_Fx0z, + SMESH_Block::ID_Fx1z, + SMESH_Block::ID_F0yz, + SMESH_Block::ID_F1yz }; + for ( int iF = 0; iF < nbSideFaces; ++iF ) + if ( !setFaceAndEdgesXYZ( sideFaceIDs[ iF ], params, z )) + return false; + + // compute coords for a new node + gp_XYZ coords; + if ( !SMESH_Block::ShellPoint( params, myShapeXYZ, coords )) + return error("Can't compute coordinates by normalized parameters"); + + SHOWYXZ("TOPFacePoint ",myShapeXYZ[ ID_TOP_FACE]); + SHOWYXZ("BOT Node "<< tBotNode.myNode->GetID(),gpXYZ(tBotNode.myNode)); + SHOWYXZ("ShellPoint ",coords); + + // create a node + node = meshDS->AddNode( coords.X(), coords.Y(), coords.Z() ); + meshDS->SetNodeInVolume( node, volumeID ); + } + } // loop on bottom nodes + } // Create volumes @@ -349,7 +540,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh { const SMDS_MeshNode* n = face->GetNode( i ); if ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) { - bot_column = myBotToColumnMap.find( n ); + TNode2ColumnMap::iterator bot_column = myBotToColumnMap.find( n ); if ( bot_column == myBotToColumnMap.end() ) return error(TCom("No nodes found above node ") << n->GetID() ); columns[ i ] = & bot_column->second; @@ -364,6 +555,10 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh AddPrisms( columns, myHelper ); } // loop on bottom mesh faces + + // clear data + myBotToColumnMap.clear(); + myBlock.Clear(); return true; } @@ -375,8 +570,8 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh //======================================================================= bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh& theMesh, - const TopoDS_Shape& theShape, - MapShapeNbElems& aResMap) + const TopoDS_Shape& theShape, + MapShapeNbElems& aResMap) { // find face contains only triangles vector < SMESH_subMesh * >meshFaces; @@ -394,8 +589,8 @@ bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh& theMesh, return false; } std::vector aVec = (*anIt).second; - int nbtri = Max(aVec[3],aVec[4]); - int nbqua = Max(aVec[5],aVec[6]); + int nbtri = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]); + int nbqua = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]); if( nbtri==0 && nbqua>0 ) { NbQFs++; } @@ -405,8 +600,8 @@ bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh& theMesh, } if(NbQFs<4) { - std::vector aResVec(17); - for(int i=0; i<17; i++) aResVec[i] = 0; + std::vector aResVec(SMDSEntity_Last); + for(int i=SMDSEntity_Node; iGetComputeError(); @@ -426,7 +621,7 @@ bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh& theMesh, MapShapeNbElemsItr anIt = aResMap.find(sm); if( anIt == aResMap.end() ) continue; std::vector aVec = (*anIt).second; - nb1d += Max(aVec[1],aVec[2]); + nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]); } } // find face opposite to base face @@ -436,8 +631,8 @@ bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh& theMesh, bool IsOpposite = true; for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) { if( Edges1.Contains(exp.Current()) ) { - IsOpposite = false; - break; + IsOpposite = false; + break; } } if(IsOpposite) { @@ -452,28 +647,29 @@ bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh& theMesh, MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] ); if( anIt == aResMap.end() ) continue; std::vector aVec = (*anIt).second; - nb2d += Max(aVec[5],aVec[6]); + nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]); } MapShapeNbElemsItr anIt = aResMap.find( meshFaces[NumBase-1] ); std::vector aVec = (*anIt).second; - bool IsQuadratic = (aVec[4]>aVec[3]) || (aVec[6]>aVec[5]); - int nb2d_face0_3 = Max(aVec[3],aVec[4]); - int nb2d_face0_4 = Max(aVec[5],aVec[6]); - int nb0d_face0 = aVec[0]; + bool IsQuadratic = (aVec[SMDSEntity_Quad_Triangle]>aVec[SMDSEntity_Triangle]) || + (aVec[SMDSEntity_Quad_Quadrangle]>aVec[SMDSEntity_Quadrangle]); + int nb2d_face0_3 = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]); + int nb2d_face0_4 = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]); + int nb0d_face0 = aVec[SMDSEntity_Node]; int nb1d_face0_int = ( nb2d_face0_3*3 + nb2d_face0_4*4 - nb1d ) / 2; - std::vector aResVec(17); - for(int i=0; i<17; i++) aResVec[i] = 0; + std::vector aResVec(SMDSEntity_Last); + for(int i=SMDSEntity_Node; iNbElements() != topSMDS->NbElements() || botSMDS->NbNodes() != topSMDS->NbNodes()) { + MESSAGE("nb elem bot " << botSMDS->NbElements() << " top " << topSMDS->NbElements()); + MESSAGE("nb node bot " << botSMDS->NbNodes() << " top " << topSMDS->NbNodes()); if ( myBlock.HasNotQuadElemOnTop() ) return error(TCom("Mesh on faces #") << botSM->GetId() <<" and #"<< topSM->GetId() << " seems different" ); @@ -651,7 +849,7 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top() // Fill myBotToColumnMap int zSize = myBlock.VerticalSize(); - TNode prevTNode; + //TNode prevTNode; TNodeNodeMap::iterator bN_tN = n2nMap.begin(); for ( ; bN_tN != n2nMap.end(); ++bN_tN ) { @@ -659,19 +857,8 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top() const SMDS_MeshNode* topNode = bN_tN->second; if ( botNode->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE ) continue; // wall columns are contained in myBlock - // compute bottom node params - TNode bN( botNode ); - if ( zSize > 2 ) { - gp_XYZ paramHint(-1,-1,-1); - if ( prevTNode.IsNeighbor( bN )) - paramHint = prevTNode.GetParams(); - if ( !myBlock.ComputeParameters( bN.GetCoords(), bN.ChangeParams(), - ID_BOT_FACE, paramHint )) - return error(TCom("Can't compute normalized parameters for node ") - << botNode->GetID() << " on the face #"<< botSM->GetId() ); - prevTNode = bN; - } // create node column + TNode bN( botNode ); TNode2ColumnMap::iterator bN_col = myBotToColumnMap.insert( make_pair ( bN, TNodeColumn() )).first; TNodeColumn & column = bN_col->second; @@ -812,10 +999,10 @@ bool StdMeshers_Prism_3D::setFaceAndEdgesXYZ( const int faceID, const gp_XYZ& pa SMESH_Block::GetFaceEdgesIDs( faceID, edgeVec ); myBlock.EdgePoint( edgeVec[ BASE ], params, myShapeXYZ[ edgeVec[ BASE ]]); - myBlock.EdgePoint( edgeVec[ TOP ], params, myShapeXYZ[ edgeVec[ TOP ]]); + myBlock.EdgePoint( edgeVec[ TOP ], params, myShapeXYZ[ edgeVec[ TOP ]]); SHOWYXZ("\nparams ", params); - SHOWYXZ("TOP is "< 0 && nbNotQuad != 2 ) - return error(COMPERR_BAD_SHAPE, - TCom("More than 2 not quadrilateral faces: ") - < 2 ) + { return error(COMPERR_BAD_INPUT_MESH, TCom("More than 2 faces with not quadrangle elements: ") < 0 && nbNotQuad != 2 ) + { + // Issue 0020843 - one of side faces is quasi-quadrilateral. + // Remove from notQuadGeomSubMesh faces meshed with regular grid + nbQuasiQuads = removeQuasiQuads( notQuadGeomSubMesh ); + nbNotQuad -= nbQuasiQuads; + if ( nbNotQuad > 0 && nbNotQuad != 2 ) + return error(COMPERR_BAD_SHAPE, + TCom("More than 2 not quadrilateral faces: ") + < 1 ); + MESSAGE("myNotQuadOnTop " << myNotQuadOnTop << " nbNotQuadMeshed " << nbNotQuadMeshed); // ---------------------------------------------------------- @@ -1110,24 +1318,47 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, // Get ordered bottom edges list< TopoDS_Edge > orderedEdges; - list< int > nbVertexInWires; + list< int > nbEInW; SMESH_Block::GetOrderedEdges( TopoDS::Face( botSM->GetSubShape().Reversed() ), - V000, orderedEdges, nbVertexInWires ); -// if ( nbVertexInWires.size() != 1 ) + V000, orderedEdges, nbEInW ); +// if ( nbEInW.size() != 1 ) // RETURN_BAD_RESULT("Wrong prism geometry"); // Get Wall faces corresponding to the ordered bottom edges list< TopoDS_Face > wallFaces; - if ( !GetWallFaces( Mesh(), shape3D, botSM->GetSubShape(), orderedEdges, wallFaces)) + if ( !GetWallFaces( Mesh(), shape3D, botSM->GetSubShape(), orderedEdges, nbEInW, wallFaces)) return error(COMPERR_BAD_SHAPE, "Can't find side faces"); + // check that the found top and bottom faces are opposite + { + for (TopExp_Explorer edge(botSM->GetSubShape(), TopAbs_EDGE); edge.More(); edge.Next()) + if ( helper->IsSubShape( edge.Current(), topSM->GetSubShape() )) + return error(notQuadGeomSubMesh.empty() ? COMPERR_BAD_INPUT_MESH : COMPERR_BAD_SHAPE, + "Non-quadrilateral faces are not opposite"); + } + + // Protect from a distorted block (test 3D_mesh_HEXA3D/B7 on 32bit platform) + // check that all wall faces have an edge common with the top face + { + list< TopoDS_Face >::iterator faceIt = wallFaces.begin(); + for ( ; faceIt != wallFaces.end(); ++faceIt ) + { + bool hasCommon = false; + for (TopExp_Explorer edge(*faceIt, TopAbs_EDGE); !hasCommon && edge.More(); edge.Next()) + if ( helper->IsSubShape( edge.Current(), topSM->GetSubShape() )) + hasCommon = true; + if ( !hasCommon ) + return error(COMPERR_BAD_SHAPE); + } + } + // Find columns of wall nodes and calculate edges' lengths // -------------------------------------------------------- myParam2ColumnMaps.clear(); myParam2ColumnMaps.resize( orderedEdges.size() ); // total nb edges - int iE, nbEdges = nbVertexInWires.front(); // nb outer edges + int iE, nbEdges = nbEInW.front(); // nb outer edges vector< double > edgeLength( nbEdges ); map< double, int > len2edgeMap; @@ -1173,11 +1404,11 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, // columns for vertices // 1 const SMDS_MeshNode* n0 = faceColumns.begin()->second.front(); - id = n0->GetPosition()->GetShapeId(); + id = n0->getshapeId(); myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward ); // 2 const SMDS_MeshNode* n1 = faceColumns.rbegin()->second.front(); - id = n1->GetPosition()->GetShapeId(); + id = n1->getshapeId(); myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward ); // SHOWYXZ("\np1 F "<second.front() )); // SHOWYXZ("p2 F "<second.front() )); @@ -1352,6 +1583,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, break; } } + //sideFace->dumpNodes( 4 ); // debug } // horizontal faces geometry { @@ -1386,11 +1618,11 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, // columns for vertices const SMDS_MeshNode* n0 = cols->begin()->second.front(); - id = n0->GetPosition()->GetShapeId(); + id = n0->getshapeId(); myShapeIndex2ColumnMap[ id ] = make_pair( cols, isForward ); const SMDS_MeshNode* n1 = cols->rbegin()->second.front(); - id = n1->GetPosition()->GetShapeId(); + id = n1->getshapeId(); myShapeIndex2ColumnMap[ id ] = make_pair( cols, !isForward ); } } @@ -1417,7 +1649,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, const TNodeColumn* StdMeshers_PrismAsBlock::GetNodeColumn(const SMDS_MeshNode* node) const { - int sID = node->GetPosition()->GetShapeId(); + int sID = node->getshapeId(); map >::const_iterator col_frw = myShapeIndex2ColumnMap.find( sID ); @@ -1431,6 +1663,89 @@ const TNodeColumn* StdMeshers_PrismAsBlock::GetNodeColumn(const SMDS_MeshNode* n 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. +//======================================================================= + +bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector & trsf) const +{ + const int zSize = VerticalSize(); + if ( zSize < 3 ) return true; + trsf.resize( zSize - 2 ); + + // Select some node columns by which we will define coordinate system of layers + + vector< const TNodeColumn* > columns; + { + const TopoDS_Shape& baseFace = Shape(ID_BOT_FACE); + list< TopoDS_Edge > orderedEdges; + list< int > nbEdgesInWires; + GetOrderedEdges( TopoDS::Face( baseFace ), TopoDS_Vertex(), orderedEdges, nbEdgesInWires ); + bool isReverse; + list< TopoDS_Edge >::iterator edgeIt = orderedEdges.begin(); + for ( int iE = 0; iE < nbEdgesInWires.front(); ++iE, ++edgeIt ) + { + if ( BRep_Tool::Degenerated( *edgeIt )) continue; + const TParam2ColumnMap* u2colMap = + GetParam2ColumnMap( myHelper->GetMeshDS()->ShapeToIndex( *edgeIt ), isReverse ); + if ( !u2colMap ) return false; + isReverse = ( edgeIt->Orientation() == TopAbs_REVERSED ); + double f = u2colMap->begin()->first, l = u2colMap->rbegin()->first; + if ( isReverse ) swap ( f, l ); + 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 ( int 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-1; ++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 ( int 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 ) + return false; + } + } + return true; +} + //================================================================================ /*! * \brief Check curve orientation of a bootom edge @@ -1448,7 +1763,7 @@ bool StdMeshers_PrismAsBlock::IsForwardEdge(SMESHDS_Mesh* meshDS, const int sideFaceID) { bool isForward = false; - if ( TAssocTool::IsClosedEdge( bottomEdge )) + if ( SMESH_MesherHelper::IsClosedEdge( bottomEdge )) { isForward = ( bottomEdge.Orientation() == TopAbs_FORWARD ); } @@ -1466,41 +1781,59 @@ bool StdMeshers_PrismAsBlock::IsForwardEdge(SMESHDS_Mesh* meshDS, } //================================================================================ - /*! - * \brief Find wall faces by bottom edges - * \param mesh - the mesh - * \param mainShape - the prism - * \param bottomFace - the bottom face - * \param bottomEdges - edges bounding the bottom face - * \param wallFaces - faces list to fill in - */ +/*! + * \brief Find wall faces by bottom edges + * \param mesh - the mesh + * \param mainShape - the prism + * \param bottomFace - the bottom face + * \param bottomEdges - edges bounding the bottom face + * \param wallFaces - faces list to fill in + */ //================================================================================ -bool StdMeshers_PrismAsBlock::GetWallFaces( SMESH_Mesh* mesh, - const TopoDS_Shape & mainShape, - const TopoDS_Shape & bottomFace, - const std::list< TopoDS_Edge >& bottomEdges, - std::list< TopoDS_Face >& wallFaces) +bool StdMeshers_PrismAsBlock::GetWallFaces( SMESH_Mesh* mesh, + const TopoDS_Shape & mainShape, + const TopoDS_Shape & bottomFace, + std::list< TopoDS_Edge >& bottomEdges, + std::list< int > & nbEInW, + std::list< TopoDS_Face >& wallFaces) { wallFaces.clear(); TopTools_IndexedMapOfShape faceMap; TopExp::MapShapes( mainShape, TopAbs_FACE, faceMap ); - list< TopoDS_Edge >::const_iterator edge = bottomEdges.begin(); - for ( ; edge != bottomEdges.end(); ++edge ) + list< TopoDS_Edge >::iterator edge = bottomEdges.begin(); + std::list< int >::iterator nbE = nbEInW.begin(); + int iE = 0; + while ( edge != bottomEdges.end() ) { - TopTools_ListIteratorOfListOfShape ancestIt = mesh->GetAncestors( *edge ); - for ( ; ancestIt.More(); ancestIt.Next() ) + ++iE; + if ( BRep_Tool::Degenerated( *edge )) + { + edge = bottomEdges.erase( edge ); + --iE; + --(*nbE); + } + else { - const TopoDS_Shape& ancestor = ancestIt.Value(); - if ( ancestor.ShapeType() == TopAbs_FACE && // face - !bottomFace.IsSame( ancestor ) && // not bottom - faceMap.FindIndex( ancestor )) // belongs to the prism + PShapeIteratorPtr fIt = myHelper->GetAncestors( *edge, *mesh, TopAbs_FACE ); + while ( fIt->more() ) { - wallFaces.push_back( TopoDS::Face( ancestor )); - break; + const TopoDS_Shape* face = fIt->next(); + if ( !bottomFace.IsSame( *face ) && // not bottom + faceMap.FindIndex( *face )) // belongs to the prism + { + wallFaces.push_back( TopoDS::Face( *face )); + break; + } } + ++edge; + } + if ( iE == *nbE ) + { + iE = 0; + ++nbE; } } return ( wallFaces.size() == bottomEdges.size() ); @@ -1726,8 +2059,6 @@ double StdMeshers_PrismAsBlock::TSideFace::GetColumns(const double U, r = 0.5; } else { -// if ( !myIsForward ) -// std::swap( col1, col2 ); double uf = col1->first; double ul = col2->first; r = ( u - uf ) / ( ul - uf ); @@ -1747,8 +2078,8 @@ double StdMeshers_PrismAsBlock::TSideFace::GetColumns(const double U, gp_Pnt StdMeshers_PrismAsBlock::TSideFace::Value(const Standard_Real U, const Standard_Real V) const { - double u; if ( !myComponents.empty() ) { + double u; TSideFace * comp = GetComponent(U,u); return comp->Value( u, V ); } @@ -1760,7 +2091,41 @@ gp_Pnt StdMeshers_PrismAsBlock::TSideFace::Value(const Standard_Real U, const SMDS_MeshNode* n2 = 0; const SMDS_MeshNode* n3 = 0; const SMDS_MeshNode* n4 = 0; - gp_XYZ pnt; + + // BEGIN issue 0020680: EDF 1252 SMESH: 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. ) + { + n1 = V < tol ? u_col1->second.front() : u_col1->second.back(); + n3 = V < tol ? u_col2->second.front() : u_col2->second.back(); + TopoDS_Edge edge; + if ( V < tol ) + { + edge = myBaseEdge; + } + else + { + TopoDS_Shape s = myHelper->GetSubShapeByNode( n1, myHelper->GetMeshDS() ); + if ( s.ShapeType() != TopAbs_EDGE ) + s = myHelper->GetSubShapeByNode( n3, myHelper->GetMeshDS() ); + if ( s.ShapeType() == TopAbs_EDGE ) + edge = TopoDS::Edge( s ); + } + if ( !edge.IsNull() ) + { + double u1 = myHelper->GetNodeU( edge, n1 ); + double u3 = myHelper->GetNodeU( edge, n3 ); + 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: EDF 1252 SMESH: Bad cell created by Radial prism in center of torus vR = getRAndNodes( & u_col1->second, V, n1, n2 ); vR = getRAndNodes( & u_col2->second, V, n3, n4 ); @@ -1774,8 +2139,9 @@ gp_Pnt StdMeshers_PrismAsBlock::TSideFace::Value(const Standard_Real U, gp_XY uv34 = uv3 * ( 1 - vR ) + uv4 * vR; gp_XY uv = uv12 * ( 1 - hR ) + uv34 * hR; - - return mySurface.Value( uv.X(), uv.Y() ); + + gp_Pnt p = mySurface.Value( uv.X(), uv.Y() ); + return p; } @@ -1836,16 +2202,9 @@ TopoDS_Edge StdMeshers_PrismAsBlock::TSideFace::GetEdge(const int iEdge) const TopoDS_Shape V2 = myHelper->GetSubShapeByNode( node, meshDS ); if ( V2.ShapeType() == TopAbs_VERTEX && !V2.IsSame( V1 )) { - TopTools_ListIteratorOfListOfShape ancestIt = - myHelper->GetMesh()->GetAncestors( V1 ); - for ( ; ancestIt.More(); ancestIt.Next() ) - { - const TopoDS_Shape & ancestor = ancestIt.Value(); - if ( ancestor.ShapeType() == TopAbs_EDGE ) - for ( TopExp_Explorer e( ancestor, TopAbs_VERTEX ); e.More(); e.Next() ) - if ( V2.IsSame( e.Current() )) - return TopoDS::Edge( ancestor ); - } + TopoDS_Shape ancestor = myHelper->GetCommonAncestor( V1, V2, *myHelper->GetMesh(), TopAbs_EDGE); + if ( !ancestor.IsNull() ) + return TopoDS::Edge( ancestor ); } return TopoDS_Edge(); } @@ -1953,6 +2312,28 @@ int StdMeshers_PrismAsBlock::TSideFace::InsertSubShapes(TBlockShapes& shapeMap) return nbInserted; } +//================================================================================ +/*! + * \brief Dump ids of nodes of sides + */ +//================================================================================ + +void StdMeshers_PrismAsBlock::TSideFace::dumpNodes(int nbNodes) const +{ +#ifdef _DEBUG_ + cout << endl << "NODES OF FACE "; SMESH_Block::DumpShapeID( myID, cout ) << endl; + THorizontalEdgeAdaptor* hSize0 = (THorizontalEdgeAdaptor*) HorizCurve(0); + cout << "Horiz side 0: "; hSize0->dumpNodes(nbNodes); cout << endl; + THorizontalEdgeAdaptor* hSize1 = (THorizontalEdgeAdaptor*) HorizCurve(1); + cout << "Horiz side 1: "; hSize1->dumpNodes(nbNodes); cout << endl; + TVerticalEdgeAdaptor* vSide0 = (TVerticalEdgeAdaptor*) VertiCurve(0); + cout << "Verti side 0: "; vSide0->dumpNodes(nbNodes); cout << endl; + TVerticalEdgeAdaptor* vSide1 = (TVerticalEdgeAdaptor*) VertiCurve(1); + cout << "Verti side 1: "; vSide1->dumpNodes(nbNodes); cout << endl; + delete hSize0; delete hSize1; delete vSide0; delete vSide1; +#endif +} + //================================================================================ /*! * \brief Creates TVerticalEdgeAdaptor @@ -1983,6 +2364,22 @@ gp_Pnt StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::Value(const Standard_Real return gpXYZ(n1) * ( 1 - r ) + gpXYZ(n2) * r; } +//================================================================================ +/*! + * \brief Dump ids of nodes + */ +//================================================================================ + +void StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::dumpNodes(int nbNodes) const +{ +#ifdef _DEBUG_ + for ( int i = 0; i < nbNodes && i < myNodeColumn->size(); ++i ) + cout << (*myNodeColumn)[i]->GetID() << " "; + if ( nbNodes < myNodeColumn->size() ) + cout << myNodeColumn->back()->GetID(); +#endif +} + //================================================================================ /*! * \brief Return coordinates for the given normalized parameter @@ -1996,6 +2393,50 @@ gp_Pnt StdMeshers_PrismAsBlock::THorizontalEdgeAdaptor::Value(const Standard_Rea return mySide->TSideFace::Value( U, myV ); } +//================================================================================ +/*! + * \brief Dump ids of first nodes and the last one + */ +//================================================================================ + +void StdMeshers_PrismAsBlock::THorizontalEdgeAdaptor::dumpNodes(int nbNodes) const +{ +#ifdef _DEBUG_ + // Not bedugged code. Last node is sometimes incorrect + const TSideFace* side = mySide; + double u = 0; + if ( mySide->IsComplex() ) + side = mySide->GetComponent(0,u); + + TParam2ColumnIt col, col2; + TParam2ColumnMap* u2cols = side->GetColumns(); + side->GetColumns( u , col, col2 ); + + int j, i = myV ? mySide->ColumnHeight()-1 : 0; + + const SMDS_MeshNode* n = 0; + const SMDS_MeshNode* lastN + = side->IsForward() ? u2cols->rbegin()->second[ i ] : u2cols->begin()->second[ i ]; + for ( j = 0; j < nbNodes && n != lastN; ++j ) + { + n = col->second[ i ]; + cout << n->GetID() << " "; + if ( side->IsForward() ) + ++col; + else + --col; + } + + // last node + u = 1; + if ( mySide->IsComplex() ) + side = mySide->GetComponent(1,u); + + side->GetColumns( u , col, col2 ); + if ( n != col->second[ i ] ) + cout << col->second[ i ]->GetID(); +#endif +} //================================================================================ /*! * \brief Return UV on pcurve for the given normalized parameter