src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx

   1 //=============================================================================
   2 // File      : NETGENPlugin_NETGEN_3D.cxx
   3 //             Moved here from SMESH_NETGEN_3D.cxx
   4 // Created   : lundi 27 Janvier 2003
   5 // Author    : Nadir BOUHAMOU (CEA)
   6 // Project   : SALOME
   7 // Copyright : CEA 2003
   8 // $Header$
   9 //=============================================================================
  10 using namespace std;
  11
  12 #include "NETGENPlugin_NETGEN_3D.hxx"
  13
  14 #include "SMESH_Gen.hxx"
  15 #include "SMESH_Mesh.hxx"
  16 #include "SMESH_ControlsDef.hxx"
  17 #include "SMESHDS_Mesh.hxx"
  18 #include "SMDS_MeshElement.hxx"
  19 #include "SMDS_MeshNode.hxx"
  20 #include "StdMeshers_Helper.hxx"
  21
  22 #include <BRep_Tool.hxx>
  23 #include <TopExp.hxx>
  24 #include <TopExp_Explorer.hxx>
  25 #include <TopoDS.hxx>
  26
  27 #include "utilities.h"
  28
  29 #include <list>
  30 #include <vector>
  31 #include <map>
  32
  33 /*
  34   Netgen include files
  35 */
  36
  37 #include "nglib.h"
  38
  39 //=============================================================================
  40 /*!
  41  *
  42  */
  43 //=============================================================================
  44
  45 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
  46                              SMESH_Gen* gen)
  47   : SMESH_3D_Algo(hypId, studyId, gen)
  48 {
  49   MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
  50   _name = "NETGEN_3D";
  51   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  52   _compatibleHypothesis.push_back("MaxElementVolume");
  53
  54   _maxElementVolume = 0.;
  55
  56   _hypMaxElementVolume = NULL;
  57 }
  58
  59 //=============================================================================
  60 /*!
  61  *
  62  */
  63 //=============================================================================
  64
  65 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
  66 {
  67   MESSAGE("NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D");
  68 }
  69
  70 //=============================================================================
  71 /*!
  72  *
  73  */
  74 //=============================================================================
  75
  76 bool NETGENPlugin_NETGEN_3D::CheckHypothesis
  77                          (SMESH_Mesh& aMesh,
  78                           const TopoDS_Shape& aShape,
  79                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
  80 {
  81   MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
  82
  83   _hypMaxElementVolume = NULL;
  84   _maxElementVolume = DBL_MAX;
  85
  86   list<const SMESHDS_Hypothesis*>::const_iterator itl;
  87   const SMESHDS_Hypothesis* theHyp;
  88
  89   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
  90   int nbHyp = hyps.size();
  91   if (!nbHyp)
  92   {
  93     aStatus = SMESH_Hypothesis::HYP_OK;
  94     //aStatus = SMESH_Hypothesis::HYP_MISSING;
  95     return true;  // can work with no hypothesis
  96   }
  97
  98   itl = hyps.begin();
  99   theHyp = (*itl); // use only the first hypothesis
 100
 101   string hypName = theHyp->GetName();
 102
 103   bool isOk = false;
 104
 105   if (hypName == "MaxElementVolume")
 106   {
 107     _hypMaxElementVolume = static_cast<const StdMeshers_MaxElementVolume*> (theHyp);
 108     ASSERT(_hypMaxElementVolume);
 109     _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
 110     isOk =true;
 111     aStatus = SMESH_Hypothesis::HYP_OK;
 112   }
 113   else
 114     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 115
 116   return isOk;
 117 }
 118
 119 //=============================================================================
 120 /*!
 121  *Here we are going to use the NETGEN mesher
 122  */
 123 //=============================================================================
 124
 125 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
 126                                      const TopoDS_Shape& aShape)
 127 {
 128   MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
 129
 130   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 131
 132   const int invalid_ID = -1;
 133
 134   SMESH::Controls::Area areaControl;
 135   SMESH::Controls::TSequenceOfXYZ nodesCoords;
 136
 137   // -------------------------------------------------------------------
 138   // get triangles on aShell and make a map of nodes to Netgen node IDs
 139   // -------------------------------------------------------------------
 140
 141   StdMeshers_Helper* myTool = new StdMeshers_Helper(aMesh);
 142   bool _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
 143
 144   typedef map< const SMDS_MeshNode*, int> TNodeToIDMap;
 145   TNodeToIDMap nodeToNetgenID;
 146   list< const SMDS_MeshElement* > triangles;
 147   list< bool >                    isReversed; // orientation of triangles
 148
 149   // for the degeneraged edge: ignore all but one node on it;
 150   // map storing ids of degen edges and vertices and their netgen id:
 151   map< int, int* > degenShapeIdToPtrNgId;
 152   map< int, int* >::iterator shId_ngId;
 153   list< int > degenNgIds;
 154
 155   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
 156   {
 157     const TopoDS_Shape& aShapeFace = exp.Current();
 158     const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
 159     if ( aSubMeshDSFace )
 160     {
 161       bool isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
 162
 163       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 164       while ( iteratorElem->more() ) // loop on elements on a face
 165       {
 166         // check element
 167         const SMDS_MeshElement* elem = iteratorElem->next();
 168         if ( !elem ||
 169              !( elem->NbNodes()==3 || ( _quadraticMesh && elem->NbNodes()==6) ) ) {
 170           INFOS( "NETGENPlugin_NETGEN_3D::Compute(), bad mesh");
 171           delete myTool; myTool = 0;
 172           return false;
 173         }
 174         // keep a triangle
 175         triangles.push_back( elem );
 176         isReversed.push_back( isRev );
 177         // put elem nodes to nodeToNetgenID map
 178         SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
 179         while ( triangleNodesIt->more() ) {
 180           const SMDS_MeshNode * node =
 181             static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 182           if(myTool->IsMedium(node))
 183             continue;
 184           nodeToNetgenID.insert( make_pair( node, invalid_ID ));
 185         }
 186 #ifdef _DEBUG_
 187         // check if a trainge is degenerated
 188         areaControl.GetPoints( elem, nodesCoords );
 189         double area = areaControl.GetValue( nodesCoords );
 190         if ( area <= DBL_MIN ) {
 191           MESSAGE( "Warning: Degenerated " << elem );
 192         }
 193 #endif
 194       }
 195       // look for degeneraged edges and vetices
 196       for (TopExp_Explorer expE(aShapeFace,TopAbs_EDGE);expE.More();expE.Next())
 197       {
 198         TopoDS_Edge aShapeEdge = TopoDS::Edge( expE.Current() );
 199         if ( BRep_Tool::Degenerated( aShapeEdge ))
 200         {
 201           degenNgIds.push_back( invalid_ID );
 202           int* ptrIdOnEdge = & degenNgIds.back();
 203           // remember edge id
 204           int edgeID = meshDS->ShapeToIndex( aShapeEdge );
 205           degenShapeIdToPtrNgId.insert( make_pair( edgeID, ptrIdOnEdge ));
 206           // remember vertex id
 207           int vertexID = meshDS->ShapeToIndex( TopExp::FirstVertex( aShapeEdge ));
 208           degenShapeIdToPtrNgId.insert( make_pair( vertexID, ptrIdOnEdge ));
 209         }
 210       }
 211     }
 212   }
 213   // ---------------------------------
 214   // Feed the Netgen with surface mesh
 215   // ---------------------------------
 216
 217   int Netgen_NbOfNodes = 0;
 218   int Netgen_param2ndOrder = 0;
 219   double Netgen_paramFine = 1.;
 220   double Netgen_paramSize = _maxElementVolume;
 221
 222   double Netgen_point[3];
 223   int Netgen_triangle[3];
 224   int Netgen_tetrahedron[4];
 225
 226   Ng_Init();
 227
 228   Ng_Mesh * Netgen_mesh = Ng_NewMesh();
 229
 230   // set nodes and remember thier netgen IDs
 231   bool isDegen = false, hasDegen = !degenShapeIdToPtrNgId.empty();
 232   TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 233   for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 234   {
 235     const SMDS_MeshNode* node = n_id->first;
 236
 237     // ignore nodes on degenerated edge
 238     if ( hasDegen ) {
 239       int shapeId = node->GetPosition()->GetShapeId();
 240       shId_ngId = degenShapeIdToPtrNgId.find( shapeId );
 241       isDegen = ( shId_ngId != degenShapeIdToPtrNgId.end() );
 242       if ( isDegen && *(shId_ngId->second) != invalid_ID ) {
 243         n_id->second = *(shId_ngId->second);
 244         continue;
 245       }
 246     }
 247     Netgen_point [ 0 ] = node->X();
 248     Netgen_point [ 1 ] = node->Y();
 249     Netgen_point [ 2 ] = node->Z();
 250     Ng_AddPoint(Netgen_mesh, Netgen_point);
 251     n_id->second = ++Netgen_NbOfNodes; // set netgen ID
 252
 253     if ( isDegen ) // all nodes on a degen edge get one netgen ID
 254       *(shId_ngId->second) = n_id->second;
 255   }
 256
 257   // set triangles
 258   list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
 259   list< bool >::iterator                 reverse = isReversed.begin();
 260   for ( ; tria != triangles.end(); ++tria, ++reverse )
 261   {
 262     int i = 0;
 263     SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
 264     while ( triangleNodesIt->more() ) {
 265       const SMDS_MeshNode * node =
 266         static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 267       if(myTool->IsMedium(node))
 268         continue;
 269       Netgen_triangle[ *reverse ? 2 - i : i ] = nodeToNetgenID[ node ];
 270       ++i;
 271     }
 272     if ( !hasDegen ||
 273          // ignore degenerated triangles, they have 2 or 3 same ids
 274          (Netgen_triangle[0] != Netgen_triangle[1] &&
 275           Netgen_triangle[0] != Netgen_triangle[2] &&
 276           Netgen_triangle[2] != Netgen_triangle[1] ))
 277       Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 278   }
 279
 280   // -------------------------
 281   // Generate the volume mesh
 282   // -------------------------
 283
 284   Ng_Meshing_Parameters Netgen_param;
 285
 286   Netgen_param.secondorder = Netgen_param2ndOrder;
 287   Netgen_param.fineness = Netgen_paramFine;
 288   Netgen_param.maxh = Netgen_paramSize;
 289
 290   Ng_Result status;
 291
 292   try {
 293     status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
 294   }
 295   catch (...) {
 296     MESSAGE("An exception has been caught during the Volume Mesh Generation ...");
 297     status = NG_VOLUME_FAILURE;
 298   }
 299
 300   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 301
 302   int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
 303
 304   MESSAGE("End of Volume Mesh Generation. status=" << status <<
 305           ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
 306           ", nb tetra: " << Netgen_NbOfTetra);
 307
 308   // -------------------------------------------------------------------
 309   // Feed back the SMESHDS with the generated Nodes and Volume Elements
 310   // -------------------------------------------------------------------
 311
 312   bool isOK = ( status == NG_OK && Netgen_NbOfTetra > 0 );
 313   if ( isOK )
 314   {
 315     // vector of nodes in which node index == netgen ID
 316     vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
 317     // insert old nodes into nodeVec
 318     for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
 319       nodeVec.at( n_id->second ) = n_id->first;
 320     }
 321     // create and insert new nodes into nodeVec
 322     int nodeIndex = Netgen_NbOfNodes + 1;
 323     int shapeID = meshDS->ShapeToIndex( aShape );
 324     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 325     {
 326       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
 327       SMDS_MeshNode * node = meshDS->AddNode(Netgen_point[0],
 328                                              Netgen_point[1],
 329                                              Netgen_point[2]);
 330       meshDS->SetNodeInVolume(node, shapeID);
 331       nodeVec.at(nodeIndex) = node;
 332     }
 333
 334     // create tetrahedrons
 335     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 336     {
 337       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 338       SMDS_MeshVolume * elt = myTool->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
 339                                                  nodeVec.at( Netgen_tetrahedron[1] ),
 340                                                  nodeVec.at( Netgen_tetrahedron[2] ),
 341                                                  nodeVec.at( Netgen_tetrahedron[3] ));
 342       meshDS->SetMeshElementOnShape(elt, shapeID );
 343     }
 344   }
 345
 346   Ng_DeleteMesh(Netgen_mesh);
 347   Ng_Exit();
 348
 349   delete myTool; myTool = 0;
 350
 351   return isOK;
 352 }
 353
 354 //=============================================================================
 355 /*!
 356  *
 357  */
 358 //=============================================================================
 359
 360 ostream & NETGENPlugin_NETGEN_3D::SaveTo(ostream & save)
 361 {
 362   return save;
 363 }
 364
 365 //=============================================================================
 366 /*!
 367  *
 368  */
 369 //=============================================================================
 370
 371 istream & NETGENPlugin_NETGEN_3D::LoadFrom(istream & load)
 372 {
 373   return load;
 374 }
 375
 376 //=============================================================================
 377 /*!
 378  *
 379  */
 380 //=============================================================================
 381
 382 ostream & operator << (ostream & save, NETGENPlugin_NETGEN_3D & hyp)
 383 {
 384   return hyp.SaveTo( save );
 385 }
 386
 387 //=============================================================================
 388 /*!
 389  *
 390  */
 391 //=============================================================================
 392
 393 istream & operator >> (istream & load, NETGENPlugin_NETGEN_3D & hyp)
 394 {
 395   return hyp.LoadFrom( load );
 396 }