src/NETGENPlugin/NETGENPlugin_NETGEN_3D.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 //=============================================================================
  23 // File      : NETGENPlugin_NETGEN_3D.cxx
  24 //             Moved here from SMESH_NETGEN_3D.cxx
  25 // Created   : lundi 27 Janvier 2003
  26 // Author    : Nadir BOUHAMOU (CEA)
  27 // Project   : SALOME
  28 //=============================================================================
  29 //
  30 #include "NETGENPlugin_NETGEN_3D.hxx"
  31
  32 #include "NETGENPlugin_Mesher.hxx"
  33
  34 #include "SMDS_MeshElement.hxx"
  35 #include "SMDS_MeshNode.hxx"
  36 #include "SMESHDS_Mesh.hxx"
  37 #include "SMESH_Comment.hxx"
  38 #include "SMESH_ControlsDef.hxx"
  39 #include "SMESH_Gen.hxx"
  40 #include "SMESH_Mesh.hxx"
  41 #include "SMESH_MesherHelper.hxx"
  42 #include "SMESH_MeshEditor.hxx"
  43 #include "StdMeshers_QuadToTriaAdaptor.hxx"
  44
  45 #include <BRepGProp.hxx>
  46 #include <BRep_Tool.hxx>
  47 #include <GProp_GProps.hxx>
  48 #include <TopExp.hxx>
  49 #include <TopExp_Explorer.hxx>
  50 #include <TopTools_ListIteratorOfListOfShape.hxx>
  51 #include <TopoDS.hxx>
  52
  53 #include <Standard_Failure.hxx>
  54 #include <Standard_ErrorHandler.hxx>
  55
  56 #include "utilities.h"
  57
  58 #include <list>
  59 #include <vector>
  60 #include <map>
  61
  62 /*
  63   Netgen include files
  64 */
  65
  66 namespace nglib {
  67 #include <nglib.h>
  68 }
  69 using namespace nglib;
  70 using namespace std;
  71
  72 //=============================================================================
  73 /*!
  74  *
  75  */
  76 //=============================================================================
  77
  78 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
  79                              SMESH_Gen* gen)
  80   : SMESH_3D_Algo(hypId, studyId, gen)
  81 {
  82   MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
  83   _name = "NETGEN_3D";
  84   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  85   _compatibleHypothesis.push_back("MaxElementVolume");
  86
  87   _maxElementVolume = 0.;
  88
  89   _hypMaxElementVolume = NULL;
  90
  91   _requireShape = false; // can work without shape
  92 }
  93
  94 //=============================================================================
  95 /*!
  96  *
  97  */
  98 //=============================================================================
  99
 100 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
 101 {
 102   MESSAGE("NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D");
 103 }
 104
 105 //=============================================================================
 106 /*!
 107  *
 108  */
 109 //=============================================================================
 110
 111 bool NETGENPlugin_NETGEN_3D::CheckHypothesis
 112                          (SMESH_Mesh& aMesh,
 113                           const TopoDS_Shape& aShape,
 114                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 115 {
 116   MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
 117
 118   _hypMaxElementVolume = NULL;
 119   _maxElementVolume = DBL_MAX;
 120
 121   list<const SMESHDS_Hypothesis*>::const_iterator itl;
 122   const SMESHDS_Hypothesis* theHyp;
 123
 124   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
 125   int nbHyp = hyps.size();
 126   if (!nbHyp)
 127   {
 128     aStatus = SMESH_Hypothesis::HYP_OK;
 129     //aStatus = SMESH_Hypothesis::HYP_MISSING;
 130     return true;  // can work with no hypothesis
 131   }
 132
 133   itl = hyps.begin();
 134   theHyp = (*itl); // use only the first hypothesis
 135
 136   string hypName = theHyp->GetName();
 137
 138   bool isOk = false;
 139
 140   if (hypName == "MaxElementVolume")
 141   {
 142     _hypMaxElementVolume = static_cast<const StdMeshers_MaxElementVolume*> (theHyp);
 143     ASSERT(_hypMaxElementVolume);
 144     _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
 145     isOk =true;
 146     aStatus = SMESH_Hypothesis::HYP_OK;
 147   }
 148   else
 149     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 150
 151   return isOk;
 152 }
 153
 154 namespace
 155 {
 156   //================================================================================
 157   /*!
 158    * \brief It correctly initializes netgen library at constructor and
 159    *        correctly finishes using netgen library at destructor
 160    */
 161   //================================================================================
 162
 163   struct TNetgenLibWrapper
 164   {
 165     Ng_Mesh* _ngMesh;
 166     TNetgenLibWrapper()
 167     {
 168       Ng_Init();
 169       _ngMesh = Ng_NewMesh();
 170     }
 171     ~TNetgenLibWrapper()
 172     {
 173       Ng_DeleteMesh( _ngMesh );
 174       Ng_Exit();
 175       NETGENPlugin_Mesher::RemoveTmpFiles();
 176     }
 177   };
 178
 179   //================================================================================
 180   /*!
 181    * \brief Find mesh faces on non-internal geom faces sharing internal edge
 182    *        some nodes of which are to be doubled to make the second border of the "crack"
 183    */
 184   //================================================================================
 185
 186   void findBorders( const set<int>&     internalShapeIds,
 187                     SMESH_MesherHelper& helper,
 188                     TIDSortedElemSet &  borderElems,
 189                     set<int> &          borderFaceIds )
 190   {
 191     SMESH_Mesh*     mesh = helper.GetMesh();
 192     SMESHDS_Mesh* meshDS = helper.GetMeshDS();
 193
 194     // loop on internal geom edges
 195     set<int>::const_iterator intShapeId = internalShapeIds.begin();
 196     for ( ; intShapeId != internalShapeIds.end(); ++intShapeId )
 197     {
 198       const TopoDS_Shape& s = meshDS->IndexToShape( *intShapeId );
 199       if ( s.ShapeType() != TopAbs_EDGE ) continue;
 200
 201       // get internal and non-internal geom faces sharing the internal edge <s>
 202       int intFace = 0;
 203       set<int>::iterator bordFace = borderFaceIds.end();
 204       TopTools_ListIteratorOfListOfShape ancIt( mesh->GetAncestors( s ));
 205       for ( ; ancIt.More(); ancIt.Next() )
 206       {
 207         if ( ancIt.Value().ShapeType() != TopAbs_FACE ) continue;
 208         int faceID = meshDS->ShapeToIndex( ancIt.Value() );
 209         if ( internalShapeIds.count( faceID ))
 210           intFace = faceID;
 211         else
 212           bordFace = borderFaceIds.insert( faceID ).first;
 213       }
 214       if ( bordFace == borderFaceIds.end() || !intFace ) continue;
 215
 216       // get all links of mesh faces on internal geom face sharing nodes on edge <s>
 217       set< SMESH_OrientedLink > links; //!< links of faces on internal geom face
 218       TIDSortedElemSet suspectFaces;   //!< mesh faces on border geom faces
 219       SMESHDS_SubMesh* intFaceSM = meshDS->MeshElements( intFace );
 220       if ( !intFaceSM || intFaceSM->NbElements() == 0 ) continue;
 221       SMESH_subMeshIteratorPtr smIt = mesh->GetSubMesh( s )->getDependsOnIterator(true,true);
 222       while ( smIt->more() )
 223       {
 224         SMESHDS_SubMesh* sm = smIt->next()->GetSubMeshDS();
 225         if ( !sm ) continue;
 226         SMDS_NodeIteratorPtr nIt = sm->GetNodes();
 227         while ( nIt->more() )
 228         {
 229           const SMDS_MeshNode* nOnEdge = nIt->next();
 230           SMDS_ElemIteratorPtr fIt = nOnEdge->GetInverseElementIterator(SMDSAbs_Face);
 231           while ( fIt->more() )
 232           {
 233             const SMDS_MeshElement* f = fIt->next();
 234             if ( intFaceSM->Contains( f ))
 235             {
 236               int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
 237               for ( int i = 0; i < nbNodes; ++i )
 238                 links.insert( SMESH_OrientedLink( f->GetNode(i), f->GetNode((i+1)%nbNodes)));
 239             }
 240             else
 241             {
 242               suspectFaces.insert( f );
 243             }
 244           }
 245         }
 246       }
 247       // <suspectFaces> having link with same orientation as mesh faces on
 248       // the internal geom face are <borderElems>.
 249       // Some of them have only one node on edge s, we collect them to decide
 250       // later by links of found <borderElems>
 251       TIDSortedElemSet posponedFaces;
 252       set< SMESH_OrientedLink > borderLinks;
 253       TIDSortedElemSet::iterator fIt = suspectFaces.begin();
 254       for ( ; fIt != suspectFaces.end(); ++fIt )
 255       {
 256         const SMDS_MeshElement* f = *fIt;
 257         bool linkFound = false, isBorder = false;
 258         list< SMESH_OrientedLink > faceLinks;
 259         int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
 260         for ( int i = 0; i < nbNodes; ++i )
 261         {
 262           SMESH_OrientedLink link( f->GetNode(i), f->GetNode((i+1)%nbNodes));
 263           faceLinks.push_back( link );
 264           if ( !linkFound )
 265           {
 266             set< SMESH_OrientedLink >::iterator foundLink = links.find( link );
 267             if ( foundLink != links.end() )
 268             {
 269               linkFound= true;
 270               isBorder = ( foundLink->_reversed == link._reversed );
 271               if ( !isBorder ) break;
 272             }
 273           }
 274         }
 275         if ( isBorder )
 276         {
 277           borderElems.insert( f );
 278           borderLinks.insert( faceLinks.begin(), faceLinks.end() );
 279         }
 280         else if ( !linkFound )
 281         {
 282           posponedFaces.insert( f );
 283         }
 284       }
 285       // decide on posponedFaces
 286       for ( fIt = posponedFaces.begin(); fIt != posponedFaces.end(); ++fIt )
 287       {
 288         const SMDS_MeshElement* f = *fIt;
 289         int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
 290         for ( int i = 0; i < nbNodes; ++i )
 291         {
 292           SMESH_OrientedLink link( f->GetNode(i), f->GetNode((i+1)%nbNodes));
 293           set< SMESH_OrientedLink >::iterator foundLink = borderLinks.find( link );
 294           if ( foundLink != borderLinks.end() )
 295           {
 296             if ( foundLink->_reversed != link._reversed )
 297               borderElems.insert( f );
 298             break;
 299           }
 300         }
 301       }
 302     }
 303   }
 304 }
 305
 306 //=============================================================================
 307 /*!
 308  *Here we are going to use the NETGEN mesher
 309  */
 310 //=============================================================================
 311
 312 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
 313                                      const TopoDS_Shape& aShape)
 314 {
 315   MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
 316
 317   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 318
 319   SMESH_MesherHelper helper(aMesh);
 320   bool _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
 321
 322   int Netgen_NbOfNodes     = 0;
 323   int Netgen_param2ndOrder = 0;
 324   double Netgen_paramFine  = 1.;
 325   double Netgen_paramSize  = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
 326
 327   double Netgen_point[3];
 328   int Netgen_triangle[3];
 329   int Netgen_tetrahedron[4];
 330
 331   TNetgenLibWrapper ngLib;
 332   Ng_Mesh * Netgen_mesh = ngLib._ngMesh;
 333
 334   // maps of 1) ordinary nodes and 2) doubled nodes on internal shapes
 335   typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 336   typedef TNodeToIDMap::value_type                     TN2ID;
 337   TNodeToIDMap nodeToNetgenID[2];
 338
 339   {
 340     const int invalid_ID = -1;
 341
 342     SMESH::Controls::Area areaControl;
 343     SMESH::Controls::TSequenceOfXYZ nodesCoords;
 344
 345     // --------------------------------------------------------------------
 346     // Issue 0020676 (StudyFiss_bugNetgen3D.hdf). Pb with internal face.
 347     // Find internal geom faces, edges and vertices.
 348     // Nodes and faces built on the found internal shapes
 349     // will be doubled in Netgen input to make two borders of the "crack".
 350     // --------------------------------------------------------------------
 351
 352     set<int> internalShapeIds;
 353     set<int> borderFaceIds; //!< non-internal geom faces sharing internal edge
 354
 355     // mesh faces on <borderFaceIds>, having nodes on internal edge that are
 356     // to be replaced by doubled nodes
 357     TIDSortedElemSet borderElems;
 358
 359     // find "internal" faces and edges
 360     TopExp_Explorer exFa(aShape,TopAbs_FACE), exEd, exVe;
 361     for ( ; exFa.More(); exFa.Next())
 362     {
 363       if ( exFa.Current().Orientation() == TopAbs_INTERNAL )
 364       {
 365         internalShapeIds.insert( meshDS->ShapeToIndex( exFa.Current() ));
 366         for ( exEd.Init( exFa.Current(), TopAbs_EDGE ); exEd.More(); exEd.Next())
 367           if ( helper.NbAncestors( exEd.Current(), aMesh, TopAbs_FACE ) > 1 )
 368             internalShapeIds.insert( meshDS->ShapeToIndex( exEd.Current() ));
 369       }
 370     }
 371     if ( !internalShapeIds.empty() )
 372     {
 373       // find internal vertices,
 374       // we consider vertex internal if it is shared by more than one internal edge
 375       TopTools_ListIteratorOfListOfShape ancIt;
 376       set<int>::iterator intShapeId = internalShapeIds.begin();
 377       for ( ; intShapeId != internalShapeIds.end(); ++intShapeId )
 378       {
 379         const TopoDS_Shape& s = meshDS->IndexToShape( *intShapeId );
 380         if ( s.ShapeType() != TopAbs_EDGE ) continue;
 381         for ( exVe.Init( s, TopAbs_VERTEX ); exVe.More(); exVe.Next())
 382         {
 383           set<int> internalEdges;
 384           for ( ancIt.Initialize( aMesh.GetAncestors( exVe.Current() ));
 385                 ancIt.More(); ancIt.Next() )
 386           {
 387             if ( ancIt.Value().ShapeType() != TopAbs_EDGE ) continue;
 388             int edgeID = meshDS->ShapeToIndex( ancIt.Value() );
 389             if ( internalShapeIds.count( edgeID ))
 390               internalEdges.insert( edgeID );
 391           }
 392           if ( internalEdges.size() > 1 )
 393             internalShapeIds.insert( meshDS->ShapeToIndex( exVe.Current() ));
 394         }
 395       }
 396       // find border shapes and mesh elements
 397       findBorders( internalShapeIds, helper, borderElems, borderFaceIds );
 398     }
 399
 400     // ---------------------------------
 401     // Feed the Netgen with surface mesh
 402     // ---------------------------------
 403
 404     TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
 405     bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
 406
 407     StdMeshers_QuadToTriaAdaptor Adaptor;
 408     if ( aMesh.NbQuadrangles() > 0 )
 409       Adaptor.Compute(aMesh,aShape);
 410
 411     for ( exFa.ReInit(); exFa.More(); exFa.Next())
 412     {
 413       const TopoDS_Shape& aShapeFace = exFa.Current();
 414       int faceID = meshDS->ShapeToIndex( aShapeFace );
 415       bool isInternalFace = internalShapeIds.count( faceID );
 416       bool isBorderFace   = borderFaceIds.count( faceID );
 417       bool isRev = false;
 418       if ( checkReverse && !isInternalFace &&
 419            helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
 420         // IsReversedSubMesh() can work wrong on strongly curved faces,
 421         // so we use it as less as possible
 422         isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
 423
 424       const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
 425       if ( !aSubMeshDSFace ) continue;
 426       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 427       while ( iteratorElem->more() ) // loop on elements on a geom face
 428       {
 429         // check mesh face
 430         const SMDS_MeshElement* elem = iteratorElem->next();
 431         if ( !elem )
 432           return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 433         vector< const SMDS_MeshElement* > trias;
 434         bool isTraingle = ( elem->NbNodes() == ( elem->IsQuadratic() ? 6 : 3 ));
 435         if ( !isTraingle )
 436         {
 437           // use adaptor to convert quadrangle face into triangles
 438           const list<const SMDS_FaceOfNodes*>* faces = Adaptor.GetTriangles(elem);
 439           if(faces==0)
 440             return error( COMPERR_BAD_INPUT_MESH,
 441                           SMESH_Comment("No triangles in adaptor for element ")<<elem->GetID());
 442           trias.assign( faces->begin(), faces->end() );
 443         }
 444         else
 445         {
 446           trias.push_back( elem );
 447         }
 448         // Add nodes of triangles and triangles them-selves to netgen mesh
 449
 450         // a triangle on internal face is added twice,
 451         // on border face, once but with doubled nodes
 452         bool isBorder = ( isBorderFace && borderElems.count( elem ));
 453         int nbDblLoops = ( isInternalFace && isTraingle || isBorder ) ? 2 : 1;
 454
 455         for ( int i = 0; i < trias.size(); ++i )
 456         {
 457           bool reverse = isRev;
 458           for ( int isDblF = isBorder; isDblF < nbDblLoops; ++isDblF, reverse = !reverse )
 459           {
 460             // add three nodes of triangle
 461             bool hasDegen = false;
 462             for ( int iN = 0; iN < 3; ++iN )
 463             {
 464               const SMDS_MeshNode* node = trias[i]->GetNode( iN );
 465               int shapeID = node->GetPosition()->GetShapeId();
 466               if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
 467                    helper.IsDegenShape( shapeID ))
 468               {
 469                 // ignore all nodes on degeneraged edge and use node on its vertex instead
 470                 TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
 471                 node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
 472                 hasDegen = true;
 473               }
 474               bool isDblN = isDblF && internalShapeIds.count( shapeID );
 475               int& ngID = nodeToNetgenID[isDblN].insert(TN2ID( node, invalid_ID )).first->second;
 476               if ( ngID == invalid_ID )
 477               {
 478                 ngID = ++Netgen_NbOfNodes;
 479                 Netgen_point [ 0 ] = node->X();
 480                 Netgen_point [ 1 ] = node->Y();
 481                 Netgen_point [ 2 ] = node->Z();
 482                 Ng_AddPoint(Netgen_mesh, Netgen_point);
 483               }
 484               Netgen_triangle[ iN ] = ngID;
 485             }
 486             // add triangle
 487             if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
 488                               Netgen_triangle[0] == Netgen_triangle[2] ||
 489                               Netgen_triangle[2] == Netgen_triangle[1] ))
 490               break;
 491             if ( reverse )
 492               swap( Netgen_triangle[1], Netgen_triangle[2] );
 493
 494             Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 495           }
 496         }
 497 #ifdef _DEBUG_
 498         // check if a trainge is degenerated
 499         areaControl.GetPoints( elem, nodesCoords );
 500         double area = areaControl.GetValue( nodesCoords );
 501         if ( area <= DBL_MIN ) {
 502           MESSAGE( "Warning: Degenerated " << elem );
 503         }
 504 #endif
 505       } // loop on elements on a face
 506     } // loop on faces of a SOLID or SHELL
 507
 508   }
 509
 510   // -------------------------
 511   // Generate the volume mesh
 512   // -------------------------
 513
 514   Ng_Meshing_Parameters Netgen_param;
 515
 516   Netgen_param.secondorder = Netgen_param2ndOrder;
 517   Netgen_param.fineness    = Netgen_paramFine;
 518   Netgen_param.maxh        = Netgen_paramSize;
 519
 520   Ng_Result status;
 521
 522   try {
 523 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 524     OCC_CATCH_SIGNALS;
 525 #endif
 526     status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
 527   }
 528   catch (Standard_Failure& exc) {
 529     error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
 530     status = NG_VOLUME_FAILURE;
 531   }
 532   catch (...) {
 533     error("Exception in Ng_GenerateVolumeMesh()");
 534     status = NG_VOLUME_FAILURE;
 535   }
 536   if ( GetComputeError()->IsOK() ) {
 537     switch ( status ) {
 538     case NG_SURFACE_INPUT_ERROR:error( status, "NG_SURFACE_INPUT_ERROR");
 539     case NG_VOLUME_FAILURE:     error( status, "NG_VOLUME_FAILURE");
 540     case NG_STL_INPUT_ERROR:    error( status, "NG_STL_INPUT_ERROR");
 541     case NG_SURFACE_FAILURE:    error( status, "NG_SURFACE_FAILURE");
 542     case NG_FILE_NOT_FOUND:     error( status, "NG_FILE_NOT_FOUND");
 543     };
 544   }
 545
 546   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 547
 548   int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
 549
 550   MESSAGE("End of Volume Mesh Generation. status=" << status <<
 551           ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
 552           ", nb tetra: " << Netgen_NbOfTetra);
 553
 554   // -------------------------------------------------------------------
 555   // Feed back the SMESHDS with the generated Nodes and Volume Elements
 556   // -------------------------------------------------------------------
 557
 558   // vector of nodes in which node index == netgen ID
 559   vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
 560   // insert old nodes into nodeVec
 561   for ( int isDbl = 0; isDbl < 2; ++isDbl )
 562   {
 563     TNodeToIDMap::iterator n_id = nodeToNetgenID[isDbl].begin();
 564     for ( ; n_id != nodeToNetgenID[isDbl].end(); ++n_id )
 565       nodeVec[ n_id->second ] = n_id->first;
 566     nodeToNetgenID[isDbl].clear();
 567   }
 568   if ( status == NG_VOLUME_FAILURE )
 569   {
 570     SMESH_ComputeErrorPtr err = NETGENPlugin_Mesher::readErrors(nodeVec);
 571     if ( err && !err->myBadElements.empty() )
 572       error( err );
 573   }
 574
 575   bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
 576   if ( isOK )
 577   {
 578     // create and insert new nodes into nodeVec
 579     int nodeIndex = Netgen_NbOfNodes + 1;
 580     int shapeID = meshDS->ShapeToIndex( aShape );
 581     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 582     {
 583       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
 584       SMDS_MeshNode * node = meshDS->AddNode(Netgen_point[0],
 585                                              Netgen_point[1],
 586                                              Netgen_point[2]);
 587       meshDS->SetNodeInVolume(node, shapeID);
 588       nodeVec.at(nodeIndex) = node;
 589     }
 590
 591     // create tetrahedrons
 592     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 593     {
 594       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 595       SMDS_MeshVolume * elt = helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
 596                                                 nodeVec.at( Netgen_tetrahedron[1] ),
 597                                                 nodeVec.at( Netgen_tetrahedron[2] ),
 598                                                 nodeVec.at( Netgen_tetrahedron[3] ));
 599       meshDS->SetMeshElementOnShape(elt, shapeID );
 600     }
 601   }
 602
 603   return (status == NG_OK);
 604 }
 605
 606 //================================================================================
 607 /*!
 608  * \brief Compute tetrahedral mesh from 2D mesh without geometry
 609  */
 610 //================================================================================
 611
 612 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
 613                                      SMESH_MesherHelper* aHelper)
 614 {
 615   MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
 616   const int invalid_ID = -1;
 617   bool _quadraticMesh = false;
 618   typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 619   TNodeToIDMap nodeToNetgenID;
 620   list< const SMDS_MeshElement* > triangles;
 621   SMESHDS_Mesh* MeshDS = aHelper->GetMeshDS();
 622
 623   SMESH_MesherHelper::MType MeshType = aHelper->IsQuadraticMesh();
 624
 625   if(MeshType == SMESH_MesherHelper::COMP)
 626     return error( COMPERR_BAD_INPUT_MESH,
 627                   SMESH_Comment("Mesh with linear and quadratic elements given."));
 628   else if (MeshType == SMESH_MesherHelper::QUADRATIC)
 629     _quadraticMesh = true;
 630
 631   StdMeshers_QuadToTriaAdaptor Adaptor;
 632   Adaptor.Compute(aMesh);
 633
 634   SMDS_FaceIteratorPtr fIt = MeshDS->facesIterator();
 635   TIDSortedElemSet sortedFaces; //  0020279: control the "random" use when using mesh algorithms
 636   while( fIt->more()) sortedFaces.insert( fIt->next() );
 637
 638   TIDSortedElemSet::iterator itFace = sortedFaces.begin(), fEnd = sortedFaces.end();
 639   for ( ; itFace != fEnd; ++itFace )
 640   {
 641     // check element
 642     const SMDS_MeshElement* elem = *itFace;
 643     if ( !elem )
 644       return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 645
 646     vector< const SMDS_MeshElement* > trias;
 647     bool isTraingle = ( elem->NbNodes() == ( elem->IsQuadratic() ? 6 : 3 ));
 648     if ( !isTraingle ) {
 649       // using adaptor
 650       const list<const SMDS_FaceOfNodes*>* faces = Adaptor.GetTriangles(elem);
 651       if(faces==0)
 652         return error( COMPERR_BAD_INPUT_MESH,
 653                       SMESH_Comment("No triangles in adaptor for element ")<<elem->GetID());
 654       trias.assign( faces->begin(), faces->end() );
 655     }
 656     else {
 657       trias.push_back( elem );
 658     }
 659     for ( int i = 0; i < trias.size(); ++i )
 660     {
 661       triangles.push_back( trias[i] );
 662       for ( int iN = 0; iN < 3; ++iN )
 663       {
 664         const SMDS_MeshNode* node = trias[i]->GetNode( iN );
 665         // put elem nodes to nodeToNetgenID map
 666         nodeToNetgenID.insert( make_pair( node, invalid_ID ));
 667       }
 668     }
 669   }
 670
 671   // ---------------------------------
 672   // Feed the Netgen with surface mesh
 673   // ---------------------------------
 674
 675   int Netgen_NbOfNodes = 0;
 676   int Netgen_param2ndOrder = 0;
 677   double Netgen_paramFine = 1.;
 678   double Netgen_paramSize = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
 679
 680   double Netgen_point[3];
 681   int Netgen_triangle[3];
 682   int Netgen_tetrahedron[4];
 683
 684   TNetgenLibWrapper ngLib;
 685   Ng_Mesh * Netgen_mesh = ngLib._ngMesh;
 686
 687   // set nodes and remember thier netgen IDs
 688
 689   TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 690   for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 691   {
 692     const SMDS_MeshNode* node = n_id->first;
 693
 694     Netgen_point [ 0 ] = node->X();
 695     Netgen_point [ 1 ] = node->Y();
 696     Netgen_point [ 2 ] = node->Z();
 697     Ng_AddPoint(Netgen_mesh, Netgen_point);
 698     n_id->second = ++Netgen_NbOfNodes; // set netgen ID
 699   }
 700
 701   // set triangles
 702   list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
 703   for ( ; tria != triangles.end(); ++tria)
 704   {
 705     int i = 0;
 706     SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
 707     while ( triangleNodesIt->more() ) {
 708       const SMDS_MeshNode * node =
 709         static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 710       if(aHelper->IsMedium(node))
 711         continue;
 712       Netgen_triangle[ i ] = nodeToNetgenID[ node ];
 713       ++i;
 714     }
 715     Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 716   }
 717
 718   // -------------------------
 719   // Generate the volume mesh
 720   // -------------------------
 721
 722   Ng_Meshing_Parameters Netgen_param;
 723
 724   Netgen_param.secondorder = Netgen_param2ndOrder;
 725   Netgen_param.fineness = Netgen_paramFine;
 726   Netgen_param.maxh = Netgen_paramSize;
 727
 728   Ng_Result status;
 729
 730   try {
 731 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 732     OCC_CATCH_SIGNALS;
 733 #endif
 734     status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
 735   }
 736   catch (Standard_Failure& exc) {
 737     error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
 738     status = NG_VOLUME_FAILURE;
 739   }
 740   catch (...) {
 741     error("Bad mesh input!!!");
 742     status = NG_VOLUME_FAILURE;
 743   }
 744   if ( GetComputeError()->IsOK() ) {
 745     error( status, "Bad mesh input!!!");
 746   }
 747
 748   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 749
 750   int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
 751
 752   MESSAGE("End of Volume Mesh Generation. status=" << status <<
 753           ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
 754           ", nb tetra: " << Netgen_NbOfTetra);
 755
 756   // -------------------------------------------------------------------
 757   // Feed back the SMESHDS with the generated Nodes and Volume Elements
 758   // -------------------------------------------------------------------
 759
 760   bool isOK = ( Netgen_NbOfTetra > 0 );// get whatever built
 761   if ( isOK )
 762   {
 763     // vector of nodes in which node index == netgen ID
 764     vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
 765     // insert old nodes into nodeVec
 766     for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
 767       nodeVec.at( n_id->second ) = n_id->first;
 768     }
 769     // create and insert new nodes into nodeVec
 770     int nodeIndex = Netgen_NbOfNodes + 1;
 771
 772     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 773     {
 774       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
 775       SMDS_MeshNode * node = aHelper->AddNode(Netgen_point[0],
 776                                               Netgen_point[1],
 777                                               Netgen_point[2]);
 778       nodeVec.at(nodeIndex) = node;
 779     }
 780
 781     // create tetrahedrons
 782     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 783     {
 784       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 785       aHelper->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
 786                           nodeVec.at( Netgen_tetrahedron[1] ),
 787                           nodeVec.at( Netgen_tetrahedron[2] ),
 788                           nodeVec.at( Netgen_tetrahedron[3] ));
 789     }
 790   }
 791
 792   return (status == NG_OK);
 793 }
 794
 795
 796 //=============================================================================
 797 /*!
 798  *
 799  */
 800 //=============================================================================
 801
 802 bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh,
 803                                       const TopoDS_Shape& aShape,
 804                                       MapShapeNbElems& aResMap)
 805 {
 806   int nbtri = 0, nbqua = 0;
 807   double fullArea = 0.0;
 808   for (TopExp_Explorer expF(aShape, TopAbs_FACE); expF.More(); expF.Next()) {
 809     TopoDS_Face F = TopoDS::Face( expF.Current() );
 810     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
 811     MapShapeNbElemsItr anIt = aResMap.find(sm);
 812     if( anIt==aResMap.end() ) {
 813       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 814       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 815       return false;
 816     }
 817     std::vector<int> aVec = (*anIt).second;
 818     nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
 819     nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 820     GProp_GProps G;
 821     BRepGProp::SurfaceProperties(F,G);
 822     double anArea = G.Mass();
 823     fullArea += anArea;
 824   }
 825
 826   // collect info from edges
 827   int nb0d_e = 0, nb1d_e = 0;
 828   bool IsQuadratic = false;
 829   bool IsFirst = true;
 830   TopTools_MapOfShape tmpMap;
 831   for (TopExp_Explorer expF(aShape, TopAbs_EDGE); expF.More(); expF.Next()) {
 832     TopoDS_Edge E = TopoDS::Edge(expF.Current());
 833     if( tmpMap.Contains(E) )
 834       continue;
 835     tmpMap.Add(E);
 836     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(expF.Current());
 837     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
 838     if( anIt==aResMap.end() ) {
 839       SMESH_ComputeErrorPtr& smError = aSubMesh->GetComputeError();
 840       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
 841                                             "Submesh can not be evaluated",this));
 842       return false;
 843     }
 844     std::vector<int> aVec = (*anIt).second;
 845     nb0d_e += aVec[SMDSEntity_Node];
 846     nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
 847     if(IsFirst) {
 848       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
 849       IsFirst = false;
 850     }
 851   }
 852   tmpMap.Clear();
 853
 854   double ELen_face = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
 855   double ELen_vol = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
 856   double ELen = Min(ELen_vol,ELen_face*2);
 857
 858   GProp_GProps G;
 859   BRepGProp::VolumeProperties(aShape,G);
 860   double aVolume = G.Mass();
 861   double tetrVol = 0.1179*ELen*ELen*ELen;
 862   double CoeffQuality = 0.9;
 863   int nbVols = int( aVolume/tetrVol/CoeffQuality );
 864   int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
 865   int nb1d_in = (nbVols*6 - nb1d_e - nb1d_f ) / 5;
 866   std::vector<int> aVec(SMDSEntity_Last);
 867   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
 868   if( IsQuadratic ) {
 869     aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
 870     aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
 871     aVec[SMDSEntity_Quad_Pyramid] = nbqua;
 872   }
 873   else {
 874     aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
 875     aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
 876     aVec[SMDSEntity_Pyramid] = nbqua;
 877   }
 878   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
 879   aResMap.insert(std::make_pair(sm,aVec));
 880
 881   return true;
 882 }