src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx

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