From 47d34023e6ec1b7501a90a081fa97b721688294b Mon Sep 17 00:00:00 2001 From: eap Date: Wed, 12 Apr 2017 16:52:32 +0300 Subject: [PATCH] 23427: [CEA 2073] No hypothesis "Viscous Layers" with Netgen 1D-2D-3D --- resources/NETGENPlugin.xml | 2 + src/NETGENPlugin/NETGENPlugin_Mesher.cxx | 201 +++++++++++++++--- src/NETGENPlugin/NETGENPlugin_Mesher.hxx | 14 +- src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.cxx | 60 +++--- src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.hxx | 16 +- 5 files changed, 221 insertions(+), 72 deletions(-) diff --git a/resources/NETGENPlugin.xml b/resources/NETGENPlugin.xml index 9e0c77e..d7d7651 100644 --- a/resources/NETGENPlugin.xml +++ b/resources/NETGENPlugin.xml @@ -152,6 +152,7 @@ group-id ="1" priority ="10" hypos ="NETGEN_Parameters, NETGEN_SimpleParameters_3D" + opt-hypos ="ViscousLayers" output ="TETRA,PYRAMID" dim ="3" support-submeshes="true"> @@ -159,6 +160,7 @@ NETGEN_2D3D=Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D) NETGEN_Parameters=Parameters() NETGEN_SimpleParameters_3D=Parameters(smeshBuilder.SIMPLE) + ViscousLayers=ViscousLayers(SetTotalThickness(),SetNumberLayers(),SetStretchFactor(),SetFaces(1),SetFaces(2),SetMethod()) diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx index be4aa07..dc3afe7 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx @@ -157,6 +157,7 @@ NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh, _progressTic(1), _totalTime(1.0), _simpleHyp(NULL), + _viscousLayersHyp(NULL), _ptrToMe(NULL) { SetDefaultParameters(); @@ -336,6 +337,17 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D* SetDefaultParameters(); } +//================================================================================ +/*! + * \brief Store a Viscous Layers hypothesis + */ +//================================================================================ + +void NETGENPlugin_Mesher::SetParameters(const StdMeshers_ViscousLayers* hyp ) +{ + _viscousLayersHyp = hyp; +} + //============================================================================= /*! * Link - a pair of integer numbers @@ -1054,13 +1066,6 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, // Find solids the geomFace bounds int solidID1 = 0, solidID2 = 0; - StdMeshers_QuadToTriaAdaptor* quadAdaptor = - dynamic_cast( proxyMesh.get() ); - if ( quadAdaptor ) - { - solidID1 = occgeom.somap.FindIndex( quadAdaptor->GetShape() ); - } - else { PShapeIteratorPtr solidIt = helper.GetAncestors( geomFace, *sm->GetFather(), TopAbs_SOLID); while ( const TopoDS_Shape * solid = solidIt->next() ) @@ -1070,6 +1075,81 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, else solidID1 = id; } } + if ( proxyMesh && proxyMesh->GetProxySubMesh( geomFace )) + { + // if a proxy sub-mesh contains temporary faces, then these faces + // should be used to mesh only one SOLID + bool hasTmp = false; + smDS = proxyMesh->GetSubMesh( geomFace ); + SMDS_ElemIteratorPtr faces = smDS->GetElements(); + while ( faces->more() ) + { + const SMDS_MeshElement* f = faces->next(); + if ( proxyMesh->IsTemporary( f )) + { + hasTmp = true; + std::vector fNodes( f->begin_nodes(), f->end_nodes() ); + std::vector vols; + if ( _mesh->GetMeshDS()->GetElementsByNodes( fNodes, vols, SMDSAbs_Volume ) == 1 ) + { + int geomID = vols[0]->getshapeId(); + const TopoDS_Shape& solid = helper.GetMeshDS()->IndexToShape( geomID ); + if ( !solid.IsNull() ) + solidID1 = occgeom.somap.FindIndex ( solid ); + solidID2 = 0; + break; + } + } + } + // exclude faces generated by NETGEN from computation of 3D mesh + const int fID = occgeom.fmap.FindIndex( geomFace ); + if ( !hasTmp ) // shrunk mesh + { + // move netgen points according to moved nodes + SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true); + while ( smIt->more() ) + { + SMESH_subMesh* sub = smIt->next(); + if ( !sub->GetSubMeshDS() ) continue; + SMDS_NodeIteratorPtr nodeIt = sub->GetSubMeshDS()->GetNodes(); + while ( nodeIt->more() ) + { + const SMDS_MeshNode* n = nodeIt->next(); + int ngID = ngNodeId( n, ngMesh, nodeNgIdMap ); + netgen::MeshPoint& ngPoint = ngMesh.Point( ngID ); + ngPoint(0) = n->X(); + ngPoint(1) = n->Y(); + ngPoint(2) = n->Z(); + } + } + // remove faces near boundary to avoid their overlapping + // with shrunk faces + for ( int i = 1; i <= ngMesh.GetNSE(); ++i ) + { + const netgen::Element2d& elem = ngMesh.SurfaceElement(i); + if ( elem.GetIndex() == fID ) + { + for ( int iN = 0; iN < elem.GetNP(); ++iN ) + if ( ngMesh[ elem[ iN ]].Type() != netgen::SURFACEPOINT ) + { + ngMesh.DeleteSurfaceElement( i ); + break; + } + } + } + } + //if ( hasTmp ) + { + faceNgID++; + ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceNgID,/*solid1=*/0,/*solid2=*/0,0 )); + for (int i = 1; i <= ngMesh.GetNSE(); ++i ) + { + const netgen::Element2d& elem = ngMesh.SurfaceElement(i); + if ( elem.GetIndex() == fID ) + const_cast< netgen::Element2d& >( elem ).SetIndex( faceNgID ); + } + } + } // Add ng face descriptors of meshed faces faceNgID++; ngMesh.AddFaceDescriptor( netgen::FaceDescriptor( faceNgID, solidID1, solidID2, 0 )); @@ -1112,8 +1192,6 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, cout << "SMESH face " << helper.GetMeshDS()->ShapeToIndex( geomFace ) << " internal="<GetSubMesh( geomFace ); SMDS_ElemIteratorPtr faces = smDS->GetElements(); while ( faces->more() ) @@ -1430,10 +1508,15 @@ namespace int ngIdCloseN; //!< ng id of closest node of the closest 2d mesh element }; - inline double dist2(const netgen::MeshPoint& p1, const netgen::MeshPoint& p2) + inline double dist2( const netgen::MeshPoint& p1, const netgen::MeshPoint& p2 ) { return gp_Pnt( NGPOINT_COORDS(p1)).SquareDistance( gp_Pnt( NGPOINT_COORDS(p2))); } + + // inline double dist2(const netgen::MeshPoint& p, const SMDS_MeshNode* n ) + // { + // return gp_Pnt( NGPOINT_COORDS(p)).SquareDistance( SMESH_NodeXYZ(n)); + // } } //================================================================================ @@ -2136,13 +2219,35 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo, if ( quadHelper && !quadHelper->GetIsQuadratic() && quadHelper->GetTLinkNodeMap().empty() ) quadHelper = 0; + int i, nbInitNod = initState._nbNodes; + if ( initState._elementsRemoved ) + { + // PAL23427. Update nodeVec to track removal of netgen free points as a result + // of removal of faces in FillNgMesh() in the case of a shrunk sub-mesh + int ngID, nodeVecSize = nodeVec.size(); + const double eps = std::numeric_limits::min(); + for ( ngID = i = 1; i < nodeVecSize; ++ngID, ++i ) + { + gp_Pnt ngPnt( NGPOINT_COORDS( ngMesh.Point( ngID ))); + gp_Pnt node ( SMESH_NodeXYZ ( nodeVec[ i ])); + if ( ngPnt.SquareDistance( node ) < eps ) + { + nodeVec[ ngID ] = nodeVec[ i ]; + } + else + { + --ngID; + } + } + nodeVec.resize( ngID ); + nbInitNod = ngID - 1; + } // ------------------------------------- // Create and insert nodes into nodeVec // ------------------------------------- nodeVec.resize( nbNod + 1 ); - int i, nbInitNod = initState._nbNodes; - for (i = nbInitNod+1; i <= nbNod; ++i ) + for ( i = nbInitNod+1; i <= nbNod; ++i ) { const netgen::MeshPoint& ngPoint = ngMesh.Point(i); SMDS_MeshNode* node = NULL; @@ -2833,37 +2938,58 @@ bool NETGENPlugin_Mesher::Compute() // generate volume mesh // --------------------- // Fill _ngMesh with nodes and faces of computed 2D submeshes - if ( !err && _isVolume && ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad )) + if ( !err && _isVolume && + ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad || _viscousLayersHyp )) { // load SMESH with computed segments and faces FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper ); + // compute prismatic boundary volumes + int nbQuad = _mesh->NbQuadrangles(); + SMESH_ProxyMesh::Ptr viscousMesh; + if ( _viscousLayersHyp ) + { + viscousMesh = _viscousLayersHyp->Compute( *_mesh, _shape ); + if ( !viscousMesh ) + return false; + } // compute pyramids on quadrangles - SMESH_ProxyMesh::Ptr proxyMesh; - if ( _mesh->NbQuadrangles() > 0 ) + vector pyramidMeshes( occgeo.somap.Extent() ); + if ( nbQuad > 0 ) for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS ) { - StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor; - proxyMesh.reset( Adaptor ); - - int nbPyrams = _mesh->NbPyramids(); - Adaptor->Compute( *_mesh, occgeo.somap(iS) ); - if ( nbPyrams != _mesh->NbPyramids() ) - { - list< SMESH_subMesh* > quadFaceSM; - for (TopExp_Explorer face(occgeo.somap(iS), TopAbs_FACE); face.More(); face.Next()) - if ( Adaptor->GetProxySubMesh( face.Current() )) - { - quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() )); - meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() ); - } - FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, proxyMesh); - } + StdMeshers_QuadToTriaAdaptor* adaptor = new StdMeshers_QuadToTriaAdaptor; + pyramidMeshes[ iS-1 ].reset( adaptor ); + bool ok = adaptor->Compute( *_mesh, occgeo.somap(iS), viscousMesh.get() ); + if ( !ok ) + return false; } + // add proxy faces to NG mesh + list< SMESH_subMesh* > viscousSM; + for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS ) + { + list< SMESH_subMesh* > quadFaceSM; + for (TopExp_Explorer face(occgeo.somap(iS), TopAbs_FACE); face.More(); face.Next()) + if ( pyramidMeshes[iS-1] && pyramidMeshes[iS-1]->GetProxySubMesh( face.Current() )) + { + quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() )); + meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() ); + } + else if ( viscousMesh && viscousMesh->GetProxySubMesh( face.Current() )) + { + viscousSM.push_back( _mesh->GetSubMesh( face.Current() )); + meshedSM[ MeshDim_2D ].remove( viscousSM.back() ); + } + if ( !quadFaceSM.empty() ) + FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, pyramidMeshes[iS-1]); + } + if ( !viscousSM.empty() ) + FillNgMesh(occgeo, *_ngMesh, nodeVec, viscousSM, &quadHelper, viscousMesh ); + // fill _ngMesh with faces of sub-meshes err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_2D ], &quadHelper)); - initState = NETGENPlugin_ngMeshInfo(_ngMesh); - //toPython( _ngMesh, "/tmp/ngPython.py"); + initState = NETGENPlugin_ngMeshInfo(_ngMesh, /*checkRemovedElems=*/true); + // toPython( _ngMesh ); } if (!err && _isVolume) { @@ -3554,8 +3680,9 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh ) */ //================================================================================ -NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh): - _copyOfLocalH(0) +NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh, + bool checkRemovedElems): + _elementsRemoved( false ), _copyOfLocalH(0) { if ( ngMesh ) { @@ -3563,6 +3690,10 @@ NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh): _nbSegments = ngMesh->GetNSeg(); _nbFaces = ngMesh->GetNSE(); _nbVolumes = ngMesh->GetNE(); + + if ( checkRemovedElems ) + for ( int i = 1; i <= ngMesh->GetNSE() && !_elementsRemoved; ++i ) + _elementsRemoved = ngMesh->SurfaceElement(i).IsDeleted(); } else { diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx index 01a387f..88e093e 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx @@ -46,14 +46,15 @@ namespace nglib { #include #include +class NETGENPlugin_Hypothesis; +class NETGENPlugin_Internals; +class NETGENPlugin_SimpleHypothesis_2D; class SMESHDS_Mesh; class SMESH_Comment; class SMESH_Mesh; class SMESH_MesherHelper; +class StdMeshers_ViscousLayers; class TopoDS_Shape; -class NETGENPlugin_Hypothesis; -class NETGENPlugin_SimpleHypothesis_2D; -class NETGENPlugin_Internals; namespace netgen { class OCCGeometry; class Mesh; @@ -66,9 +67,10 @@ namespace netgen { struct NETGENPlugin_ngMeshInfo { - int _nbNodes, _nbSegments, _nbFaces, _nbVolumes; + int _nbNodes, _nbSegments, _nbFaces, _nbVolumes; + bool _elementsRemoved; // case where netgen can remove free nodes char* _copyOfLocalH; - NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh=0); + NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh=0, bool checkRemovedElems=false ); void transferLocalH( netgen::Mesh* fromMesh, netgen::Mesh* toMesh ); void restoreLocalH ( netgen::Mesh* ngMesh); }; @@ -120,6 +122,7 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher void SetParameters(const NETGENPlugin_Hypothesis* hyp); void SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp); + void SetParameters(const StdMeshers_ViscousLayers* hyp ); void SetViscousLayers2DAssigned(bool isAssigned) { _isViscousLayers2D = isAssigned; } static void SetLocalSize( netgen::OCCGeometry& occgeo, netgen::Mesh& ngMesh ); @@ -210,6 +213,7 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher volatile double _totalTime; const NETGENPlugin_SimpleHypothesis_2D * _simpleHyp; + const StdMeshers_ViscousLayers* _viscousLayersHyp; // a pointer to NETGENPlugin_Mesher* field of the holder, that will be // nullified at destruction of this diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.cxx index d2381f1..3b5c529 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.cxx @@ -32,10 +32,12 @@ #include "NETGENPlugin_SimpleHypothesis_3D.hxx" #include "NETGENPlugin_Mesher.hxx" +#include +#include #include #include -#include -#include +#include + #include #include @@ -62,6 +64,7 @@ NETGENPlugin_NETGEN_2D3D::NETGENPlugin_NETGEN_2D3D(int hypId, int studyId, _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type _compatibleHypothesis.push_back("NETGEN_Parameters"); _compatibleHypothesis.push_back("NETGEN_SimpleParameters_3D"); + _compatibleHypothesis.push_back( StdMeshers_ViscousLayers::GetHypType() ); _requireDiscreteBoundary = false; _onlyUnaryInput = false; _hypothesis = NULL; @@ -70,7 +73,7 @@ NETGENPlugin_NETGEN_2D3D::NETGENPlugin_NETGEN_2D3D(int hypId, int studyId, //============================================================================= /*! - * + * */ //============================================================================= @@ -81,39 +84,45 @@ NETGENPlugin_NETGEN_2D3D::~NETGENPlugin_NETGEN_2D3D() //============================================================================= /*! - * + * */ //============================================================================= -bool NETGENPlugin_NETGEN_2D3D::CheckHypothesis - (SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape, - SMESH_Hypothesis::Hypothesis_Status& aStatus) +bool NETGENPlugin_NETGEN_2D3D::CheckHypothesis (SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + Hypothesis_Status& aStatus) { - _hypothesis = NULL; - _mesher = NULL; + _hypothesis = NULL; + _viscousLayersHyp = NULL; + _mesher = NULL; - const list& hyps = GetUsedHypothesis(aMesh, aShape); - int nbHyp = hyps.size(); - if (!nbHyp) + const list& hyps = GetUsedHypothesis(aMesh, aShape, /*noAux=*/false); + if ( hyps.empty() ) { aStatus = SMESH_Hypothesis::HYP_OK; return true; // can work with no hypothesis } - const SMESHDS_Hypothesis* theHyp = hyps.front(); // use only the first hypothesis - - string hypName = theHyp->GetName(); - - if ( find( _compatibleHypothesis.begin(), _compatibleHypothesis.end(), - hypName ) != _compatibleHypothesis.end() ) + list::const_iterator h = hyps.begin(); + for ( ; h != hyps.end(); ++h ) { - _hypothesis = theHyp; - aStatus = SMESH_Hypothesis::HYP_OK; - } - else - { - aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE; + const SMESHDS_Hypothesis* aHyp = *h; + std::string hypName = aHyp->GetName(); + + if ( std::find( _compatibleHypothesis.begin(), _compatibleHypothesis.end(), + hypName ) != _compatibleHypothesis.end() ) + { + if ( hypName == StdMeshers_ViscousLayers::GetHypType() ) + _viscousLayersHyp = dynamic_cast( aHyp ); + else + _hypothesis = aHyp; + aStatus = SMESH_Hypothesis::HYP_OK; + } + else + { + aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE; + break; + } } return aStatus == SMESH_Hypothesis::HYP_OK; @@ -133,6 +142,7 @@ bool NETGENPlugin_NETGEN_2D3D::Compute(SMESH_Mesh& aMesh, NETGENPlugin_Mesher mesher(&aMesh, aShape, true); mesher.SetParameters(dynamic_cast(_hypothesis)); mesher.SetParameters(dynamic_cast(_hypothesis)); + mesher.SetParameters(_viscousLayersHyp); mesher.SetSelfPointer( &_mesher ); return mesher.Compute(); } diff --git a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.hxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.hxx index 23ca46f..0728470 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.hxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D3D.hxx @@ -35,6 +35,7 @@ #include class NETGENPlugin_Mesher; +class StdMeshers_ViscousLayers; class NETGENPLUGIN_EXPORT NETGENPlugin_NETGEN_2D3D: public SMESH_3D_Algo { @@ -42,11 +43,11 @@ public: NETGENPlugin_NETGEN_2D3D(int hypId, int studyId, SMESH_Gen* gen); virtual ~NETGENPlugin_NETGEN_2D3D(); - virtual bool CheckHypothesis(SMESH_Mesh& aMesh, + virtual bool CheckHypothesis(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape, - SMESH_Hypothesis::Hypothesis_Status& aStatus); + Hypothesis_Status& aStatus); - virtual bool Compute(SMESH_Mesh& aMesh, + virtual bool Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape); virtual void CancelCompute(); @@ -54,13 +55,14 @@ public: virtual double GetProgress() const; - virtual bool Evaluate(SMESH_Mesh& aMesh, + virtual bool Evaluate(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape, - MapShapeNbElems& aResMap); + MapShapeNbElems& aResMap); protected: - const SMESHDS_Hypothesis* _hypothesis; - NETGENPlugin_Mesher * _mesher; + const SMESHDS_Hypothesis* _hypothesis; + const StdMeshers_ViscousLayers* _viscousLayersHyp; + NETGENPlugin_Mesher * _mesher; }; #endif