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