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