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