src/NETGENPlugin/NETGENPlugin_Mesher.cxx

   1 //  NETGENPlugin : C++ implementation
   2 //
   3 //  Copyright (C) 2006  OPEN CASCADE, CEA/DEN, EDF R&D
   4 //
   5 //  This library is free software; you can redistribute it and/or
   6 //  modify it under the terms of the GNU Lesser General Public
   7 //  License as published by the Free Software Foundation; either
   8 //  version 2.1 of the License.
   9 //
  10 //  This library is distributed in the hope that it will be useful,
  11 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 //  Lesser General Public License for more details.
  14 //
  15 //  You should have received a copy of the GNU Lesser General Public
  16 //  License along with this library; if not, write to the Free Software
  17 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 //
  19 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 //
  21 //
  22 // File      : NETGENPlugin_Mesher.cxx
  23 // Author    : Michael Sazonov (OCN)
  24 // Date      : 31/03/2006
  25 // Project   : SALOME
  26 // $Header$
  27 //=============================================================================
  28
  29 #include "NETGENPlugin_Mesher.hxx"
  30 #include "NETGENPlugin_Hypothesis_2D.hxx"
  31
  32 #include <SMESH_Mesh.hxx>
  33 #include <SMESH_Comment.hxx>
  34 #include <SMESH_ComputeError.hxx>
  35 #include <SMESH_subMesh.hxx>
  36 #include <SMESHDS_Mesh.hxx>
  37 #include <SMDS_MeshElement.hxx>
  38 #include <SMDS_MeshNode.hxx>
  39 #include <utilities.h>
  40
  41 #include <vector>
  42
  43 #include <BRep_Tool.hxx>
  44 #include <TopExp.hxx>
  45 #include <TopExp_Explorer.hxx>
  46 #include <TopoDS.hxx>
  47 #include <NCollection_Map.hxx>
  48 #include <OSD_Path.hxx>
  49 #include <OSD_File.hxx>
  50 #include <TCollection_AsciiString.hxx>
  51
  52 // Netgen include files
  53 namespace nglib {
  54 #include <nglib.h>
  55 }
  56 #define OCCGEOMETRY
  57 #include <occgeom.hpp>
  58 #include <meshing.hpp>
  59 //#include <ngexception.hpp>
  60 namespace netgen {
  61   extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
  62   extern MeshingParameters mparam;
  63 }
  64
  65 using namespace std;
  66
  67 //=============================================================================
  68 /*!
  69  *
  70  */
  71 //=============================================================================
  72
  73 NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh,
  74                                           const TopoDS_Shape& aShape,
  75                                           const bool isVolume)
  76   : _mesh  (mesh),
  77     _shape   (aShape),
  78     _isVolume(isVolume),
  79     _optimize(true)
  80 {
  81 #ifdef WNT
  82   netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
  83 #else
  84   netgen::MeshingParameters& mparams = netgen::mparam;
  85 #endif
  86   // Initialize global NETGEN parameters by default values:
  87   // maximal mesh edge size
  88   mparams.maxh = NETGENPlugin_Hypothesis::GetDefaultMaxSize();
  89   // minimal number of segments per edge
  90   mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
  91   // rate of growth of size between elements
  92   mparams.grading = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
  93   // safety factor for curvatures (elements per radius)
  94   mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
  95   // create elements of second order
  96   mparams.secondorder = NETGENPlugin_Hypothesis::GetDefaultSecondOrder() ? 1 : 0;
  97   // quad-dominated surface meshing
  98   if (_isVolume)
  99     mparams.quad = 0;
 100   else
 101     mparams.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed() ? 1 : 0;
 102 }
 103
 104 //=============================================================================
 105 /*!
 106  * Pass parameters to NETGEN
 107  */
 108 //=============================================================================
 109 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
 110 {
 111   if (hyp)
 112   {
 113 #ifdef WNT
 114     netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
 115 #else
 116     netgen::MeshingParameters& mparams = netgen::mparam;
 117 #endif
 118     // Initialize global NETGEN parameters:
 119     // maximal mesh segment size
 120     mparams.maxh = hyp->GetMaxSize();
 121     // minimal number of segments per edge
 122     mparams.segmentsperedge = hyp->GetNbSegPerEdge();
 123     // rate of growth of size between elements
 124     mparams.grading = hyp->GetGrowthRate();
 125     // safety factor for curvatures (elements per radius)
 126     mparams.curvaturesafety = hyp->GetNbSegPerRadius();
 127     // create elements of second order
 128     mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0;
 129     // quad-dominated surface meshing
 130     // only triangles are allowed for volumic mesh
 131     if (!_isVolume)
 132       mparams.quad = static_cast<const NETGENPlugin_Hypothesis_2D*>
 133         (hyp)->GetQuadAllowed() ? 1 : 0;
 134     _optimize = hyp->GetOptimize();
 135   }
 136 }
 137
 138 //=============================================================================
 139 /*!
 140  *  Link - a pair of integer numbers
 141  */
 142 //=============================================================================
 143 struct Link
 144 {
 145   int n1, n2;
 146   Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
 147   Link() : n1(0), n2(0) {}
 148 };
 149
 150 int HashCode(const Link& aLink, int aLimit)
 151 {
 152   return HashCode(aLink.n1 + aLink.n2, aLimit);
 153 }
 154
 155 Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
 156 {
 157   return (aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ||
 158           aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1);
 159 }
 160
 161 //================================================================================
 162 /*!
 163  * \brief Initialize netgen::OCCGeometry with OCCT shape
 164  */
 165 //================================================================================
 166
 167 void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry& occgeo,
 168                                              const TopoDS_Shape&  shape)
 169 {
 170   occgeo.shape = shape;
 171   occgeo.changed = 1;
 172   occgeo.BuildFMap();
 173
 174   BRepTools::Clean (shape);
 175   BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (shape, 0.01, true);
 176   Bnd_Box bb;
 177   BRepBndLib::Add (shape, bb);
 178   double x1,y1,z1,x2,y2,z2;
 179   bb.Get (x1,y1,z1,x2,y2,z2);
 180   MESSAGE("shape bounding box:\n" <<
 181           "(" << x1 << " " << y1 << " " << z1 << ") " <<
 182           "(" << x2 << " " << y2 << " " << z2 << ")");
 183   netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
 184   netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
 185   occgeo.boundingbox = netgen::Box<3> (p1,p2);
 186 }
 187
 188 //=============================================================================
 189 /*!
 190  * Here we are going to use the NETGEN mesher
 191  */
 192 //=============================================================================
 193 bool NETGENPlugin_Mesher::Compute()
 194 {
 195 #ifdef WNT
 196   netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
 197 #else
 198   netgen::MeshingParameters& mparams = netgen::mparam;
 199 #endif
 200   MESSAGE("Compute with:\n"
 201           " max size = " << mparams.maxh << "\n"
 202           " segments per edge = " << mparams.segmentsperedge);
 203   MESSAGE("\n"
 204           " growth rate = " << mparams.grading << "\n"
 205           " elements per radius = " << mparams.curvaturesafety << "\n"
 206           " second order = " << mparams.secondorder << "\n"
 207           " quad allowed = " << mparams.quad);
 208
 209   SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
 210   nglib::Ng_Init();
 211
 212   // -------------------------
 213   // Prepare OCC geometry
 214   // -------------------------
 215
 216   netgen::OCCGeometry occgeo;
 217   PrepareOCCgeometry( occgeo, _shape );
 218
 219   // -------------------------
 220   // Generate the mesh
 221   // -------------------------
 222
 223   netgen::Mesh *ngMesh = NULL;
 224   // we start always with ANALYSE,
 225   // but end depending on _optimize and _isVolume
 226   int startWith = netgen::MESHCONST_ANALYSE;
 227   int endWith = (_optimize
 228                  ? (_isVolume ? netgen::MESHCONST_OPTVOLUME : netgen::MESHCONST_OPTSURFACE)
 229                  : netgen::MESHCONST_MESHSURFACE);
 230   char *optstr = 0;
 231
 232   int err = 0;
 233   SMESH_Comment comment;
 234   try
 235   {
 236     err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
 237     if (err) comment << "Error in netgen::OCCGenerateMesh()";
 238     if (!err && !_optimize)
 239     {
 240       // we have got surface mesh only, so generate volume mesh
 241       startWith = endWith = netgen::MESHCONST_MESHVOLUME;
 242       err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
 243       if (err) comment << "Error in netgen::OCCGenerateMesh()";
 244     }
 245     if (!err && mparams.secondorder > 0)
 246     {
 247       netgen::OCCRefinementSurfaces ref (occgeo);
 248       ref.MakeSecondOrder (*ngMesh);
 249     }
 250   }
 251   catch (netgen::NgException exc)
 252   {
 253     error->myName = err = COMPERR_ALGO_FAILED;
 254     comment << exc.What();
 255   }
 256
 257   int nbNod = ngMesh->GetNP();
 258   int nbSeg = ngMesh->GetNSeg();
 259   int nbFac = ngMesh->GetNSE();
 260   int nbVol = ngMesh->GetNE();
 261
 262   MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
 263           ", nb nodes: " << nbNod <<
 264           ", nb segments: " << nbSeg <<
 265           ", nb faces: " << nbFac <<
 266           ", nb volumes: " << nbVol);
 267
 268   // -----------------------------------------------------------
 269   // Feed back the SMESHDS with the generated Nodes and Elements
 270   // -----------------------------------------------------------
 271
 272   SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
 273   bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
 274   if ( true /*isOK*/ ) // get whatever built
 275   {
 276     // vector of nodes in which node index == netgen ID
 277     vector< SMDS_MeshNode* > nodeVec ( nbNod + 1 );
 278     // map of nodes assigned to submeshes
 279     NCollection_Map<int> pindMap;
 280     // create and insert nodes into nodeVec
 281     int i;
 282     for (i = 1; i <= nbNod /*&& isOK*/; ++i )
 283     {
 284       const netgen::MeshPoint& ngPoint = ngMesh->Point(i);
 285       SMDS_MeshNode* node = NULL;
 286       bool newNodeOnVertex = false;
 287       TopoDS_Vertex aVert;
 288       if (i <= occgeo.vmap.Extent())
 289       {
 290         // point on vertex
 291         aVert = TopoDS::Vertex(occgeo.vmap(i));
 292         SMESHDS_SubMesh * submesh = meshDS->MeshElements(aVert);
 293         if (submesh)
 294         {
 295           SMDS_NodeIteratorPtr it = submesh->GetNodes();
 296           if (it->more())
 297           {
 298             node = const_cast<SMDS_MeshNode*> (it->next());
 299             pindMap.Add(i);
 300           }
 301         }
 302         if (!node)
 303           newNodeOnVertex = true;
 304       }
 305       if (!node)
 306         node = meshDS->AddNode(ngPoint.X(), ngPoint.Y(), ngPoint.Z());
 307       if (!node)
 308       {
 309         MESSAGE("Cannot create a mesh node");
 310         if ( !comment.size() ) comment << "Cannot create a mesh node";
 311         nbSeg = nbFac = nbVol = isOK = 0;
 312         break;
 313       }
 314       nodeVec.at(i) = node;
 315       if (newNodeOnVertex)
 316       {
 317         // point on vertex
 318         meshDS->SetNodeOnVertex(node, aVert);
 319         pindMap.Add(i);
 320       }
 321     }
 322
 323     // create mesh segments along geometric edges
 324     NCollection_Map<Link> linkMap;
 325     for (i = 1; i <= nbSeg/* && isOK*/; ++i )
 326     {
 327       const netgen::Segment& seg = ngMesh->LineSegment(i);
 328       Link link(seg.p1, seg.p2);
 329       if (linkMap.Contains(link))
 330         continue;
 331       linkMap.Add(link);
 332       TopoDS_Edge aEdge;
 333       int pinds[3] = { seg.p1, seg.p2, seg.pmid };
 334       int nbp = 0;
 335       double param2 = 0;
 336       for (int j=0; j < 3; ++j)
 337       {
 338         int pind = pinds[j];
 339         if (pind <= 0) continue;
 340         ++nbp;
 341         double param;
 342         if (j < 2)
 343         {
 344           if (aEdge.IsNull())
 345           {
 346             int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
 347             if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
 348               aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
 349           }
 350           param = seg.epgeominfo[j].dist;
 351           param2 += param;
 352         }
 353         else
 354           param = param2 * 0.5;
 355         if (pindMap.Contains(pind))
 356           continue;
 357         if (!aEdge.IsNull())
 358         {
 359           meshDS->SetNodeOnEdge(nodeVec.at(pind), aEdge, param);
 360           pindMap.Add(pind);
 361         }
 362       }
 363       SMDS_MeshEdge* edge;
 364       if (nbp < 3) // second order ?
 365         edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]));
 366       else
 367         edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]),
 368                                 nodeVec.at(pinds[2]));
 369       if (!edge)
 370       {
 371         if ( !comment.size() ) comment << "Cannot create a mesh edge";
 372         MESSAGE("Cannot create a mesh edge");
 373         nbSeg = nbFac = nbVol = isOK = 0;
 374         break;
 375       }
 376       if (!aEdge.IsNull())
 377         meshDS->SetMeshElementOnShape(edge, aEdge);
 378     }
 379
 380     // create mesh faces along geometric faces
 381     for (i = 1; i <= nbFac/* && isOK*/; ++i )
 382     {
 383       const netgen::Element2d& elem = ngMesh->SurfaceElement(i);
 384       int aGeomFaceInd = elem.GetIndex();
 385       TopoDS_Face aFace;
 386       if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
 387         aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
 388       vector<SMDS_MeshNode*> nodes;
 389       for (int j=1; j <= elem.GetNP(); ++j)
 390       {
 391         int pind = elem.PNum(j);
 392         SMDS_MeshNode* node = nodeVec.at(pind);
 393         nodes.push_back(node);
 394         if (pindMap.Contains(pind))
 395           continue;
 396         if (!aFace.IsNull())
 397         {
 398           const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
 399           meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
 400           pindMap.Add(pind);
 401         }
 402       }
 403       SMDS_MeshFace* face = NULL;
 404       switch (elem.GetType())
 405       {
 406       case netgen::TRIG:
 407         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
 408         break;
 409       case netgen::QUAD:
 410         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
 411         break;
 412       case netgen::TRIG6:
 413         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
 414         break;
 415       case netgen::QUAD8:
 416         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
 417                                nodes[4],nodes[7],nodes[5],nodes[6]);
 418         break;
 419       default:
 420         MESSAGE("NETGEN created a face of unexpected type, ignoring");
 421         continue;
 422       }
 423       if (!face)
 424       {
 425         if ( !comment.size() ) comment << "Cannot create a mesh face";
 426         MESSAGE("Cannot create a mesh face");
 427         nbSeg = nbFac = nbVol = isOK = 0;
 428         break;
 429       }
 430       if (!aFace.IsNull())
 431         meshDS->SetMeshElementOnShape(face, aFace);
 432     }
 433
 434     // create tetrahedra
 435     for (i = 1; i <= nbVol/* && isOK*/; ++i)
 436     {
 437       const netgen::Element& elem = ngMesh->VolumeElement(i);
 438       int aSolidInd = elem.GetIndex();
 439       TopoDS_Solid aSolid;
 440       if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
 441         aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
 442       vector<SMDS_MeshNode*> nodes;
 443       for (int j=1; j <= elem.GetNP(); ++j)
 444       {
 445         int pind = elem.PNum(j);
 446         SMDS_MeshNode* node = nodeVec.at(pind);
 447         nodes.push_back(node);
 448         if (pindMap.Contains(pind))
 449           continue;
 450         if (!aSolid.IsNull())
 451         {
 452           // point in solid
 453           meshDS->SetNodeInVolume(node, aSolid);
 454           pindMap.Add(pind);
 455         }
 456       }
 457       SMDS_MeshVolume* vol = NULL;
 458       switch (elem.GetType())
 459       {
 460       case netgen::TET:
 461         vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]);
 462         break;
 463       case netgen::TET10:
 464         vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3],
 465                                 nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]);
 466         break;
 467       default:
 468         MESSAGE("NETGEN created a volume of unexpected type, ignoring");
 469         continue;
 470       }
 471       if (!vol)
 472       {
 473         if ( !comment.size() ) comment << "Cannot create a mesh volume";
 474         MESSAGE("Cannot create a mesh volume");
 475         nbSeg = nbFac = nbVol = isOK = 0;
 476         break;
 477       }
 478       if (!aSolid.IsNull())
 479         meshDS->SetMeshElementOnShape(vol, aSolid);
 480     }
 481   }
 482
 483   if ( error->IsOK() && ( !isOK || comment.size() > 0 ))
 484     error->myName = COMPERR_ALGO_FAILED;
 485   if ( !comment.empty() )
 486     error->myComment = comment;
 487
 488   // set bad compute error to subshapes of all failed subshapes shapes
 489   if ( !error->IsOK() && err )
 490   {
 491     for (int i = 1; i <= occgeo.fmap.Extent(); i++) {
 492       int status = occgeo.facemeshstatus[i-1];
 493       if (status == 1 ) continue;
 494       if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.fmap( i ))) {
 495         SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 496         if ( !smError || smError->IsOK() ) {
 497           if ( status == -1 )
 498             smError.reset( new SMESH_ComputeError( error->myName, error->myComment ));
 499           else
 500             smError.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, "Ignored" ));
 501         }
 502       }
 503     }
 504   }
 505
 506   nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);
 507   nglib::Ng_Exit();
 508
 509   RemoveTmpFiles();
 510
 511   return error->IsOK();
 512 }
 513
 514 //================================================================================
 515 /*!
 516  * \brief Remove "test.out" and "problemfaces" files in current directory
 517  */
 518 //================================================================================
 519
 520 void NETGENPlugin_Mesher::RemoveTmpFiles()
 521 {
 522   TCollection_AsciiString str("test.out");
 523   OSD_Path path1( str );
 524   OSD_File file1( path1 );
 525   file1.Remove();
 526   str = "problemfaces";
 527   OSD_Path path2( str );
 528   OSD_File file2( path2 );
 529   file2.Remove();
 530 }