src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 // File      : NETGENPlugin_NETGEN_2D_ONLY.cxx
  21 // Author    : Edward AGAPOV (OCC)
  22 // Project   : SALOME
  23 //
  24 #include "NETGENPlugin_NETGEN_2D_ONLY.hxx"
  25
  26 #include "NETGENPlugin_Mesher.hxx"
  27 #include "NETGENPlugin_Hypothesis_2D.hxx"
  28
  29 #include "SMDS_MeshElement.hxx"
  30 #include "SMDS_MeshNode.hxx"
  31 #include "SMESHDS_Mesh.hxx"
  32 #include "SMESH_Comment.hxx"
  33 #include "SMESH_Gen.hxx"
  34 #include "SMESH_Mesh.hxx"
  35 #include "SMESH_MesherHelper.hxx"
  36 #include "StdMeshers_FaceSide.hxx"
  37 #include "StdMeshers_MaxElementArea.hxx"
  38 #include "StdMeshers_LengthFromEdges.hxx"
  39 #include "StdMeshers_QuadranglePreference.hxx"
  40
  41 #include <Precision.hxx>
  42 #include <Standard_ErrorHandler.hxx>
  43 #include <Standard_Failure.hxx>
  44
  45 #include "utilities.h"
  46
  47 #include <list>
  48 #include <vector>
  49 #include <limits>
  50
  51 /*
  52   Netgen include files
  53 */
  54 namespace nglib {
  55 #include <nglib.h>
  56 }
  57 #ifndef OCCGEOMETRY
  58 #define OCCGEOMETRY
  59 #endif
  60 #include <occgeom.hpp>
  61 #include <meshing.hpp>
  62 //#include <meshtype.hpp>
  63 namespace netgen {
  64   extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
  65   extern MeshingParameters mparam;
  66 }
  67
  68 using namespace std;
  69 using namespace netgen;
  70 using namespace nglib;
  71
  72 //#define DUMP_SEGMENTS
  73
  74 //=============================================================================
  75 /*!
  76  *
  77  */
  78 //=============================================================================
  79
  80 NETGENPlugin_NETGEN_2D_ONLY::NETGENPlugin_NETGEN_2D_ONLY(int hypId, int studyId,
  81                                                          SMESH_Gen* gen)
  82   : SMESH_2D_Algo(hypId, studyId, gen)
  83 {
  84   MESSAGE("NETGENPlugin_NETGEN_2D_ONLY::NETGENPlugin_NETGEN_2D_ONLY");
  85   _name = "NETGEN_2D_ONLY";
  86
  87   _shapeType = (1 << TopAbs_FACE);// 1 bit /shape type
  88
  89   _compatibleHypothesis.push_back("MaxElementArea");
  90   _compatibleHypothesis.push_back("LengthFromEdges");
  91   _compatibleHypothesis.push_back("QuadranglePreference");
  92   _compatibleHypothesis.push_back("NETGEN_Parameters_2D");
  93
  94   _hypMaxElementArea = 0;
  95   _hypLengthFromEdges = 0;
  96   _hypQuadranglePreference = 0;
  97   _hypParameters = 0;
  98 }
  99
 100 //=============================================================================
 101 /*!
 102  *
 103  */
 104 //=============================================================================
 105
 106 NETGENPlugin_NETGEN_2D_ONLY::~NETGENPlugin_NETGEN_2D_ONLY()
 107 {
 108   MESSAGE("NETGENPlugin_NETGEN_2D_ONLY::~NETGENPlugin_NETGEN_2D_ONLY");
 109 }
 110
 111 //=============================================================================
 112 /*!
 113  *
 114  */
 115 //=============================================================================
 116
 117 bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh&         aMesh,
 118                                                    const TopoDS_Shape& aShape,
 119                                                    Hypothesis_Status&  aStatus)
 120 {
 121   _hypMaxElementArea = 0;
 122   _hypLengthFromEdges = 0;
 123   _hypQuadranglePreference = 0;
 124
 125   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape, false);
 126
 127   if (hyps.empty())
 128   {
 129     aStatus = HYP_OK; //SMESH_Hypothesis::HYP_MISSING;
 130     return true;  // (PAL13464) can work with no hypothesis, LengthFromEdges is default one
 131   }
 132
 133   aStatus = HYP_MISSING;
 134
 135   list<const SMESHDS_Hypothesis*>::const_iterator ith;
 136   for (ith = hyps.begin(); ith != hyps.end(); ++ith )
 137   {
 138     const SMESHDS_Hypothesis* hyp = (*ith);
 139
 140     string hypName = hyp->GetName();
 141
 142     if      ( hypName == "MaxElementArea")
 143       _hypMaxElementArea = static_cast<const StdMeshers_MaxElementArea*> (hyp);
 144     else if ( hypName == "LengthFromEdges" )
 145       _hypLengthFromEdges = static_cast<const StdMeshers_LengthFromEdges*> (hyp);
 146     else if ( hypName == "QuadranglePreference" )
 147       _hypQuadranglePreference = static_cast<const StdMeshers_QuadranglePreference*>(hyp);
 148     else if ( hypName == "NETGEN_Parameters_2D" )
 149       _hypParameters = static_cast<const NETGENPlugin_Hypothesis_2D*>(hyp);
 150     else {
 151       aStatus = HYP_INCOMPATIBLE;
 152       return false;
 153     }
 154   }
 155
 156   int nbHyps = bool(_hypMaxElementArea) + bool(_hypLengthFromEdges) + bool(_hypParameters );
 157   if ( nbHyps > 1 )
 158     aStatus = HYP_CONCURENT;
 159   else if ( nbHyps == 1)
 160     aStatus = HYP_OK;
 161
 162   return ( aStatus == HYP_OK );
 163 }
 164
 165 //================================================================================
 166 /*!
 167  * \brief Fill netgen mesh with segments
 168   * \retval SMESH_ComputeErrorPtr - error description
 169  */
 170 //================================================================================
 171
 172 static TError addSegmentsToMesh(netgen::Mesh&                    ngMesh,
 173                                 OCCGeometry&                     geom,
 174                                 const TSideVector&               wires,
 175                                 SMESH_MesherHelper&              helper,
 176                                 vector< const SMDS_MeshNode* > & nodeVec)
 177 {
 178   // ----------------------------
 179   // Check wires and count nodes
 180   // ----------------------------
 181   int nbNodes = 0;
 182   double totalLength = 0;
 183   for ( int iW = 0; iW < wires.size(); ++iW )
 184   {
 185     StdMeshers_FaceSidePtr wire = wires[ iW ];
 186     if ( wire->MissVertexNode() )
 187     {
 188       // Commented for issue 0020960. It worked for the case, let's wait for case where it doesn't.
 189       // It seems that there is no reason for this limitation
 190 //       return TError
 191 //         (new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH, "Missing nodes on vertices"));
 192       if (getenv("USER") && string("eap")==getenv("USER"))
 193         cout << "Warning: NETGENPlugin_NETGEN_2D_ONLY : try to work with missing nodes on vertices"<<endl;
 194     }
 195     const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
 196     if ( uvPtVec.size() != wire->NbPoints() )
 197       return TError
 198         (new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,
 199                                 SMESH_Comment("Unexpected nb of points on wire ") << iW
 200                                 << ": " << uvPtVec.size()<<" != "<<wire->NbPoints()));
 201     nbNodes += wire->NbPoints();
 202     totalLength += wire->Length();
 203   }
 204   nodeVec.reserve( nbNodes );
 205
 206   // -----------------
 207   // Fill netgen mesh
 208   // -----------------
 209
 210 //   netgen::Box<3> bb = geom.GetBoundingBox();
 211 //   bb.Increase (bb.Diam()/10);
 212 //   ngMesh.SetLocalH (bb.PMin(), bb.PMax(), 0.5); // set grading
 213
 214   // map for nodes on vertices since they can be shared between wires
 215   // ( issue 0020676, face_int_box.brep)
 216   map<const SMDS_MeshNode*, int > node2ngID;
 217
 218   const int faceID = 1, solidID = 0;
 219   if ( ngMesh.GetNFD() < 1 )
 220     ngMesh.AddFaceDescriptor (FaceDescriptor(faceID, solidID, solidID, 0));
 221
 222   for ( int iW = 0; iW < wires.size(); ++iW )
 223   {
 224     StdMeshers_FaceSidePtr wire = wires[ iW ];
 225     const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
 226     const int nbSegments = wire->NbPoints() - 1;
 227
 228     // compute length of every segment
 229     vector<double> segLen( nbSegments );
 230     for ( int i = 0; i < nbSegments; ++i )
 231       segLen[i] = SMESH_TNodeXYZ( uvPtVec[ i ].node ).Distance( uvPtVec[ i+1 ].node );
 232
 233     int edgeID = 1, posID = -2;
 234     bool isInternalWire = false;
 235     for ( int i = 0; i < nbSegments; ++i ) // loop on segments
 236     {
 237       // Add the first point of a segment
 238       const SMDS_MeshNode * n = uvPtVec[ i ].node;
 239       const int posShapeID = n->getshapeId();
 240       bool onVertex = ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX );
 241
 242       // skip nodes on degenerated edges
 243       if ( helper.IsDegenShape( posShapeID ) &&
 244            helper.IsDegenShape( uvPtVec[ i+1 ].node->getshapeId() ))
 245         continue;
 246
 247       int ngID1 = ngMesh.GetNP() + 1, ngID2 = ngID1+1;
 248       if ( onVertex )
 249         ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second;
 250       if ( ngID1 > ngMesh.GetNP() )
 251       {
 252         MeshPoint mp( Point<3> (n->X(), n->Y(), n->Z()) );
 253         ngMesh.AddPoint ( mp, 1, EDGEPOINT );
 254         nodeVec.push_back( n );
 255       }
 256       else
 257       {
 258         ngID2 = ngMesh.GetNP() + 1;
 259         if ( i > 0 ) // prev segment belongs to same wire
 260         {
 261           Segment& prevSeg = ngMesh.LineSegment( ngMesh.GetNSeg() );
 262           prevSeg[1] = ngID1;
 263         }
 264       }
 265
 266       // Add the segment
 267       Segment seg;
 268
 269       seg[0] = ngID1;  // ng node id
 270       seg[1] = ngID2;  // ng node id
 271
 272       seg.edgenr = ngMesh.GetNSeg() + 1;// segment id
 273       seg.si = faceID;                  // = geom.fmap.FindIndex (face);
 274
 275       for ( int iEnd = 0; iEnd < 2; ++iEnd)
 276       {
 277         const UVPtStruct& pnt = uvPtVec[ i + iEnd ];
 278
 279         seg.epgeominfo[ iEnd ].dist = pnt.param; // param on curve
 280         seg.epgeominfo[ iEnd ].u    = pnt.u;
 281         seg.epgeominfo[ iEnd ].v    = pnt.v;
 282
 283         // find out edge id and node parameter on edge
 284         onVertex = ( pnt.node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX );
 285         if ( onVertex || posShapeID != posID )
 286         {
 287           // get edge id
 288           double normParam = pnt.normParam;
 289           if ( onVertex )
 290             normParam = 0.5 * ( uvPtVec[ i ].normParam + uvPtVec[ i+1 ].normParam );
 291           const TopoDS_Edge& edge = wire->Edge( wire->EdgeIndex( normParam ));
 292           edgeID = geom.emap.FindIndex( edge );
 293           posID  = posShapeID;
 294           isInternalWire = ( edge.Orientation() == TopAbs_INTERNAL );
 295           // if ( onVertex ) // param on curve is different on each of two edges
 296           //   seg.epgeominfo[ iEnd ].dist = helper.GetNodeU( edge, pnt.node );
 297         }
 298         seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
 299       }
 300
 301       ngMesh.AddSegment (seg);
 302       {
 303         // restrict size of elements near the segment
 304         SMESH_TNodeXYZ np1( n ), np2( uvPtVec[ i+1 ].node );
 305         // get an average size of adjacent segments to avoid sharp change of
 306         // element size (regression on issue 0020452, note 0010898)
 307         int iPrev = SMESH_MesherHelper::WrapIndex( i-1, nbSegments );
 308         int iNext = SMESH_MesherHelper::WrapIndex( i+1, nbSegments );
 309         double avgH = ( segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ]) / 3;
 310
 311         NETGENPlugin_Mesher::RestrictLocalSize( ngMesh, 0.5*(np1+np2), avgH );
 312       }
 313 #ifdef DUMP_SEGMENTS
 314         cout << "Segment: " << seg.edgenr << endl
 315            << "\tp1: " << seg[0] << endl
 316            << "\tp2: " << seg[1] << endl
 317            << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
 318            << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
 319            << "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
 320            << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
 321            << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl
 322            << "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
 323 #endif
 324       if ( isInternalWire )
 325       {
 326         swap (seg[0], seg[1]);
 327         swap( seg.epgeominfo[0], seg.epgeominfo[1] );
 328         seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
 329         ngMesh.AddSegment (seg);
 330 #ifdef DUMP_SEGMENTS
 331         cout << "Segment: " << seg.edgenr << endl << "\tis REVRESE of the previous one" << endl;
 332 #endif
 333       }
 334     } // loop on segments on a wire
 335
 336     // close chain of segments
 337     if ( nbSegments > 0 )
 338     {
 339       Segment& lastSeg = ngMesh.LineSegment( ngMesh.GetNSeg() - int( isInternalWire));
 340       const SMDS_MeshNode * lastNode = uvPtVec.back().node;
 341       lastSeg[1] = node2ngID.insert( make_pair( lastNode, lastSeg[1] )).first->second;
 342       if ( lastSeg[1] > ngMesh.GetNP() )
 343       {
 344         MeshPoint mp( Point<3> (lastNode->X(), lastNode->Y(), lastNode->Z()) );
 345         ngMesh.AddPoint ( mp, 1, EDGEPOINT );
 346         nodeVec.push_back( lastNode );
 347       }
 348       if ( isInternalWire )
 349       {
 350         Segment& realLastSeg = ngMesh.LineSegment( ngMesh.GetNSeg() );
 351         realLastSeg[0] = lastSeg[1];
 352       }
 353     }
 354
 355   } // loop on wires of a face
 356
 357   // add a segment instead of internal vertex
 358   NETGENPlugin_Internals intShapes( *helper.GetMesh(), helper.GetSubShape(), /*is3D=*/false );
 359   NETGENPlugin_Mesher::addIntVerticesInFaces( geom, ngMesh, nodeVec, intShapes );
 360
 361   ngMesh.CalcSurfacesOfNode();
 362
 363   return TError();
 364 }
 365
 366 //=============================================================================
 367 /*!
 368  *Here we are going to use the NETGEN mesher
 369  */
 370 //=============================================================================
 371
 372 bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
 373                                           const TopoDS_Shape& aShape)
 374 {
 375 #ifdef WITH_SMESH_CANCEL_COMPUTE
 376   netgen::multithread.terminate = 0;
 377 #endif
 378   MESSAGE("NETGENPlugin_NETGEN_2D_ONLY::Compute()");
 379
 380   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 381   int faceID = meshDS->ShapeToIndex( aShape );
 382
 383   SMESH_MesherHelper helper(aMesh);
 384   _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
 385   helper.SetElementsOnShape( true );
 386   const bool ignoreMediumNodes = _quadraticMesh;
 387
 388   // ------------------------
 389   // get all edges of a face
 390   // ------------------------
 391   const TopoDS_Face F = TopoDS::Face( aShape.Oriented( TopAbs_FORWARD ));
 392   TError problem;
 393   TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, aMesh, ignoreMediumNodes, problem );
 394   if ( problem && !problem->IsOK() )
 395     return error( problem );
 396   int nbWires = wires.size();
 397   if ( nbWires == 0 )
 398     return error( "Problem in StdMeshers_FaceSide::GetFaceWires()");
 399   if ( wires[0]->NbSegments() < 3 ) // ex: a circle with 2 segments
 400     return error(COMPERR_BAD_INPUT_MESH,
 401                  SMESH_Comment("Too few segments: ")<<wires[0]->NbSegments());
 402
 403   // --------------------
 404   // compute edge length
 405   // --------------------
 406
 407   NETGENPlugin_Mesher aMesher( &aMesh, aShape, /*isVolume=*/false);
 408   netgen::OCCGeometry occgeo;
 409   aMesher.PrepareOCCgeometry( occgeo, F, aMesh );
 410   occgeo.fmap.Clear(); // face can be reversed, which is wrong in this case (issue 19978)
 411   occgeo.fmap.Add( F );
 412
 413   if ( _hypParameters )
 414   {
 415     aMesher.SetParameters(_hypParameters);
 416   }
 417   else
 418   {
 419     double edgeLength = 0;
 420     if (_hypLengthFromEdges || (!_hypLengthFromEdges && !_hypMaxElementArea))
 421     {
 422       int nbSegments = 0;
 423       for ( int iW = 0; iW < nbWires; ++iW )
 424       {
 425         edgeLength += wires[ iW ]->Length();
 426         nbSegments += wires[ iW ]->NbSegments();
 427       }
 428       if ( nbSegments )
 429         edgeLength /= nbSegments;
 430     }
 431     if ( _hypMaxElementArea )
 432     {
 433       double maxArea = _hypMaxElementArea->GetMaxArea();
 434       edgeLength = sqrt(2. * maxArea/sqrt(3.0));
 435     }
 436     if ( edgeLength < DBL_MIN )
 437       edgeLength = occgeo.GetBoundingBox().Diam();
 438
 439     netgen::mparam.maxh = edgeLength;
 440     netgen::mparam.minh = aMesher.GetDefaultMinSize( aShape, netgen::mparam.maxh );
 441     netgen::mparam.quad = _hypQuadranglePreference ? 1 : 0;
 442     netgen::mparam.grading = 0.7; // very coarse mesh by default
 443   }
 444 #ifdef NETGEN_NEW
 445   occgeo.face_maxh = netgen::mparam.maxh;
 446 #endif
 447
 448   // -------------------------
 449   // Make input netgen mesh
 450   // -------------------------
 451
 452   NETGENPlugin_NetgenLibWrapper ngLib;
 453   netgen::Mesh * ngMesh = (netgen::Mesh*) ngLib._ngMesh;
 454
 455   Box<3> bb = occgeo.GetBoundingBox();
 456   bb.Increase (bb.Diam()/10);
 457   ngMesh->SetLocalH (bb.PMin(), bb.PMax(), netgen::mparam.grading);
 458   ngMesh->SetGlobalH (netgen::mparam.maxh);
 459
 460   vector< const SMDS_MeshNode* > nodeVec;
 461   problem = addSegmentsToMesh( *ngMesh, occgeo, wires, helper, nodeVec );
 462   if ( problem && !problem->IsOK() )
 463     return error( problem );
 464
 465   // -------------------------
 466   // Generate surface mesh
 467   // -------------------------
 468
 469   char *optstr = 0;
 470   int startWith = MESHCONST_MESHSURFACE;
 471   int endWith   = MESHCONST_OPTSURFACE;
 472   int err = 1;
 473
 474   try {
 475 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 476     OCC_CATCH_SIGNALS;
 477 #endif
 478     err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
 479 #ifdef WITH_SMESH_CANCEL_COMPUTE
 480     if(netgen::multithread.terminate)
 481       return false;
 482 #endif
 483     if ( err )
 484       error(SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task);
 485   }
 486   catch (Standard_Failure& ex)
 487   {
 488     SMESH_Comment str("Exception in  netgen::OCCGenerateMesh()");
 489     str << " at " << netgen::multithread.task
 490         << ": " << ex.DynamicType()->Name();
 491     if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
 492       str << ": " << ex.GetMessageString();
 493     error(str);
 494   }
 495   catch (...) {
 496     SMESH_Comment str("Exception in  netgen::OCCGenerateMesh()");
 497     str << " at " << netgen::multithread.task;
 498     error(str);
 499   }
 500
 501   // ----------------------------------------------------
 502   // Fill the SMESHDS with the generated nodes and faces
 503   // ----------------------------------------------------
 504
 505   int nbNodes = ngMesh->GetNP();
 506   int nbFaces = ngMesh->GetNSE();
 507
 508   int nbInputNodes = nodeVec.size();
 509   nodeVec.resize( nbNodes, 0 );
 510
 511   // add nodes
 512   for ( int i = nbInputNodes + 1; i <= nbNodes; ++i )
 513   {
 514     const MeshPoint& ngPoint = ngMesh->Point(i);
 515 #ifdef NETGEN_NEW
 516     SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));
 517 #else
 518     SMDS_MeshNode * node = meshDS->AddNode(ngPoint.X(), ngPoint.Y(), ngPoint.Z());
 519 #endif
 520     nodeVec[ i-1 ] = node;
 521   }
 522
 523   // create faces
 524   bool reverse = ( aShape.Orientation() == TopAbs_REVERSED );
 525   int i,j;
 526   for ( i = 1; i <= nbFaces ; ++i )
 527   {
 528     const Element2d& elem = ngMesh->SurfaceElement(i);
 529     vector<const SMDS_MeshNode*> nodes( elem.GetNP() );
 530     for (j=1; j <= elem.GetNP(); ++j)
 531     {
 532       int pind = elem.PNum(j);
 533       if ( pind-1 < 0 )
 534         break;
 535       const SMDS_MeshNode* node = nodeVec.at(pind-1);
 536       if ( reverse )
 537         nodes[ nodes.size()-j ] = node;
 538       else
 539         nodes[ j-1 ] = node;
 540       if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE )
 541       {
 542         const PointGeomInfo& pgi = elem.GeomInfoPi(j);
 543         meshDS->SetNodeOnFace((SMDS_MeshNode*)node, faceID, pgi.u, pgi.v);
 544       }
 545     }
 546     if ( j > elem.GetNP() )
 547     {
 548       SMDS_MeshFace* face = 0;
 549       if ( elem.GetType() == TRIG )
 550         face = helper.AddFace(nodes[0],nodes[1],nodes[2]);
 551       else
 552         face = helper.AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
 553     }
 554   }
 555
 556   return !err;
 557 }
 558
 559 #ifdef WITH_SMESH_CANCEL_COMPUTE
 560 void NETGENPlugin_NETGEN_2D_ONLY::CancelCompute()
 561 {
 562   SMESH_Algo::CancelCompute();
 563   netgen::multithread.terminate = 1;
 564 }
 565 #endif
 566
 567 //=============================================================================
 568 /*!
 569  *
 570  */
 571 //=============================================================================
 572
 573 bool NETGENPlugin_NETGEN_2D_ONLY::Evaluate(SMESH_Mesh& aMesh,
 574                                            const TopoDS_Shape& aShape,
 575                                            MapShapeNbElems& aResMap)
 576 {
 577   TopoDS_Face F = TopoDS::Face(aShape);
 578   if(F.IsNull())
 579     return false;
 580
 581   // collect info from edges
 582   int nb0d = 0, nb1d = 0;
 583   bool IsQuadratic = false;
 584   bool IsFirst = true;
 585   double fullLen = 0.0;
 586   TopTools_MapOfShape tmpMap;
 587   for (TopExp_Explorer exp(F, TopAbs_EDGE); exp.More(); exp.Next()) {
 588     TopoDS_Edge E = TopoDS::Edge(exp.Current());
 589     if( tmpMap.Contains(E) )
 590       continue;
 591     tmpMap.Add(E);
 592     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
 593     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
 594     if( anIt==aResMap.end() ) {
 595       SMESH_subMesh *sm = aMesh.GetSubMesh(F);
 596       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 597       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 598       return false;
 599     }
 600     std::vector<int> aVec = (*anIt).second;
 601     nb0d += aVec[SMDSEntity_Node];
 602     nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
 603     double aLen = SMESH_Algo::EdgeLength(E);
 604     fullLen += aLen;
 605     if(IsFirst) {
 606       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
 607       IsFirst = false;
 608     }
 609   }
 610   tmpMap.Clear();
 611
 612   // compute edge length
 613   double ELen = 0;
 614   if (_hypLengthFromEdges || !_hypLengthFromEdges && !_hypMaxElementArea) {
 615     if ( nb1d > 0 )
 616       ELen = fullLen / nb1d;
 617   }
 618   if ( _hypMaxElementArea ) {
 619     double maxArea = _hypMaxElementArea->GetMaxArea();
 620     ELen = sqrt(2. * maxArea/sqrt(3.0));
 621   }
 622   GProp_GProps G;
 623   BRepGProp::SurfaceProperties(F,G);
 624   double anArea = G.Mass();
 625
 626   const int hugeNb = numeric_limits<int>::max()/10;
 627   if ( anArea / hugeNb > ELen*ELen )
 628   {
 629     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
 630     SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 631     smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated.\nToo small element length",this));
 632     return false;
 633   }
 634   int nbFaces = (int) ( anArea / ( ELen*ELen*sqrt(3.) / 4 ) );
 635   int nbNodes = (int) ( ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
 636   std::vector<int> aVec(SMDSEntity_Last);
 637   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
 638   if( IsQuadratic ) {
 639     aVec[SMDSEntity_Node] = nbNodes;
 640     aVec[SMDSEntity_Quad_Triangle] = nbFaces;
 641   }
 642   else {
 643     aVec[SMDSEntity_Node] = nbNodes;
 644     aVec[SMDSEntity_Triangle] = nbFaces;
 645   }
 646   SMESH_subMesh *sm = aMesh.GetSubMesh(F);
 647   aResMap.insert(std::make_pair(sm,aVec));
 648
 649   return true;
 650 }