src/StdMeshers/StdMeshers_MEFISTO_2D.cxx

   1 //  SMESH SMESH : implementaion of SMESH idl descriptions
   2 //
   3 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  21 //
  22 //
  23 //
  24 //  File   : StdMeshers_MEFISTO_2D.cxx
  25 //           Moved here from SMESH_MEFISTO_2D.cxx
  26 //  Author : Paul RASCLE, EDF
  27 //  Module : SMESH
  28 //  $Header$
  29
  30 #include "StdMeshers_MEFISTO_2D.hxx"
  31 #include "SMESH_Gen.hxx"
  32 #include "SMESH_Mesh.hxx"
  33 #include "SMESH_subMesh.hxx"
  34
  35 #include "StdMeshers_MaxElementArea.hxx"
  36 #include "StdMeshers_LengthFromEdges.hxx"
  37
  38 #include "Rn.h"
  39 #include "aptrte.h"
  40
  41 #include "SMDS_MeshElement.hxx"
  42 #include "SMDS_MeshNode.hxx"
  43 #include "SMDS_EdgePosition.hxx"
  44 #include "SMDS_FacePosition.hxx"
  45
  46 #include "utilities.h"
  47
  48 #include <TopoDS_Face.hxx>
  49 #include <TopoDS_Edge.hxx>
  50 #include <TopoDS_Shape.hxx>
  51 #include <Geom_Surface.hxx>
  52 #include <GeomAdaptor_Curve.hxx>
  53 #include <Geom2d_Curve.hxx>
  54 #include <gp_Pnt2d.hxx>
  55 #include <BRep_Tool.hxx>
  56 #include <BRepTools.hxx>
  57 #include <BRepTools_WireExplorer.hxx>
  58 #include <GCPnts_AbscissaPoint.hxx>
  59 #include <GCPnts_UniformAbscissa.hxx>
  60 #include <TColStd_ListIteratorOfListOfInteger.hxx>
  61 #include <TopTools_ListIteratorOfListOfShape.hxx>
  62 #include <TopTools_ListOfShape.hxx>
  63
  64 #include <string>
  65
  66 using namespace std;
  67 //#include <algorithm>
  68
  69 //=============================================================================
  70 /*!
  71  *
  72  */
  73 //=============================================================================
  74
  75 StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D(int hypId, int studyId,
  76         SMESH_Gen * gen):SMESH_2D_Algo(hypId, studyId, gen)
  77 {
  78         MESSAGE("StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D");
  79         _name = "MEFISTO_2D";
  80 //   _shapeType = TopAbs_FACE;
  81         _shapeType = (1 << TopAbs_FACE);
  82         _compatibleHypothesis.push_back("MaxElementArea");
  83         _compatibleHypothesis.push_back("LengthFromEdges");
  84
  85         _edgeLength = 0;
  86         _maxElementArea = 0;
  87         _hypMaxElementArea = NULL;
  88         _hypLengthFromEdges = NULL;
  89 }
  90
  91 //=============================================================================
  92 /*!
  93  *
  94  */
  95 //=============================================================================
  96
  97 StdMeshers_MEFISTO_2D::~StdMeshers_MEFISTO_2D()
  98 {
  99         MESSAGE("StdMeshers_MEFISTO_2D::~StdMeshers_MEFISTO_2D");
 100 }
 101
 102 //=============================================================================
 103 /*!
 104  *
 105  */
 106 //=============================================================================
 107
 108 bool StdMeshers_MEFISTO_2D::CheckHypothesis
 109                          (SMESH_Mesh& aMesh,
 110                           const TopoDS_Shape& aShape,
 111                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 112 {
 113         //MESSAGE("StdMeshers_MEFISTO_2D::CheckHypothesis");
 114
 115         _hypMaxElementArea = NULL;
 116         _hypLengthFromEdges = NULL;
 117
 118         list <const SMESHDS_Hypothesis * >::const_iterator itl;
 119         const SMESHDS_Hypothesis *theHyp;
 120
 121         const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
 122         int nbHyp = hyps.size();
 123         if (!nbHyp)
 124         {
 125           aStatus = SMESH_Hypothesis::HYP_MISSING;
 126           return false;  // can't work with no hypothesis
 127         }
 128
 129         itl = hyps.begin();
 130         theHyp = (*itl); // use only the first hypothesis
 131
 132         string hypName = theHyp->GetName();
 133         //int hypId = theHyp->GetID();
 134         //SCRUTE(hypName);
 135
 136         bool isOk = false;
 137
 138         if (hypName == "MaxElementArea")
 139         {
 140                 _hypMaxElementArea = static_cast<const StdMeshers_MaxElementArea *>(theHyp);
 141                 ASSERT(_hypMaxElementArea);
 142                 _maxElementArea = _hypMaxElementArea->GetMaxArea();
 143                 _edgeLength = 0;
 144                 isOk = true;
 145                 aStatus = SMESH_Hypothesis::HYP_OK;
 146         }
 147
 148         else if (hypName == "LengthFromEdges")
 149         {
 150                 _hypLengthFromEdges = static_cast<const StdMeshers_LengthFromEdges *>(theHyp);
 151                 ASSERT(_hypLengthFromEdges);
 152                 _edgeLength = 0;
 153                 _maxElementArea = 0;
 154                 isOk = true;
 155                 aStatus = SMESH_Hypothesis::HYP_OK;
 156         }
 157         else
 158           aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 159
 160         if (isOk)
 161         {
 162                 isOk = false;
 163                 if (_maxElementArea > 0)
 164                 {
 165 //                      _edgeLength = 2 * sqrt(_maxElementArea);        // triangles : minorant
 166                         _edgeLength = 2 * sqrt(_maxElementArea/sqrt(3.0));
 167                         isOk = true;
 168                 }
 169                 else
 170                         isOk = (_hypLengthFromEdges != NULL);   // **** check mode
 171                 if (!isOk)
 172                   aStatus = SMESH_Hypothesis::HYP_BAD_PARAMETER;
 173         }
 174
 175         //SCRUTE(_edgeLength);
 176         //SCRUTE(_maxElementArea);
 177         return isOk;
 178 }
 179
 180 //=============================================================================
 181 /*!
 182  *
 183  */
 184 //=============================================================================
 185
 186 bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape)
 187 {
 188         MESSAGE("StdMeshers_MEFISTO_2D::Compute");
 189
 190         if (_hypLengthFromEdges)
 191                 _edgeLength = ComputeEdgeElementLength(aMesh, aShape);
 192
 193         bool isOk = false;
 194         //const SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 195         //SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
 196
 197         const TopoDS_Face & FF = TopoDS::Face(aShape);
 198         bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
 199         TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
 200
 201         Z nblf;                                         //nombre de lignes fermees (enveloppe en tete)
 202         Z *nudslf = NULL;                       //numero du dernier sommet de chaque ligne fermee
 203         R2 *uvslf = NULL;
 204         Z nbpti = 0;                            //nombre points internes futurs sommets de la triangulation
 205         R2 *uvpti = NULL;
 206
 207         Z nbst;
 208         R2 *uvst = NULL;
 209         Z nbt;
 210         Z *nust = NULL;
 211         Z ierr = 0;
 212
 213         Z nutysu = 1;                           // 1: il existe un fonction areteideale_()
 214         // Z  nutysu=0;              // 0: on utilise aretmx
 215         R aretmx = _edgeLength;         // longueur max aretes future triangulation
 216
 217         nblf = NumberOfWires(F);
 218
 219         nudslf = new Z[1 + nblf];
 220         nudslf[0] = 0;
 221         int iw = 1;
 222         int nbpnt = 0;
 223
 224         myOuterWire = BRepTools::OuterWire(F);
 225         nbpnt += NumberOfPoints(aMesh, myOuterWire);
 226         if ( nbpnt < 3 ) // ex: a circle with 2 segments
 227           return false;
 228         nudslf[iw++] = nbpnt;
 229
 230         for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
 231         {
 232                 const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
 233                 if (!myOuterWire.IsSame(W))
 234                 {
 235                         nbpnt += NumberOfPoints(aMesh, W);
 236                         nudslf[iw++] = nbpnt;
 237                 }
 238         }
 239
 240         // avoid passing same uv points for a vertex common to 2 wires
 241         TopTools_IndexedDataMapOfShapeListOfShape VWMap;
 242         if ( iw - 1 > 1 ) // nbofWires > 1
 243           TopExp::MapShapesAndAncestors( F , TopAbs_VERTEX, TopAbs_WIRE, VWMap );
 244
 245         uvslf = new R2[nudslf[nblf]];
 246         int m = 0;
 247
 248         double scalex, scaley;
 249         ComputeScaleOnFace(aMesh, F, scalex, scaley);
 250
 251         map<int, const SMDS_MeshNode*> mefistoToDS;     // correspondence mefisto index--> points IDNodes
 252         if ( !LoadPoints(aMesh, F, myOuterWire, uvslf, m,
 253                          mefistoToDS, scalex, scaley, VWMap))
 254           return false;
 255
 256         for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
 257         {
 258                 const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
 259                 if (!myOuterWire.IsSame(W))
 260                 {
 261                         if (! LoadPoints(aMesh, F, W, uvslf, m,
 262                                          mefistoToDS, scalex, scaley, VWMap ))
 263                           return false;
 264                 }
 265         }
 266
 267         uvst = NULL;
 268         nust = NULL;
 269         aptrte(nutysu, aretmx,
 270                 nblf, nudslf, uvslf, nbpti, uvpti, nbst, uvst, nbt, nust, ierr);
 271
 272         if (ierr == 0)
 273           {
 274             MESSAGE("... End Triangulation Generated Triangle Number " << nbt);
 275             MESSAGE("                                    Node Number " << nbst);
 276             StoreResult(aMesh, nbst, uvst, nbt, nust, F,
 277                         faceIsForward, mefistoToDS, scalex, scaley);
 278             isOk = true;
 279           }
 280         else
 281         {
 282                 MESSAGE("Error in Triangulation");
 283                 isOk = false;
 284         }
 285         if (nudslf != NULL)
 286                 delete[]nudslf;
 287         if (uvslf != NULL)
 288                 delete[]uvslf;
 289         if (uvst != NULL)
 290                 delete[]uvst;
 291         if (nust != NULL)
 292                 delete[]nust;
 293         return isOk;
 294 }
 295
 296 //=======================================================================
 297 //function : fixOverlappedLinkUV
 298 //purpose  : prevent failure due to overlapped adjacent links
 299 //=======================================================================
 300
 301 static bool fixOverlappedLinkUV( R2& uv0, const R2& uv1, const R2& uv2 )
 302 {
 303   gp_XY v1( uv0.x - uv1.x, uv0.y - uv1.y );
 304   gp_XY v2( uv2.x - uv1.x, uv2.y - uv1.y );
 305
 306   double tol2 = DBL_MIN * DBL_MIN;
 307   double sqMod1 = v1.SquareModulus();
 308   if ( sqMod1 <= tol2 ) return false;
 309   double sqMod2 = v2.SquareModulus();
 310   if ( sqMod2 <= tol2 ) return false;
 311
 312   double dot = v1*v2;
 313
 314   // check sinus >= 1.e-3
 315   const double minSin = 1.e-3;
 316   if ( dot > 0 && 1 - dot * dot / ( sqMod1 * sqMod2 ) < minSin * minSin ) {
 317     MESSAGE(" ___ FIX UV ____" << uv0.x << " " << uv0.y);
 318     v1.SetCoord( -v1.Y(), v1.X() );
 319     double delta = sqrt( sqMod1 ) * minSin;
 320     if ( v1.X() < 0 )
 321       uv0.x -= delta;
 322     else
 323       uv0.x += delta;
 324     if ( v1.Y() < 0 )
 325       uv0.y -= delta;
 326     else
 327       uv0.y += delta;
 328 // #ifdef _DEBUG_
 329 //     MESSAGE(" -> " << uv0.x << " " << uv0.y << " ");
 330 //     MESSAGE("v1( " << v1.X() << " " << v1.Y() << " ) " <<
 331 //       "v2( " << v2.X() << " " << v2.Y() << " ) ");
 332 //    MESSAGE("SIN: " << sqrt(1 - dot * dot / (sqMod1 * sqMod2)));
 333 //     v1.SetCoord( uv0.x - uv1.x, uv0.y - uv1.y );
 334 //     v2.SetCoord( uv2.x - uv1.x, uv2.y - uv1.y );
 335 //     gp_XY v3( uv2.x - uv0.x, uv2.y - uv0.y );
 336 //     sqMod1 = v1.SquareModulus();
 337 //     sqMod2 = v2.SquareModulus();
 338 //     dot = v1*v2;
 339 //     double sin = sqrt(1 - dot * dot / (sqMod1 * sqMod2));
 340 //     MESSAGE("NEW SIN: " << sin);
 341 // #endif
 342     return true;
 343   }
 344   return false;
 345 }
 346
 347 //=======================================================================
 348 //function : fixCommonVertexUV
 349 //purpose  :
 350 //=======================================================================
 351
 352 static bool fixCommonVertexUV (gp_Pnt2d &           theUV,
 353                                const TopoDS_Vertex& theV,
 354                                const TopoDS_Wire&   theW,
 355                                const TopoDS_Wire&   theOW,
 356                                const TopoDS_Face&   theF,
 357                                const TopTools_IndexedDataMapOfShapeListOfShape & theVWMap,
 358                                SMESH_Mesh &         theMesh)
 359 {
 360   if( theW.IsSame( theOW ) ||
 361       !theVWMap.Contains( theV )) return false;
 362
 363   // check if there is another wire sharing theV
 364   const TopTools_ListOfShape& WList = theVWMap.FindFromKey( theV );
 365   TopTools_ListIteratorOfListOfShape aWIt;
 366   for ( aWIt.Initialize( WList ); aWIt.More(); aWIt.Next() )
 367     if ( !theW.IsSame( aWIt.Value() ))
 368       break;
 369   if ( !aWIt.More() ) return false;
 370
 371   TopTools_ListOfShape EList;
 372   list< double >       UList;
 373
 374   // find edges of theW sharing theV
 375   // and find 2d normal to them at theV
 376   gp_Vec2d N(0.,0.);
 377   TopoDS_Iterator itE( theW );
 378   for (  ; itE.More(); itE.Next() )
 379   {
 380     const TopoDS_Edge& E = TopoDS::Edge( itE.Value() );
 381     TopoDS_Iterator itV( E );
 382     for ( ; itV.More(); itV.Next() )
 383     {
 384       const TopoDS_Vertex & V = TopoDS::Vertex( itV.Value() );
 385       if ( !V.IsSame( theV ))
 386         continue;
 387       EList.Append( E );
 388       Standard_Real u = BRep_Tool::Parameter( V, E );
 389       UList.push_back( u );
 390       double f, l;
 391       Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, theF, f, l);
 392       gp_Vec2d d1;
 393       gp_Pnt2d p;
 394       C2d->D1( u, p, d1 );
 395       gp_Vec2d n( d1.Y(), -d1.X() );
 396       if ( E.Orientation() == TopAbs_REVERSED )
 397         n.Reverse();
 398       N += n.Normalized();
 399     }
 400   }
 401
 402   // define step size by which to move theUV
 403
 404   gp_Pnt2d nextUV; // uv of next node on edge, most distant of the four
 405   double maxDist = -DBL_MAX;
 406   TopTools_ListIteratorOfListOfShape aEIt (EList);
 407   list< double >::iterator aUIt = UList.begin();
 408   for ( ; aEIt.More(); aEIt.Next(), aUIt++ )
 409   {
 410     const TopoDS_Edge& E = TopoDS::Edge( aEIt.Value() );
 411     double f, l;
 412     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, theF, f, l);
 413
 414     double umin = DBL_MAX, umax = -DBL_MAX;
 415     SMDS_NodeIteratorPtr nIt = theMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
 416     if ( !nIt->more() ) // no nodes on edge, only on vertices
 417     {
 418       umin = l;
 419       umax = f;
 420     }
 421     else
 422     {
 423       while ( nIt->more() ) {
 424         const SMDS_MeshNode* node = nIt->next();
 425         const SMDS_EdgePosition* epos =
 426           static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
 427         double u = epos->GetUParameter();
 428         if ( u < umin )
 429           umin = u;
 430         if ( u > umax )
 431           umax = u;
 432       }
 433     }
 434     bool isFirstCommon = ( *aUIt == f );
 435     gp_Pnt2d uv = C2d->Value( isFirstCommon ? umin : umax );
 436     double dist = theUV.SquareDistance( uv );
 437     if ( dist > maxDist ) {
 438       maxDist = dist;
 439       nextUV  = uv;
 440     }
 441   }
 442   R2 uv0, uv1, uv2;
 443   uv0.x = theUV.X();    uv0.y = theUV.Y();
 444   uv1.x = nextUV.X();   uv1.y = nextUV.Y();
 445   uv2.x = uv0.x;        uv2.y = uv0.y;
 446   if ( fixOverlappedLinkUV( uv0, uv1, uv2 ))
 447   {
 448     double step = theUV.Distance( gp_Pnt2d( uv0.x, uv0.y ));
 449
 450     // move theUV along the normal by the step
 451
 452     N *= step;
 453
 454     MESSAGE("--fixCommonVertexUV move(" << theUV.X() << " " << theUV.Y()
 455             << ") by (" << N.X() << " " << N.Y() << ")"
 456             << endl << "--- MAX DIST " << maxDist);
 457
 458     theUV.SetXY( theUV.XY() + N.XY() );
 459
 460     return true;
 461   }
 462   return false;
 463 }
 464
 465 //=============================================================================
 466 /*!
 467  *
 468  */
 469 //=============================================================================
 470
 471 bool StdMeshers_MEFISTO_2D::LoadPoints(SMESH_Mesh &        aMesh,
 472                                        const TopoDS_Face & FF,
 473                                        const TopoDS_Wire & WW,
 474                                        R2 *                uvslf,
 475                                        int &               m,
 476                                        map<int, const SMDS_MeshNode*>&mefistoToDS,
 477                                        double              scalex,
 478                                        double              scaley,
 479                                        const TopTools_IndexedDataMapOfShapeListOfShape& VWMap)
 480 {
 481   //  MESSAGE("StdMeshers_MEFISTO_2D::LoadPoints");
 482
 483   TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
 484
 485   list< int > mOnVertex;
 486
 487   gp_XY scale( scalex, scaley );
 488
 489   TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
 490   BRepTools_WireExplorer wexp(W, F);
 491   for ( wexp.Init(W, F); wexp.More(); wexp.Next() )
 492   {
 493     const TopoDS_Edge & E = wexp.Current();
 494     bool isForward = (E.Orientation() == TopAbs_FORWARD);
 495
 496     // --- IDNodes of first and last Vertex
 497
 498     TopoDS_Vertex VFirst, VLast, V;
 499     TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
 500     V = isForward ? VFirst : VLast;
 501
 502     ASSERT(!V.IsNull());
 503     SMDS_NodeIteratorPtr lid= aMesh.GetSubMesh(V)->GetSubMeshDS()->GetNodes();
 504     if ( !lid->more() ) {
 505       MESSAGE (" NO NODE BUILT ON VERTEX ");
 506       return false;
 507     }
 508     const SMDS_MeshNode* idFirst = lid->next();
 509
 510     // --- edge internal IDNodes (relies on good order storage, not checked)
 511
 512     map<double, const SMDS_MeshNode*> params;
 513     const SMDS_MeshNode * node;
 514
 515     SMDS_NodeIteratorPtr nodeIt= aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
 516     while ( nodeIt->more() )
 517     {
 518       node = nodeIt->next();
 519       const SMDS_EdgePosition* epos =
 520         static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
 521       params[ epos->GetUParameter() ] = node;
 522     }
 523     int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
 524     if ( nbPoints != params.size())
 525     {
 526       MESSAGE( "BAD NODE ON EDGE POSITIONS" );
 527       return false;
 528     }
 529
 530     // --- load 2D values into MEFISTO structure,
 531     //     add IDNodes in mefistoToDS map
 532
 533     double f, l, uFirst, u;
 534     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
 535     uFirst = isForward ? f : l;
 536
 537     // vertex node
 538     gp_Pnt2d p = C2d->Value( uFirst ).XY().Multiplied( scale );
 539     if ( fixCommonVertexUV( p, V, W, myOuterWire, F, VWMap, aMesh ))
 540       myNodesOnCommonV.push_back( idFirst );
 541     mOnVertex.push_back( m );
 542     uvslf[m].x = p.X();
 543     uvslf[m].y = p.Y();
 544     mefistoToDS[m + 1] = idFirst;
 545     //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
 546     //MESSAGE("__ f "<<uFirst<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
 547     m++;
 548
 549     // internal nodes
 550     map<double, const SMDS_MeshNode*>::iterator u_n = params.begin();
 551     map<double, const SMDS_MeshNode*>::reverse_iterator u_n_rev = params.rbegin();
 552     for ( int i = 0; i < nbPoints; ++i )
 553     {
 554       if ( isForward ) {
 555         u    = u_n->first;
 556         node = u_n->second;
 557         ++u_n;
 558       } else {
 559         u    = u_n_rev->first;
 560         node = u_n_rev->second;
 561         ++u_n_rev;
 562       }
 563       gp_Pnt2d p = C2d->Value( u ).XY().Multiplied( scale );
 564       uvslf[m].x = p.X();
 565       uvslf[m].y = p.Y();
 566       mefistoToDS[m + 1] = node;
 567       //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
 568       //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
 569       m++;
 570     }
 571   } // for  wexp
 572
 573   // prevent failure on overlapped adjacent links,
 574   // check only links ending in vertex nodes
 575   int mFirst = mOnVertex.front(), mLast = m - 1;
 576   list< int >::iterator mIt = mOnVertex.begin();
 577   for ( ; mIt != mOnVertex.end(); ++mIt ) {
 578     int i = *mIt;
 579     int iB = i - 1, iA = i + 1; // indices Before and After
 580     if ( iB < mFirst ) iB = mLast;
 581     if ( iA > mLast )  iA = mFirst;
 582     fixOverlappedLinkUV (uvslf[ iB ], uvslf[ i ], uvslf[ iA ]);
 583   }
 584
 585   return true;
 586 }
 587
 588 //=============================================================================
 589 /*!
 590  *
 591  */
 592 //=============================================================================
 593
 594 void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh,
 595         const TopoDS_Face & aFace, double &scalex, double &scaley)
 596 {
 597   //MESSAGE("StdMeshers_MEFISTO_2D::ComputeScaleOnFace");
 598   TopoDS_Face F = TopoDS::Face(aFace.Oriented(TopAbs_FORWARD));
 599   TopoDS_Wire W = BRepTools::OuterWire(F);
 600
 601   double xmin = 1.e300;         // min & max of face 2D parametric coord.
 602   double xmax = -1.e300;
 603   double ymin = 1.e300;
 604   double ymax = -1.e300;
 605   int nbp = 23;
 606   scalex = 1;
 607   scaley = 1;
 608
 609   TopExp_Explorer wexp(W, TopAbs_EDGE);
 610   for ( ; wexp.More(); wexp.Next())
 611   {
 612     const TopoDS_Edge & E = TopoDS::Edge( wexp.Current() );
 613     double f, l;
 614     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
 615     if ( C2d.IsNull() ) continue;
 616     double du = (l - f) / double (nbp);
 617     for (int i = 0; i <= nbp; i++)
 618     {
 619       double param = f + double (i) * du;
 620       gp_Pnt2d p = C2d->Value(param);
 621       if (p.X() < xmin)
 622         xmin = p.X();
 623       if (p.X() > xmax)
 624         xmax = p.X();
 625       if (p.Y() < ymin)
 626         ymin = p.Y();
 627       if (p.Y() > ymax)
 628         ymax = p.Y();
 629       //    MESSAGE(" "<< f<<" "<<l<<" "<<param<<" "<<xmin<<" "<<xmax<<" "<<ymin<<" "<<ymax);
 630     }
 631   }
 632   //   SCRUTE(xmin);
 633   //   SCRUTE(xmax);
 634   //   SCRUTE(ymin);
 635   //   SCRUTE(ymax);
 636   double xmoy = (xmax + xmin) / 2.;
 637   double ymoy = (ymax + ymin) / 2.;
 638   double xsize = xmax - xmin;
 639   double ysize = ymax - ymin;
 640
 641   Handle(Geom_Surface) S = BRep_Tool::Surface(F);       // 3D surface
 642
 643   double length_x = 0;
 644   double length_y = 0;
 645   gp_Pnt PX0 = S->Value(xmin, ymoy);
 646   gp_Pnt PY0 = S->Value(xmoy, ymin);
 647   double dx = xsize / double (nbp);
 648   double dy = ysize / double (nbp);
 649   for (int i = 1; i <= nbp; i++)
 650   {
 651     double x = xmin + double (i) * dx;
 652     gp_Pnt PX = S->Value(x, ymoy);
 653     double y = ymin + double (i) * dy;
 654     gp_Pnt PY = S->Value(xmoy, y);
 655     length_x += PX.Distance(PX0);
 656     length_y += PY.Distance(PY0);
 657     PX0 = PX;
 658     PY0 = PY;
 659   }
 660   scalex = length_x / xsize;
 661   scaley = length_y / ysize;
 662 //   SCRUTE(xsize);
 663 //   SCRUTE(ysize);
 664   double xyratio = xsize*scalex/(ysize*scaley);
 665   const double maxratio = 1.e2;
 666   //SCRUTE(xyratio);
 667   if (xyratio > maxratio) {
 668     SCRUTE( scaley );
 669     scaley *= xyratio / maxratio;
 670     SCRUTE( scaley );
 671   }
 672   else if (xyratio < 1./maxratio) {
 673     SCRUTE( scalex );
 674     scalex *= 1 / xyratio / maxratio;
 675     SCRUTE( scalex );
 676   }
 677   ASSERT(scalex);
 678   ASSERT(scaley);
 679 }
 680
 681 //=============================================================================
 682 /*!
 683  *
 684  */
 685 //=============================================================================
 686
 687 void StdMeshers_MEFISTO_2D::StoreResult(SMESH_Mesh & aMesh,
 688                                         Z nbst, R2 * uvst, Z nbt, Z * nust,
 689                                         const TopoDS_Face & F, bool faceIsForward,
 690                                         map<int, const SMDS_MeshNode*>&mefistoToDS,
 691                                         double scalex, double scaley)
 692 {
 693   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 694   int faceID = meshDS->ShapeToIndex( F );
 695
 696   Z n, m;
 697   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 698
 699   for (n = 0; n < nbst; n++)
 700   {
 701     if (mefistoToDS.find(n + 1) == mefistoToDS.end())
 702     {
 703       double u = uvst[n][0] / scalex;
 704       double v = uvst[n][1] / scaley;
 705       gp_Pnt P = S->Value(u, v);
 706
 707       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 708       meshDS->SetNodeOnFace(node, faceID, u, v);
 709
 710       //MESSAGE(P.X()<<" "<<P.Y()<<" "<<P.Z());
 711       mefistoToDS[n + 1] = node;
 712       //MESSAGE("NEW: "<<n<<" "<<mefistoToDS[n+1]);
 713     }
 714   }
 715
 716   m = 0;
 717   //int mt = 0;
 718
 719   //SCRUTE(faceIsForward);
 720   for (n = 1; n <= nbt; n++)
 721   {
 722     int inode1 = nust[m++];
 723     int inode2 = nust[m++];
 724     int inode3 = nust[m++];
 725
 726     const SMDS_MeshNode *n1, *n2, *n3;
 727     n1 = mefistoToDS[inode1];
 728     n2 = mefistoToDS[inode2];
 729     n3 = mefistoToDS[inode3];
 730     //MESSAGE("-- "<<inode1<<" "<<inode2<<" "<<inode3);
 731
 732     // triangle points must be in trigonometric order if face is Forward
 733     // else they must be put clockwise
 734
 735     bool triangleIsWellOriented = faceIsForward;
 736
 737     SMDS_MeshElement * elt;
 738     if (triangleIsWellOriented)
 739       elt = meshDS->AddFace(n1, n2, n3);
 740     else
 741       elt = meshDS->AddFace(n1, n3, n2);
 742
 743     meshDS->SetMeshElementOnShape(elt, faceID);
 744     m++;
 745   }
 746
 747   // remove bad elements build on vertices shared by wires
 748
 749   list<const SMDS_MeshNode*>::iterator itN = myNodesOnCommonV.begin();
 750   for ( ; itN != myNodesOnCommonV.end(); itN++ )
 751   {
 752     const SMDS_MeshNode* node = *itN;
 753     SMDS_ElemIteratorPtr invElemIt = node->GetInverseElementIterator();
 754     while ( invElemIt->more() )
 755     {
 756       const SMDS_MeshElement* elem = invElemIt->next();
 757       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
 758       int nbSame = 0;
 759       while ( itN->more() )
 760         if ( itN->next() == node)
 761           nbSame++;
 762       if (nbSame > 1) {
 763         MESSAGE( "RM bad element " << elem->GetID());
 764         meshDS->RemoveElement( elem );
 765       }
 766     }
 767   }
 768
 769 }
 770
 771 //=============================================================================
 772 /*!
 773  *
 774  */
 775 //=============================================================================
 776
 777 double StdMeshers_MEFISTO_2D::ComputeEdgeElementLength(SMESH_Mesh & aMesh,
 778         const TopoDS_Shape & aShape)
 779 {
 780         //MESSAGE("StdMeshers_MEFISTO_2D::ComputeEdgeElementLength");
 781         // **** a mettre dans SMESH_2D_Algo ?
 782
 783         const TopoDS_Face & FF = TopoDS::Face(aShape);
 784         //bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
 785         TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
 786
 787         double meanElementLength = 100;
 788         double wireLength = 0;
 789         int wireElementsNumber = 0;
 790         for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
 791         {
 792                 const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
 793                 for (TopExp_Explorer expe(W, TopAbs_EDGE); expe.More(); expe.Next())
 794                 {
 795                         const TopoDS_Edge & E = TopoDS::Edge(expe.Current());
 796                         int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
 797                         double length = EdgeLength(E);
 798                         wireLength += length;
 799                         wireElementsNumber += nb;
 800                 }
 801         }
 802         if (wireElementsNumber)
 803                 meanElementLength = wireLength / wireElementsNumber;
 804         //SCRUTE(meanElementLength);
 805         return meanElementLength;
 806 }
 807
 808 //=============================================================================
 809 /*!
 810  *
 811  */
 812 //=============================================================================
 813
 814 ostream & StdMeshers_MEFISTO_2D::SaveTo(ostream & save)
 815 {
 816   return save;
 817 }
 818
 819 //=============================================================================
 820 /*!
 821  *
 822  */
 823 //=============================================================================
 824
 825 istream & StdMeshers_MEFISTO_2D::LoadFrom(istream & load)
 826 {
 827   return load;
 828 }
 829
 830 //=============================================================================
 831 /*!
 832  *
 833  */
 834 //=============================================================================
 835
 836 ostream & operator <<(ostream & save, StdMeshers_MEFISTO_2D & hyp)
 837 {
 838   return hyp.SaveTo( save );
 839 }
 840
 841 //=============================================================================
 842 /*!
 843  *
 844  */
 845 //=============================================================================
 846
 847 istream & operator >>(istream & load, StdMeshers_MEFISTO_2D & hyp)
 848 {
 849   return hyp.LoadFrom( load );
 850 }