Salome HOME
Merge multi-study removal branch.
[plugins/netgenplugin.git] / src / NETGENPlugin / NETGENPlugin_Mesher.cxx
index 6076844cdc510469767d94e1a4df26c4e065b77b..73e8c6d89387e14aa11aab0cc125c5ecad149264 100644 (file)
@@ -82,16 +82,16 @@ namespace netgen {
   extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
 #endif
   //extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh);
-#ifdef NETGEN_V5 && WIN32
+#if defined(NETGEN_V5) && defined(WIN32)
   DLL_HEADER 
 #endif
   extern MeshingParameters mparam;
-#ifdef NETGEN_V5 && WIN32
+#if defined(NETGEN_V5) && defined(WIN32)
   DLL_HEADER 
 #endif
   extern volatile multithreadt multithread;
 
-#ifdef NETGEN_V5 && WIN32
+#if defined(NETGEN_V5) && defined(WIN32)
   DLL_HEADER 
 #endif
   extern bool merge_solids;
@@ -157,6 +157,7 @@ NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh*         mesh,
     _progressTic(1),
     _totalTime(1.0),
     _simpleHyp(NULL),
+    _viscousLayersHyp(NULL),
     _ptrToMe(NULL)
 {
   SetDefaultParameters();
@@ -297,28 +298,23 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
     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();
+    NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin();
+    for ( ; it != localSizes.end() ; it++)
     {
-      const NETGENPlugin_Hypothesis::TLocalSize   localSizes = hyp->GetLocalSizesAndEntries();
-      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);
+      std::string entry = (*it).first;
+      double        val = (*it).second;
+      // --
+      GEOM::GEOM_Object_var aGeomObj;
+      SALOMEDS::SObject_var aSObj = SMESH_Gen_i::getStudyServant()->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);
     }
   }
 }
@@ -336,6 +332,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
@@ -808,7 +815,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 = " <<minh << endl;
   }
   if ( minh > 0.5 * maxSize )
@@ -926,7 +933,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
 
         // get all nodes from connected <edges>
         const bool isQuad = smDS->IsQuadratic();
-        StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, isQuad );
+        StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, isQuad, &helper );
         const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
         if ( points.empty() )
           return false; // invalid node params?
@@ -1054,13 +1061,6 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
 
       // Find solids the geomFace bounds
       int solidID1 = 0, solidID2 = 0;
-      StdMeshers_QuadToTriaAdaptor* quadAdaptor =
-        dynamic_cast<StdMeshers_QuadToTriaAdaptor*>( 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 +1070,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<const SMDS_MeshNode*> fNodes( f->begin_nodes(), f->end_nodes() );
+            std::vector<const SMDS_MeshElement*> 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 +1187,6 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
       cout << "SMESH face " << helper.GetMeshDS()->ShapeToIndex( geomFace )
            << " internal="<<isInternalFace << endl;
 #endif
-      if ( proxyMesh )
-        smDS = proxyMesh->GetSubMesh( geomFace );
 
       SMDS_ElemIteratorPtr faces = smDS->GetElements();
       while ( faces->more() )
@@ -1430,10 +1503,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));
+  // }
 }
 
 //================================================================================
@@ -1641,7 +1719,7 @@ void NETGENPlugin_Mesher::AddIntVerticesInSolids(const netgen::OCCGeometry&
   ofstream py(DUMP_TRIANGLES_SCRIPT);
   py << "import SMESH"<< endl
      << "from salome.smesh import smeshBuilder"<<endl
-     << "smesh = smeshBuilder.New(salome.myStudy)"<<endl
+     << "smesh = smeshBuilder.New()"<<endl
      << "m = smesh.Mesh(name='triangles')" << endl;
 #endif
   if ((int) nodeVec.size() < ngMesh.GetNP() )
@@ -2136,13 +2214,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<double>::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 +2603,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
   // -------------------------
@@ -2791,7 +2889,7 @@ bool NETGENPlugin_Mesher::Compute()
           helper.SetSubShape( F );
           TSideVector wires =
             StdMeshers_FaceSide::GetFaceWires( F, *_mesh, /*skipMediumNodes=*/true,
-                                               error, viscousMesh );
+                                               error, &helper, viscousMesh );
           error = AddSegmentsToMesh( *_ngMesh, occgeo, wires, helper, nodeVec );
 
           if ( !error ) error = SMESH_ComputeError::New();
@@ -2835,37 +2933,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<SMESH_ProxyMesh::Ptr> 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)
     {
@@ -3492,9 +3611,9 @@ 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
+          << "smesh = smeshBuilder.New()" << endl
           << "mesh = smesh.Mesh()" << endl << endl;
 
   using namespace netgen;
@@ -3556,8 +3675,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 +3685,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
   {
@@ -4100,10 +4224,10 @@ void NETGENPlugin_NetgenLibWrapper::removeOutputFile()
     }
     string    tmpDir = SALOMEDS_Tool::GetDirFromPath ( _outputFileName );
     string aFileName = SALOMEDS_Tool::GetNameFromPath( _outputFileName ) + ".out";
-    SALOMEDS::ListOfFileNames_var aFiles = new SALOMEDS::ListOfFileNames;
-    aFiles->length(1);
-    aFiles[0] = aFileName.c_str();
+    SALOMEDS_Tool::ListOfFiles aFiles;
+    aFiles.reserve(1);
+    aFiles.push_back(aFileName.c_str());
 
-    SALOMEDS_Tool::RemoveTemporaryFiles( tmpDir.c_str(), aFiles.in(), true );
+    SALOMEDS_Tool::RemoveTemporaryFiles( tmpDir.c_str(), aFiles, true );
   }
 }