X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=242900a749d0877a3c4f9e25ace56a0c2fdf43eb;hp=35973eb99f8a336ef77d6f1f1cfa402d6bd6976d;hb=88141f757b048eaa5aae0be49faaf274448bbcaf;hpb=004925bca40d46c25cf9cc837193335163c2d44f diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index 35973eb99..242900a74 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -31,8 +31,10 @@ #include "SMDS_FacePosition.hxx" #include "SMDS_IteratorOnIterators.hxx" #include "SMDS_VolumeTool.hxx" +#include "SMESHDS_Mesh.hxx" #include "SMESH_Block.hxx" #include "SMESH_HypoFilter.hxx" +#include "SMESH_Mesh.hxx" #include "SMESH_MeshAlgos.hxx" #include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" @@ -113,11 +115,33 @@ SMESH_MesherHelper::~SMESH_MesherHelper() } } +//================================================================================ +/*! + * \brief Return SMESH_Gen + */ +//================================================================================ + +SMESH_Gen* SMESH_MesherHelper::GetGen() const +{ + return GetMesh()->GetGen(); +} + +//================================================================================ +/*! + * \brief Return mesh DS + */ +//================================================================================ + +SMESHDS_Mesh* SMESH_MesherHelper::GetMeshDS() const +{ + return GetMesh()->GetMeshDS(); +} + //======================================================================= //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) @@ -125,7 +149,6 @@ 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(); @@ -138,9 +161,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 ))) @@ -280,6 +300,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) double u2 = uv1.Coord(1); myPar1[0] = Min( u1, u2 ); myPar2[0] = Max( u1, u2 ); + myParIndex |= U_periodic; } else { @@ -289,6 +310,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) double v2 = uv1.Coord(2); myPar1[1] = Min( v1, v2 ); myPar2[1] = Max( v1, v2 ); + myParIndex |= V_periodic; } } else //if ( !isSeam ) @@ -347,6 +369,37 @@ 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 +//======================================================================= + +int SMESH_MesherHelper::ShapeToIndex( const TopoDS_Shape& S ) const +{ + return GetMeshDS()->ShapeToIndex( S ); +} + //======================================================================= //function : GetNodeUVneedInFaceNode //purpose : Check if inFaceNode argument is necessary for call GetNodeUV(F,..) @@ -771,7 +824,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, MESSAGE( "SMESH_MesherHelper::CheckNodeUV() failed to project" ); return false; } - Quantity_Parameter U,V; + Standard_Real U,V; projector.LowerDistanceParameters(U,V); uv.SetCoord( U,V ); surfPnt = surface->Value( U, V ); @@ -1077,7 +1130,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 ); @@ -1853,7 +1906,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*/ ) { @@ -2126,13 +2179,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 ); @@ -2532,6 +2602,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 ) { @@ -2599,18 +2670,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 @@ -2783,8 +2859,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; @@ -2792,12 +2869,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; @@ -2822,12 +2913,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; } @@ -2836,6 +2935,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 ) @@ -2856,7 +2963,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 */ //================================================================================ @@ -2894,7 +3001,7 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace) TopoDS_Shape s0 = GetSubShapeByNode( nn[0], GetMeshDS() ); TopoDS_Shape s1 = GetSubShapeByNode( nn[1], GetMeshDS() ); TopoDS_Shape E = GetCommonAncestor( s0, s1, *myMesh, TopAbs_EDGE ); - if ( !E.IsNull() && !s0.IsSame( s1 )) + if ( !E.IsNull() && !s0.IsSame( s1 ) && E.Orientation() != TopAbs_INTERNAL ) { // is E seam edge? int nb = 0; @@ -2908,10 +3015,16 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace) bool ok = true; double u0 = GetNodeU( TopoDS::Edge( E ), nn[0], nn[1], &ok ); double u1 = GetNodeU( TopoDS::Edge( E ), nn[1], nn[0], &ok ); - // check that the 2 nodes are connected with a segment (IPAL53055) - if ( SMESHDS_SubMesh* sm = GetMeshDS()->MeshElements( E )) - if ( sm->NbElements() > 0 && !GetMeshDS()->FindEdge( nn[0], nn[1] )) - ok = false; + if ( ok ) + { + // check that the 2 nodes are connected with a segment (IPAL53055) + ok = false; + const SMDS_MeshElement* seg; + if ( SMESHDS_SubMesh* sm = GetMeshDS()->MeshElements( E )) + if (( sm->NbElements() > 0 ) && + ( seg = GetMeshDS()->FindEdge( nn[0], nn[1] ))) + ok = sm->Contains( seg ); + } if ( ok ) { isReversed = ( u0 > u1 ); @@ -3039,7 +3152,7 @@ TopAbs_Orientation SMESH_MesherHelper::GetSubShapeOri(const TopoDS_Shape& shape, //======================================================================= //function : IsSubShape -//purpose : +//purpose : //======================================================================= bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, @@ -3137,7 +3250,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! */ //================================================================================ @@ -3347,6 +3460,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 { //======================================================================= @@ -3360,11 +3510,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() ); } } @@ -3377,25 +3532,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)); } //======================================================================= @@ -3673,7 +3835,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 @@ -3687,7 +3849,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 ); @@ -3733,7 +3895,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 ); @@ -3861,7 +4023,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 @@ -3981,9 +4143,9 @@ namespace { // Structures used by FixQuadraticElements() { // 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 ) + 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())); @@ -3999,12 +4161,12 @@ namespace { // Structures used by FixQuadraticElements() return true; } return false; - + } //================================================================================ /*! - * \brief Find pairs of continues faces + * \brief Find pairs of continues faces */ //================================================================================ @@ -4216,7 +4378,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 ) { @@ -4592,8 +4754,8 @@ 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 - std::auto_ptr< SMESH_ElementSearcher > faceSearcher + // search to check if a volume is close to a concave face + SMESHUtils::Deleter< SMESH_ElementSearcher > faceSearcher ( SMESH_MeshAlgos::GetElementSearcher( *theHelper.GetMeshDS(), faceIter )); // classifier @@ -4695,7 +4857,7 @@ 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.5 )); + isDistorted = ( Abs( hMedium ) > Abs( hVol * 0.75 )); } } } @@ -4739,6 +4901,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; @@ -4776,9 +4939,11 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, } // fix nodes on geom faces #ifdef _DEBUG_ - int nbfaces = faces.Extent(); /*avoid "unused varianbles": */ nbfaces++, nbfaces--; + int nbfaces = nbSolids; + 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() ); @@ -4979,7 +5144,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 @@ -5117,6 +5282,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, "uv2: "<Move( move, /*sum=*/false, /*is2dFixed=*/true ); @@ -5134,11 +5300,11 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, // 4. Move nodes // ------------- - TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa; - const SMDS_MeshElement *biQuadQua, *triQuadHex; + TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa, biQuadPenta; const bool toFixCentralNodes = ( myMesh->NbBiQuadQuadrangles() + myMesh->NbBiQuadTriangles() + - myMesh->NbTriQuadraticHexas() ); + myMesh->NbTriQuadraticHexas() + + myMesh->NbBiQuadPrisms()); double distXYZ[4]; faceHlp.ToFixNodeParameters( true ); @@ -5166,7 +5332,6 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, // collect bi-quadratic elements if ( toFixCentralNodes ) { - biQuadQua = triQuadHex = 0; SMDS_ElemIteratorPtr eIt = pLink->_mediumNode->GetInverseElementIterator(); while ( eIt->more() ) { @@ -5175,6 +5340,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:; } } @@ -5325,6 +5491,20 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, nCenterCoords.X(), nCenterCoords.Y(), nCenterCoords.Z()); } } + // treat tri-quadratic hexahedra + { + SMDS_VolumeTool volExp; + TIDSortedElemSet::iterator pentIt = biQuadPenta.begin(); + for ( ; pentIt != biQuadPenta.end(); ++pentIt ) + { + MESSAGE("---"); + volExp.Set( *pentIt, /*ignoreCentralNodes=*/false ); + } + } +#ifdef _DEBUG_ + // avoid warning: defined but not used operator<<() + SMESH_Comment() << *links.begin() << *faces.begin(); +#endif } //================================================================================