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