From 1a836a25cc85e8f09fda5e375c1d460aa4a148fd Mon Sep 17 00:00:00 2001 From: eap Date: Mon, 25 Sep 2017 20:48:41 +0300 Subject: [PATCH] refs 254#: 7.5.5. Mesh generation for semi-analytical boundaries Add "Chordal Error" parameter to NETGEN --- .../NETGENPLUGIN/input/netgen_2d_3d_hypo.doc | 19 +- idl/NETGENPlugin_Algorithm.idl | 5 + src/GUI/NETGENPluginGUI_HypothesisCreator.cxx | 80 ++- src/GUI/NETGENPluginGUI_HypothesisCreator.h | 17 +- src/GUI/NETGENPlugin_msg_en.ts | 4 + src/NETGENPlugin/NETGENPluginBuilder.py | 8 + src/NETGENPlugin/NETGENPlugin_Hypothesis.cxx | 80 ++- src/NETGENPlugin/NETGENPlugin_Hypothesis.hxx | 8 + .../NETGENPlugin_Hypothesis_i.cxx | 30 ++ .../NETGENPlugin_Hypothesis_i.hxx | 33 +- src/NETGENPlugin/NETGENPlugin_Mesher.cxx | 479 ++++++++++++++---- src/NETGENPlugin/NETGENPlugin_Mesher.hxx | 17 +- .../NETGENPlugin_NETGEN_2D_ONLY.cxx | 2 + 13 files changed, 630 insertions(+), 152 deletions(-) diff --git a/doc/salome/gui/NETGENPLUGIN/input/netgen_2d_3d_hypo.doc b/doc/salome/gui/NETGENPLUGIN/input/netgen_2d_3d_hypo.doc index 4ef8e06..1113607 100644 --- a/doc/salome/gui/NETGENPLUGIN/input/netgen_2d_3d_hypo.doc +++ b/doc/salome/gui/NETGENPLUGIN/input/netgen_2d_3d_hypo.doc @@ -31,13 +31,19 @@ parameters below. You can select \a Custom to define them manually. - Growth rate - allows to define how much the linear dimensions of two adjacent cells can differ (e.g. 0.3 means 30%). - Nb. Segs per Edge - allows to define the minimum number of -mesh segments in which edges will be split. This parameter is used -only if Limit Size by Surface Curvature is checked. +mesh segments in which edges will be split. Size of elements computed using +this value is trimmed between Min Size and Max Size +bounds. This parameter is used only if Limit Size by Surface + Curvature is checked. - Nb Segs per Radius - allows to define the size of mesh segments and mesh faces in which curved edges and surfaces will -be split. This value divided by a radius of curvature gives an element -size at a given point. This parameter is used only if Limit Size by - Surface Curvature is checked. +be split. A radius of local curvature divided by this value gives an element +size at a given point. Element size computed this way is then trimmed +between Min Size and Max Size bounds. This parameter is +used only if Limit Size by Surface Curvature is checked. +- Chordal Error - allows to define the maximum distance between +the generated 2D element and the surface. Size of elements computed using +this criterion is trimmed between Min Size and Max Size bounds. - Limit Size by Surface Curvature - if this box is checked in, then size of mesh segments and mesh faces on curved edges and surfaces is defined using value of Nb Segs per Radius parameter, and @@ -87,6 +93,9 @@ section.
"25 25 0 25 25 200 0.3" means that along the line between points (25, 25, 0) and (25, 25, 200) size of elements should be 0.3. +Sizes defined in the file are trimmed between Min Size and Max Size +bounds. + \image html netgen2d3d_simple.png NETGEN 2D simple parameters and NETGEN 3D simple diff --git a/idl/NETGENPlugin_Algorithm.idl b/idl/NETGENPlugin_Algorithm.idl index 5a0b252..71b84f5 100644 --- a/idl/NETGENPlugin_Algorithm.idl +++ b/idl/NETGENPlugin_Algorithm.idl @@ -92,6 +92,11 @@ module NETGENPlugin void SetNbSegPerEdge(in double value); double GetNbSegPerEdge(); + void SetChordalErrorEnabled(in boolean value); + boolean GetChordalErrorEnabled(); + void SetChordalError(in double value); + double GetChordalError(); + void SetNbSegPerRadius(in double value); double GetNbSegPerRadius(); diff --git a/src/GUI/NETGENPluginGUI_HypothesisCreator.cxx b/src/GUI/NETGENPluginGUI_HypothesisCreator.cxx index 2db7d1e..e207e26 100644 --- a/src/GUI/NETGENPluginGUI_HypothesisCreator.cxx +++ b/src/GUI/NETGENPluginGUI_HypothesisCreator.cxx @@ -177,7 +177,7 @@ QFrame* NETGENPluginGUI_HypothesisCreator::buildFrame() myFineness = new QComboBox( GroupC1 ); QStringList types; types << tr( "NETGEN_VERYCOARSE" ) << tr( "NETGEN_COARSE" ) << tr( "NETGEN_MODERATE" ) << - tr( "NETGEN_FINE" ) << tr( "NETGEN_VERYFINE" ) << tr( "NETGEN_CUSTOM" ); + tr( "NETGEN_FINE" ) << tr( "NETGEN_VERYFINE" ) << tr( "NETGEN_CUSTOM" ); myFineness->addItems( types ); aGroupLayout->addWidget( myFineness, row, 1 ); connect( myFineness, SIGNAL( activated( int ) ), this, SLOT( onFinenessChanged() ) ); @@ -208,6 +208,19 @@ QFrame* NETGENPluginGUI_HypothesisCreator::buildFrame() row++; } + myChordalErrorEnabled = 0; + myChordalError = 0; + if ( myIs2D || !myIsONLY ) + { + myChordalErrorEnabled = new QCheckBox( tr( "NETGEN_CHORDAL_ERROR" ), GroupC1 ); + aGroupLayout->addWidget( myChordalErrorEnabled, row, 0 ); + myChordalError = new SMESHGUI_SpinBox( GroupC1 ); + myChordalError->RangeStepAndValidator( COORD_MIN, COORD_MAX, .1, "length_precision" ); + aGroupLayout->addWidget( myChordalError, row, 1 ); + connect( myChordalErrorEnabled, SIGNAL( stateChanged(int)), SLOT( onChordalErrorEnabled())); + row++; + } + mySurfaceCurvature = 0; if ( myIs2D || !myIsONLY ) { @@ -337,6 +350,15 @@ void NETGENPluginGUI_HypothesisCreator::retrieveParams() const else myNbSegPerRadius->setText( data.myNbSegPerRadiusVar ); } + if ( myChordalError ) + { + myChordalErrorEnabled->setChecked( data.myChordalErrorEnabled && data.myChordalError > 0 ); + if(data.myChordalErrorVar.isEmpty()) + myChordalError->setValue( data.myChordalError > 0 ? data.myChordalError : 0.1 ); + else + myChordalError->setText( data.myChordalErrorVar ); + myChordalError->setEnabled( myChordalErrorEnabled->isChecked() ); + } if (myAllowQuadrangles) myAllowQuadrangles->setChecked( data.myAllowQuadrangles ); @@ -408,26 +430,28 @@ bool NETGENPluginGUI_HypothesisCreator::readParamsFromHypo( NetgenHypothesisData NETGENPlugin::NETGENPlugin_Hypothesis_var h = NETGENPlugin::NETGENPlugin_Hypothesis::_narrow( initParamsHypothesis() ); - //HypothesisData* data = SMESH::GetHypothesisData( hypType() ); h_data.myName = isCreation() ? hypName() : ""; - h_data.myMaxSize = h->GetMaxSize(); - h_data.myMaxSizeVar = getVariableName("SetMaxSize"); + h_data.myMaxSize = h->GetMaxSize(); + h_data.myMaxSizeVar = getVariableName("SetMaxSize"); h_data.mySecondOrder = h->GetSecondOrder(); - h_data.myOptimize = h->GetOptimize(); - - h_data.myFineness = (int) h->GetFineness(); - h_data.myGrowthRate = h->GetGrowthRate(); - h_data.myGrowthRateVar = getVariableName("SetGrowthRate"); - h_data.myNbSegPerEdge = h->GetNbSegPerEdge(); - h_data.myNbSegPerEdgeVar = getVariableName("SetNbSegPerEdge"); - h_data.myNbSegPerRadius = h->GetNbSegPerRadius(); - h_data.myNbSegPerRadiusVar = getVariableName("SetNbSegPerRadius"); - h_data.myMinSize = h->GetMinSize(); - h_data.myMinSizeVar = getVariableName("SetMinSize"); - h_data.mySurfaceCurvature = h->GetUseSurfaceCurvature(); - h_data.myFuseEdges = h->GetFuseEdges(); - h_data.myMeshSizeFile = h->GetMeshSizeFile(); + h_data.myOptimize = h->GetOptimize(); + + h_data.myFineness = (int) h->GetFineness(); + h_data.myGrowthRate = h->GetGrowthRate(); + h_data.myGrowthRateVar = getVariableName("SetGrowthRate"); + h_data.myNbSegPerEdge = h->GetNbSegPerEdge(); + h_data.myNbSegPerEdgeVar = getVariableName("SetNbSegPerEdge"); + h_data.myNbSegPerRadius = h->GetNbSegPerRadius(); + h_data.myNbSegPerRadiusVar = getVariableName("SetNbSegPerRadius"); + h_data.myChordalError = h->GetChordalError(); + h_data.myChordalErrorVar = getVariableName("SetChordalError"); + h_data.myChordalErrorEnabled = h->GetChordalErrorEnabled(); + h_data.myMinSize = h->GetMinSize(); + h_data.myMinSizeVar = getVariableName("SetMinSize"); + h_data.mySurfaceCurvature = h->GetUseSurfaceCurvature(); + h_data.myFuseEdges = h->GetFuseEdges(); + h_data.myMeshSizeFile = h->GetMeshSizeFile(); //if ( myIs2D ) { @@ -486,6 +510,9 @@ bool NETGENPluginGUI_HypothesisCreator::storeParamsToHypo( const NetgenHypothesi h->SetVarParameter ( h_data.myNbSegPerRadiusVar.toLatin1().constData(), "SetNbSegPerRadius"); h->SetNbSegPerRadius( h_data.myNbSegPerRadius ); } + h->SetVarParameter ( h_data.myChordalErrorVar.toLatin1().constData(), "SetChordalError"); + h->SetChordalError ( h_data.myChordalError ); + h->SetChordalErrorEnabled( h_data.myChordalErrorEnabled ); h->SetVarParameter ( h_data.myMinSizeVar.toLatin1().constData(), "SetMinSize"); h->SetMinSize ( h_data.myMinSize ); h->SetUseSurfaceCurvature( h_data.mySurfaceCurvature ); @@ -551,8 +578,13 @@ bool NETGENPluginGUI_HypothesisCreator::readParamsFromWidgets( NetgenHypothesisD h_data.myNbSegPerEdgeVar = myNbSegPerEdge->text(); if ( myNbSegPerRadius ) h_data.myNbSegPerRadiusVar = myNbSegPerRadius->text(); + if ( myChordalError ) + { + h_data.myChordalErrorVar = myChordalError->text(); + h_data.myChordalError = myChordalError->value(); + h_data.myChordalErrorEnabled = myChordalError->isEnabled(); + } - if ( myAllowQuadrangles ) h_data.myAllowQuadrangles = myAllowQuadrangles->isChecked(); @@ -577,6 +609,11 @@ bool NETGENPluginGUI_HypothesisCreator::readParamsFromWidgets( NetgenHypothesisD return true; } +void NETGENPluginGUI_HypothesisCreator::onChordalErrorEnabled() +{ + myChordalError->setEnabled( myChordalErrorEnabled->isChecked() ); +} + void NETGENPluginGUI_HypothesisCreator::onSurfaceCurvatureChanged() { bool isSurfaceCurvature = (mySurfaceCurvature ? mySurfaceCurvature->isChecked() : true); @@ -587,6 +624,11 @@ void NETGENPluginGUI_HypothesisCreator::onSurfaceCurvatureChanged() myNbSegPerEdge->setEnabled(isCustom && isSurfaceCurvature); if ( myNbSegPerRadius ) myNbSegPerRadius->setEnabled(isCustom && isSurfaceCurvature); + // if ( myChordalError ) + // { + // myChordalError->setEnabled( isSurfaceCurvature ); + // myChordalErrorEnabled->setEnabled( isSurfaceCurvature ); + // } } void NETGENPluginGUI_HypothesisCreator::onFinenessChanged() diff --git a/src/GUI/NETGENPluginGUI_HypothesisCreator.h b/src/GUI/NETGENPluginGUI_HypothesisCreator.h index 64a1616..e471ba1 100644 --- a/src/GUI/NETGENPluginGUI_HypothesisCreator.h +++ b/src/GUI/NETGENPluginGUI_HypothesisCreator.h @@ -44,11 +44,11 @@ class QTableWidget; typedef struct { - double myMaxSize, myMinSize, myGrowthRate, myNbSegPerEdge, myNbSegPerRadius; - int myFineness; - bool mySecondOrder, myAllowQuadrangles, myOptimize, mySurfaceCurvature, myFuseEdges; - QString myName, myMeshSizeFile; - QString myMaxSizeVar, myMinSizeVar, myGrowthRateVar, myNbSegPerEdgeVar, myNbSegPerRadiusVar; + double myMaxSize, myMinSize, myGrowthRate, myNbSegPerEdge, myNbSegPerRadius, myChordalError; + int myFineness; + bool mySecondOrder, myAllowQuadrangles, myOptimize, mySurfaceCurvature, myFuseEdges, myChordalErrorEnabled; + QString myName, myMeshSizeFile; + QString myMaxSizeVar, myMinSizeVar, myGrowthRateVar, myNbSegPerEdgeVar, myNbSegPerRadiusVar, myChordalErrorVar; } NetgenHypothesisData; /*! @@ -76,6 +76,7 @@ protected: protected slots: virtual void onFinenessChanged(); + virtual void onChordalErrorEnabled(); virtual void onSurfaceCurvatureChanged(); virtual void onAddLocalSizeOnVertex(); virtual void onAddLocalSizeOnEdge(); @@ -102,12 +103,14 @@ private: SMESHGUI_SpinBox* myGrowthRate; SMESHGUI_SpinBox* myNbSegPerEdge; SMESHGUI_SpinBox* myNbSegPerRadius; + QCheckBox* myChordalErrorEnabled; + SMESHGUI_SpinBox* myChordalError; QCheckBox* myAllowQuadrangles; QCheckBox* mySurfaceCurvature; QCheckBox* myFuseEdges; - bool myIs2D; - bool myIsONLY; + bool myIs2D; // 2D or 3D + bool myIsONLY; // one dim or several QLineEdit* myMeshSizeFile; QTableWidget* myLocalSizeTable; diff --git a/src/GUI/NETGENPlugin_msg_en.ts b/src/GUI/NETGENPlugin_msg_en.ts index 478cc25..d73590d 100644 --- a/src/GUI/NETGENPlugin_msg_en.ts +++ b/src/GUI/NETGENPlugin_msg_en.ts @@ -91,6 +91,10 @@ NETGEN_SEG_PER_RADIUS Nb. Segs per Radius + + NETGEN_CHORDAL_ERROR + Chordal Error + NETGEN_SURFACE_CURVATURE Limit Size by Surface Curvature diff --git a/src/NETGENPlugin/NETGENPluginBuilder.py b/src/NETGENPlugin/NETGENPluginBuilder.py index d5f9460..7a2c7fa 100644 --- a/src/NETGENPlugin/NETGENPluginBuilder.py +++ b/src/NETGENPlugin/NETGENPluginBuilder.py @@ -221,6 +221,14 @@ class NETGEN_1D2D3D_Algorithm(NETGEN_Algorithm): if self.Parameters(): self.params.SetNbSegPerRadius(theVal) pass + ## Sets @c ChordalError parameter + # @param theVal new value of the @c ChordalError parameter + def SetChordalError(self, theVal): + if self.Parameters(): + self.params.SetChordalError(theVal) + self.params.SetChordalErrorEnabled( theVal > 0 ) + pass + ## Sets @c QuadAllowed flag # @param toAllow new value of the @c QuadAllowed parameter (@c True by default) def SetQuadAllowed(self, toAllow=True): diff --git a/src/NETGENPlugin/NETGENPlugin_Hypothesis.cxx b/src/NETGENPlugin/NETGENPlugin_Hypothesis.cxx index fe4c2ee..9bd9ec8 100644 --- a/src/NETGENPlugin/NETGENPlugin_Hypothesis.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Hypothesis.cxx @@ -43,18 +43,20 @@ using namespace std; NETGENPlugin_Hypothesis::NETGENPlugin_Hypothesis (int hypId, int studyId, SMESH_Gen * gen) : SMESH_Hypothesis(hypId, studyId, gen), - _maxSize (GetDefaultMaxSize()), - _minSize (0), - _growthRate (GetDefaultGrowthRate()), - _nbSegPerEdge (GetDefaultNbSegPerEdge()), - _nbSegPerRadius (GetDefaultNbSegPerRadius()), - _fineness (GetDefaultFineness()), - _secondOrder (GetDefaultSecondOrder()), - _optimize (GetDefaultOptimize()), - _localSize (GetDefaultLocalSize()), - _quadAllowed (GetDefaultQuadAllowed()), - _surfaceCurvature(GetDefaultSurfaceCurvature()), - _fuseEdges (GetDefaultFuseEdges()) + _maxSize (GetDefaultMaxSize()), + _minSize (0), + _growthRate (GetDefaultGrowthRate()), + _nbSegPerEdge (GetDefaultNbSegPerEdge()), + _nbSegPerRadius (GetDefaultNbSegPerRadius()), + _fineness (GetDefaultFineness()), + _chordalErrorEnabled(GetDefaultChordalError() > 0), + _chordalError (GetDefaultChordalError() ), + _secondOrder (GetDefaultSecondOrder()), + _optimize (GetDefaultOptimize()), + _localSize (GetDefaultLocalSize()), + _quadAllowed (GetDefaultQuadAllowed()), + _surfaceCurvature (GetDefaultSurfaceCurvature()), + _fuseEdges (GetDefaultFuseEdges()) { _name = "NETGEN_Parameters"; _param_algo_dim = 3; @@ -208,6 +210,34 @@ void NETGENPlugin_Hypothesis::SetNbSegPerRadius(double theVal) } } +//============================================================================= +/*! + * + */ +//============================================================================= +void NETGENPlugin_Hypothesis::SetChordalErrorEnabled(bool theVal) +{ + if (theVal != _chordalErrorEnabled) + { + _chordalErrorEnabled = theVal; + NotifySubMeshesHypothesisModification(); + } +} + +//============================================================================= +/*! + * + */ +//============================================================================= +void NETGENPlugin_Hypothesis::SetChordalError(double theVal) +{ + if (theVal != _chordalError) + { + _chordalError = theVal; + NotifySubMeshesHypothesisModification(); + } +} + //============================================================================= /*! * @@ -351,8 +381,8 @@ ostream & NETGENPlugin_Hypothesis::SaveTo(ostream & save) if (it_sm != _localSize.end()) { save << " " << "__LOCALSIZE_BEGIN__"; for ( ; it_sm != _localSize.end(); ++it_sm ) { - save << " " << it_sm->first - << " " << it_sm->second << "%#"; // "%#" is a mark of value end + save << " " << it_sm->first + << " " << it_sm->second << "%#"; // "%#" is a mark of value end } save << " " << "__LOCALSIZE_END__"; } @@ -363,6 +393,8 @@ ostream & NETGENPlugin_Hypothesis::SaveTo(ostream & save) save << " " << _meshSizeFile.size() << " " << _meshSizeFile; + save << " " << ( _chordalErrorEnabled ? _chordalError : 0. ); + return save; } @@ -476,12 +508,19 @@ istream & NETGENPlugin_Hypothesis::LoadFrom(istream & load) load.get( &_meshSizeFile[0], is+1 ); } + isOK = static_cast(load >> val); + if (isOK) + _chordalError = val; + else + load.clear(ios::badbit | load.rdstate()); + _chordalErrorEnabled = ( _chordalError > 0 ); + return load; } //============================================================================= /*! - * + * */ //============================================================================= ostream & operator <<(ostream & save, NETGENPlugin_Hypothesis & hyp) @@ -491,7 +530,7 @@ ostream & operator <<(ostream & save, NETGENPlugin_Hypothesis & hyp) //============================================================================= /*! - * + * */ //============================================================================= istream & operator >>(istream & load, NETGENPlugin_Hypothesis & hyp) @@ -584,6 +623,15 @@ double NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius() { return 2; } +//============================================================================= +/*! + * + */ +//============================================================================= +double NETGENPlugin_Hypothesis::GetDefaultChordalError() +{ + return -1; // disabled by default +} //============================================================================= /*! diff --git a/src/NETGENPlugin/NETGENPlugin_Hypothesis.hxx b/src/NETGENPlugin/NETGENPlugin_Hypothesis.hxx index 7d6da32..ec2fbb5 100644 --- a/src/NETGENPlugin/NETGENPlugin_Hypothesis.hxx +++ b/src/NETGENPlugin/NETGENPlugin_Hypothesis.hxx @@ -83,6 +83,11 @@ public: void SetNbSegPerRadius(double theVal); double GetNbSegPerRadius() const { return _nbSegPerRadius; } + void SetChordalErrorEnabled(bool value); + double GetChordalErrorEnabled() const { return _chordalErrorEnabled; } + void SetChordalError(double value); + double GetChordalError() const { return _chordalError; } + typedef std::map TLocalSize; static TLocalSize GetDefaultLocalSize() { return TLocalSize(); } void SetLocalSizeOnEntry(const std::string& entry, double localSize); @@ -109,6 +114,7 @@ public: static double GetDefaultGrowthRate(); static double GetDefaultNbSegPerEdge(); static double GetDefaultNbSegPerRadius(); + static double GetDefaultChordalError(); static bool GetDefaultSecondOrder(); static bool GetDefaultOptimize(); static bool GetDefaultQuadAllowed(); @@ -141,6 +147,8 @@ private: double _nbSegPerEdge; double _nbSegPerRadius; Fineness _fineness; + bool _chordalErrorEnabled; + double _chordalError; bool _secondOrder; bool _optimize; TLocalSize _localSize; diff --git a/src/NETGENPlugin/NETGENPlugin_Hypothesis_i.cxx b/src/NETGENPlugin/NETGENPlugin_Hypothesis_i.cxx index d60d666..11461d6 100644 --- a/src/NETGENPlugin/NETGENPlugin_Hypothesis_i.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Hypothesis_i.cxx @@ -307,6 +307,36 @@ CORBA::Double NETGENPlugin_Hypothesis_i::GetNbSegPerRadius() //============================================================================= +void NETGENPlugin_Hypothesis_i::SetChordalErrorEnabled(CORBA::Boolean theValue) +{ + if ( isToSetParameter( GetChordalErrorEnabled(), theValue, METH_SetChordalErrorEnabled )) + { + this->GetImpl()->SetChordalErrorEnabled(theValue); + SMESH::TPythonDump() << _this() << ".SetChordalErrorEnabled( " << theValue << " )"; + } +} + +CORBA::Boolean NETGENPlugin_Hypothesis_i::GetChordalErrorEnabled() +{ + return GetImpl()->GetChordalErrorEnabled(); +} + +void NETGENPlugin_Hypothesis_i::SetChordalError(CORBA::Double theValue) +{ + if ( isToSetParameter( GetChordalError(), theValue, METH_SetChordalError )) + { + this->GetImpl()->SetChordalError(theValue); + SMESH::TPythonDump() << _this() << ".SetChordalError( " << SMESH::TVar(theValue) << " )"; + } +} + +CORBA::Double NETGENPlugin_Hypothesis_i::GetChordalError() +{ + return GetImpl()->GetChordalError(); +} + +//============================================================================= + void NETGENPlugin_Hypothesis_i::SetLocalSizeOnShape(GEOM::GEOM_Object_ptr GeomObj, CORBA::Double localSize) throw (SALOME::SALOME_Exception) diff --git a/src/NETGENPlugin/NETGENPlugin_Hypothesis_i.hxx b/src/NETGENPlugin/NETGENPlugin_Hypothesis_i.hxx index 12ffa7b..e9b908d 100644 --- a/src/NETGENPlugin/NETGENPlugin_Hypothesis_i.hxx +++ b/src/NETGENPlugin/NETGENPlugin_Hypothesis_i.hxx @@ -79,6 +79,11 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Hypothesis_i: void SetNbSegPerRadius(CORBA::Double theVal); CORBA::Double GetNbSegPerRadius(); + void SetChordalErrorEnabled(CORBA::Boolean value); + CORBA::Boolean GetChordalErrorEnabled(); + void SetChordalError(CORBA::Double value); + CORBA::Double GetChordalError(); + void SetLocalSizeOnShape(GEOM::GEOM_Object_ptr GeomObj, CORBA::Double localSize) throw (SALOME::SALOME_Exception); void SetLocalSizeOnEntry(const char* entry, CORBA::Double localSize); @@ -109,19 +114,21 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Hypothesis_i: // to remember whether a parameter is already set (issue 0021364) enum SettingMethod { - METH_SetMaxSize = 1, - METH_SetMinSize = 2, - METH_SetSecondOrder = 4, - METH_SetOptimize = 8, - METH_SetFineness = 16, - METH_SetGrowthRate = 32, - METH_SetNbSegPerEdge = 64, - METH_SetNbSegPerRadius = 128, - METH_SetLocalSizeOnEntry = 256, - METH_SetQuadAllowed = METH_SetLocalSizeOnEntry * 2, - METH_SetSurfaceCurvature = METH_SetQuadAllowed * 2, - METH_SetFuseEdges = METH_SetSurfaceCurvature * 2, - METH_LAST = METH_SetFuseEdges + METH_SetMaxSize = 1, + METH_SetMinSize = 2, + METH_SetSecondOrder = 4, + METH_SetOptimize = 8, + METH_SetFineness = 16, + METH_SetGrowthRate = 32, + METH_SetNbSegPerEdge = 64, + METH_SetNbSegPerRadius = 128, + METH_SetLocalSizeOnEntry = 256, + METH_SetQuadAllowed = METH_SetLocalSizeOnEntry * 2, + METH_SetSurfaceCurvature = METH_SetQuadAllowed * 2, + METH_SetFuseEdges = METH_SetSurfaceCurvature * 2, + METH_SetChordalErrorEnabled = METH_SetFuseEdges * 2, + METH_SetChordalError = METH_SetChordalErrorEnabled * 2, + METH_LAST = METH_SetFuseEdges }; int mySetMethodFlags; diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx index ecf18ad..1fa59e4 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx @@ -52,21 +52,29 @@ #include +#include #include +#include +#include +#include #include #include +#include #include +#include #include #include #include #include #include +#include #include #include #include #include #include #include +#include // Netgen include files #ifndef OCCGEOMETRY @@ -151,12 +159,14 @@ NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh, _optimize(true), _fineness(NETGENPlugin_Hypothesis::GetDefaultFineness()), _isViscousLayers2D(false), + _chordalError(-1), // means disabled _ngMesh(NULL), _occgeom(NULL), _curShapeIndex(-1), _progressTic(1), _totalTime(1.0), _simpleHyp(NULL), + _viscousLayersHyp(NULL), _ptrToMe(NULL) { SetDefaultParameters(); @@ -292,32 +302,36 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp) _fineness = hyp->GetFineness(); mparams.uselocalh = hyp->GetSurfaceCurvature(); netgen::merge_solids = hyp->GetFuseEdges(); + _chordalError = hyp->GetChordalErrorEnabled() ? hyp->GetChordalError() : -1.; _simpleHyp = NULL; // mesh size file mparams.meshsizefilename= hyp->GetMeshSizeFile().empty() ? 0 : hyp->GetMeshSizeFile().c_str(); - SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen(); - CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager"); - SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject); - SALOMEDS::Study_var myStudy = aStudyMgr->GetStudyByID(hyp->GetStudyId()); - if ( !myStudy->_is_nil() ) + const NETGENPlugin_Hypothesis::TLocalSize& localSizes = hyp->GetLocalSizesAndEntries(); + if ( !localSizes.empty() ) { - const NETGENPlugin_Hypothesis::TLocalSize localSizes = hyp->GetLocalSizesAndEntries(); - NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin(); - for ( ; it != localSizes.end() ; it++) + SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen(); + CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager"); + SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject); + SALOMEDS::Study_var myStudy = aStudyMgr->GetStudyByID(hyp->GetStudyId()); + if ( !myStudy->_is_nil() ) { - std::string entry = (*it).first; - double val = (*it).second; - // -- - GEOM::GEOM_Object_var aGeomObj; - SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() ); - if ( !aSObj->_is_nil() ) { - CORBA::Object_var obj = aSObj->GetObject(); - aGeomObj = GEOM::GEOM_Object::_narrow(obj); - aSObj->UnRegister(); + NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin(); + for ( ; it != localSizes.end() ; it++) + { + std::string entry = (*it).first; + double val = (*it).second; + // -- + GEOM::GEOM_Object_var aGeomObj; + SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() ); + if ( !aSObj->_is_nil() ) { + CORBA::Object_var obj = aSObj->GetObject(); + aGeomObj = GEOM::GEOM_Object::_narrow(obj); + aSObj->UnRegister(); + } + TopoDS_Shape S = smeshGen_i->GeomObjectToShape( aGeomObj.in() ); + ::SetLocalSize(S, val); } - TopoDS_Shape S = smeshGen_i->GeomObjectToShape( aGeomObj.in() ); - ::SetLocalSize(S, val); } } } @@ -336,6 +350,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 @@ -585,7 +610,8 @@ namespace void setLocalSize(const TopoDS_Edge& edge, double size, - netgen::Mesh& mesh) + netgen::Mesh& mesh, + const bool overrideMinH = true) { if ( size <= std::numeric_limits::min() ) return; @@ -596,7 +622,7 @@ namespace TopoDS_Iterator vIt( edge ); if ( !vIt.More() ) return; gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vIt.Value() )); - NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size ); + NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size, overrideMinH ); } else { @@ -606,15 +632,29 @@ namespace { Standard_Real u = u1 + delta*i; gp_Pnt p = curve->Value(u); - NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size ); + NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size, overrideMinH ); netgen::Point3d pi(p.X(), p.Y(), p.Z()); double resultSize = mesh.GetH(pi); if ( resultSize - size > 0.1*size ) // netgen does restriction iff oldH/newH > 1.2 (localh.cpp:136) - NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), resultSize/1.201 ); + NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), resultSize/1.201, overrideMinH ); } } } + + //================================================================================ + /*! + * \brief Return triangle size for a given chordalError and radius of curvature + */ + //================================================================================ + + double elemSizeForChordalError( double chordalError, double radius ) + { + if ( 2 * radius < chordalError ) + return 1.5 * radius; + return Sqrt( 3 ) * Sqrt( chordalError * ( 2 * radius - chordalError )); + } + } // namespace //================================================================================ @@ -624,16 +664,19 @@ namespace //================================================================================ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo, - netgen::Mesh& ngMesh ) + netgen::Mesh& ngMesh) { - for(std::map::const_iterator it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++) + // edges + std::map::const_iterator it; + for( it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++) { int key = (*it).first; double hi = (*it).second; const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key); setLocalSize( TopoDS::Edge(shape), hi, ngMesh ); } - for(std::map::const_iterator it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++) + // vertices + for(it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++) { int key = (*it).first; double hi = (*it).second; @@ -641,7 +684,8 @@ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo, gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex(shape) ); NETGENPlugin_Mesher::RestrictLocalSize( ngMesh, p.XYZ(), hi ); } - for(map::const_iterator it=FaceId2LocalSize.begin(); it!=FaceId2LocalSize.end(); it++) + // faces + for(it=FaceId2LocalSize.begin(); it!=FaceId2LocalSize.end(); it++) { int key = (*it).first; double val = (*it).second; @@ -659,7 +703,8 @@ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo, ShapesWithControlPoints.insert( key ); } } - for(map::const_iterator it=SolidId2LocalSize.begin(); it!=SolidId2LocalSize.end(); it++) + //solids + for(it=SolidId2LocalSize.begin(); it!=SolidId2LocalSize.end(); it++) { int key = (*it).first; double val = (*it).second; @@ -676,6 +721,146 @@ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo, for ( size_t i = 0; i < ControlPoints.size(); ++i ) NETGENPlugin_Mesher::RestrictLocalSize( ngMesh, ControlPoints[i].XYZ(), ControlPoints[i].Size() ); } + return; +} + +//================================================================================ +/*! + * \brief Restrict local size to achieve a required _chordalError + */ +//================================================================================ + +void NETGENPlugin_Mesher::SetLocalSizeForChordalError( netgen::OCCGeometry& occgeo, + netgen::Mesh& ngMesh) +{ + if ( _chordalError <= 0. ) + return; + + TopLoc_Location loc; + BRepLProp_SLProps surfProp( 2, 1e-6 ); + const double sizeCoef = 0.95; + + // find non-planar FACEs with non-constant curvature + std::vector fInd; + for ( int i = 1; i <= occgeo.fmap.Extent(); ++i ) + { + const TopoDS_Face& face = TopoDS::Face( occgeo.fmap( i )); + BRepAdaptor_Surface surfAd( face, false ); + switch ( surfAd.GetType() ) + { + case GeomAbs_Plane: + continue; + case GeomAbs_Cylinder: + case GeomAbs_Sphere: + case GeomAbs_Torus: // constant curvature + { + surfProp.SetSurface( surfAd ); + surfProp.SetParameters( 0, 0 ); + double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() )); + double size = elemSizeForChordalError( _chordalError, 1 / maxCurv ); + occgeo.SetFaceMaxH( i, size * sizeCoef ); + // limit size one edges + TopTools_MapOfShape edgeMap; + for ( TopExp_Explorer eExp( face, TopAbs_EDGE ); eExp.More(); eExp.Next() ) + if ( edgeMap.Add( eExp.Current() )) + setLocalSize( TopoDS::Edge( eExp.Current() ), size, ngMesh, /*overrideMinH=*/false ); + break; + } + default: + Handle(Geom_Surface) surf = BRep_Tool::Surface( face, loc ); + if ( GeomLib_IsPlanarSurface( surf ).IsPlanar() ) + continue; + fInd.push_back( i ); + } + } + // set local size + if ( !fInd.empty() ) + { + BRep_Builder b; + TopoDS_Compound allFacesComp; + b.MakeCompound( allFacesComp ); + for ( size_t i = 0; i < fInd.size(); ++i ) + b.Add( allFacesComp, occgeo.fmap( fInd[i] )); + + // copy the shape to avoid spoiling its triangulation + TopoDS_Shape allFacesCompCopy = BRepBuilderAPI_Copy( allFacesComp ); + + // create triangulation with desired chordal error + BRepMesh_IncrementalMesh( allFacesCompCopy, + _chordalError, + /*isRelative = */Standard_False, + /*theAngDeflection = */ 0.5, + /*isInParallel = */Standard_True); + + // loop on FACEs + for ( TopExp_Explorer fExp( allFacesCompCopy, TopAbs_FACE ); fExp.More(); fExp.Next() ) + { + const TopoDS_Face& face = TopoDS::Face( fExp.Current() ); + Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation ( face, loc ); + if ( triangulation.IsNull() ) continue; + + BRepAdaptor_Surface surf( face, false ); + surfProp.SetSurface( surf ); + + gp_XY uv[3]; + gp_XYZ p[3]; + double size[3]; + for ( int i = 1; i <= triangulation->NbTriangles(); ++i ) + { + Standard_Integer n1,n2,n3; + triangulation->Triangles()(i).Get( n1,n2,n3 ); + p [0] = triangulation->Nodes()(n1).Transformed(loc).XYZ(); + p [1] = triangulation->Nodes()(n2).Transformed(loc).XYZ(); + p [2] = triangulation->Nodes()(n3).Transformed(loc).XYZ(); + uv[0] = triangulation->UVNodes()(n1).XY(); + uv[1] = triangulation->UVNodes()(n2).XY(); + uv[2] = triangulation->UVNodes()(n3).XY(); + surfProp.SetParameters( uv[0].X(), uv[0].Y() ); + if ( !surfProp.IsCurvatureDefined() ) + break; + + for ( int n = 0; n < 3; ++n ) // get size at triangle nodes + { + surfProp.SetParameters( uv[n].X(), uv[n].Y() ); + double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() )); + size[n] = elemSizeForChordalError( _chordalError, 1 / maxCurv ); + } + for ( int n1 = 0; n1 < 3; ++n1 ) // limit size along each triangle edge + { + int n2 = ( n1 + 1 ) % 3; + double minSize = size[n1], maxSize = size[n2]; + if ( size[n1] > size[n2] ) + minSize = size[n2], maxSize = size[n1]; + + if ( maxSize / minSize < 1.2 ) // netgen ignores size difference < 1.2 + { + ngMesh.RestrictLocalHLine ( netgen::Point3d( p[n1].X(), p[n1].Y(), p[n1].Z() ), + netgen::Point3d( p[n2].X(), p[n2].Y(), p[n2].Z() ), + sizeCoef * minSize ); + } + else + { + gp_XY uvVec( uv[n2] - uv[n1] ); + double len = ( p[n1] - p[n2] ).Modulus(); + int nb = int( len / minSize ) + 1; + for ( int j = 0; j <= nb; ++j ) + { + double r = double( j ) / nb; + gp_XY uvj = uv[n1] + r * uvVec; + + surfProp.SetParameters( uvj.X(), uvj.Y() ); + double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() )); + double h = elemSizeForChordalError( _chordalError, 1 / maxCurv ); + + const gp_Pnt& pj = surfProp.Value(); + netgen::Point3d ngP( pj.X(), pj.Y(), pj.Z()); + ngMesh.RestrictLocalH( ngP, h * sizeCoef ); + } + } + } + } + } + } } //================================================================================ @@ -808,7 +993,7 @@ double NETGENPlugin_Mesher::GetDefaultMinSize(const TopoDS_Shape& geom, } else { - minh = 3 * sqrt( minh ); // triangulation for visualization is rather fine + minh = sqrt( minh ); // triangulation for visualization is rather fine //cout << "TRIANGULATION minh = " < 0.5 * maxSize ) @@ -926,7 +1111,8 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry& occgeom, // get all nodes from connected const bool isQuad = smDS->IsQuadratic(); - StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, isQuad ); + //StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, isQuad, &helper ); -- master + StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, isQuad ); // -- V8_2_BR const vector& points = fSide.GetUVPtStruct(); if ( points.empty() ) return false; // invalid node params? @@ -1054,13 +1240,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 +1249,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 +1366,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 +1682,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 +2393,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; @@ -2503,8 +2782,6 @@ bool NETGENPlugin_Mesher::Compute() SMESH_MesherHelper quadHelper( *_mesh ); quadHelper.SetIsQuadratic( mparams.secondorder ); - static string debugFile = "/tmp/ngMesh.py"; /* to call toPython( _ngMesh, debugFile ) - while debugging netgen */ // ------------------------- // Prepare OCC geometry // ------------------------- @@ -2621,6 +2898,7 @@ bool NETGENPlugin_Mesher::Compute() { // Local size on shapes SetLocalSize( occgeo, *_ngMesh ); + SetLocalSizeForChordalError( occgeo, *_ngMesh ); } // Precompute internal edges (issue 0020676) in order to @@ -2791,7 +3069,8 @@ bool NETGENPlugin_Mesher::Compute() helper.SetSubShape( F ); TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *_mesh, /*skipMediumNodes=*/true, - error, viscousMesh ); + error, viscousMesh ); // -- V8_2_BR + // error, &helper, viscousMesh ); -- master error = AddSegmentsToMesh( *_ngMesh, occgeo, wires, helper, nodeVec ); if ( !error ) error = SMESH_ComputeError::New(); @@ -2835,37 +3114,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() ) + 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() )) { - 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); + 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) { @@ -3166,25 +3466,25 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap) sm->GetComputeError().reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED )); return false; } - if ( _simpleHyp ) - { - // Pass 1D simple parameters to NETGEN - // -------------------------------- - int nbSeg = _simpleHyp->GetNumberOfSegments(); - double segSize = _simpleHyp->GetLocalLength(); - for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE ) - { - const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE)); - if ( nbSeg ) - segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 ); - setLocalSize( e, segSize, *ngMesh ); - } - } - else // if ( ! _simpleHyp ) - { - // Local size on shapes - SetLocalSize( occgeo, *ngMesh ); - } + // if ( _simpleHyp ) + // { + // // Pass 1D simple parameters to NETGEN + // // -------------------------------- + // int nbSeg = _simpleHyp->GetNumberOfSegments(); + // double segSize = _simpleHyp->GetLocalLength(); + // for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE ) + // { + // const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE)); + // if ( nbSeg ) + // segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 ); + // setLocalSize( e, segSize, *ngMesh ); + // } + // } + // else // if ( ! _simpleHyp ) + // { + // // Local size on shapes + // SetLocalSize( occgeo, *ngMesh ); + // } // calculate total nb of segments and length of edges double fullLen = 0.0; int fullNbSeg = 0; @@ -3492,7 +3792,7 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh ) ofstream outfile( pyFile, ios::out ); if ( !outfile ) return; - outfile << "import SMESH" << endl + outfile << "import salome, SMESH" << endl << "from salome.smesh import smeshBuilder" << endl << "smesh = smeshBuilder.New(salome.myStudy)" << endl << "mesh = smesh.Mesh()" << endl << endl; @@ -3556,8 +3856,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 ) { @@ -3565,6 +3866,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..41d9aa5 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,7 +122,10 @@ 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; } + + void SetLocalSizeForChordalError( netgen::OCCGeometry& occgeo, netgen::Mesh& ngMesh ); static void SetLocalSize( netgen::OCCGeometry& occgeo, netgen::Mesh& ngMesh ); bool Compute(); @@ -201,6 +206,7 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher bool _optimize; int _fineness; bool _isViscousLayers2D; + double _chordalError; netgen::Mesh* _ngMesh; netgen::OCCGeometry* _occgeom; @@ -210,6 +216,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_2D_ONLY.cxx b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx index b4c7919..907c93c 100644 --- a/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx +++ b/src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx @@ -319,6 +319,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, // set local size defined on shapes aMesher.SetLocalSize( occgeoComm, *ngMeshes[0] ); + aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes[0] ); try { ngMeshes[0]->LoadLocalMeshSize( mparam.meshsizefilename ); } catch (NgException & ex) { @@ -460,6 +461,7 @@ bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh& aMesh, bb.Increase (bb.Diam()/10); ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading); aMesher.SetLocalSize( occgeom, *ngMesh ); + aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMesh ); try { ngMesh->LoadLocalMeshSize( mparam.meshsizefilename ); } catch (NgException & ex) { -- 2.39.2