src/NETGENPlugin/NETGENPlugin_Mesher.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 //  NETGENPlugin : C++ implementation
  23 // File      : NETGENPlugin_Mesher.cxx
  24 // Author    : Michael Sazonov (OCN)
  25 // Date      : 31/03/2006
  26 // Project   : SALOME
  27 //=============================================================================
  28 //
  29 #include "NETGENPlugin_Mesher.hxx"
  30 #include "NETGENPlugin_Hypothesis_2D.hxx"
  31 #include "NETGENPlugin_SimpleHypothesis_3D.hxx"
  32
  33 #include <SMESH_Mesh.hxx>
  34 #include <SMESH_Comment.hxx>
  35 #include <SMESH_ComputeError.hxx>
  36 #include <SMESH_subMesh.hxx>
  37 #include <SMESH_MesherHelper.hxx>
  38 #include <SMESHDS_Mesh.hxx>
  39 #include <SMDS_MeshElement.hxx>
  40 #include <SMDS_MeshNode.hxx>
  41 #include <utilities.h>
  42
  43 #include <vector>
  44
  45 #include <BRep_Tool.hxx>
  46 #include <TopExp.hxx>
  47 #include <TopExp_Explorer.hxx>
  48 #include <TopoDS.hxx>
  49 #include <NCollection_Map.hxx>
  50 #include <OSD_Path.hxx>
  51 #include <OSD_File.hxx>
  52 #include <TCollection_AsciiString.hxx>
  53 #include <TopTools_ListIteratorOfListOfShape.hxx>
  54 #include <Standard_ErrorHandler.hxx>
  55
  56 // Netgen include files
  57 namespace nglib {
  58 #include <nglib.h>
  59 }
  60 #define OCCGEOMETRY
  61 #include <occgeom.hpp>
  62 #include <meshing.hpp>
  63 //#include <ngexception.hpp>
  64 namespace netgen {
  65   extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
  66   extern MeshingParameters mparam;
  67 }
  68
  69 using namespace std;
  70
  71 //=============================================================================
  72 /*!
  73  *
  74  */
  75 //=============================================================================
  76
  77 NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh,
  78                                           const TopoDS_Shape& aShape,
  79                                           const bool isVolume)
  80   : _mesh    (mesh),
  81     _shape   (aShape),
  82     _isVolume(isVolume),
  83     _optimize(true),
  84     _simpleHyp(NULL)
  85 {
  86   defaultParameters();
  87 }
  88
  89 //================================================================================
  90 /*!
  91  * \brief Initialize global NETGEN parameters with default values
  92  */
  93 //================================================================================
  94
  95 void NETGENPlugin_Mesher::defaultParameters()
  96 {
  97 #ifdef WNT
  98   netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
  99 #else
 100   netgen::MeshingParameters& mparams = netgen::mparam;
 101 #endif
 102   // maximal mesh edge size
 103   mparams.maxh = NETGENPlugin_Hypothesis::GetDefaultMaxSize();
 104   // minimal number of segments per edge
 105   mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
 106   // rate of growth of size between elements
 107   mparams.grading = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
 108   // safety factor for curvatures (elements per radius)
 109   mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
 110   // create elements of second order
 111   mparams.secondorder = NETGENPlugin_Hypothesis::GetDefaultSecondOrder() ? 1 : 0;
 112   // quad-dominated surface meshing
 113   if (_isVolume)
 114     mparams.quad = 0;
 115   else
 116     mparams.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed() ? 1 : 0;
 117 }
 118
 119 //=============================================================================
 120 /*!
 121  * Pass parameters to NETGEN
 122  */
 123 //=============================================================================
 124 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
 125 {
 126   if (hyp)
 127   {
 128 #ifdef WNT
 129     netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
 130 #else
 131     netgen::MeshingParameters& mparams = netgen::mparam;
 132 #endif
 133     // Initialize global NETGEN parameters:
 134     // maximal mesh segment size
 135     mparams.maxh = hyp->GetMaxSize();
 136     // minimal number of segments per edge
 137     mparams.segmentsperedge = hyp->GetNbSegPerEdge();
 138     // rate of growth of size between elements
 139     mparams.grading = hyp->GetGrowthRate();
 140     // safety factor for curvatures (elements per radius)
 141     mparams.curvaturesafety = hyp->GetNbSegPerRadius();
 142     // create elements of second order
 143     mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0;
 144     // quad-dominated surface meshing
 145     // only triangles are allowed for volumic mesh
 146     if (!_isVolume)
 147       mparams.quad = static_cast<const NETGENPlugin_Hypothesis_2D*>
 148         (hyp)->GetQuadAllowed() ? 1 : 0;
 149     _optimize = hyp->GetOptimize();
 150     _simpleHyp = NULL;
 151   }
 152 }
 153
 154 //=============================================================================
 155 /*!
 156  * Pass simple parameters to NETGEN
 157  */
 158 //=============================================================================
 159
 160 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp)
 161 {
 162   _simpleHyp = hyp;
 163   if ( _simpleHyp )
 164     defaultParameters();
 165 }
 166
 167 //=============================================================================
 168 /*!
 169  *  Link - a pair of integer numbers
 170  */
 171 //=============================================================================
 172 struct Link
 173 {
 174   int n1, n2;
 175   Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
 176   Link() : n1(0), n2(0) {}
 177 };
 178
 179 int HashCode(const Link& aLink, int aLimit)
 180 {
 181   return HashCode(aLink.n1 + aLink.n2, aLimit);
 182 }
 183
 184 Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
 185 {
 186   return (aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ||
 187           aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1);
 188 }
 189
 190 //================================================================================
 191 /*!
 192  * \brief Initialize netgen::OCCGeometry with OCCT shape
 193  */
 194 //================================================================================
 195
 196 void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry&     occgeo,
 197                                              const TopoDS_Shape&      shape,
 198                                              SMESH_Mesh&              mesh,
 199                                              list< SMESH_subMesh* > * meshedSM)
 200 {
 201   BRepTools::Clean (shape);
 202   try {
 203 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 204     OCC_CATCH_SIGNALS;
 205 #endif
 206     BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (shape, 0.01, true);
 207   } catch (Standard_Failure) {
 208   }
 209   Bnd_Box bb;
 210   BRepBndLib::Add (shape, bb);
 211   double x1,y1,z1,x2,y2,z2;
 212   bb.Get (x1,y1,z1,x2,y2,z2);
 213   MESSAGE("shape bounding box:\n" <<
 214           "(" << x1 << " " << y1 << " " << z1 << ") " <<
 215           "(" << x2 << " " << y2 << " " << z2 << ")");
 216   netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
 217   netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
 218   occgeo.boundingbox = netgen::Box<3> (p1,p2);
 219
 220   occgeo.shape = shape;
 221   occgeo.changed = 1;
 222   //occgeo.BuildFMap();
 223
 224   // fill maps of shapes of occgeo with not yet meshed subshapes
 225
 226   // get root submeshes
 227   list< SMESH_subMesh* > rootSM;
 228   if ( SMESH_subMesh* sm = mesh.GetSubMeshContaining( shape )) {
 229     rootSM.push_back( sm );
 230   }
 231   else {
 232     for ( TopoDS_Iterator it( shape ); it.More(); it.Next() )
 233       rootSM.push_back( mesh.GetSubMesh( it.Value() ));
 234   }
 235
 236   // add subshapes of empty submeshes
 237   list< SMESH_subMesh* >::iterator rootIt = rootSM.begin(), rootEnd = rootSM.end();
 238   for ( ; rootIt != rootEnd; ++rootIt ) {
 239     SMESH_subMesh * root = *rootIt;
 240     SMESH_subMeshIteratorPtr smIt = root->getDependsOnIterator(/*includeSelf=*/true,
 241                                                                /*complexShapeFirst=*/true);
 242     // to find a right orientation of subshapes (PAL20462)
 243     TopTools_IndexedMapOfShape subShapes;
 244     TopExp::MapShapes(root->GetSubShape(), subShapes);
 245     while ( smIt->more() ) {
 246       SMESH_subMesh* sm = smIt->next();
 247       if ( sm->IsEmpty() ) {
 248         TopoDS_Shape shape = sm->GetSubShape();
 249         if ( shape.ShapeType() != TopAbs_VERTEX )
 250           shape = subShapes( subShapes.FindIndex( shape ));// - shape->index->oriented shape
 251         switch ( shape.ShapeType() ) {
 252         case TopAbs_FACE  : occgeo.fmap.Add( shape ); break;
 253         case TopAbs_EDGE  : occgeo.emap.Add( shape ); break;
 254         case TopAbs_VERTEX: occgeo.vmap.Add( shape ); break;
 255         case TopAbs_SOLID :occgeo.somap.Add( shape ); break;
 256         default:;
 257         }
 258       }
 259       // collect submeshes of meshed shapes
 260       else if (meshedSM) {
 261         meshedSM->push_back( sm );
 262       }
 263     }
 264   }
 265   occgeo.facemeshstatus.SetSize (occgeo.fmap.Extent());
 266   occgeo.facemeshstatus = 0;
 267
 268 }
 269
 270 //================================================================================
 271 /*!
 272  * \brief return id of netgen point corresponding to SMDS node
 273  */
 274 //================================================================================
 275 typedef map< const SMDS_MeshNode*, int > TNode2IdMap;
 276
 277 static int ngNodeId( const SMDS_MeshNode* node,
 278                      netgen::Mesh&        ngMesh,
 279                      TNode2IdMap&         nodeNgIdMap)
 280 {
 281   int newNgId = ngMesh.GetNP() + 1;
 282
 283   pair< TNode2IdMap::iterator, bool > it_isNew = nodeNgIdMap.insert( make_pair( node, newNgId ));
 284
 285   if ( it_isNew.second ) {
 286     netgen::MeshPoint p( netgen::Point<3> (node->X(), node->Y(), node->Z()) );
 287     ngMesh.AddPoint( p );
 288   }
 289   return it_isNew.first->second;
 290 }
 291
 292 //================================================================================
 293 /*!
 294  * \brief fill ngMesh with nodes and elements of computed submeshes
 295  */
 296 //================================================================================
 297
 298 bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry&           occgeom,
 299                                      netgen::Mesh&                  ngMesh,
 300                                      vector<SMDS_MeshNode*>&        nodeVec,
 301                                      const list< SMESH_subMesh* > & meshedSM)
 302 {
 303   TNode2IdMap nodeNgIdMap;
 304
 305   TopTools_MapOfShape visitedShapes;
 306
 307   SMESH_MesherHelper helper (*_mesh);
 308
 309   int faceID = occgeom.fmap.Extent();
 310   _faceDescriptors.clear();
 311
 312   list< SMESH_subMesh* >::const_iterator smIt, smEnd = meshedSM.end();
 313   for ( smIt = meshedSM.begin(); smIt != smEnd; ++smIt )
 314   {
 315     SMESH_subMesh* sm = *smIt;
 316     if ( !visitedShapes.Add( sm->GetSubShape() ))
 317       continue;
 318
 319     SMESHDS_SubMesh * smDS = sm->GetSubMeshDS();
 320
 321     switch ( sm->GetSubShape().ShapeType() )
 322     {
 323     case TopAbs_EDGE: { // EDGE
 324       // ----------------------
 325       const TopoDS_Edge& geomEdge  = TopoDS::Edge( sm->GetSubShape() );
 326
 327       // Add ng segments for each not meshed face the edge bounds
 328       TopTools_MapOfShape visitedAncestors;
 329       const TopTools_ListOfShape& ancestors = _mesh->GetAncestors( geomEdge );
 330       TopTools_ListIteratorOfListOfShape ancestorIt ( ancestors );
 331       for ( ; ancestorIt.More(); ancestorIt.Next() )
 332       {
 333         const TopoDS_Shape & ans = ancestorIt.Value();
 334         if ( ans.ShapeType() != TopAbs_FACE || !visitedAncestors.Add( ans ))
 335           continue;
 336         const TopoDS_Face& face = TopoDS::Face( ans );
 337
 338         int faceID = occgeom.fmap.FindIndex( face );
 339         if ( faceID < 1 )
 340           continue; // meshed face
 341
 342         // find out orientation of geomEdge within face
 343         bool isForwad = false;
 344         for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() ) {
 345           if ( geomEdge.IsSame( exp.Current() )) {
 346             isForwad = ( exp.Current().Orientation() == geomEdge.Orientation() );
 347             break;
 348           }
 349         }
 350         bool isQuad = smDS->GetElements()->next()->IsQuadratic();
 351
 352         // get all nodes from geomEdge
 353         StdMeshers_FaceSide fSide( face, geomEdge, _mesh, isForwad, isQuad );
 354         const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
 355         int i, nbSeg = fSide.NbSegments();
 356
 357         double otherSeamParam = 0;
 358         helper.SetSubShape( face );
 359         bool isSeam = helper.IsRealSeam( geomEdge );
 360         if ( isSeam )
 361           otherSeamParam =
 362             helper.GetOtherParam( helper.GetPeriodicIndex() == 1 ? points[0].u : points[0].v );
 363
 364         // add segments
 365
 366         int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );
 367
 368         for ( i = 0; i < nbSeg; ++i )
 369         {
 370           const UVPtStruct& p1 = points[ i ];
 371           const UVPtStruct& p2 = points[ i+1 ];
 372
 373           netgen::Segment seg;
 374           // ng node ids
 375           seg.p1 = prevNgId;
 376           seg.p2 = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
 377           // node param on curve
 378           seg.epgeominfo[ 0 ].dist = p1.param;
 379           seg.epgeominfo[ 1 ].dist = p2.param;
 380           // uv on face
 381           seg.epgeominfo[ 0 ].u = p1.u;
 382           seg.epgeominfo[ 0 ].v = p1.v;
 383           seg.epgeominfo[ 1 ].u = p2.u;
 384           seg.epgeominfo[ 1 ].v = p2.v;
 385
 386           //seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
 387           seg.si = faceID;                   // = geom.fmap.FindIndex (face);
 388           seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
 389           ngMesh.AddSegment (seg);
 390
 391           if ( isSeam )
 392           {
 393             if ( helper.GetPeriodicIndex() == 1 ) {
 394               seg.epgeominfo[ 0 ].u = otherSeamParam;
 395               seg.epgeominfo[ 1 ].u = otherSeamParam;
 396               swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
 397             } else {
 398               seg.epgeominfo[ 0 ].v = otherSeamParam;
 399               seg.epgeominfo[ 1 ].v = otherSeamParam;
 400               swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
 401             }
 402             swap (seg.p1, seg.p2);
 403             swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
 404             seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
 405             ngMesh.AddSegment (seg);
 406           }
 407         }
 408       } // loop on geomEdge ancestors
 409
 410       break;
 411     } // case TopAbs_EDGE
 412
 413     case TopAbs_FACE: { // FACE
 414       // ----------------------
 415       const TopoDS_Face& geomFace  = TopoDS::Face( sm->GetSubShape() );
 416       helper.SetSubShape( geomFace );
 417
 418       // Find solids the geomFace bounds
 419       int solidID1 = 0, solidID2 = 0;
 420       const TopTools_ListOfShape& ancestors = _mesh->GetAncestors( geomFace );
 421       TopTools_ListIteratorOfListOfShape ancestorIt ( ancestors );
 422       for ( ; ancestorIt.More(); ancestorIt.Next() )
 423       {
 424         const TopoDS_Shape & solid = ancestorIt.Value();
 425         if ( solid.ShapeType() == TopAbs_SOLID  ) {
 426           int id = occgeom.somap.FindIndex ( solid );
 427           if ( solidID1 && id != solidID1 ) solidID2 = id;
 428           else                              solidID1 = id;
 429         }
 430       }
 431       faceID++;
 432       _faceDescriptors[ faceID ].first  = solidID1;
 433       _faceDescriptors[ faceID ].second = solidID2;
 434
 435       // Orient the face correctly in solidID1 (issue 0020206)
 436       bool reverse = false;
 437       if ( solidID1 ) {
 438         TopoDS_Shape solid = occgeom.somap( solidID1 );
 439         for ( TopExp_Explorer f( solid, TopAbs_FACE ); f.More(); f.Next() ) {
 440           if ( geomFace.IsSame( f.Current() )) {
 441             reverse = SMESH_Algo::IsReversedSubMesh( TopoDS::Face( f.Current()), helper.GetMeshDS() );
 442             break;
 443           }
 444         }
 445       }
 446
 447       // Add surface elements
 448       SMDS_ElemIteratorPtr faces = smDS->GetElements();
 449       while ( faces->more() ) {
 450
 451         const SMDS_MeshElement* f = faces->next();
 452         if ( f->NbNodes() % 3 != 0 ) { // not triangle
 453           for ( ancestorIt.Initialize(ancestors); ancestorIt.More(); ancestorIt.Next() )
 454             if ( ancestorIt.Value().ShapeType() == TopAbs_SOLID  ) {
 455               sm = _mesh->GetSubMesh( ancestorIt.Value() );
 456               break;
 457             }
 458           SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 459           smError.reset( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"Not triangle submesh"));
 460           smError->myBadElements.push_back( f );
 461           return false;
 462         }
 463
 464         netgen::Element2d tri(3);
 465         tri.SetIndex ( faceID );
 466
 467         for ( int i = 0; i < 3; ++i ) {
 468           const SMDS_MeshNode* node = f->GetNode( i ), * inFaceNode=0;
 469           if ( helper.IsSeamShape( node->GetPosition()->GetShapeId() ))
 470             if ( helper.IsSeamShape( f->GetNodeWrap( i+1 )->GetPosition()->GetShapeId() ))
 471               inFaceNode = f->GetNodeWrap( i-1 );
 472             else
 473               inFaceNode = f->GetNodeWrap( i+1 );
 474
 475           gp_XY uv = helper.GetNodeUV( geomFace, node, inFaceNode );
 476           if ( reverse ) {
 477             tri.GeomInfoPi(3-i).u = uv.X();
 478             tri.GeomInfoPi(3-i).v = uv.Y();
 479             tri.PNum      (3-i) = ngNodeId( node, ngMesh, nodeNgIdMap );
 480           } else {
 481             tri.GeomInfoPi(i+1).u = uv.X();
 482             tri.GeomInfoPi(i+1).v = uv.Y();
 483             tri.PNum      (i+1) = ngNodeId( node, ngMesh, nodeNgIdMap );
 484           }
 485         }
 486
 487         ngMesh.AddSurfaceElement (tri);
 488
 489       }
 490       break;
 491     } //
 492
 493     case TopAbs_VERTEX: { // VERTEX
 494       // --------------------------
 495       SMDS_NodeIteratorPtr nodeIt = smDS->GetNodes();
 496       if ( nodeIt->more() )
 497         ngNodeId( nodeIt->next(), ngMesh, nodeNgIdMap );
 498       break;
 499     }
 500     default:;
 501     } // switch
 502   } // loop on submeshes
 503
 504   // fill nodeVec
 505   nodeVec.resize( ngMesh.GetNP() + 1 );
 506   TNode2IdMap::iterator node_NgId, nodeNgIdEnd = nodeNgIdMap.end();
 507   for ( node_NgId = nodeNgIdMap.begin(); node_NgId != nodeNgIdEnd; ++node_NgId)
 508     nodeVec[ node_NgId->second ] = (SMDS_MeshNode*) node_NgId->first;
 509
 510   return true;
 511 }
 512
 513 //=============================================================================
 514 /*!
 515  * Here we are going to use the NETGEN mesher
 516  */
 517 //=============================================================================
 518 bool NETGENPlugin_Mesher::Compute()
 519 {
 520 #ifdef WNT
 521   netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
 522 #else
 523   netgen::MeshingParameters& mparams = netgen::mparam;
 524 #endif
 525   MESSAGE("Compute with:\n"
 526           " max size = " << mparams.maxh << "\n"
 527           " segments per edge = " << mparams.segmentsperedge);
 528   MESSAGE("\n"
 529           " growth rate = " << mparams.grading << "\n"
 530           " elements per radius = " << mparams.curvaturesafety << "\n"
 531           " second order = " << mparams.secondorder << "\n"
 532           " quad allowed = " << mparams.quad);
 533
 534   SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
 535   nglib::Ng_Init();
 536
 537   // -------------------------
 538   // Prepare OCC geometry
 539   // -------------------------
 540
 541   netgen::OCCGeometry occgeo;
 542   list< SMESH_subMesh* > meshedSM;
 543   PrepareOCCgeometry( occgeo, _shape, *_mesh, &meshedSM );
 544
 545   // -------------------------
 546   // Generate the mesh
 547   // -------------------------
 548
 549   netgen::Mesh *ngMesh = NULL;
 550
 551   SMESH_Comment comment;
 552   int err = 0;
 553   int nbInitNod = 0;
 554   int nbInitSeg = 0;
 555   int nbInitFac = 0;
 556   // vector of nodes in which node index == netgen ID
 557   vector< SMDS_MeshNode* > nodeVec;
 558   try
 559   {
 560     // ----------------
 561     // compute 1D mesh
 562     // ----------------
 563     // pass 1D simple parameters to NETGEN
 564     if ( _simpleHyp ) {
 565       if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) {
 566         // nb of segments
 567         mparams.segmentsperedge = nbSeg + 0.1;
 568         mparams.maxh = occgeo.boundingbox.Diam();
 569         mparams.grading = 0.01;
 570       }
 571       else {
 572         // segment length
 573         mparams.segmentsperedge = 1;
 574         mparams.maxh = _simpleHyp->GetLocalLength();
 575       }
 576     }
 577     // let netgen create ngMesh and calculate element size on not meshed shapes
 578     char *optstr = 0;
 579     int startWith = netgen::MESHCONST_ANALYSE;
 580     int endWith   = netgen::MESHCONST_ANALYSE;
 581     err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
 582     if (err) comment << "Error in netgen::OCCGenerateMesh() at MESHCONST_ANALYSE step";
 583
 584     // fill ngMesh with nodes and elements of computed submeshes
 585     err = ! fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM);
 586     nbInitNod = ngMesh->GetNP();
 587     nbInitSeg = ngMesh->GetNSeg();
 588     nbInitFac = ngMesh->GetNSE();
 589
 590     // compute mesh
 591     if (!err)
 592     {
 593       startWith = endWith = netgen::MESHCONST_MESHEDGES;
 594       err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
 595       if (err) comment << "Error in netgen::OCCGenerateMesh() at 1D mesh generation";
 596     }
 597     // ---------------------
 598     // compute surface mesh
 599     // ---------------------
 600     if (!err)
 601     {
 602       // pass 2D simple parameters to NETGEN
 603       if ( _simpleHyp ) {
 604         if ( double area = _simpleHyp->GetMaxElementArea() ) {
 605           // face area
 606           mparams.maxh = sqrt(2. * area/sqrt(3.0));
 607           mparams.grading = 0.4; // moderate size growth
 608         }
 609         else {
 610           // length from edges
 611           double length = 0;
 612           for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() )
 613             length += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() ));
 614           if ( ngMesh->GetNSeg() )
 615             mparams.maxh = length / ngMesh->GetNSeg();
 616           else
 617             mparams.maxh = 1000;
 618           mparams.grading = 0.2; // slow size growth
 619         }
 620         mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
 621         ngMesh->SetGlobalH (mparams.maxh);
 622         netgen::Box<3> bb = occgeo.GetBoundingBox();
 623         bb.Increase (bb.Diam()/20);
 624         ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading);
 625       }
 626       // let netgen compute 2D mesh
 627       startWith = netgen::MESHCONST_MESHSURFACE;
 628       endWith = _optimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
 629       err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
 630       if (err) comment << "Error in netgen::OCCGenerateMesh() at surface mesh generation";
 631     }
 632     // ---------------------
 633     // generate volume mesh
 634     // ---------------------
 635     if (!err && _isVolume)
 636     {
 637       // add ng face descriptors of meshed faces
 638       std::map< int, std::pair<int,int> >::iterator fId_soIds = _faceDescriptors.begin();
 639       for ( ; fId_soIds != _faceDescriptors.end(); ++fId_soIds ) {
 640         int faceID   = fId_soIds->first;
 641         int solidID1 = fId_soIds->second.first;
 642         int solidID2 = fId_soIds->second.second;
 643         ngMesh->AddFaceDescriptor (netgen::FaceDescriptor(faceID, solidID1, solidID2, 0));
 644       }
 645       // pass 3D simple parameters to NETGEN
 646       const NETGENPlugin_SimpleHypothesis_3D* simple3d =
 647         dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
 648       if ( simple3d ) {
 649         if ( double vol = simple3d->GetMaxElementVolume() ) {
 650           // max volume
 651           mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
 652           mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
 653         }
 654         else {
 655           // length from faces
 656           mparams.maxh = ngMesh->AverageH();
 657         }
 658 //      netgen::ARRAY<double> maxhdom;
 659 //      maxhdom.SetSize (occgeo.NrSolids());
 660 //      maxhdom = mparams.maxh;
 661 //      ngMesh->SetMaxHDomain (maxhdom);
 662         ngMesh->SetGlobalH (mparams.maxh);
 663         mparams.grading = 0.4;
 664         ngMesh->CalcLocalH();
 665       }
 666       // let netgen compute 3D mesh
 667       startWith = netgen::MESHCONST_MESHVOLUME;
 668       endWith = _optimize ? netgen::MESHCONST_OPTVOLUME : netgen::MESHCONST_MESHVOLUME;
 669       err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
 670       if (err) comment << "Error in netgen::OCCGenerateMesh()";
 671     }
 672     if (!err && mparams.secondorder > 0)
 673     {
 674       netgen::OCCRefinementSurfaces ref (occgeo);
 675       ref.MakeSecondOrder (*ngMesh);
 676     }
 677   }
 678   catch (netgen::NgException exc)
 679   {
 680     error->myName = err = COMPERR_ALGO_FAILED;
 681     comment << exc.What();
 682   }
 683
 684   int nbNod = ngMesh->GetNP();
 685   int nbSeg = ngMesh->GetNSeg();
 686   int nbFac = ngMesh->GetNSE();
 687   int nbVol = ngMesh->GetNE();
 688
 689   MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
 690           ", nb nodes: " << nbNod <<
 691           ", nb segments: " << nbSeg <<
 692           ", nb faces: " << nbFac <<
 693           ", nb volumes: " << nbVol);
 694
 695   // -----------------------------------------------------------
 696   // Feed back the SMESHDS with the generated Nodes and Elements
 697   // -----------------------------------------------------------
 698
 699   SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
 700   bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
 701   if ( true /*isOK*/ ) // get whatever built
 702   {
 703     // map of nodes assigned to submeshes
 704     NCollection_Map<int> pindMap;
 705     // create and insert nodes into nodeVec
 706     nodeVec.resize( nbNod + 1 );
 707     int i;
 708     for (i = nbInitNod+1; i <= nbNod /*&& isOK*/; ++i )
 709     {
 710       const netgen::MeshPoint& ngPoint = ngMesh->Point(i);
 711       SMDS_MeshNode* node = NULL;
 712       bool newNodeOnVertex = false;
 713       TopoDS_Vertex aVert;
 714       if (i-nbInitNod <= occgeo.vmap.Extent())
 715       {
 716         // point on vertex
 717         aVert = TopoDS::Vertex(occgeo.vmap(i-nbInitNod));
 718         SMESHDS_SubMesh * submesh = meshDS->MeshElements(aVert);
 719         if (submesh)
 720         {
 721           SMDS_NodeIteratorPtr it = submesh->GetNodes();
 722           if (it->more())
 723           {
 724             node = const_cast<SMDS_MeshNode*> (it->next());
 725             pindMap.Add(i);
 726           }
 727         }
 728         if (!node)
 729           newNodeOnVertex = true;
 730       }
 731       if (!node)
 732         node = meshDS->AddNode(ngPoint.X(), ngPoint.Y(), ngPoint.Z());
 733       if (!node)
 734       {
 735         MESSAGE("Cannot create a mesh node");
 736         if ( !comment.size() ) comment << "Cannot create a mesh node";
 737         nbSeg = nbFac = nbVol = isOK = 0;
 738         break;
 739       }
 740       nodeVec.at(i) = node;
 741       if (newNodeOnVertex)
 742       {
 743         // point on vertex
 744         meshDS->SetNodeOnVertex(node, aVert);
 745         pindMap.Add(i);
 746       }
 747     }
 748
 749     // create mesh segments along geometric edges
 750     NCollection_Map<Link> linkMap;
 751     for (i = nbInitSeg+1; i <= nbSeg/* && isOK*/; ++i )
 752     {
 753       const netgen::Segment& seg = ngMesh->LineSegment(i);
 754       Link link(seg.p1, seg.p2);
 755       if (linkMap.Contains(link))
 756         continue;
 757       linkMap.Add(link);
 758       TopoDS_Edge aEdge;
 759       int pinds[3] = { seg.p1, seg.p2, seg.pmid };
 760       int nbp = 0;
 761       double param2 = 0;
 762       for (int j=0; j < 3; ++j)
 763       {
 764         int pind = pinds[j];
 765         if (pind <= 0) continue;
 766         ++nbp;
 767         double param;
 768         if (j < 2)
 769         {
 770           if (aEdge.IsNull())
 771           {
 772             int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
 773             if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
 774               aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
 775           }
 776           param = seg.epgeominfo[j].dist;
 777           param2 += param;
 778         }
 779         else
 780           param = param2 * 0.5;
 781         if (pind <= nbInitNod || pindMap.Contains(pind))
 782           continue;
 783         if (!aEdge.IsNull())
 784         {
 785           meshDS->SetNodeOnEdge(nodeVec.at(pind), aEdge, param);
 786           pindMap.Add(pind);
 787         }
 788       }
 789       SMDS_MeshEdge* edge;
 790       if (nbp < 3) // second order ?
 791         edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]));
 792       else
 793         edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]),
 794                                 nodeVec.at(pinds[2]));
 795       if (!edge)
 796       {
 797         if ( !comment.size() ) comment << "Cannot create a mesh edge";
 798         MESSAGE("Cannot create a mesh edge");
 799         nbSeg = nbFac = nbVol = isOK = 0;
 800         break;
 801       }
 802       if (!aEdge.IsNull())
 803         meshDS->SetMeshElementOnShape(edge, aEdge);
 804     }
 805
 806     // create mesh faces along geometric faces
 807     for (i = nbInitFac+1; i <= nbFac/* && isOK*/; ++i )
 808     {
 809       const netgen::Element2d& elem = ngMesh->SurfaceElement(i);
 810       int aGeomFaceInd = elem.GetIndex();
 811       TopoDS_Face aFace;
 812       if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
 813         aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
 814       vector<SMDS_MeshNode*> nodes;
 815       for (int j=1; j <= elem.GetNP(); ++j)
 816       {
 817         int pind = elem.PNum(j);
 818         SMDS_MeshNode* node = nodeVec.at(pind);
 819         nodes.push_back(node);
 820         if (pind <= nbInitNod || pindMap.Contains(pind))
 821           continue;
 822         if (!aFace.IsNull())
 823         {
 824           const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
 825           meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
 826           pindMap.Add(pind);
 827         }
 828       }
 829       SMDS_MeshFace* face = NULL;
 830       switch (elem.GetType())
 831       {
 832       case netgen::TRIG:
 833         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
 834         break;
 835       case netgen::QUAD:
 836         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
 837         break;
 838       case netgen::TRIG6:
 839         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
 840         break;
 841       case netgen::QUAD8:
 842         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
 843                                nodes[4],nodes[7],nodes[5],nodes[6]);
 844         break;
 845       default:
 846         MESSAGE("NETGEN created a face of unexpected type, ignoring");
 847         continue;
 848       }
 849       if (!face)
 850       {
 851         if ( !comment.size() ) comment << "Cannot create a mesh face";
 852         MESSAGE("Cannot create a mesh face");
 853         nbSeg = nbFac = nbVol = isOK = 0;
 854         break;
 855       }
 856       if (!aFace.IsNull())
 857         meshDS->SetMeshElementOnShape(face, aFace);
 858     }
 859
 860     // create tetrahedra
 861     for (i = 1; i <= nbVol/* && isOK*/; ++i)
 862     {
 863       const netgen::Element& elem = ngMesh->VolumeElement(i);
 864       int aSolidInd = elem.GetIndex();
 865       TopoDS_Solid aSolid;
 866       if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
 867         aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
 868       vector<SMDS_MeshNode*> nodes;
 869       for (int j=1; j <= elem.GetNP(); ++j)
 870       {
 871         int pind = elem.PNum(j);
 872         SMDS_MeshNode* node = nodeVec.at(pind);
 873         nodes.push_back(node);
 874         if (pind <= nbInitNod || pindMap.Contains(pind))
 875           continue;
 876         if (!aSolid.IsNull())
 877         {
 878           // point in solid
 879           meshDS->SetNodeInVolume(node, aSolid);
 880           pindMap.Add(pind);
 881         }
 882       }
 883       SMDS_MeshVolume* vol = NULL;
 884       switch (elem.GetType())
 885       {
 886       case netgen::TET:
 887         vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]);
 888         break;
 889       case netgen::TET10:
 890         vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3],
 891                                 nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]);
 892         break;
 893       default:
 894         MESSAGE("NETGEN created a volume of unexpected type, ignoring");
 895         continue;
 896       }
 897       if (!vol)
 898       {
 899         if ( !comment.size() ) comment << "Cannot create a mesh volume";
 900         MESSAGE("Cannot create a mesh volume");
 901         nbSeg = nbFac = nbVol = isOK = 0;
 902         break;
 903       }
 904       if (!aSolid.IsNull())
 905         meshDS->SetMeshElementOnShape(vol, aSolid);
 906     }
 907   }
 908
 909   if ( error->IsOK() && ( !isOK || comment.size() > 0 ))
 910     error->myName = COMPERR_ALGO_FAILED;
 911   if ( !comment.empty() )
 912     error->myComment = comment;
 913
 914   // set bad compute error to subshapes of all failed subshapes shapes
 915   if ( !error->IsOK() && err )
 916   {
 917     for (int i = 1; i <= occgeo.fmap.Extent(); i++) {
 918       int status = occgeo.facemeshstatus[i-1];
 919       if (status == 1 ) continue;
 920       if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.fmap( i ))) {
 921         SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 922         if ( !smError || smError->IsOK() ) {
 923           if ( status == -1 )
 924             smError.reset( new SMESH_ComputeError( error->myName, error->myComment ));
 925           else
 926             smError.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, "Ignored" ));
 927         }
 928       }
 929     }
 930   }
 931
 932   nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);
 933   nglib::Ng_Exit();
 934
 935   RemoveTmpFiles();
 936
 937   return error->IsOK();
 938 }
 939
 940 //================================================================================
 941 /*!
 942  * \brief Remove "test.out" and "problemfaces" files in current directory
 943  */
 944 //================================================================================
 945
 946 void NETGENPlugin_Mesher::RemoveTmpFiles()
 947 {
 948   TCollection_AsciiString str("test.out");
 949   OSD_Path path1( str );
 950   OSD_File file1( path1 );
 951   file1.Remove();
 952   str = "problemfaces";
 953   OSD_Path path2( str );
 954   OSD_File file2( path2 );
 955   file2.Remove();
 956 }