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