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 <BRep_Tool.hxx>
  46 #include <GProp_GProps.hxx>
  47 #include <BRepGProp.hxx>
  48 #include <TopExp.hxx>
  49 #include <TopExp_Explorer.hxx>
  50 #include <TopoDS.hxx>
  51
  52 #include <Standard_Failure.hxx>
  53 #include <Standard_ErrorHandler.hxx>
  54
  55 #include "utilities.h"
  56
  57 #include <list>
  58 #include <vector>
  59 #include <map>
  60
  61 /*
  62   Netgen include files
  63 */
  64
  65 namespace nglib {
  66 #include <nglib.h>
  67 }
  68 using namespace nglib;
  69 using namespace std;
  70
  71 //=============================================================================
  72 /*!
  73  *
  74  */
  75 //=============================================================================
  76
  77 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
  78                              SMESH_Gen* gen)
  79   : SMESH_3D_Algo(hypId, studyId, gen)
  80 {
  81   MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
  82   _name = "NETGEN_3D";
  83   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  84   _compatibleHypothesis.push_back("MaxElementVolume");
  85
  86   _maxElementVolume = 0.;
  87
  88   _hypMaxElementVolume = NULL;
  89
  90   _requireShape = false; // can work without shape
  91 }
  92
  93 //=============================================================================
  94 /*!
  95  *
  96  */
  97 //=============================================================================
  98
  99 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
 100 {
 101   MESSAGE("NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D");
 102 }
 103
 104 //=============================================================================
 105 /*!
 106  *
 107  */
 108 //=============================================================================
 109
 110 bool NETGENPlugin_NETGEN_3D::CheckHypothesis
 111                          (SMESH_Mesh& aMesh,
 112                           const TopoDS_Shape& aShape,
 113                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 114 {
 115   MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
 116
 117   _hypMaxElementVolume = NULL;
 118   _maxElementVolume = DBL_MAX;
 119
 120   list<const SMESHDS_Hypothesis*>::const_iterator itl;
 121   const SMESHDS_Hypothesis* theHyp;
 122
 123   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
 124   int nbHyp = hyps.size();
 125   if (!nbHyp)
 126   {
 127     aStatus = SMESH_Hypothesis::HYP_OK;
 128     //aStatus = SMESH_Hypothesis::HYP_MISSING;
 129     return true;  // can work with no hypothesis
 130   }
 131
 132   itl = hyps.begin();
 133   theHyp = (*itl); // use only the first hypothesis
 134
 135   string hypName = theHyp->GetName();
 136
 137   bool isOk = false;
 138
 139   if (hypName == "MaxElementVolume")
 140   {
 141     _hypMaxElementVolume = static_cast<const StdMeshers_MaxElementVolume*> (theHyp);
 142     ASSERT(_hypMaxElementVolume);
 143     _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
 144     isOk =true;
 145     aStatus = SMESH_Hypothesis::HYP_OK;
 146   }
 147   else
 148     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 149
 150   return isOk;
 151 }
 152
 153 //=============================================================================
 154 /*!
 155  *Here we are going to use the NETGEN mesher
 156  */
 157 //=============================================================================
 158
 159 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
 160                                      const TopoDS_Shape& aShape)
 161 {
 162   MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
 163
 164   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 165
 166   const int invalid_ID = -1;
 167
 168   SMESH::Controls::Area areaControl;
 169   SMESH::Controls::TSequenceOfXYZ nodesCoords;
 170
 171   // -------------------------------------------------------------------
 172   // get triangles on aShell and make a map of nodes to Netgen node IDs
 173   // -------------------------------------------------------------------
 174
 175   SMESH_MesherHelper helper(aMesh);
 176   SMESH_MesherHelper* myTool = &helper;
 177   bool _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
 178
 179   typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 180   TNodeToIDMap nodeToNetgenID;
 181   list< const SMDS_MeshElement* > triangles;
 182   list< bool >                    isReversed; // orientation of triangles
 183
 184   TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
 185   bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
 186
 187   // for the degeneraged edge: ignore all but one node on it;
 188   // map storing ids of degen edges and vertices and their netgen id:
 189   map< int, int* > degenShapeIdToPtrNgId;
 190   map< int, int* >::iterator shId_ngId;
 191   list< int > degenNgIds;
 192
 193   StdMeshers_QuadToTriaAdaptor Adaptor;
 194   Adaptor.Compute(aMesh,aShape);
 195
 196   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
 197   {
 198     const TopoDS_Shape& aShapeFace = exp.Current();
 199     const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
 200     if ( aSubMeshDSFace )
 201     {
 202       bool isRev = false;
 203       if ( checkReverse && helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
 204         // IsReversedSubMesh() can work wrong on strongly curved faces,
 205         // so we use it as less as possible
 206         isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
 207
 208       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 209       while ( iteratorElem->more() ) // loop on elements on a face
 210       {
 211         // check element
 212         const SMDS_MeshElement* elem = iteratorElem->next();
 213         if ( !elem )
 214           return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 215         bool isTraingle = ( elem->NbNodes() == ( elem->IsQuadratic() ? 6 : 3 ));
 216         if ( !isTraingle ) {
 217           // using adaptor
 218           const list<const SMDS_FaceOfNodes*>* faces = Adaptor.GetTriangles(elem);
 219           if(faces==0) {
 220             return error( COMPERR_BAD_INPUT_MESH,
 221                           SMESH_Comment("No triangles in adaptor for element ")<<elem->GetID());
 222           }
 223           list<const SMDS_FaceOfNodes*>::const_iterator itf = faces->begin();
 224           for(; itf!=faces->end(); itf++ ) {
 225             triangles.push_back( (*itf) );
 226             isReversed.push_back( isRev );
 227             // put triange's nodes to nodeToNetgenID map
 228             SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator();
 229             while ( triangleNodesIt->more() ) {
 230               const SMDS_MeshNode * node =
 231                 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 232               if(myTool->IsMedium(node))
 233                 continue;
 234               nodeToNetgenID.insert( make_pair( node, invalid_ID ));
 235             }
 236           }
 237         }
 238         else {
 239           // keep a triangle
 240           triangles.push_back( elem );
 241           isReversed.push_back( isRev );
 242           // put elem nodes to nodeToNetgenID map
 243           SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
 244           while ( triangleNodesIt->more() ) {
 245             const SMDS_MeshNode * node =
 246               static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 247             if(myTool->IsMedium(node))
 248               continue;
 249             nodeToNetgenID.insert( make_pair( node, invalid_ID ));
 250           }
 251         }
 252 #ifdef _DEBUG_
 253         // check if a trainge is degenerated
 254         areaControl.GetPoints( elem, nodesCoords );
 255         double area = areaControl.GetValue( nodesCoords );
 256         if ( area <= DBL_MIN ) {
 257           MESSAGE( "Warning: Degenerated " << elem );
 258         }
 259 #endif
 260       }
 261       // look for degeneraged edges and vetices
 262       for (TopExp_Explorer expE(aShapeFace,TopAbs_EDGE);expE.More();expE.Next())
 263       {
 264         TopoDS_Edge aShapeEdge = TopoDS::Edge( expE.Current() );
 265         if ( BRep_Tool::Degenerated( aShapeEdge ))
 266         {
 267           degenNgIds.push_back( invalid_ID );
 268           int* ptrIdOnEdge = & degenNgIds.back();
 269           // remember edge id
 270           int edgeID = meshDS->ShapeToIndex( aShapeEdge );
 271           degenShapeIdToPtrNgId.insert( make_pair( edgeID, ptrIdOnEdge ));
 272           // remember vertex id
 273           int vertexID = meshDS->ShapeToIndex( TopExp::FirstVertex( aShapeEdge ));
 274           degenShapeIdToPtrNgId.insert( make_pair( vertexID, ptrIdOnEdge ));
 275         }
 276       }
 277     }
 278   }
 279   // ---------------------------------
 280   // Feed the Netgen with surface mesh
 281   // ---------------------------------
 282
 283   int Netgen_NbOfNodes = 0;
 284   int Netgen_param2ndOrder = 0;
 285   double Netgen_paramFine = 1.;
 286   double Netgen_paramSize = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
 287
 288   double Netgen_point[3];
 289   int Netgen_triangle[3];
 290   int Netgen_tetrahedron[4];
 291
 292   Ng_Init();
 293
 294   Ng_Mesh * Netgen_mesh = Ng_NewMesh();
 295
 296   // set nodes and remember thier netgen IDs
 297   bool isDegen = false, hasDegen = !degenShapeIdToPtrNgId.empty();
 298   TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 299   for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 300   {
 301     const SMDS_MeshNode* node = n_id->first;
 302
 303     // ignore nodes on degenerated edge
 304     if ( hasDegen ) {
 305       int shapeId = node->GetPosition()->GetShapeId();
 306       shId_ngId = degenShapeIdToPtrNgId.find( shapeId );
 307       isDegen = ( shId_ngId != degenShapeIdToPtrNgId.end() );
 308       if ( isDegen && *(shId_ngId->second) != invalid_ID ) {
 309         n_id->second = *(shId_ngId->second);
 310         continue;
 311       }
 312     }
 313     Netgen_point [ 0 ] = node->X();
 314     Netgen_point [ 1 ] = node->Y();
 315     Netgen_point [ 2 ] = node->Z();
 316     Ng_AddPoint(Netgen_mesh, Netgen_point);
 317     n_id->second = ++Netgen_NbOfNodes; // set netgen ID
 318
 319     if ( isDegen ) // all nodes on a degen edge get one netgen ID
 320       *(shId_ngId->second) = n_id->second;
 321   }
 322
 323   // set triangles
 324   list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
 325   list< bool >::iterator                 reverse = isReversed.begin();
 326   for ( ; tria != triangles.end(); ++tria, ++reverse )
 327   {
 328     int i = 0;
 329     SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
 330     while ( triangleNodesIt->more() ) {
 331       const SMDS_MeshNode * node =
 332         static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 333       if(myTool->IsMedium(node))
 334         continue;
 335       Netgen_triangle[ *reverse ? 2 - i : i ] = nodeToNetgenID[ node ];
 336       ++i;
 337     }
 338     if ( !hasDegen ||
 339          // ignore degenerated triangles, they have 2 or 3 same ids
 340          (Netgen_triangle[0] != Netgen_triangle[1] &&
 341           Netgen_triangle[0] != Netgen_triangle[2] &&
 342           Netgen_triangle[2] != Netgen_triangle[1] ))
 343     {
 344       Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 345     }
 346   }
 347
 348   // -------------------------
 349   // Generate the volume mesh
 350   // -------------------------
 351
 352   Ng_Meshing_Parameters Netgen_param;
 353
 354   Netgen_param.secondorder = Netgen_param2ndOrder;
 355   Netgen_param.fineness = Netgen_paramFine;
 356   Netgen_param.maxh = Netgen_paramSize;
 357
 358   Ng_Result status;
 359
 360   try {
 361 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 362     OCC_CATCH_SIGNALS;
 363 #endif
 364     status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
 365   }
 366   catch (Standard_Failure& exc) {
 367     error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
 368     status = NG_VOLUME_FAILURE;
 369   }
 370   catch (...) {
 371     error("Exception in Ng_GenerateVolumeMesh()");
 372     status = NG_VOLUME_FAILURE;
 373   }
 374   if ( GetComputeError()->IsOK() ) {
 375     switch ( status ) {
 376     case NG_SURFACE_INPUT_ERROR:error( status, "NG_SURFACE_INPUT_ERROR");
 377     case NG_VOLUME_FAILURE:     error( status, "NG_VOLUME_FAILURE");
 378     case NG_STL_INPUT_ERROR:    error( status, "NG_STL_INPUT_ERROR");
 379     case NG_SURFACE_FAILURE:    error( status, "NG_SURFACE_FAILURE");
 380     case NG_FILE_NOT_FOUND:     error( status, "NG_FILE_NOT_FOUND");
 381     };
 382   }
 383
 384   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 385
 386   int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
 387
 388   MESSAGE("End of Volume Mesh Generation. status=" << status <<
 389           ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
 390           ", nb tetra: " << Netgen_NbOfTetra);
 391
 392   // -------------------------------------------------------------------
 393   // Feed back the SMESHDS with the generated Nodes and Volume Elements
 394   // -------------------------------------------------------------------
 395
 396   bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
 397   if ( isOK )
 398   {
 399     // vector of nodes in which node index == netgen ID
 400     vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
 401     // insert old nodes into nodeVec
 402     for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
 403       nodeVec.at( n_id->second ) = n_id->first;
 404     }
 405     // create and insert new nodes into nodeVec
 406     int nodeIndex = Netgen_NbOfNodes + 1;
 407     int shapeID = meshDS->ShapeToIndex( aShape );
 408     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 409     {
 410       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
 411       SMDS_MeshNode * node = meshDS->AddNode(Netgen_point[0],
 412                                              Netgen_point[1],
 413                                              Netgen_point[2]);
 414       meshDS->SetNodeInVolume(node, shapeID);
 415       nodeVec.at(nodeIndex) = node;
 416     }
 417
 418     // create tetrahedrons
 419     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 420     {
 421       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 422       SMDS_MeshVolume * elt = myTool->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
 423                                                  nodeVec.at( Netgen_tetrahedron[1] ),
 424                                                  nodeVec.at( Netgen_tetrahedron[2] ),
 425                                                  nodeVec.at( Netgen_tetrahedron[3] ));
 426       meshDS->SetMeshElementOnShape(elt, shapeID );
 427     }
 428   }
 429
 430   Ng_DeleteMesh(Netgen_mesh);
 431   Ng_Exit();
 432
 433   NETGENPlugin_Mesher::RemoveTmpFiles();
 434
 435   return (status == NG_OK);
 436 }
 437
 438 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
 439                                      SMESH_MesherHelper* aHelper)
 440 {
 441   MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
 442   const int invalid_ID = -1;
 443   bool _quadraticMesh = false;
 444   typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 445   TNodeToIDMap nodeToNetgenID;
 446   list< const SMDS_MeshElement* > triangles;
 447   SMESHDS_Mesh* MeshDS = aHelper->GetMeshDS();
 448
 449   SMESH_MesherHelper::MType MeshType = aHelper->IsQuadraticMesh();
 450
 451   if(MeshType == SMESH_MesherHelper::COMP)
 452     return error( COMPERR_BAD_INPUT_MESH,
 453                   SMESH_Comment("Mesh with linear and quadratic elements given."));
 454   else if (MeshType == SMESH_MesherHelper::QUADRATIC)
 455     _quadraticMesh = true;
 456
 457   StdMeshers_QuadToTriaAdaptor Adaptor;
 458   Adaptor.Compute(aMesh);
 459
 460   SMDS_FaceIteratorPtr fIt = MeshDS->facesIterator();
 461   TIDSortedElemSet sortedFaces; //  0020279: control the "random" use when using mesh algorithms
 462   while( fIt->more()) sortedFaces.insert( fIt->next() );
 463
 464   TIDSortedElemSet::iterator itFace = sortedFaces.begin(), fEnd = sortedFaces.end();
 465   for ( ; itFace != fEnd; ++itFace ) {
 466     // check element
 467     const SMDS_MeshElement* elem = *itFace;
 468     if ( !elem )
 469       return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 470     bool isTraingle = ( elem->NbNodes() == ( elem->IsQuadratic() ? 6 : 3 ));
 471     if ( !isTraingle ) {
 472       // using adaptor
 473       const list<const SMDS_FaceOfNodes*>* faces = Adaptor.GetTriangles(elem);
 474       if(faces==0) {
 475         continue; // Issue 0020682. There already can be 3d mesh
 476       }
 477       list<const SMDS_FaceOfNodes*>::const_iterator itf = faces->begin();
 478       for(; itf!=faces->end(); itf++ ) {
 479         triangles.push_back( (*itf) );
 480         // put triange's nodes to nodeToNetgenID map
 481         SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator();
 482         while ( triangleNodesIt->more() ) {
 483           const SMDS_MeshNode * node =
 484             static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 485           if(aHelper->IsMedium(node))
 486             continue;
 487           nodeToNetgenID.insert( make_pair( node, invalid_ID ));
 488         }
 489       }
 490     }
 491     else {
 492       // keep a triangle
 493       triangles.push_back( elem );
 494       // put elem nodes to nodeToNetgenID map
 495       SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
 496       while ( triangleNodesIt->more() ) {
 497         const SMDS_MeshNode * node =
 498           static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 499         if(aHelper->IsMedium(node))
 500           continue;
 501         nodeToNetgenID.insert( make_pair( node, invalid_ID ));
 502       }
 503     }
 504   }
 505
 506   // ---------------------------------
 507   // Feed the Netgen with surface mesh
 508   // ---------------------------------
 509
 510   int Netgen_NbOfNodes = 0;
 511   int Netgen_param2ndOrder = 0;
 512   double Netgen_paramFine = 1.;
 513   double Netgen_paramSize = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
 514
 515   double Netgen_point[3];
 516   int Netgen_triangle[3];
 517   int Netgen_tetrahedron[4];
 518
 519   Ng_Init();
 520
 521   Ng_Mesh * Netgen_mesh = Ng_NewMesh();
 522
 523     // set nodes and remember thier netgen IDs
 524
 525   TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 526   for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 527   {
 528     const SMDS_MeshNode* node = n_id->first;
 529
 530     Netgen_point [ 0 ] = node->X();
 531     Netgen_point [ 1 ] = node->Y();
 532     Netgen_point [ 2 ] = node->Z();
 533     Ng_AddPoint(Netgen_mesh, Netgen_point);
 534     n_id->second = ++Netgen_NbOfNodes; // set netgen ID
 535   }
 536
 537   // set triangles
 538   list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
 539   for ( ; tria != triangles.end(); ++tria)
 540   {
 541     int i = 0;
 542     SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
 543     while ( triangleNodesIt->more() ) {
 544       const SMDS_MeshNode * node =
 545         static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 546       if(aHelper->IsMedium(node))
 547         continue;
 548       Netgen_triangle[ i ] = nodeToNetgenID[ node ];
 549       ++i;
 550     }
 551     Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 552   }
 553
 554   // -------------------------
 555   // Generate the volume mesh
 556   // -------------------------
 557
 558   Ng_Meshing_Parameters Netgen_param;
 559
 560   Netgen_param.secondorder = Netgen_param2ndOrder;
 561   Netgen_param.fineness = Netgen_paramFine;
 562   Netgen_param.maxh = Netgen_paramSize;
 563
 564   Ng_Result status;
 565
 566   try {
 567 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 568     OCC_CATCH_SIGNALS;
 569 #endif
 570     status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
 571   }
 572   catch (Standard_Failure& exc) {
 573     error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
 574     status = NG_VOLUME_FAILURE;
 575   }
 576   catch (...) {
 577     error("Bad mesh input!!!");
 578     status = NG_VOLUME_FAILURE;
 579   }
 580   if ( GetComputeError()->IsOK() ) {
 581     error( status, "Bad mesh input!!!");
 582   }
 583
 584   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 585
 586   int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
 587
 588   MESSAGE("End of Volume Mesh Generation. status=" << status <<
 589           ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
 590           ", nb tetra: " << Netgen_NbOfTetra);
 591
 592   // -------------------------------------------------------------------
 593   // Feed back the SMESHDS with the generated Nodes and Volume Elements
 594   // -------------------------------------------------------------------
 595
 596   bool isOK = ( Netgen_NbOfTetra > 0 );// get whatever built
 597   if ( isOK )
 598   {
 599     // vector of nodes in which node index == netgen ID
 600     vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
 601     // insert old nodes into nodeVec
 602     for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
 603       nodeVec.at( n_id->second ) = n_id->first;
 604     }
 605     // create and insert new nodes into nodeVec
 606     int nodeIndex = Netgen_NbOfNodes + 1;
 607
 608     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 609     {
 610       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
 611       SMDS_MeshNode * node = aHelper->AddNode(Netgen_point[0],
 612                                               Netgen_point[1],
 613                                               Netgen_point[2]);
 614       nodeVec.at(nodeIndex) = node;
 615     }
 616
 617     // create tetrahedrons
 618     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 619     {
 620       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 621       aHelper->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
 622                           nodeVec.at( Netgen_tetrahedron[1] ),
 623                           nodeVec.at( Netgen_tetrahedron[2] ),
 624                           nodeVec.at( Netgen_tetrahedron[3] ));
 625     }
 626   }
 627
 628   Ng_DeleteMesh(Netgen_mesh);
 629   Ng_Exit();
 630
 631   NETGENPlugin_Mesher::RemoveTmpFiles();
 632
 633   return (status == NG_OK);
 634 }
 635
 636
 637 //=============================================================================
 638 /*!
 639  *
 640  */
 641 //=============================================================================
 642
 643 bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh,
 644                                       const TopoDS_Shape& aShape,
 645                                       MapShapeNbElems& aResMap)
 646 {
 647   int nbtri = 0, nbqua = 0;
 648   double fullArea = 0.0;
 649   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
 650     TopoDS_Face F = TopoDS::Face( exp.Current() );
 651     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
 652     MapShapeNbElemsItr anIt = aResMap.find(sm);
 653     if( anIt==aResMap.end() ) {
 654       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 655       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 656       return false;
 657     }
 658     std::vector<int> aVec = (*anIt).second;
 659     nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
 660     nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 661     GProp_GProps G;
 662     BRepGProp::SurfaceProperties(F,G);
 663     double anArea = G.Mass();
 664     fullArea += anArea;
 665   }
 666
 667   // collect info from edges
 668   int nb0d_e = 0, nb1d_e = 0;
 669   bool IsQuadratic = false;
 670   bool IsFirst = true;
 671   TopTools_MapOfShape tmpMap;
 672   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
 673     TopoDS_Edge E = TopoDS::Edge(exp.Current());
 674     if( tmpMap.Contains(E) )
 675       continue;
 676     tmpMap.Add(E);
 677     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
 678     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
 679     if( anIt==aResMap.end() ) {
 680       SMESH_ComputeErrorPtr& smError = aSubMesh->GetComputeError();
 681       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
 682                                             "Submesh can not be evaluated",this));
 683       return false;
 684     }
 685     std::vector<int> aVec = (*anIt).second;
 686     nb0d_e += aVec[SMDSEntity_Node];
 687     nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
 688     if(IsFirst) {
 689       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
 690       IsFirst = false;
 691     }
 692   }
 693   tmpMap.Clear();
 694
 695   double ELen_face = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
 696   double ELen_vol = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
 697   double ELen = Min(ELen_vol,ELen_face*2);
 698
 699   GProp_GProps G;
 700   BRepGProp::VolumeProperties(aShape,G);
 701   double aVolume = G.Mass();
 702   double tetrVol = 0.1179*ELen*ELen*ELen;
 703   double CoeffQuality = 0.9;
 704   int nbVols = int( aVolume/tetrVol/CoeffQuality );
 705   int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
 706   int nb1d_in = (nbVols*6 - nb1d_e - nb1d_f ) / 5;
 707   std::vector<int> aVec(SMDSEntity_Last);
 708   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
 709   if( IsQuadratic ) {
 710     aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
 711     aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
 712     aVec[SMDSEntity_Quad_Pyramid] = nbqua;
 713   }
 714   else {
 715     aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
 716     aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
 717     aVec[SMDSEntity_Pyramid] = nbqua;
 718   }
 719   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
 720   aResMap.insert(std::make_pair(sm,aVec));
 721
 722   return true;
 723 }