From 28cdc5e829db988c6d017224d4f939ce02f40186 Mon Sep 17 00:00:00 2001 From: eap Date: Wed, 29 Dec 2010 18:55:42 +0000 Subject: [PATCH] 0020832: EDF 1359 SMESH : Automatic meshing of boundary layers for StdMeshers_QuadToTriaAdaptor to work after StdMeshers_ViscousLayers --- src/StdMeshers/StdMeshers_ProxyMesh.cxx | 141 ++++- src/StdMeshers/StdMeshers_ProxyMesh.hxx | 39 +- .../StdMeshers_QuadToTriaAdaptor.cxx | 550 +++++++++--------- .../StdMeshers_QuadToTriaAdaptor.hxx | 39 +- src/StdMeshers/StdMeshers_ViscousLayers.cxx | 18 +- 5 files changed, 476 insertions(+), 311 deletions(-) diff --git a/src/StdMeshers/StdMeshers_ProxyMesh.cxx b/src/StdMeshers/StdMeshers_ProxyMesh.cxx index 9d37b2364..ace63005a 100644 --- a/src/StdMeshers/StdMeshers_ProxyMesh.cxx +++ b/src/StdMeshers/StdMeshers_ProxyMesh.cxx @@ -27,7 +27,8 @@ #include "SMESH_MesherHelper.hxx" #include -#include +#include +#include //================================================================================ /*! @@ -52,6 +53,9 @@ StdMeshers_ProxyMesh::StdMeshers_ProxyMesh(vector& co for ( unsigned i = 0; i < components.size(); ++i ) { StdMeshers_ProxyMesh* m = components[i].get(); + if ( !m ) continue; + + takeTmpElemsInMesh( m ); if ( !_mesh ) _mesh = m->_mesh; if ( _allowedTypes.empty() ) _allowedTypes = m->_allowedTypes; @@ -69,6 +73,7 @@ StdMeshers_ProxyMesh::StdMeshers_ProxyMesh(vector& co elems.insert( m->_subMeshes[j]->_elements.begin(), m->_subMeshes[j]->_elements.end()); _subMeshes[j]->_elements.assign( elems.begin(), elems.end() ); + m->_subMeshes[j]->_elements.clear(); if ( !_subMeshes[j]->_n2n ) _subMeshes[j]->_n2n = m->_subMeshes[j]->_n2n, m->_subMeshes[j]->_n2n = 0; @@ -88,7 +93,7 @@ StdMeshers_ProxyMesh::StdMeshers_ProxyMesh(vector& co //================================================================================ /*! - * \brief Destructor deletes proxy submeshes + * \brief Destructor deletes proxy submeshes and tmp elemens */ //================================================================================ @@ -97,6 +102,11 @@ StdMeshers_ProxyMesh::~StdMeshers_ProxyMesh() for ( unsigned i = 0; i < _subMeshes.size(); ++i ) delete _subMeshes[i]; _subMeshes.clear(); + + set< const SMDS_MeshElement* >::iterator i = _elemsInMesh.begin(); + for ( ; i != _elemsInMesh.end(); ++i ) + getMeshDS()->RemoveFreeElement( *i, 0 ); + _elemsInMesh.clear(); } //================================================================================ @@ -186,6 +196,7 @@ namespace const SMDS_ElemIteratorPtr& elemIterator) :_iter(elemIterator), _curElem(0), _okTypes(okTypes) { + next(); } virtual bool more() { @@ -198,7 +209,7 @@ namespace while ( _iter->more() && !_curElem ) { _curElem = _iter->next(); - if ( find( _okTypes.begin(), _okTypes.end(), _curElem->GetEntityType()) != _okTypes.end()) + if ( find( _okTypes.begin(), _okTypes.end(), _curElem->GetEntityType()) == _okTypes.end()) _curElem = 0; } return res; @@ -219,10 +230,12 @@ SMDS_ElemIteratorPtr StdMeshers_ProxyMesh::GetFaces(const TopoDS_Shape& shape) c _subContainer.RemoveAllSubmeshes(); - TopExp_Explorer fExp( shape, TopAbs_FACE ); - for ( ; fExp.More(); fExp.Next() ) - if ( const SMESHDS_SubMesh* sm = GetSubMesh( fExp.Current())) + TopTools_IndexedMapOfShape FF; + TopExp::MapShapes( shape, TopAbs_FACE, FF ); + for ( int i = 1; i <= FF.Extent(); ++i ) + if ( const SMESHDS_SubMesh* sm = GetSubMesh( FF(i))) _subContainer.AddSubMesh( sm ); + return _subContainer.SMESHDS_SubMesh::GetElements(); } @@ -247,7 +260,7 @@ SMDS_ElemIteratorPtr StdMeshers_ProxyMesh::GetFaces() const return getMeshDS()->elementsIterator(SMDSAbs_Face); // if _allowedTypes is empty, only elements from _subMeshes are returned,... - SMDS_ElemIteratorPtr proxyIter = _subContainer.GetElements(); + SMDS_ElemIteratorPtr proxyIter = _subContainer.SMESHDS_SubMesh::GetElements(); if ( _allowedTypes.empty() || NbFaces() == _mesh->NbFaces() ) return proxyIter; @@ -271,7 +284,38 @@ SMDS_ElemIteratorPtr StdMeshers_ProxyMesh::GetFaces() const int StdMeshers_ProxyMesh::NbFaces() const { - // TODO + int nb = 0; + if ( _mesh->HasShapeToMesh() ) + { + TopTools_IndexedMapOfShape FF; + TopExp::MapShapes( _mesh->GetShapeToMesh(), TopAbs_FACE, FF ); + for ( int i = 1; i <= FF.Extent(); ++i ) + if ( const SMESHDS_SubMesh* sm = GetSubMesh( FF(i))) + nb += sm->NbElements(); + } + else + { + if ( _subMeshes.empty() ) + return getMeshDS()->NbFaces(); + + for ( unsigned i = 0; i < _subMeshes.size(); ++i ) + if ( _subMeshes[i] ) + nb += _subMeshes[i]->NbElements(); + + // if _allowedTypes is empty, only elements from _subMeshes are returned, + // else elements filtered using allowedTypes are additionally returned + if ( !_allowedTypes.empty() ) + { + for ( int t = SMDSEntity_Triangle; t <= SMDSEntity_Quad_Quadrangle; ++t ) + { + bool allowed = + ( find( _allowedTypes.begin(), _allowedTypes.end(), t ) != _allowedTypes.end() ); + if ( allowed ) + nb += getMeshDS()->GetMeshInfo().NbEntities( SMDSAbs_EntityType( t )); + } + } + } + return nb; } //================================================================================ @@ -317,9 +361,83 @@ StdMeshers_ProxyMesh::SubMesh* StdMeshers_ProxyMesh::findProxySubMesh(int shapeI */ //================================================================================ -const SMESHDS_Mesh* StdMeshers_ProxyMesh::getMeshDS() const +SMESHDS_Mesh* StdMeshers_ProxyMesh::getMeshDS() const { - return _mesh ? _mesh->GetMeshDS() : 0; + return (SMESHDS_Mesh*)( _mesh ? _mesh->GetMeshDS() : 0 ); +} + +//================================================================================ +/*! + * \brief Move proxy sub-mesh from other proxy mesh to this, returns true if sub-mesh found + */ +//================================================================================ + +bool StdMeshers_ProxyMesh::takeProxySubMesh( const TopoDS_Shape& shape, + StdMeshers_ProxyMesh* proxyMesh ) +{ + if ( proxyMesh && proxyMesh->_mesh == _mesh ) + { + int iS = shapeIndex( shape ); + if ( SubMesh* sm = proxyMesh->findProxySubMesh( iS )) + { + if ( iS >= int(_subMeshes.size()) ) + _subMeshes.resize( iS + 1, 0 ); + _subMeshes[iS] = sm; + proxyMesh->_subMeshes[iS] = 0; + return true; + } + } + return false; +} + +//================================================================================ +/*! + * \brief Move tmp elements residing the _mesh from other proxy mesh to this + */ +//================================================================================ + +void StdMeshers_ProxyMesh::takeTmpElemsInMesh( StdMeshers_ProxyMesh* proxyMesh ) +{ + if ( proxyMesh ) + { + _elemsInMesh.insert( proxyMesh->_elemsInMesh.begin(), + proxyMesh->_elemsInMesh.end()); + proxyMesh->_elemsInMesh.clear(); + } +} + +//================================================================================ +/*! + * \brief Removes tmp faces from the _mesh + */ +//================================================================================ + +void StdMeshers_ProxyMesh::removeTmpElement( const SMDS_MeshElement* face ) +{ + if ( face && face->GetID() > 0 ) + { + set< const SMDS_MeshElement* >::iterator i = _elemsInMesh.find( face ); + if ( i != _elemsInMesh.end() ) + { + getMeshDS()->RemoveFreeElement( face, 0 ); + _elemsInMesh.erase( i ); + } + } + else + { + delete face; + } +} + +//================================================================================ +/*! + * \brief Stores tmp element residing the _mesh + */ +//================================================================================ + +void StdMeshers_ProxyMesh::storeTmpElement( const SMDS_MeshElement* face ) +{ + _elemsInMesh.insert( face ); } //================================================================================ @@ -345,7 +463,7 @@ const SMDS_MeshNode* StdMeshers_ProxyMesh::SubMesh::GetProxyNode( const SMDS_Mes void StdMeshers_ProxyMesh::SubMesh::Clear() { for ( unsigned i = 0; i < _elements.size(); ++i ) - if ( _elements[i]->GetID() > 0 ) + if ( _elements[i]->GetID() < 0 ) delete _elements[i]; _elements.clear(); } @@ -383,3 +501,4 @@ void StdMeshers_ProxyMesh::SubMesh::AddElement(const SMDS_MeshElement * e) { _elements.push_back( e ); } + diff --git a/src/StdMeshers/StdMeshers_ProxyMesh.hxx b/src/StdMeshers/StdMeshers_ProxyMesh.hxx index dd24b28ee..094609647 100644 --- a/src/StdMeshers/StdMeshers_ProxyMesh.hxx +++ b/src/StdMeshers/StdMeshers_ProxyMesh.hxx @@ -65,10 +65,18 @@ public: virtual SMDS_ElemIteratorPtr GetElements() const; virtual void Clear(); + template< class ITERATOR > + void ChangeElements( ITERATOR it, ITERATOR end ) + { + // change SubMesh contents without deleting tmp faces + // for which the caller is responsible + _elements.clear(); + while ( it != end ) _elements.push_back( *it++ ); + } SubMesh():_n2n(0) {} ~SubMesh() { Clear(); } - protected: + private: std::vector _elements; TN2NMap* _n2n; friend class StdMeshers_ProxyMesh; @@ -78,6 +86,7 @@ public: StdMeshers_ProxyMesh(); StdMeshers_ProxyMesh(std::vector& components); + StdMeshers_ProxyMesh(const SMESH_Mesh& mesh) { _mesh = &mesh; } virtual ~StdMeshers_ProxyMesh(); // Returns the submesh of a face; it can be a proxy sub-mesh @@ -106,7 +115,9 @@ public: void setMesh(const SMESH_Mesh& mesh) { _mesh = &mesh; } - const SMESHDS_Mesh* getMeshDS() const; + const SMESH_Mesh* getMesh() const { return _mesh; } + + SMESHDS_Mesh* getMeshDS() const; int shapeIndex(const TopoDS_Shape& shape) const; @@ -114,24 +125,38 @@ public: SubMesh* findProxySubMesh(int shapeIndex=0) const; // returns a proxy sub-mesh; it is created if not yet exists - SubMesh* getProxySubMesh(int shapeIndex=0); + SubMesh* getProxySubMesh(int shapeIndex); // returns a proxy sub-mesh; it is created if not yet exists SubMesh* getProxySubMesh(const TopoDS_Shape& shape=TopoDS_Shape()); + // move proxy sub-mesh from other proxy mesh to this, returns true if sub-mesh found + bool takeProxySubMesh( const TopoDS_Shape& shape, StdMeshers_ProxyMesh* proxyMesh ); - private: + // move tmp elements residing the _mesh from other proxy mesh to this + void takeTmpElemsInMesh( StdMeshers_ProxyMesh* proxyMesh ); - const SMESH_Mesh* _mesh; + // removes tmp faces from the _mesh + void removeTmpElement( const SMDS_MeshElement* face ); - // proxy sub-meshes; index in vector == shapeIndex(shape) - std::vector< SubMesh* > _subMeshes; + // stores tmp element residing the _mesh + void storeTmpElement( const SMDS_MeshElement* face ); // types of elements needed to implement NbFaces() and GetFaces(); // if _allowedTypes is empty, only elements from _subMeshes are returned, // else elements of _mesh filtered using allowedTypes are additionally returned std::vector< SMDSAbs_EntityType> _allowedTypes; + private: + + const SMESH_Mesh* _mesh; + + // proxy sub-meshes; index in vector == shapeIndex(shape) + std::vector< SubMesh* > _subMeshes; + + // tmp elements residing the _mesh, to be deleted at destruction + std::set< const SMDS_MeshElement* > _elemsInMesh; + // Complex submesh used to iterate over elements in other sub-meshes mutable SubMesh _subContainer; }; diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx index 8e2718dd0..c87f83b57 100644 --- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx +++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx @@ -44,61 +44,13 @@ using namespace std; -enum EQuadNature { NOT_QUAD, QUAD, DEGEN_QUAD }; +enum EQuadNature { NOT_QUAD, QUAD, DEGEN_QUAD, PYRAM_APEX = 4, TRIA_APEX = 0 }; // std-like iterator used to get coordinates of nodes of mesh element typedef SMDS_StdIterator< SMESH_MeshEditor::TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator; namespace { - const int Q2TAbs_TmpTriangle = SMDSEntity_Last + 1; - - //================================================================================ - /*! - * \brief Temporary face. It's key feature is that it adds itself to an apex node - * as inverse element, so that tmp triangles of a piramid can be easily found - */ - //================================================================================ - - class STDMESHERS_EXPORT Q2TAdaptor_Triangle : public SMDS_MeshFace - { - const SMDS_MeshNode* _nodes[3]; - public: - Q2TAdaptor_Triangle(const SMDS_MeshNode* apexNode, - const SMDS_MeshNode* node2, - const SMDS_MeshNode* node3) - { - _nodes[0]=0; ChangeApex(apexNode); - _nodes[1]=node2; - _nodes[2]=node3; - } - ~Q2TAdaptor_Triangle() { MarkAsRemoved(); } - void ChangeApex(const SMDS_MeshNode* node) - { - MarkAsRemoved(); - _nodes[0]=node; - const_cast(node)->AddInverseElement(this); - } - void MarkAsRemoved() - { - if ( _nodes[0] ) - const_cast(_nodes[0])->RemoveInverseElement(this), _nodes[0] = 0; - } - bool IsRemoved() const { return !_nodes[0]; } - virtual int NbNodes() const { return 3; } - virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nodes[ind]; } - virtual SMDSAbs_EntityType GetEntityType() const - { - return SMDSAbs_EntityType( Q2TAbs_TmpTriangle ); - } - virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const - { - if ( type == SMDSAbs_Node ) - return SMDS_ElemIteratorPtr( new SMDS_NodeArrayElemIterator( _nodes, _nodes+3 )); - throw SALOME_Exception(LOCALIZED("Not implemented")); - } - }; - //================================================================================ /*! * \brief Return true if two nodes of triangles are equal @@ -111,75 +63,6 @@ namespace ( F1->GetNode(1)==F2->GetNode(2) && F1->GetNode(2)==F2->GetNode(1) ) || ( F1->GetNode(1)==F2->GetNode(1) && F1->GetNode(2)==F2->GetNode(2) ); } - - //================================================================================ - /*! - * \brief Merge the two pyramids (i.e. fuse their apex) and others already merged with them - */ - //================================================================================ - - void MergePiramids( const SMDS_MeshElement* PrmI, - const SMDS_MeshElement* PrmJ, - SMESHDS_Mesh* meshDS, - set & nodesToMove) - { - const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove - int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume ); - SMESH_MeshEditor::TNodeXYZ Pj( Nrem ); - - // an apex node to make common to all merged pyramids - SMDS_MeshNode* CommonNode = const_cast(PrmI->GetNode(4)); - if ( CommonNode == Nrem ) return; // already merged - int nbI = CommonNode->NbInverseElements( SMDSAbs_Volume ); - SMESH_MeshEditor::TNodeXYZ Pi( CommonNode ); - gp_XYZ Pnew = ( nbI*Pi + nbJ*Pj ) / (nbI+nbJ); - CommonNode->setXYZ( Pnew.X(), Pnew.Y(), Pnew.Z() ); - - nodesToMove.insert( CommonNode ); - nodesToMove.erase ( Nrem ); - - // find and remove coincided faces of merged pyramids - SMDS_ElemIteratorPtr triItI = CommonNode->GetInverseElementIterator(SMDSAbs_Face); - while ( triItI->more() ) - { - const SMDS_MeshElement* FI = triItI->next(); - const SMDS_MeshElement* FJEqual = 0; - SMDS_ElemIteratorPtr triItJ = Nrem->GetInverseElementIterator(SMDSAbs_Face); - while ( !FJEqual && triItJ->more() ) - { - const SMDS_MeshElement* FJ = triItJ->next(); - if ( EqualTriangles( FJ, FI )) - FJEqual = FJ; - } - if ( FJEqual ) - { - ((Q2TAdaptor_Triangle*) FI)->MarkAsRemoved(); - ((Q2TAdaptor_Triangle*) FJEqual)->MarkAsRemoved(); - } - } - - // set the common apex node to pyramids and triangles merged with J - SMDS_ElemIteratorPtr itJ = Nrem->GetInverseElementIterator(); - while ( itJ->more() ) - { - const SMDS_MeshElement* elem = itJ->next(); - if ( elem->GetType() == SMDSAbs_Volume ) // pyramid - { - vector< const SMDS_MeshNode* > nodes( elem->begin_nodes(), elem->end_nodes() ); - nodes[4] = CommonNode; - meshDS->ChangeElementNodes( elem, &nodes[0], nodes.size()); - } - else if ( elem->GetEntityType() == Q2TAbs_TmpTriangle ) // tmp triangle - { - ((Q2TAdaptor_Triangle*) elem )->ChangeApex( CommonNode ); - } - } - ASSERT( Nrem->NbInverseElements() == 0 ); - meshDS->RemoveFreeNode( Nrem, - meshDS->MeshElements( Nrem->GetPosition()->GetShapeId()), - /*fromGroups=*/false); - } - //================================================================================ /*! * \brief Return true if two adjacent pyramids are too close one to another @@ -286,44 +169,110 @@ namespace return tooClose; } +} - //================================================================================ - /*! - * \brief Merges adjacent pyramids - */ - //================================================================================ +//================================================================================ +/*! + * \brief Merge the two pyramids (i.e. fuse their apex) and others already merged with them + */ +//================================================================================ - void MergeAdjacent(const SMDS_MeshElement* PrmI, - SMESH_Mesh& mesh, - set& nodesToMove) +void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* PrmI, + const SMDS_MeshElement* PrmJ, + set & nodesToMove) +{ + const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove + int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume ); + SMESH_MeshEditor::TNodeXYZ Pj( Nrem ); + + // an apex node to make common to all merged pyramids + SMDS_MeshNode* CommonNode = const_cast(PrmI->GetNode(4)); + if ( CommonNode == Nrem ) return; // already merged + int nbI = CommonNode->NbInverseElements( SMDSAbs_Volume ); + SMESH_MeshEditor::TNodeXYZ Pi( CommonNode ); + gp_XYZ Pnew = ( nbI*Pi + nbJ*Pj ) / (nbI+nbJ); + CommonNode->setXYZ( Pnew.X(), Pnew.Y(), Pnew.Z() ); + + nodesToMove.insert( CommonNode ); + nodesToMove.erase ( Nrem ); + + typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator; + TStdElemIterator itEnd; + + // find and remove coincided faces of merged pyramids + vector< const SMDS_MeshElement* > inverseElems + // copy inverse elements to avoid iteration on changing conainer + ( TStdElemIterator( CommonNode->GetInverseElementIterator(SMDSAbs_Face)), itEnd); + for ( unsigned i = 0; i < inverseElems.size(); ++i ) { - TIDSortedElemSet adjacentPyrams, mergedPyrams; - for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI + const SMDS_MeshElement* FI = inverseElems[i]; + const SMDS_MeshElement* FJEqual = 0; + SMDS_ElemIteratorPtr triItJ = Nrem->GetInverseElementIterator(SMDSAbs_Face); + while ( !FJEqual && triItJ->more() ) { - const SMDS_MeshNode* n = PrmI->GetNode(k); - SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); - while ( vIt->more() ) - { - const SMDS_MeshElement* PrmJ = vIt->next(); - if ( PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second ) - continue; - if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, mesh.HasShapeToMesh() )) - { - MergePiramids( PrmI, PrmJ, mesh.GetMeshDS(), nodesToMove ); - mergedPyrams.insert( PrmJ ); - } - } + const SMDS_MeshElement* FJ = triItJ->next(); + if ( EqualTriangles( FJ, FI )) + FJEqual = FJ; } - if ( !mergedPyrams.empty() ) + if ( FJEqual ) { - TIDSortedElemSet::iterator prm; -// for (prm = mergedPyrams.begin(); prm != mergedPyrams.end(); ++prm) -// MergeAdjacent( *prm, mesh, nodesToMove ); + removeTmpElement( FI ); + removeTmpElement( FJEqual ); + myRemovedTrias.insert( FI ); + myRemovedTrias.insert( FJEqual ); + } + } - for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm) - MergeAdjacent( *prm, mesh, nodesToMove ); + // set the common apex node to pyramids and triangles merged with J + inverseElems.assign( TStdElemIterator( Nrem->GetInverseElementIterator()), itEnd ); + for ( unsigned i = 0; i < inverseElems.size(); ++i ) + { + const SMDS_MeshElement* elem = inverseElems[i]; + vector< const SMDS_MeshNode* > nodes( elem->begin_nodes(), elem->end_nodes() ); + nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode; + getMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size()); + } + ASSERT( Nrem->NbInverseElements() == 0 ); + getMeshDS()->RemoveFreeNode( Nrem, + getMeshDS()->MeshElements( Nrem->GetPosition()->GetShapeId()), + /*fromGroups=*/false); +} + +//================================================================================ +/*! + * \brief Merges adjacent pyramids + */ +//================================================================================ + +void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement* PrmI, + set& nodesToMove) +{ + TIDSortedElemSet adjacentPyrams; + bool mergedPyrams = false; + for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI + { + const SMDS_MeshNode* n = PrmI->GetNode(k); + SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); + while ( vIt->more() ) + { + const SMDS_MeshElement* PrmJ = vIt->next(); + if ( PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second ) + continue; + if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, getMesh()->HasShapeToMesh() )) + { + MergePiramids( PrmI, PrmJ, nodesToMove ); + mergedPyrams = true; + // container of inverse elements can change + vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); + } } } + if ( mergedPyrams ) + { + TIDSortedElemSet::iterator prm; + for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm) + MergeAdjacent( *prm, nodesToMove ); + } } //================================================================================ @@ -333,7 +282,7 @@ namespace //================================================================================ StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor(): - myElemSearcher(0), myNbTriangles(0) + myElemSearcher(0)//, myNbTriangles(0) { } @@ -345,17 +294,7 @@ StdMeshers_QuadToTriaAdaptor::StdMeshers_QuadToTriaAdaptor(): StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor() { - // delete temporary faces - TQuad2Trias::iterator f_f = myResMap.begin(), ffEnd = myResMap.end(); - for ( ; f_f != ffEnd; ++f_f ) - { - TTriaList& fList = f_f->second; - TTriaList::iterator f = fList.begin(), fEnd = fList.end(); - for ( ; f != fEnd; ++f ) - delete *f; - } - myResMap.clear(); - + // temporary faces are deleted by ~StdMeshers_ProxyMesh() if ( myElemSearcher ) delete myElemSearcher; myElemSearcher=0; } @@ -582,7 +521,7 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement* face { if( face->NbCornerNodes() != 4 ) { - myNbTriangles += int( face->NbCornerNodes() == 3 ); + //myNbTriangles += int( face->NbCornerNodes() == 3 ); return NOT_QUAD; } @@ -694,22 +633,38 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement* face //purpose : //======================================================================= -bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape) +bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + StdMeshers_ProxyMesh* aProxyMesh) { - myResMap.clear(); - myPyramids.clear(); - myNbTriangles = 0; - myShape = aShape; + StdMeshers_ProxyMesh::setMesh( aMesh ); + //myResMap.clear(); + vector myPyramids; + //myNbTriangles = 0; + //myShape = aShape; SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); SMESH_MesherHelper helper(aMesh); helper.IsQuadraticSubMesh(aShape); helper.SetElementsOnShape( true ); + if ( myElemSearcher ) delete myElemSearcher; + if ( aProxyMesh ) + myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher( aProxyMesh->GetFaces(aShape)); + else + myElemSearcher = SMESH_MeshEditor(&aMesh).GetElementSearcher(); + + const SMESHDS_SubMesh * aSubMeshDSFace; + for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) { const TopoDS_Shape& aShapeFace = exp.Current(); - const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace ); + if ( aProxyMesh ) + aSubMeshDSFace = aProxyMesh->GetSubMesh( aShapeFace ); + else + aSubMeshDSFace = meshDS->MeshElements( aShapeFace ); + vector trias, quads; + bool hasNewTrias = false; if ( aSubMeshDSFace ) { bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS ); @@ -725,88 +680,123 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape vector FNodes(5); gp_Pnt PC; gp_Vec VNorm; - int stat = Preparation(face, PN, VN, FNodes, PC, VNorm); - if(stat==0) - continue; - - if(stat==2) + int stat = Preparation(face, PN, VN, FNodes, PC, VNorm); + switch ( stat ) { - // degenerate face - // add triangles to result map - SMDS_MeshFace* NewFace; - if(!isRev) - NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[1], FNodes[2] ); - else - NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[2], FNodes[1] ); - TTriaList aList( 1, NewFace ); - myResMap.insert(make_pair(face,aList)); - continue; - } + case NOT_QUAD: - if(!isRev) VNorm.Reverse(); - double xc = 0., yc = 0., zc = 0.; - int i = 1; - for(; i<=4; i++) { - gp_Pnt Pbest; - if(!isRev) - Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i).Reversed()); - else - Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i)); - xc += Pbest.X(); - yc += Pbest.Y(); - zc += Pbest.Z(); - } - gp_Pnt PCbest(xc/4., yc/4., zc/4.); + trias.push_back( face ); + break; - // check PCbest - double height = PCbest.Distance(PC); - if(height<1.e-6) { - // create new PCbest using a bit shift along VNorm - PCbest = PC.XYZ() + VNorm.XYZ() * 0.001; - } - else { - // check possible intersection with other faces - gp_Pnt Pint; - bool check = CheckIntersection(PCbest, PC, Pint, aMesh, aShape, face); - if(check) { - //cout<<"--PC("<AddFace( FNodes[0], FNodes[1], FNodes[2] ); + else + NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] ); + storeTmpElement( NewFace ); + trias.push_back ( NewFace ); + quads.push_back( face ); + hasNewTrias = true; + break; } - else { - gp_Vec VB(PC,PCbest); - gp_Pnt PCbestTmp = PC.XYZ() + VB.XYZ() * 3.0; - check = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, face); - if(check) { - double dist = PC.Distance(Pint)/3.; - if(distValue(i), PN->Value(i+1), PC, VN->Value(i).Reversed()); + else + Pbest = FindBestPoint(PN->Value(i), PN->Value(i+1), PC, VN->Value(i)); + xc += Pbest.X(); + yc += Pbest.Y(); + zc += Pbest.Z(); + } + gp_Pnt PCbest(xc/4., yc/4., zc/4.); + + // check PCbest + double height = PCbest.Distance(PC); + if(height<1.e-6) { + // create new PCbest using a bit shift along VNorm + PCbest = PC.XYZ() + VNorm.XYZ() * 0.001; + } + else { + // check possible intersection with other faces + gp_Pnt Pint; + bool check = CheckIntersection(PCbest, PC, Pint, aMesh, aShape, face); + if(check) { + //cout<<"--PC("<AddFace( NewNode, FNodes[i], FNodes[i+1] )); + storeTmpElement( trias.back() ); + } + // create a pyramid + if ( isRev ) swap( FNodes[1], FNodes[3]); + SMDS_MeshVolume* aPyram = + helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); + myPyramids.push_back(aPyram); + + quads.push_back( face ); + hasNewTrias = true; + break; + + } // case QUAD: - // add triangles to result map - TTriaList& triaList = myResMap.insert( make_pair( face, TTriaList() ))->second; - for(i=0; i<4; i++) - triaList.push_back( new Q2TAdaptor_Triangle( NewNode, FNodes[i], FNodes[i+1] )); + } // switch ( stat ) - // create a pyramid - if ( isRev ) swap( FNodes[1], FNodes[3]); - SMDS_MeshVolume* aPyram = - helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); - myPyramids.push_back(aPyram); + bool sourceSubMeshIsProxy = false; + if ( aProxyMesh ) + { + // move proxy sub-mesh from other proxy mesh to this + sourceSubMeshIsProxy = takeProxySubMesh( aShapeFace, aProxyMesh ); + // move tmp elements added in mesh + takeTmpElemsInMesh( aProxyMesh ); + } + if ( hasNewTrias ) + { + StdMeshers_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh( aShapeFace ); + prxSubMesh->ChangeElements( trias.begin(), trias.end() ); + // delete tmp quadrangles removed from aProxyMesh + if ( sourceSubMeshIsProxy ) + for ( unsigned i = 0; i < quads.size(); ++i ) + removeTmpElement( quads[i] ); + } } // end loop on elements on a face submesh } } // end for(TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) { - return Compute2ndPart(aMesh); + return Compute2ndPart(aMesh, myPyramids); } //================================================================================ @@ -817,8 +807,13 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) { - myResMap.clear(); - myPyramids.clear(); + StdMeshers_ProxyMesh::setMesh( aMesh ); + StdMeshers_ProxyMesh::_allowedTypes.push_back( SMDSEntity_Triangle ); + StdMeshers_ProxyMesh::_allowedTypes.push_back( SMDSEntity_Quad_Triangle ); + if ( aMesh.NbQuadrangles() < 1 ) + return false; + + vector myPyramids; SMESH_MesherHelper helper(aMesh); helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh()); helper.SetElementsOnShape( true ); @@ -828,13 +823,13 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) SMESH_ElementSearcher* searcher = const_cast(myElemSearcher); SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); + StdMeshers_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh(); SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(/*idInceasingOrder=*/true); while( fIt->more()) { const SMDS_MeshElement* face = fIt->next(); if ( !face ) continue; - //cout<GetID() = "<GetID()<Value(1).Distance(PN->Value(2)) + PN->Value(2).Distance(PN->Value(3)); // far points in VNorm direction gp_Pnt Ptmp1 = PC.XYZ() + VNorm.XYZ() * tmp * 1.e6; @@ -916,14 +910,16 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) } } if(!IsRev) - NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[1], FNodes[2] ); + NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] ); else - NewFace = new Q2TAdaptor_Triangle( FNodes[0], FNodes[2], FNodes[1] ); - aList.push_back(NewFace); - myResMap.insert(make_pair(face,aList)); + NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] ); + storeTmpElement( NewFace ); + prxSubMesh->AddElement( NewFace ); continue; } + // Case of non-degenerated quadrangle + // Find pyramid peak gp_XYZ PCbest(0., 0., 0.); // pyramid peak @@ -940,7 +936,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) PCbest = PC.XYZ() + VNorm.XYZ() * 0.001; height = PC.Distance(PCbest); } - //cout<<" PCbest("<second; for(i=0; i<4; i++) { SMDS_MeshFace* NewFace; if(isRev) - NewFace = new Q2TAdaptor_Triangle( NewNode, FNodes[i], FNodes[i+1] ); + NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ); else - NewFace = new Q2TAdaptor_Triangle( NewNode, FNodes[i+1], FNodes[i] ); - aList.push_back(NewFace); + NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] ); + storeTmpElement( NewFace ); + prxSubMesh->AddElement( NewFace ); } // create a pyramid SMDS_MeshVolume* aPyram; @@ -1011,7 +1006,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) } } // end loop on all faces - return Compute2ndPart(aMesh); + return Compute2ndPart(aMesh, myPyramids); } //================================================================================ @@ -1020,7 +1015,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) */ //================================================================================ -bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) +bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh, + const vector& myPyramids) { if(myPyramids.empty()) return true; @@ -1039,7 +1035,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) for ( i = 0; i < myPyramids.size(); ++i ) { const SMDS_MeshElement* PrmI = myPyramids[i]; - MergeAdjacent( PrmI, aMesh, nodesToMove ); + MergeAdjacent( PrmI, nodesToMove ); } // iterate on all pyramids @@ -1116,7 +1112,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) if(nbc>0) { // Merge the two pyramids and others already merged with them - MergePiramids( PrmI, PrmJ, meshDS, nodesToMove ); + MergePiramids( PrmI, PrmJ, nodesToMove ); } else { // nbc==0 @@ -1149,8 +1145,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) nodesToMove.insert( aNode2 ); } // fix intersections that could appear after apex movement - MergeAdjacent( PrmI, aMesh, nodesToMove ); - MergeAdjacent( PrmJ, aMesh, nodesToMove ); + MergeAdjacent( PrmI, nodesToMove ); + MergeAdjacent( PrmJ, nodesToMove ); } // end if(hasInt) } // loop on suspectPyrams @@ -1165,26 +1161,28 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) meshDS->MoveNode( *n, (*n)->X(), (*n)->Y(), (*n)->Z() ); } - // rebind triangles of pyramids sharing the same base quadrangle to the first - // entrance of the base quadrangle - TQuad2Trias::iterator q2t = myResMap.begin(), q2tPrev = q2t; - for ( ++q2t; q2t != myResMap.end(); ++q2t, ++q2tPrev ) - { - if ( q2t->first == q2tPrev->first ) - q2tPrev->second.splice( q2tPrev->second.end(), q2t->second ); - } - // delete removed triangles and count resulting nb of triangles - for ( q2t = myResMap.begin(); q2t != myResMap.end(); ++q2t ) + // erase removed triangles from the proxy mesh + if ( !myRemovedTrias.empty() ) { - TTriaList & trias = q2t->second; - for ( TTriaList::iterator tri = trias.begin(); tri != trias.end(); ) - if ( ((const Q2TAdaptor_Triangle*) *tri)->IsRemoved() ) - delete *tri, trias.erase( tri++ ); - else - tri++, myNbTriangles++; + for ( int i = 0; i <= meshDS->MaxShapeIndex(); ++i ) + if ( StdMeshers_ProxyMesh::SubMesh* sm = findProxySubMesh(i)) + { + vector faces; + faces.reserve( sm->NbElements() ); + SMDS_ElemIteratorPtr fIt = sm->GetElements(); + while ( fIt->more() ) + { + const SMDS_MeshElement* tria = fIt->next(); + set::iterator rmTria = myRemovedTrias.find( tria ); + if ( rmTria != myRemovedTrias.end() ) + myRemovedTrias.erase( rmTria ); + else + faces.push_back( tria ); + } + sm->ChangeElements( faces.begin(), faces.end() ); + } } - myPyramids.clear(); // no more needed myDegNodes.clear(); delete myElemSearcher; @@ -1199,8 +1197,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh) */ //================================================================================ -const list* StdMeshers_QuadToTriaAdaptor::GetTriangles (const SMDS_MeshElement* aQuad) -{ - TQuad2Trias::iterator it = myResMap.find(aQuad); - return ( it != myResMap.end() ? & it->second : 0 ); -} +// const list* StdMeshers_QuadToTriaAdaptor::GetTriangles (const SMDS_MeshElement* aQuad) +// { +// TQuad2Trias::iterator it = myResMap.find(aQuad); +// return ( it != myResMap.end() ? & it->second : 0 ); +// } diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx index ceccb2542..e3c8e51c0 100644 --- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx +++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx @@ -26,6 +26,8 @@ #include "SMESH_StdMeshers.hxx" +#include "StdMeshers_ProxyMesh.hxx" + class SMESH_Mesh; class SMESH_ElementSearcher; class SMDS_MeshElement; @@ -37,7 +39,7 @@ class gp_Pnt; class gp_Vec; -#include +#include #include #include @@ -46,25 +48,27 @@ class gp_Vec; /*! * \brief "Transforms" quadrilateral faces into triangular ones by creation of pyramids */ -class STDMESHERS_EXPORT StdMeshers_QuadToTriaAdaptor +class STDMESHERS_EXPORT StdMeshers_QuadToTriaAdaptor : public StdMeshers_ProxyMesh { public: StdMeshers_QuadToTriaAdaptor(); ~StdMeshers_QuadToTriaAdaptor(); - bool Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape); + bool Compute(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + StdMeshers_ProxyMesh* aProxyMesh=0); bool Compute(SMESH_Mesh& aMesh); - const std::list* GetTriangles(const SMDS_MeshElement* aFace); + //const std::list* GetTriangles(const SMDS_MeshElement* aFace); /*! * \brief Return sum of generated and already present triangles */ - int TotalNbOfTriangles() const { return myNbTriangles; } + //int TotalNbOfTriangles() const { return myNbTriangles; } - TopoDS_Shape GetShape() const { return myShape; } + //TopoDS_Shape GetShape() const { return myShape; } protected: @@ -82,22 +86,31 @@ protected: const TopoDS_Shape& aShape, const SMDS_MeshElement* NotCheckedFace); - bool Compute2ndPart(SMESH_Mesh& aMesh); + bool Compute2ndPart(SMESH_Mesh& aMesh, + const std::vector& pyramids); + + + void MergePiramids( const SMDS_MeshElement* PrmI, + const SMDS_MeshElement* PrmJ, + set & nodesToMove); + void MergeAdjacent(const SMDS_MeshElement* PrmI, + set& nodesToMove); + //typedef std::list TTriaList; + //typedef std::multimap TQuad2Trias; - typedef std::list TTriaList; - typedef std::multimap TQuad2Trias; + //TQuad2Trias myResMap; + //std::vector myPyramids; - TQuad2Trias myResMap; - std::vector myPyramids; + std::set myRemovedTrias; std::list< const SMDS_MeshNode* > myDegNodes; const SMESH_ElementSearcher* myElemSearcher; - int myNbTriangles; + //int myNbTriangles; - TopoDS_Shape myShape; + //TopoDS_Shape myShape; }; #endif diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 5eecf1a81..0a0da267a 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -1174,8 +1174,8 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) if ( data._edges[i]->IsOnEdge()) for ( int j = 0; j < 2; ++j ) { - //if ( !data._edges[i]->_2neibors ) - //break; // _LayerEdge is shared by two _SolidData's + if ( data._edges[i]->_nodes.back()->NbInverseElements(SMDSAbs_Volume) > 0 ) + break; // _LayerEdge is shared by two _SolidData's const SMDS_MeshNode* & n = data._edges[i]->_2neibors->_nodes[j]; if (( n2e = data._n2eMap.find( n )) == data._n2eMap.end() ) return error("_LayerEdge not found by src node", &data); @@ -2927,8 +2927,8 @@ bool _ViscousBuilder::refine(_SolidData& data) bool _ViscousBuilder::shrink() { - // make map of (ids of faces to shrink mesh on) to (_SolidData containing _LayerEdge's - // inflated along face or edge) + // make map of (ids of FACEs to shrink mesh on) to (_SolidData containing _LayerEdge's + // inflated along FACE or EDGE) map< TGeomID, _SolidData* > f2sdMap; for ( unsigned i = 0 ; i < _sdVec.size(); ++i ) { @@ -2960,6 +2960,16 @@ bool _ViscousBuilder::shrink() helper.SetSubShape(F); + // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid + // usage of mesh faces made in addBoundaryElements() by the 3D algo or + // by StdMeshers_QuadToTriaAdaptor + { + StdMeshers_ProxyMesh::SubMesh* proxySub = + data._proxyMesh->getFaceSubM( F, /*create=*/true); + SMDS_ElemIteratorPtr fIt = smDS->GetElements(); + while ( fIt->more() ) + proxySub->AddElement( fIt->next() ); + } // =========================== // Prepare data for shrinking // =========================== -- 2.39.2