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