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