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