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
  21 #include <TopExp.hxx>
  22
  23 #include "utilities.h"
  24
  25 #include <list>
  26 #include <vector>
  27 #include <map>
  28
  29 /*
  30   Netgen include files
  31 */
  32
  33 #include "nglib.h"
  34
  35 //=============================================================================
  36 /*!
  37  *
  38  */
  39 //=============================================================================
  40
  41 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
  42                              SMESH_Gen* gen)
  43   : SMESH_3D_Algo(hypId, studyId, gen)
  44 {
  45   MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
  46   _name = "NETGEN_3D";
  47 //   _shapeType = TopAbs_SOLID;
  48   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  49 //   MESSAGE("_shapeType octal " << oct << _shapeType);
  50   _compatibleHypothesis.push_back("MaxElementVolume");
  51
  52   _maxElementVolume = 0.;
  53
  54   _hypMaxElementVolume = NULL;
  55 }
  56
  57 //=============================================================================
  58 /*!
  59  *
  60  */
  61 //=============================================================================
  62
  63 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
  64 {
  65   MESSAGE("NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D");
  66 }
  67
  68 //=============================================================================
  69 /*!
  70  *
  71  */
  72 //=============================================================================
  73
  74 bool NETGENPlugin_NETGEN_3D::CheckHypothesis
  75                          (SMESH_Mesh& aMesh,
  76                           const TopoDS_Shape& aShape,
  77                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
  78 {
  79   MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
  80
  81   _hypMaxElementVolume = NULL;
  82   _maxElementVolume = DBL_MAX;
  83
  84   list<const SMESHDS_Hypothesis*>::const_iterator itl;
  85   const SMESHDS_Hypothesis* theHyp;
  86
  87   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
  88   int nbHyp = hyps.size();
  89   if (!nbHyp)
  90   {
  91     aStatus = SMESH_Hypothesis::HYP_OK;
  92     //aStatus = SMESH_Hypothesis::HYP_MISSING;
  93     return true;  // can work with no hypothesis
  94   }
  95
  96   itl = hyps.begin();
  97   theHyp = (*itl); // use only the first hypothesis
  98
  99   string hypName = theHyp->GetName();
 100
 101   bool isOk = false;
 102
 103   if (hypName == "MaxElementVolume")
 104   {
 105     _hypMaxElementVolume = static_cast<const StdMeshers_MaxElementVolume*> (theHyp);
 106     ASSERT(_hypMaxElementVolume);
 107     _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
 108     isOk =true;
 109     aStatus = SMESH_Hypothesis::HYP_OK;
 110   }
 111   else
 112     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 113
 114   return isOk;
 115 }
 116
 117 //=============================================================================
 118 /*!
 119  *Here we are going to use the NETGEN mesher
 120  */
 121 //=============================================================================
 122
 123 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
 124                                      const TopoDS_Shape& aShape)
 125 {
 126   MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
 127
 128   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 129
 130   // -------------------------------------------------------------------
 131   // get triangles on aShell and make a map of nodes to Netgen node IDs
 132   // -------------------------------------------------------------------
 133
 134   typedef map< const SMDS_MeshNode*, int> TNodeToIDMap;
 135   TNodeToIDMap nodeToNetgenID;
 136   list< const SMDS_MeshElement* > triangles;
 137   list< bool >                    isReversed; // orientation of triangles
 138
 139   SMESH::Controls::Area areaControl;
 140   SMESH::Controls::TSequenceOfXYZ nodesCoords;
 141
 142   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
 143   {
 144     TopoDS_Shape aShapeFace = exp.Current();
 145     int faceID = meshDS->ShapeToIndex( aShapeFace );
 146     TopoDS_Shape aMeshedFace = meshDS->IndexToShape( faceID );
 147     const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( faceID );
 148     if ( aSubMeshDSFace )
 149     {
 150       bool isRev = ( aShapeFace.Orientation() != aMeshedFace.Orientation() );
 151       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 152       while ( iteratorElem->more() ) // loop on elements on a face
 153       {
 154         // check element
 155         const SMDS_MeshElement* elem = iteratorElem->next();
 156         if ( !elem || elem->NbNodes() != 3 ) {
 157           INFOS( "NETGENPlugin_NETGEN_3D::Compute(), bad mesh");
 158           return false;
 159         }
 160         // keep a triangle
 161         triangles.push_back( elem );
 162         isReversed.push_back( isRev );
 163         // put elem nodes to nodeToNetgenID map
 164         SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
 165         while ( triangleNodesIt->more() ) {
 166           const SMDS_MeshNode * node =
 167             static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 168           nodeToNetgenID.insert( make_pair( node, 0 ));
 169         }
 170 #ifdef _DEBUG_
 171         // check if a trainge is degenerated
 172         areaControl.GetPoints( elem, nodesCoords );
 173         double area = areaControl.GetValue( nodesCoords );
 174         if ( area <= DBL_MIN ) {
 175           MESSAGE( "Warning: Degenerated " << elem );
 176         }
 177 #endif
 178       }
 179     }
 180   }
 181   // ---------------------------------
 182   // Feed the Netgen with surface mesh
 183   // ---------------------------------
 184
 185   int Netgen_NbOfNodes = nodeToNetgenID.size();
 186   int Netgen_param2ndOrder = 0;
 187   double Netgen_paramFine = 1.;
 188   double Netgen_paramSize = _maxElementVolume;
 189
 190   double Netgen_point[3];
 191   int Netgen_triangle[3];
 192   int Netgen_tetrahedron[4];
 193
 194   Ng_Init();
 195
 196   Ng_Mesh * Netgen_mesh = Ng_NewMesh();
 197
 198   // set nodes and remember thier netgen IDs
 199   TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 200   for ( int id = 0; n_id != nodeToNetgenID.end(); ++n_id )
 201   {
 202     const SMDS_MeshNode* node = n_id->first;
 203     Netgen_point [ 0 ] = node->X();
 204     Netgen_point [ 1 ] = node->Y();
 205     Netgen_point [ 2 ] = node->Z();
 206     n_id->second = ++id;
 207     Ng_AddPoint(Netgen_mesh, Netgen_point);
 208   }
 209   // set triangles
 210   list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
 211   list< bool >::iterator                 reverse = isReversed.begin();
 212   for ( ; tria != triangles.end(); ++tria, ++reverse )
 213   {
 214     int i = 0;
 215     SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
 216     while ( triangleNodesIt->more() ) {
 217       const SMDS_MeshNode * node =
 218         static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 219       Netgen_triangle[ *reverse ? 2 - i : i ] = nodeToNetgenID[ node ];
 220       ++i;
 221     }
 222     Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 223   }
 224
 225   // -------------------------
 226   // Generate the volume mesh
 227   // -------------------------
 228
 229   Ng_Meshing_Parameters Netgen_param;
 230
 231   Netgen_param.secondorder = Netgen_param2ndOrder;
 232   Netgen_param.fineness = Netgen_paramFine;
 233   Netgen_param.maxh = Netgen_paramSize;
 234
 235   Ng_Result status;
 236
 237   try {
 238     status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
 239   } catch (...) {
 240     MESSAGE("An exception has been caught during the Volume Mesh Generation ...");
 241     status = NG_VOLUME_FAILURE;
 242   }
 243
 244   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 245
 246   int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
 247
 248   MESSAGE("End of Volume Mesh Generation. status=" << status <<
 249           ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
 250           ", nb tetra: " << Netgen_NbOfTetra);
 251
 252   // -------------------------------------------------------------------
 253   // Feed back the SMESHDS with the generated Nodes and Volume Elements
 254   // -------------------------------------------------------------------
 255
 256   bool isOK = ( status == NG_OK && Netgen_NbOfTetra > 0 );
 257   if ( isOK )
 258   {
 259     // vector of nodes in which node index == netgen ID
 260     vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
 261     // insert old nodes into nodeVec
 262     for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id )
 263       nodeVec[ n_id->second ] = n_id->first;
 264     // create and insert new nodes into nodeVec
 265     int nodeIndex = Netgen_NbOfNodes + 1;
 266     int shapeID = meshDS->ShapeToIndex( aShape );
 267     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 268     {
 269       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
 270       SMDS_MeshNode * node = meshDS->AddNode(Netgen_point[0],
 271                                              Netgen_point[1],
 272                                              Netgen_point[2]);
 273       meshDS->SetNodeInVolume(node, shapeID);
 274       nodeVec[nodeIndex] = node;
 275     }
 276
 277     // create tetrahedrons
 278     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 279     {
 280       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 281       SMDS_MeshVolume * elt = meshDS->AddVolume (nodeVec[ Netgen_tetrahedron[0] ],
 282                                                  nodeVec[ Netgen_tetrahedron[1] ],
 283                                                  nodeVec[ Netgen_tetrahedron[2] ],
 284                                                  nodeVec[ Netgen_tetrahedron[3] ]);
 285       meshDS->SetMeshElementOnShape(elt, shapeID );
 286     }
 287   }
 288
 289   Ng_DeleteMesh(Netgen_mesh);
 290   Ng_Exit();
 291
 292   return isOK;
 293 }
 294
 295 //=============================================================================
 296 /*!
 297  *
 298  */
 299 //=============================================================================
 300
 301 ostream & NETGENPlugin_NETGEN_3D::SaveTo(ostream & save)
 302 {
 303   return save;
 304 }
 305
 306 //=============================================================================
 307 /*!
 308  *
 309  */
 310 //=============================================================================
 311
 312 istream & NETGENPlugin_NETGEN_3D::LoadFrom(istream & load)
 313 {
 314   return load;
 315 }
 316
 317 //=============================================================================
 318 /*!
 319  *
 320  */
 321 //=============================================================================
 322
 323 ostream & operator << (ostream & save, NETGENPlugin_NETGEN_3D & hyp)
 324 {
 325   return hyp.SaveTo( save );
 326 }
 327
 328 //=============================================================================
 329 /*!
 330  *
 331  */
 332 //=============================================================================
 333
 334 istream & operator >> (istream & load, NETGENPlugin_NETGEN_3D & hyp)
 335 {
 336   return hyp.LoadFrom( load );
 337 }