X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=3b5483b13e49c1a72d5677277d0a49067b6a063d;hp=c6c34ed47fda1f8f105c340b2135671ff3f2432c;hb=refs%2Ftags%2FV9_7_0a1;hpb=a1920ff31054e2c882bd94d4f3c04abe53980ce0 diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index c6c34ed47..3b5483b13 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2021 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 @@ -36,6 +36,7 @@ #include "SMESH_HypoFilter.hxx" #include "SMESH_Mesh.hxx" #include "SMESH_MeshAlgos.hxx" +#include "SMESH_MeshEditor.hxx" #include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" @@ -73,7 +74,7 @@ using namespace std; namespace { - inline SMESH_TNodeXYZ XYZ(const SMDS_MeshNode* n) { return SMESH_TNodeXYZ(n); } + inline SMESH_NodeXYZ XYZ(const SMDS_MeshNode* n) { return SMESH_NodeXYZ(n); } enum { U_periodic = 1, V_periodic = 2 }; } @@ -139,9 +140,9 @@ SMESHDS_Mesh* SMESH_MesherHelper::GetMeshDS() const //======================================================================= //function : IsQuadraticSubMesh -//purpose : Check submesh for given shape: if all elements on this shape +//purpose : Check sub-meshes of a given shape: if all elements on sub-shapes // are quadratic, quadratic elements will be created. -// Also fill myTLinkNodeMap +// Fill myTLinkNodeMap //======================================================================= bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) @@ -149,10 +150,7 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) SMESHDS_Mesh* meshDS = GetMeshDS(); // we can create quadratic elements only if all elements // created on sub-shapes of given shape are quadratic - // also we have to fill myTLinkNodeMap myCreateQuadratic = true; - mySeamShapeIds.clear(); - myDegenShapeIds.clear(); TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE ); if ( aSh.ShapeType()==TopAbs_COMPOUND ) { @@ -162,9 +160,6 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) } SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge ); - - //int nbOldLinks = myTLinkNodeMap.size(); - if ( !myMesh->HasShapeToMesh() ) { if (( myCreateQuadratic = myMesh->NbFaces( ORDER_QUADRATIC ))) @@ -373,6 +368,27 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) } } +//======================================================================= +/*! + * \brief Copy shape information from another helper. Used to improve performance + * since SetSubShape() can be time consuming if there are many edges + */ +//======================================================================= + +void SMESH_MesherHelper::CopySubShapeInfo(const SMESH_MesherHelper& other) +{ + this->myShape = other.myShape; + this->myShapeID = other.myShapeID; + this->myDegenShapeIds = other.myDegenShapeIds; + this->mySeamShapeIds = other.mySeamShapeIds; + this->myPar1[0] = other.myPar1[0]; + this->myPar1[1] = other.myPar1[1]; + this->myPar2[0] = other.myPar2[0]; + this->myPar2[1] = other.myPar2[1]; + this->myParIndex = other.myParIndex; + this->myFace2Surface = other.myFace2Surface; +} + //======================================================================= //function : ShapeToIndex //purpose : Convert a shape to its index in the SMESHDS_Mesh @@ -626,22 +642,22 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, { gp_Pnt2d uv( Precision::Infinite(), Precision::Infinite() ); - const SMDS_PositionPtr Pos = n->GetPosition(); + SMDS_PositionPtr pos = n->GetPosition(); bool uvOK = false; - if ( Pos->GetTypeOfPosition() == SMDS_TOP_FACE ) + if ( pos->GetTypeOfPosition() == SMDS_TOP_FACE ) { // node has position on face - const SMDS_FacePosition* fpos = static_cast( Pos ); + SMDS_FacePositionPtr fpos = pos; uv.SetCoord( fpos->GetUParameter(), fpos->GetVParameter() ); if ( check ) uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 2.*getFaceMaxTol( F )); // 2. from 22830 } - else if ( Pos->GetTypeOfPosition() == SMDS_TOP_EDGE ) + else if ( pos->GetTypeOfPosition() == SMDS_TOP_EDGE ) { // node has position on EDGE => it is needed to find // corresponding EDGE from FACE, get pcurve for this // EDGE and retrieve value from this pcurve - const SMDS_EdgePosition* epos = static_cast( Pos ); + SMDS_EdgePositionPtr epos = pos; const int edgeID = n->getshapeId(); const TopoDS_Edge& E = TopoDS::Edge( GetMeshDS()->IndexToShape( edgeID )); double f, l, u = epos->GetUParameter(); @@ -682,7 +698,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, uv = newUV; } } - else if ( Pos->GetTypeOfPosition() == SMDS_TOP_VERTEX ) + else if ( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX ) { if ( int vertexID = n->getshapeId() ) { const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID)); @@ -723,7 +739,8 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, if ( !C2d.IsNull() ) { double u = ( V == IthVertex( 0, edge )) ? f : l; uv = C2d->Value( u ); - uvOK = true; + gp_Pnt p = GetSurface( F )->Value( uv ); + uvOK = ( p.Distance( BRep_Tool::Pnt( V )) < getFaceMaxTol( F )); break; } } @@ -787,33 +804,20 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, // check that uv is correct TopLoc_Location loc; Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc ); - gp_Pnt nodePnt = XYZ( n ), surfPnt(0,0,0); + SMESH_NodeXYZ nXYZ( n ); + gp_Pnt nodePnt = nXYZ, surfPnt(0,0,0); double dist = 0; if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() ); if ( infinit || (dist = nodePnt.Distance( surfPnt = surface->Value( uv.X(), uv.Y() ))) > tol ) { setPosOnShapeValidity( shapeID, false ); - if ( !infinit && distXYZ ) { - surfPnt.Transform( loc ); - distXYZ[0] = dist; - distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z(); - } // uv incorrect, project the node to surface - GeomAPI_ProjectPointOnSurf& projector = GetProjector( F, loc, tol ); - projector.Perform( nodePnt ); - if ( !projector.IsDone() || projector.NbPoints() < 1 ) - { - MESSAGE( "SMESH_MesherHelper::CheckNodeUV() failed to project" ); - return false; - } - Quantity_Parameter U,V; - projector.LowerDistanceParameters(U,V); - uv.SetCoord( U,V ); - surfPnt = surface->Value( U, V ); - dist = nodePnt.Distance( surfPnt ); + Handle(ShapeAnalysis_Surface) sprojector = GetSurface( F ); + uv = sprojector->ValueOfUV( nXYZ, tol ).XY(); + surfPnt = sprojector->Value( uv ); + dist = surfPnt.Distance( nXYZ ); if ( distXYZ ) { - surfPnt.Transform( loc ); distXYZ[0] = dist; distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z(); } @@ -825,7 +829,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, // store the fixed UV on the face if ( myShape.IsSame(F) && shapeID == myShapeID && myFixNodeParameters ) const_cast(n)->SetPosition - ( SMDS_PositionPtr( new SMDS_FacePosition( U, V ))); + ( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() ))); } else if ( myShape.IsSame(F) && uv.Modulus() > numeric_limits::min() ) { @@ -837,7 +841,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, //======================================================================= //function : GetProjector -//purpose : Return projector intitialized by given face without location, which is returned +//purpose : Return projector initialized by given face without location, which is returned //======================================================================= GeomAPI_ProjectPointOnSurf& SMESH_MesherHelper::GetProjector(const TopoDS_Face& F, @@ -1014,8 +1018,7 @@ double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge& E, const SMDS_PositionPtr pos = n->GetPosition(); if ( pos->GetTypeOfPosition()==SMDS_TOP_EDGE ) { - const SMDS_EdgePosition* epos = static_cast( pos ); - param = epos->GetUParameter(); + param = pos->GetParameters()[0]; } else if( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX ) { @@ -1032,6 +1035,16 @@ double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge& E, int vertexID = n->getshapeId(); const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID)); param = BRep_Tool::Parameter( V, E ); + + if ( inEdgeNode ) + { + BRepAdaptor_Curve curve( E ); + if ( curve.IsPeriodic() ) + { + double uInEdge = GetNodeU( E, inEdgeNode ); + param += ShapeAnalysis::AdjustByPeriod( param, uInEdge, curve.Period() ); + } + } } } if ( check ) @@ -1113,7 +1126,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, MESSAGE( "SMESH_MesherHelper::CheckNodeU() failed to project" ); return false; } - Quantity_Parameter U = projector->LowerDistanceParameter(); + Standard_Real U = projector->LowerDistanceParameter(); u = double( U ); curvPnt = curve->Value( u ); dist = nodePnt.Distance( curvPnt ); @@ -1659,30 +1672,8 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, // get positions of the given nodes on shapes if ( pos.second == TopAbs_FACE ) { - F = TopoDS::Face(meshDS->IndexToShape( faceID = pos.first )); + F = TopoDS::Face( meshDS->IndexToShape( faceID = pos.first )); uv[0] = GetNodeUV(F,n1,n2, force3d ? 0 : &uvOK[0]); - if (( !force3d ) && - ( HasDegeneratedEdges() || GetSurface( F )->HasSingularities( 1e-7 ))) - { - // IPAL52850 (degen VERTEX not at singularity) - // project middle point to a surface - SMESH_TNodeXYZ p1( n1 ), p2( n2 ); - gp_Pnt pMid = 0.5 * ( p1 + p2 ); - Handle(ShapeAnalysis_Surface) projector = GetSurface( F ); - gp_Pnt2d uvMid; - if ( uvOK[0] ) - uvMid = projector->NextValueOfUV( uv[0], pMid, BRep_Tool::Tolerance( F )); - else - uvMid = projector->ValueOfUV( pMid, getFaceMaxTol( F )); - if ( projector->Gap() * projector->Gap() < ( p1 - p2 ).SquareModulus() / 4 ) - { - gp_Pnt pProj = projector->Value( uvMid ); - n12 = meshDS->AddNode( pProj.X(), pProj.Y(), pProj.Z() ); - meshDS->SetNodeOnFace( n12, faceID, uvMid.X(), uvMid.Y() ); - myTLinkNodeMap.insert( make_pair ( link, n12 )); - return n12; - } - } uv[1] = GetNodeUV(F,n2,n1, force3d ? 0 : &uvOK[1]); } else if ( pos.second == TopAbs_EDGE ) @@ -1715,26 +1706,43 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, // nodes, else - medium between corresponding 3d points if( ! F.IsNull() ) { - //if ( uvOK[0] && uvOK[1] ) + if ( IsDegenShape( n1->getshapeId() )) { + if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 )); + else uv[0].SetCoord( 2, uv[1].Coord( 2 )); + } + else if ( IsDegenShape( n2->getshapeId() )) { + if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 )); + else uv[1].SetCoord( 2, uv[0].Coord( 2 )); + } + TopLoc_Location loc; + Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc); + gp_XY UV = GetMiddleUV( S, uv[0], uv[1] ); + gp_Pnt P = S->Value( UV.X(), UV.Y() ).Transformed(loc); + + SMESH_TNodeXYZ p1( n1 ), p2( n2 ); + gp_Pnt pMid = 0.5 * ( p1 + p2 ); + double distMid = pMid.SquareDistance( P ); + double dist12 = ( p1 - p2 ).SquareModulus(); + Handle(ShapeAnalysis_Surface) surfInfo = GetSurface( F ); + if ( distMid > dist12 || + HasDegeneratedEdges() || + surfInfo->HasSingularities( 1e-7 ) ) { - if ( IsDegenShape( n1->getshapeId() )) { - if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 )); - else uv[0].SetCoord( 2, uv[1].Coord( 2 )); - } - else if ( IsDegenShape( n2->getshapeId() )) { - if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 )); - else uv[1].SetCoord( 2, uv[0].Coord( 2 )); - } - TopLoc_Location loc; - Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc); - gp_XY UV = GetMiddleUV( S, uv[0], uv[1] ); - gp_Pnt P = S->Value( UV.X(), UV.Y() ).Transformed(loc); - n12 = meshDS->AddNode(P.X(), P.Y(), P.Z()); - // if ( mySetElemOnShape ) node is not elem! - meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y()); - myTLinkNodeMap.insert(make_pair(link,n12)); - return n12; + // IPAL52850 (degen VERTEX not at singularity) + // project middle point to a surface + gp_Pnt2d uvMid; + if ( uvOK[0] ) + uvMid = surfInfo->NextValueOfUV( uv[0], pMid, BRep_Tool::Tolerance( F )); + else + uvMid = surfInfo->ValueOfUV( pMid, getFaceMaxTol( F )); + if ( surfInfo->Gap() * surfInfo->Gap() < distMid ) + P = surfInfo->Value( uvMid ); } + n12 = meshDS->AddNode(P.X(), P.Y(), P.Z()); + // if ( mySetElemOnShape ) node is not elem! + meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y()); + myTLinkNodeMap.insert(make_pair(link,n12)); + return n12; } else if ( !E.IsNull() ) { @@ -1889,7 +1897,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_ if ( !bestEdge.IsNull() ) { - // move n12 to position of a successfull projection + // move n12 to position of a successful projection //double tol = BRep_Tool::Tolerance(edges[ iOkEdge ]); if ( !force3d /*&& distMiddleProj > 2*tol*/ ) { @@ -1916,7 +1924,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_ //purpose : Creates a node //======================================================================= -SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID, +SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, smIdType ID, double u, double v) { SMESHDS_Mesh * meshDS = GetMeshDS(); @@ -1945,7 +1953,7 @@ SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID, SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, - const int id, + const smIdType id, const bool force3d) { SMESHDS_Mesh * meshDS = GetMeshDS(); @@ -1979,7 +1987,7 @@ SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1, SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshNode* n3, - const int id, + const smIdType id, const bool force3d) { SMESHDS_Mesh * meshDS = GetMeshDS(); @@ -2029,7 +2037,7 @@ SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshNode* n3, const SMDS_MeshNode* n4, - const int id, + const smIdType id, const bool force3d) { SMESHDS_Mesh * meshDS = GetMeshDS(); @@ -2093,7 +2101,7 @@ SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1, //======================================================================= SMDS_MeshFace* SMESH_MesherHelper::AddPolygonalFace (const vector& nodes, - const int id, + const smIdType id, const bool force3d) { SMESHDS_Mesh * meshDS = GetMeshDS(); @@ -2139,7 +2147,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, const SMDS_MeshNode* n4, const SMDS_MeshNode* n5, const SMDS_MeshNode* n6, - const int id, + const smIdType id, const bool force3d) { SMESHDS_Mesh * meshDS = GetMeshDS(); @@ -2162,13 +2170,30 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, const SMDS_MeshNode* n14 = GetMediumNode( n1, n4, force3d, TopAbs_SOLID ); const SMDS_MeshNode* n25 = GetMediumNode( n2, n5, force3d, TopAbs_SOLID ); const SMDS_MeshNode* n36 = GetMediumNode( n3, n6, force3d, TopAbs_SOLID ); + if ( myCreateBiQuadratic ) + { + const SMDS_MeshNode* n1245 = GetCentralNode( n1,n2,n4,n5,n12,n25,n45,n14,force3d ); + const SMDS_MeshNode* n1346 = GetCentralNode( n1,n3,n4,n6,n31,n36,n64,n14,force3d ); + const SMDS_MeshNode* n2356 = GetCentralNode( n2,n3,n6,n5,n23,n36,n56,n25,force3d ); - if(id) - elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, - n12, n23, n31, n45, n56, n64, n14, n25, n36, id); + if(id) + elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, + n12, n23, n31, n45, n56, n64, n14, n25, n36, + n1245, n2356, n1346, id); + else + elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, + n12, n23, n31, n45, n56, n64, n14, n25, n36, + n1245, n2356, n1346); + } else - elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, - n12, n23, n31, n45, n56, n64, n14, n25, n36); + { + if(id) + elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, + n12, n23, n31, n45, n56, n64, n14, n25, n36, id); + else + elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, + n12, n23, n31, n45, n56, n64, n14, n25, n36); + } } if ( mySetElemOnShape && myShapeID > 0 ) meshDS->SetMeshElementOnShape( elem, myShapeID ); @@ -2185,7 +2210,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshNode* n3, const SMDS_MeshNode* n4, - const int id, + const smIdType id, const bool force3d) { SMESHDS_Mesh * meshDS = GetMeshDS(); @@ -2226,7 +2251,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, const SMDS_MeshNode* n3, const SMDS_MeshNode* n4, const SMDS_MeshNode* n5, - const int id, + const smIdType id, const bool force3d) { SMDS_MeshVolume* elem = 0; @@ -2276,7 +2301,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, const SMDS_MeshNode* n6, const SMDS_MeshNode* n7, const SMDS_MeshNode* n8, - const int id, + const smIdType id, const bool force3d) { SMESHDS_Mesh * meshDS = GetMeshDS(); @@ -2395,8 +2420,8 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, const SMDS_MeshNode* n10, const SMDS_MeshNode* n11, const SMDS_MeshNode* n12, - const int id, - bool force3d) + const smIdType id, + bool /*force3d*/) { SMESHDS_Mesh * meshDS = GetMeshDS(); SMDS_MeshVolume* elem = 0; @@ -2417,7 +2442,7 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, SMDS_MeshVolume* SMESH_MesherHelper::AddPolyhedralVolume (const std::vector& nodes, const std::vector& quantities, - const int id, + const smIdType id, const bool force3d) { SMESHDS_Mesh * meshDS = GetMeshDS(); @@ -2568,6 +2593,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 } // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them + const SMDS_MeshNode* prevEndNodes[2] = { 0, 0 }; edge = theBaseSide.begin(); for ( int iE = 0; edge != theBaseSide.end(); ++edge, ++iE ) { @@ -2635,18 +2661,23 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 const double prevPar = theParam2ColumnMap.empty() ? 0 : theParam2ColumnMap.rbegin()->first; for ( u_n = sortedBaseNN.begin(); u_n != sortedBaseNN.end(); u_n++ ) { + if ( u_n->second == prevEndNodes[0] || + u_n->second == prevEndNodes[1] ) + continue; double par = prevPar + coeff * ( u_n->first - f ); TParam2ColumnMap::iterator u2nn = theParam2ColumnMap.insert( theParam2ColumnMap.end(), make_pair( par, TNodeColumn())); u2nn->second.push_back( u_n->second ); } + prevEndNodes[0] = sortedBaseNN.begin()->second; + prevEndNodes[1] = sortedBaseNN.rbegin()->second; } if ( theParam2ColumnMap.size() < 2 ) return false; } // nb rows of nodes - size_t prevNbRows = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here + size_t prevNbRows = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here size_t expectNbRows = faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 ); // to be added // fill theParam2ColumnMap column by column by passing from nodes on @@ -2819,8 +2850,9 @@ bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM ) //purpose : Return true if 2D mesh on FACE is ditorted //======================================================================= -bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, - bool checkUV) +bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, + bool checkUV, + SMESH_MesherHelper* faceHelper) { if ( !faceSM || faceSM->GetSubShape().ShapeType() != TopAbs_FACE ) return false; @@ -2828,12 +2860,26 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, bool haveBadFaces = false; SMESH_MesherHelper helper( *faceSM->GetFather() ); + if ( faceHelper ) + helper.CopySubShapeInfo( *faceHelper ); helper.SetSubShape( faceSM->GetSubShape() ); const TopoDS_Face& F = TopoDS::Face( faceSM->GetSubShape() ); SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( F ); if ( !smDS || smDS->NbElements() == 0 ) return false; + bool subIdsValid = true; // shape ID of nodes is OK + if ( helper.HasSeam() ) + { + // check if nodes are bound to seam edges + SMESH_subMeshIteratorPtr smIt = faceSM->getDependsOnIterator(/*includeSelf=*/false); + while ( smIt->more() && subIdsValid ) + { + SMESH_subMesh* sm = smIt->next(); + if ( helper.IsSeamShape( sm->GetId() ) && sm->IsEmpty() ) + subIdsValid = false; + } + } SMDS_ElemIteratorPtr faceIt = smDS->GetElements(); double prevArea = 0; vector< const SMDS_MeshNode* > nodes; @@ -2849,7 +2895,7 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, for ( size_t i = 0; i < nodes.size(); ++n, ++i ) nodes[ i ] = *n; - // avoid elems on degenarate shapes as UV on them can be wrong + // avoid elems on degenerate shapes as UV on them can be wrong if ( helper.HasDegeneratedEdges() ) { bool isOnDegen = false; @@ -2858,12 +2904,20 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, if ( isOnDegen ) continue; } - // prepare to getting UVs + // prepare for getting UVs const SMDS_MeshNode* inFaceNode = 0; if ( helper.HasSeam() ) { for ( size_t i = 0; ( i < nodes.size() && !inFaceNode ); ++i ) if ( !helper.IsSeamShape( nodes[ i ]->getshapeId() )) + { inFaceNode = nodes[ i ]; + if ( !subIdsValid ) + { + gp_XY uv = helper.GetNodeUV( F, inFaceNode ); + if ( helper.IsOnSeam( uv )) + inFaceNode = NULL; + } + } if ( !inFaceNode ) continue; } @@ -2872,6 +2926,14 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, for ( size_t i = 0; i < nodes.size(); ++i ) uv[ i ] = helper.GetNodeUV( F, nodes[ i ], inFaceNode, toCheckUV ); + if ( !subIdsValid ) // fix uv on seam + { + gp_XY uvInFace = helper.GetNodeUV( F, inFaceNode ); + for ( size_t i = 0; i < uv.size(); ++i ) + if ( helper.IsOnSeam( uv[i] )) + uv[i] = helper.getUVOnSeam( uv[i], uvInFace ).XY(); + } + // compare orientation of triangles double faceArea = 0; for ( int iT = 0, nbT = nodes.size()-2; iT < nbT; ++iT ) @@ -2892,7 +2954,7 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, * \brief Find out elements orientation on a geometrical face * \param theFace - The face correctly oriented in the shape being meshed * \retval bool - true if the face normal and the normal of first element - * in the correspoding submesh point in different directions + * in the corresponding submesh point in different directions */ //================================================================================ @@ -2910,7 +2972,7 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace) if ( !aSubMeshDSFace ) return isReversed; - // find an element on a bounday of theFace + // find an element on a boundary of theFace SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); const SMDS_MeshNode* nn[2]; while ( iteratorElem->more() ) // loop on elements on theFace @@ -3081,7 +3143,7 @@ TopAbs_Orientation SMESH_MesherHelper::GetSubShapeOri(const TopoDS_Shape& shape, //======================================================================= //function : IsSubShape -//purpose : +//purpose : //======================================================================= bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, @@ -3179,7 +3241,7 @@ double SMESH_MesherHelper::getFaceMaxTol( const TopoDS_Shape& face ) const * of the FACE normal * \return double - the angle (between -Pi and Pi), negative if the angle is concave, * 1e100 in case of failure - * \waring Care about order of the EDGEs and their orientation to be as they are + * \warning Care about order of the EDGEs and their orientation to be as they are * within the FACE! Don't pass degenerated EDGEs neither! */ //================================================================================ @@ -3351,9 +3413,9 @@ TopoDS_Shape SMESH_MesherHelper::GetShapeOfHypothesis( const SMESHDS_Hypothesis SMESH_MesherHelper:: MType SMESH_MesherHelper::IsQuadraticMesh() { - int NbAllEdgsAndFaces=0; - int NbQuadFacesAndEdgs=0; - int NbFacesAndEdges=0; + smIdType NbAllEdgsAndFaces=0; + smIdType NbQuadFacesAndEdgs=0; + smIdType NbFacesAndEdges=0; //All faces and edges NbAllEdgsAndFaces = myMesh->NbEdges() + myMesh->NbFaces(); if ( NbAllEdgsAndFaces == 0 ) @@ -3389,6 +3451,43 @@ double SMESH_MesherHelper::GetOtherParam(const double param) const return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i]; } +//======================================================================= +//function : NbRealSeam +//purpose : Return a number of real seam edges in the shape set through +// IsQuadraticSubMesh() or SetSubShape(). A real seam edge encounters twice in a wire +//======================================================================= + +size_t SMESH_MesherHelper::NbRealSeam() const +{ + size_t nb = 0; + + std::set< int >::const_iterator id = mySeamShapeIds.begin(); + for ( ; id != mySeamShapeIds.end(); ++id ) + if ( *id < 0 ) ++nb; + else break; + + return nb; +} + +//======================================================================= +//function : IsOnSeam +//purpose : Check if UV is on seam. Return 0 if not, 1 for U seam, 2 for V seam +//======================================================================= + +int SMESH_MesherHelper::IsOnSeam(const gp_XY& uv) const +{ + for ( int i = U_periodic; i <= V_periodic ; ++i ) + if ( myParIndex & i ) + { + double p = uv.Coord( i ); + double tol = ( myPar2[i-1] - myPar1[i-1] ) / 100.; + if ( Abs( p - myPar1[i-1] ) < tol || + Abs( p - myPar2[i-1] ) < tol ) + return i; + } + return 0; +} + namespace { //======================================================================= @@ -3402,11 +3501,16 @@ namespace { TopTools_ListIteratorOfListOfShape _ancIter; TopAbs_ShapeEnum _type; TopTools_MapOfShape _encountered; - TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type) + TopTools_IndexedMapOfShape _allowed; + TAncestorsIterator( const TopTools_ListOfShape& ancestors, + TopAbs_ShapeEnum type, + const TopoDS_Shape* container/* = 0*/) : _ancIter( ancestors ), _type( type ) { + if ( container && !container->IsNull() ) + TopExp::MapShapes( *container, type, _allowed); if ( _ancIter.More() ) { - if ( _ancIter.Value().ShapeType() != _type ) next(); + if ( !isCurrentAllowed() ) next(); else _encountered.Add( _ancIter.Value() ); } } @@ -3419,25 +3523,32 @@ namespace { const TopoDS_Shape* s = _ancIter.More() ? & _ancIter.Value() : 0; if ( _ancIter.More() ) for ( _ancIter.Next(); _ancIter.More(); _ancIter.Next()) - if ( _ancIter.Value().ShapeType() == _type && _encountered.Add( _ancIter.Value() )) + if ( isCurrentAllowed() && _encountered.Add( _ancIter.Value() )) break; return s; } + bool isCurrentAllowed() + { + return (( _ancIter.Value().ShapeType() == _type ) && + ( _allowed.IsEmpty() || _allowed.Contains( _ancIter.Value() ))); + } }; } // namespace //======================================================================= /*! - * \brief Return iterator on ancestors of the given type + * \brief Return iterator on ancestors of the given type, included into a container shape */ //======================================================================= PShapeIteratorPtr SMESH_MesherHelper::GetAncestors(const TopoDS_Shape& shape, const SMESH_Mesh& mesh, - TopAbs_ShapeEnum ancestorType) + TopAbs_ShapeEnum ancestorType, + const TopoDS_Shape* container) { - return PShapeIteratorPtr( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType)); + return PShapeIteratorPtr + ( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType, container)); } //======================================================================= @@ -3680,7 +3791,7 @@ namespace { // Structures used by FixQuadraticElements() //================================================================================ /*! - * \brief Construct QFace from QLinks + * \brief Construct QFace from QLinks */ //================================================================================ @@ -3697,7 +3808,7 @@ namespace { // Structures used by FixQuadraticElements() gp_Vec v1( XYZ( l1->node2()), XYZ( l1->node1())); gp_Vec v2( XYZ( l2->node1()), XYZ( l2->node2())); if ( l1->node1() != l2->node1() && l1->node2() != l2->node2() ) - v1.Reverse(); + v1.Reverse(); _normal += v1 ^ v2; } double normSqSize = _normal.SquareMagnitude(); @@ -3708,6 +3819,8 @@ namespace { // Structures used by FixQuadraticElements() #ifdef _DEBUG_ _face = face; +#else + (void)face; // unused in release mode #endif } //================================================================================ @@ -3715,7 +3828,7 @@ namespace { // Structures used by FixQuadraticElements() * \brief Make up a chain of links * \param iSide - link to add first * \param chain - chain to fill in - * \param pos - postion of medium nodes the links should have + * \param pos - position of medium nodes the links should have * \param error - out, specifies what is wrong * \retval bool - false if valid chain can't be built; "valid" means that links * of the chain belongs to rectangles bounding hexahedrons @@ -3729,7 +3842,7 @@ namespace { // Structures used by FixQuadraticElements() if ( _sideIsAdded[ iSide ]) // already in chain return true; - if ( _sides.size() != 4 ) { // triangle - visit all my continous faces + if ( _sides.size() != 4 ) { // triangle - visit all my continuous faces MSGBEG( *this ); TLinkSet links; list< const QFace* > faces( 1, this ); @@ -3775,7 +3888,7 @@ namespace { // Structures used by FixQuadraticElements() if ( link->MediumPos() >= pos ) { int nbLinkFaces = link->_faces.size(); if ( nbLinkFaces == 4 || (/*nbLinkFaces < 4 && */link->OnBoundary())) { - // hexahedral mesh or boundary quadrangles - goto a continous face + // hexahedral mesh or boundary quadrangles - goto a continuous face if ( const QFace* f = link->GetContinuesFace( this )) if ( f->_sides.size() == 4 ) return f->GetLinkChain( *chLink, chain, pos, error ); @@ -3903,7 +4016,7 @@ namespace { // Structures used by FixQuadraticElements() * \brief Move medium node of theLink according to its distance from boundary * \param theLink - link to fix * \param theRefVec - movement of boundary - * \param theLinks - all adjacent links of continous triangles + * \param theLinks - all adjacent links of continuous triangles * \param theFaceHelper - helper is not used so far * \param thePrevLen - distance from the boundary * \param theStep - number of steps till movement propagation limit @@ -4019,34 +4132,34 @@ namespace { // Structures used by FixQuadraticElements() */ //================================================================================ - bool QFace::IsSpoiled(const QLink* bentLink ) const - { - // code is valid for convex faces only - gp_XYZ gc(0,0,0); - for ( TIDSortedNodeSet::const_iterator n = begin(); n!=end(); ++n) - gc += XYZ( *n ) / size(); - for (unsigned i = 0; i < _sides.size(); ++i ) - { - if ( _sides[i] == bentLink ) continue; - gp_Vec linkNorm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2())); - gp_Vec vecOut( gc, _sides[i]->MiddlePnt() ); - if ( linkNorm * vecOut < 0 ) - linkNorm.Reverse(); - double mag2 = linkNorm.SquareMagnitude(); - if ( mag2 > numeric_limits::min() ) - linkNorm /= sqrt( mag2 ); - gp_Vec vecBent ( _sides[i]->MiddlePnt(), bentLink->MediumPnt()); - gp_Vec vecStraight( _sides[i]->MiddlePnt(), bentLink->MiddlePnt()); - if ( vecBent * linkNorm > -0.1*vecStraight.Magnitude() ) - return true; - } - return false; - - } + // bool QFace::IsSpoiled(const QLink* bentLink ) const + // { + // // code is valid for convex faces only + // gp_XYZ gc(0,0,0); + // for ( TIDSortedNodeSet::const_iterator n = begin(); n != end(); ++n ) + // gc += XYZ( *n ) / double( size() ); + // for ( size_t i = 0; i < _sides.size(); ++i ) + // { + // if ( _sides[i] == bentLink ) continue; + // gp_Vec linkNorm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2())); + // gp_Vec vecOut( gc, _sides[i]->MiddlePnt() ); + // if ( linkNorm * vecOut < 0 ) + // linkNorm.Reverse(); + // double mag2 = linkNorm.SquareMagnitude(); + // if ( mag2 > numeric_limits::min() ) + // linkNorm /= sqrt( mag2 ); + // gp_Vec vecBent ( _sides[i]->MiddlePnt(), bentLink->MediumPnt()); + // gp_Vec vecStraight( _sides[i]->MiddlePnt(), bentLink->MiddlePnt()); + // if ( vecBent * linkNorm > -0.1*vecStraight.Magnitude() ) + // return true; + // } + // return false; + + // } //================================================================================ /*! - * \brief Find pairs of continues faces + * \brief Find pairs of continues faces */ //================================================================================ @@ -4057,7 +4170,7 @@ namespace { // Structures used by FixQuadraticElements() // | Between _faces of link x2 two vertical faces are continues // x1----x2-----x3 and two horizontal faces are continues. We set vertical faces // | to _faces[0] and _faces[1] and horizontal faces to - // v2 | v3 _faces[2] and _faces[3] (or vise versa). + // v2 | v3 _faces[2] and _faces[3] (or vice versa). // x4 if ( _faces.empty() ) @@ -4258,7 +4371,7 @@ namespace { // Structures used by FixQuadraticElements() vector< TChain> & resultChains, SMDS_TypeOfPosition pos ) { - // put links in the set and evalute number of result chains by number of boundary links + // put links in the set and evaluate number of result chains by number of boundary links TLinkSet linkSet; size_t nbBndLinks = 0; for ( TChain::iterator lnk = allLinks.begin(); lnk != allLinks.end(); ++lnk ) { @@ -4442,8 +4555,16 @@ namespace { // Structures used by FixQuadraticElements() TopoDS_Shape shape = theHelper.GetSubShape().Oriented( TopAbs_FORWARD ); if ( shape.IsNull() ) return; - if ( !theError ) theError = SMESH_ComputeError::New(); - + if ( !dynamic_cast( theError.get() )) + { + if ( !theError ) + theError.reset( new SMESH_BadInputElements( meshDS )); + else + theError.reset( new SMESH_BadInputElements( meshDS, + theError->myName, + theError->myComment, + theError->myAlgo)); + } gp_XYZ faceNorm; if ( shape.ShapeType() == TopAbs_FACE ) // 2D @@ -4494,7 +4615,7 @@ namespace { // Structures used by FixQuadraticElements() if ( curvNorm * D2 > 0 ) continue; // convex edge } - catch ( Standard_Failure ) + catch ( Standard_Failure& ) { continue; } @@ -4526,6 +4647,7 @@ namespace { // Structures used by FixQuadraticElements() const SMDS_MeshElement* f = faceIt->next(); if ( !faceSM->Contains( f ) || f->NbNodes() < 6 || // check quadratic triangles only + f->NbNodes() > 7 || !checkedFaces.insert( f ).second ) continue; @@ -4561,13 +4683,13 @@ namespace { // Structures used by FixQuadraticElements() gp_XYZ pMid3D = 0.5 * ( pN0 + SMESH_TNodeXYZ( nOnEdge[1] )); meshDS->MoveNode( n, pMid3D.X(), pMid3D.Y(), pMid3D.Z() ); MSG( "move OUT of face " << n ); - theError->myBadElements.push_back( f ); + static_cast( theError.get() )->add( f ); } } } } } - if ( !theError->myBadElements.empty() ) + if ( theError->HasBadElems() ) theError->myName = EDITERR_NO_MEDIUM_ON_GEOM; return; @@ -4608,7 +4730,7 @@ namespace { // Structures used by FixQuadraticElements() if ( concaveU || concaveV ) concaveFaces.push_back( face ); } - catch ( Standard_Failure ) + catch ( Standard_Failure& ) { concaveFaces.push_back( face ); } @@ -4634,7 +4756,7 @@ namespace { // Structures used by FixQuadraticElements() < const SMDS_MeshElement*, vector< SMDS_ElemIteratorPtr > > TIterOnIter; SMDS_ElemIteratorPtr faceIter( new TIterOnIter( faceIterVec )); - // a seacher to check if a volume is close to a concave face + // search to check if a volume is close to a concave face SMESHUtils::Deleter< SMESH_ElementSearcher > faceSearcher ( SMESH_MeshAlgos::GetElementSearcher( *theHelper.GetMeshDS(), faceIter )); @@ -4737,7 +4859,12 @@ namespace { // Structures used by FixQuadraticElements() gp_Pnt pMedium = SMESH_TNodeXYZ( linkIt->second ); double hMedium = faceNorm * gp_Vec( pOnFace0, pMedium ).XYZ(); double hVol = faceNorm * gp_Vec( pOnFace0, pInSolid ).XYZ(); - isDistorted = ( Abs( hMedium ) > Abs( hVol * 0.75 )); + if ( Abs( hMedium ) > Abs( hVol * 0.75 )) + { + SMESH_TNodeXYZ pI( nOnFace[i]), pJ( nOnFace[j]); + double angle = gp_Vec( pI, pMedium ).Angle( gp_Vec( pI, pJ )); + isDistorted = ( angle > M_PI / 20 ); + } } } } @@ -4755,13 +4882,13 @@ namespace { // Structures used by FixQuadraticElements() MSG( "move OUT of solid " << nMedium ); } } - theError->myBadElements.push_back( vol ); + static_cast( theError.get() )->add( vol ); } } // loop on volumes sharing a node on FACE } // loop on nodes on FACE } // loop on FACEs of a SOLID - if ( !theError->myBadElements.empty() ) + if ( theError->HasBadElems() ) theError->myName = EDITERR_NO_MEDIUM_ON_GEOM; } // 3D case } @@ -4773,7 +4900,7 @@ namespace { // Structures used by FixQuadraticElements() * \brief Move medium nodes of faces and volumes to fix distorted elements * \param error - container of fixed distorted elements * \param volumeOnly - to fix nodes on faces or not, if the shape is solid - * + * * Issue 0020307: EDF 992 SMESH : Linea/Quadratic with Medium Node on Geometry */ //======================================================================= @@ -4781,6 +4908,7 @@ namespace { // Structures used by FixQuadraticElements() void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, bool volumeOnly) { + //MESSAGE("FixQuadraticElements " << volumeOnly); // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion if ( getenv("NO_FixQuadraticElements") ) return; @@ -4822,6 +4950,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, nbfaces = faces.Extent(); /*avoid "unused varianbles": */ nbfaces++, nbfaces--; #endif for ( TopTools_MapIteratorOfMapOfShape fIt( faces ); fIt.More(); fIt.Next() ) { + MESSAGE("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key())); MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key())); SMESH_MesherHelper h(*myMesh); h.SetSubShape( fIt.Key() ); @@ -5022,7 +5151,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, MSG("Internal chain - ignore"); continue; } - // mesure chain length and compute link position along the chain + // measure chain length and compute link position along the chain double chainLen = 0; vector< double > linkPos; TChain savedChain; // backup @@ -5133,7 +5262,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, try { gp_Vec x = x01.Normalized() + x12.Normalized(); trsf.SetTransformation( gp_Ax3( gp::Origin(), link1->Normal(), x), gp_Ax3() ); - } catch ( Standard_Failure ) { + } catch ( Standard_Failure& ) { trsf.Invert(); } move.Transform(trsf); @@ -5146,8 +5275,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, gp_XY newUV = ApplyIn2D( s, oldUV, gp_XY( move.X(),move.Y()), gp_XY_Added ); gp_Pnt newPnt = s->Value( newUV.X(), newUV.Y()); move = gp_Vec( XYZ((*link1)->_mediumNode), newPnt.Transformed(loc) ); - if ( SMDS_FacePosition* nPos = - dynamic_cast< SMDS_FacePosition* >((*link1)->_mediumNode->GetPosition())) + if ( SMDS_FacePositionPtr nPos = (*link1)->_mediumNode->GetPosition()) nPos->SetParameters( newUV.X(), newUV.Y() ); #ifdef _DEBUG_ if ( (XYZ((*link1)->node1()) - XYZ((*link1)->node2())).SquareModulus() < @@ -5178,10 +5306,11 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, // 4. Move nodes // ------------- - TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa; + TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa, biQuadPenta; const bool toFixCentralNodes = ( myMesh->NbBiQuadQuadrangles() + myMesh->NbBiQuadTriangles() + - myMesh->NbTriQuadraticHexas() ); + myMesh->NbTriQuadraticHexas() + + myMesh->NbBiQuadPrisms()); double distXYZ[4]; faceHlp.ToFixNodeParameters( true ); @@ -5217,6 +5346,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, case SMDSEntity_BiQuad_Quadrangle: biQuadQuas.insert( e ); break; case SMDSEntity_BiQuad_Triangle: biQuadTris.insert( e ); break; case SMDSEntity_TriQuad_Hexa: triQuadHexa.insert( e ); break; + case SMDSEntity_BiQuad_Penta: biQuadPenta.insert( e ); break; default:; } } @@ -5367,10 +5497,22 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, nCenterCoords.X(), nCenterCoords.Y(), nCenterCoords.Z()); } } -#ifdef _DEBUG_ - // avoid warning: defined but not used operator<<() - SMESH_Comment() << *links.begin() << *faces.begin(); -#endif + // treat tri-quadratic hexahedra + { + SMDS_VolumeTool volExp; + TIDSortedElemSet::iterator pentIt = biQuadPenta.begin(); + for ( ; pentIt != biQuadPenta.end(); ++pentIt ) + { + MESSAGE("---"); + volExp.Set( *pentIt, /*ignoreCentralNodes=*/false ); + } + } + + if ( false ) + // avoid warning: defined but not used operator<<() + SMESH_Comment() << *links.begin() << *faces.begin(); + + return; } //================================================================================