src/SMESH/SMESH_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   : SMESH_MEFISTO_2D.cxx
  25 //  Author : Paul RASCLE, EDF
  26 //  Module : SMESH
  27 //  $Header$
  28
  29 using namespace std;
  30 using namespace std;
  31 #include "SMESH_MEFISTO_2D.hxx"
  32 #include "SMESH_Gen.hxx"
  33 #include "SMESH_Mesh.hxx"
  34
  35 #include "SMESH_MaxElementArea.hxx"
  36 #include "SMESH_LengthFromEdges.hxx"
  37
  38 #include "Rn.h"
  39 #include "aptrte.h"
  40
  41 #include "SMESHDS_ListOfPtrHypothesis.hxx"
  42 #include "SMESHDS_ListIteratorOfListOfPtrHypothesis.hxx"
  43 #include "SMDS_MeshElement.hxx"
  44 #include "SMDS_MeshNode.hxx"
  45 #include "SMDS_EdgePosition.hxx"
  46 #include "SMDS_FacePosition.hxx"
  47
  48 #include "utilities.h"
  49
  50 #include <TopoDS_Face.hxx>
  51 #include <TopoDS_Edge.hxx>
  52 #include <TopoDS_Shape.hxx>
  53 #include <Geom_Surface.hxx>
  54 #include <GeomAdaptor_Curve.hxx>
  55 #include <Geom2d_Curve.hxx>
  56 #include <gp_Pnt2d.hxx>
  57 #include <BRep_Tool.hxx>
  58 #include <BRepTools.hxx>
  59 #include <BRepTools_WireExplorer.hxx>
  60 #include <GCPnts_AbscissaPoint.hxx>
  61 #include <GCPnts_UniformAbscissa.hxx>
  62 #include <TColStd_ListIteratorOfListOfInteger.hxx>
  63
  64 #include <string>
  65 #include <algorithm>
  66
  67 //=============================================================================
  68 /*!
  69  *
  70  */
  71 //=============================================================================
  72
  73 SMESH_MEFISTO_2D::SMESH_MEFISTO_2D(int hypId, int studyId, SMESH_Gen* gen)
  74   : SMESH_2D_Algo(hypId, studyId, gen)
  75 {
  76   MESSAGE("SMESH_MEFISTO_2D::SMESH_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 SMESH_MEFISTO_2D::~SMESH_MEFISTO_2D()
  96 {
  97   MESSAGE("SMESH_MEFISTO_2D::~SMESH_MEFISTO_2D");
  98 }
  99
 100 //=============================================================================
 101 /*!
 102  *
 103  */
 104 //=============================================================================
 105
 106 bool SMESH_MEFISTO_2D::CheckHypothesis(SMESH_Mesh& aMesh,
 107                                const TopoDS_Shape& aShape)
 108 {
 109   //MESSAGE("SMESH_MEFISTO_2D::CheckHypothesis");
 110
 111   _hypMaxElementArea = NULL;
 112   _hypLengthFromEdges = NULL;
 113
 114   list<SMESHDS_Hypothesis*>::const_iterator itl;
 115   SMESHDS_Hypothesis* theHyp;
 116
 117   const list<SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
 118   int nbHyp = hyps.size();
 119   if (nbHyp != 1) return false;  // only one compatible hypothesis allowed
 120
 121   itl = hyps.begin();
 122   theHyp = (*itl);
 123
 124   string hypName = theHyp->GetName();
 125   int hypId = theHyp->GetID();
 126   //SCRUTE(hypName);
 127
 128   bool isOk = false;
 129
 130   if (hypName == "MaxElementArea")
 131     {
 132       _hypMaxElementArea = dynamic_cast<SMESH_MaxElementArea*> (theHyp);
 133       ASSERT(_hypMaxElementArea);
 134       _maxElementArea = _hypMaxElementArea->GetMaxArea();
 135       _edgeLength = 0;
 136       isOk =true;
 137     }
 138
 139   if (hypName == "LengthFromEdges")
 140     {
 141       _hypLengthFromEdges = dynamic_cast<SMESH_LengthFromEdges*> (theHyp);
 142       ASSERT(_hypLengthFromEdges);
 143       _edgeLength = 0;
 144       _maxElementArea = 0;
 145       isOk =true;
 146     }
 147
 148
 149   if (isOk)
 150     {
 151       isOk = false;
 152       if (_maxElementArea > 0)
 153         {
 154           _edgeLength = 2*sqrt(_maxElementArea);  // triangles : minorant
 155           isOk = true;
 156         }
 157       else isOk = (_hypLengthFromEdges != NULL); // **** check mode
 158     }
 159
 160   //SCRUTE(_edgeLength);
 161   //SCRUTE(_maxElementArea);
 162   return isOk;
 163 }
 164
 165
 166 //=============================================================================
 167 /*!
 168  *
 169  */
 170 //=============================================================================
 171
 172 bool SMESH_MEFISTO_2D::Compute(SMESH_Mesh& aMesh,
 173                        const TopoDS_Shape& aShape)
 174 {
 175   MESSAGE("SMESH_MEFISTO_2D::Compute");
 176
 177   if (_hypLengthFromEdges)
 178     _edgeLength = ComputeEdgeElementLength(aMesh, aShape);
 179
 180   bool isOk = false;
 181   const Handle(SMESHDS_Mesh)& meshDS = aMesh.GetMeshDS();
 182   SMESH_subMesh* theSubMesh = aMesh.GetSubMesh(aShape);
 183
 184   const TopoDS_Face& FF = TopoDS::Face(aShape);
 185   bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
 186   TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
 187
 188   Z   nblf;        //nombre de lignes fermees (enveloppe en tete)
 189   Z  *nudslf=NULL; //numero du dernier sommet de chaque ligne fermee
 190   R2 *uvslf=NULL;
 191   Z   nbpti=0;     //nombre points internes futurs sommets de la triangulation
 192   R2 *uvpti=NULL;
 193
 194   Z   nbst;
 195   R2 *uvst=NULL;
 196   Z   nbt;
 197   Z  *nust=NULL;
 198   Z  ierr=0;
 199
 200   Z  nutysu=1;                 // 1: il existe un fonction areteideale_()
 201   // Z  nutysu=0;              // 0: on utilise aretmx
 202   R  aretmx=_edgeLength;       // longueur max aretes future triangulation
 203   //SCRUTE(aretmx);
 204
 205   nblf =  NumberOfWires(F);
 206   //SCRUTE(nblf);
 207
 208   nudslf = new Z[1+nblf];
 209   nudslf[0] = 0;
 210   int iw = 1;
 211   int nbpnt = 0;
 212
 213   const TopoDS_Wire OW1 = BRepTools::OuterWire(F);
 214   nbpnt += NumberOfPoints (aMesh, OW1);
 215   nudslf [iw++] = nbpnt;
 216   //SCRUTE(nbpnt);
 217
 218   for (TopExp_Explorer exp (F, TopAbs_WIRE); exp.More(); exp.Next())
 219     {
 220       const TopoDS_Wire& W = TopoDS::Wire(exp.Current());
 221       if (!OW1.IsSame(W))
 222         {
 223           nbpnt += NumberOfPoints (aMesh, W);
 224           nudslf [iw++] = nbpnt;
 225           //SCRUTE(nbpnt);
 226         }
 227     }
 228
 229   uvslf = new R2[nudslf[nblf]];
 230   //SCRUTE(nudslf[nblf]);
 231   int m = 0;
 232
 233   map<int,int> mefistoToDS;  // correspondence mefisto index--> points IDNodes
 234   TopoDS_Wire OW = BRepTools::OuterWire(F);
 235   LoadPoints (aMesh, F, OW, uvslf, m, mefistoToDS);
 236   //SCRUTE(m);
 237
 238   for (TopExp_Explorer exp (F, TopAbs_WIRE); exp.More(); exp.Next())
 239     {
 240       const TopoDS_Wire& W = TopoDS::Wire(exp.Current());
 241       if (!OW.IsSame(W))
 242         {
 243           LoadPoints (aMesh, F, W, uvslf, m, mefistoToDS);
 244           //SCRUTE(m);
 245         }
 246     }
 247 //   SCRUTE(nudslf[nblf]);
 248 //   for (int i=0; i<=nblf; i++)
 249 //     {
 250 //       MESSAGE(" -+- " <<i<< " "<< nudslf[i]);
 251 //     }
 252 //   for (int i=0; i<nudslf[nblf]; i++)
 253 //     {
 254 //       MESSAGE(" -+- " <<i<< " "<< uvslf[i]);
 255 //     }
 256 //   SCRUTE(nutysu);
 257 //   SCRUTE(aretmx);
 258 //   SCRUTE(nblf);
 259
 260   MESSAGE("MEFISTO triangulation ...");
 261   uvst = NULL;
 262   nust = NULL;
 263   aptrte( nutysu, aretmx,
 264           nblf, nudslf, uvslf,
 265           nbpti, uvpti,
 266           nbst, uvst, nbt, nust,
 267           ierr );
 268
 269   if( ierr == 0 )
 270     {
 271       MESSAGE("... End Triangulation");
 272       //SCRUTE(nbst);
 273       //SCRUTE(nbt);
 274       StoreResult (aMesh, nbst, uvst, nbt, nust, F,
 275                    faceIsForward, mefistoToDS);
 276       isOk = true;
 277     }
 278   else
 279     {
 280       MESSAGE("Error in Triangulation");
 281       isOk = false;
 282     }
 283   if (nudslf != NULL) delete [] nudslf;
 284   if (uvslf  != NULL) delete [] uvslf;
 285   if (uvst   != NULL) delete [] uvst;
 286   if (nust   != NULL) delete [] nust;
 287   return isOk;
 288 }
 289
 290 //=============================================================================
 291 /*!
 292  *
 293  */
 294 //=============================================================================
 295
 296 void SMESH_MEFISTO_2D::LoadPoints(SMESH_Mesh& aMesh,
 297                                   const TopoDS_Face& FF,
 298                                   const TopoDS_Wire& WW,
 299                                   R2* uvslf,
 300                                   int& m,
 301                                   map<int,int>& mefistoToDS)
 302 {
 303   MESSAGE("SMESH_MEFISTO_2D::LoadPoints");
 304
 305   Handle (SMDS_Mesh) meshDS = aMesh.GetMeshDS();
 306
 307   double scalex;
 308   double scaley;
 309   TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
 310   ComputeScaleOnFace(aMesh, F, scalex, scaley);
 311
 312   TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
 313   BRepTools_WireExplorer wexp(W,F);
 314   for (wexp.Init(W,F);wexp.More(); wexp.Next())
 315     {
 316       const TopoDS_Edge& E = wexp.Current();
 317
 318       // --- IDNodes of first and last Vertex
 319
 320       TopoDS_Vertex VFirst, VLast;
 321       TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
 322
 323       ASSERT(!VFirst.IsNull());
 324       SMESH_subMesh* firstSubMesh = aMesh.GetSubMesh(VFirst);
 325       const TColStd_ListOfInteger& lidf
 326         = firstSubMesh->GetSubMeshDS()->GetIDNodes();
 327       int idFirst= lidf.First();
 328 //       SCRUTE(idFirst);
 329
 330       ASSERT(!VLast.IsNull());
 331       SMESH_subMesh* lastSubMesh = aMesh.GetSubMesh(VLast);
 332       const TColStd_ListOfInteger& lidl
 333         = lastSubMesh->GetSubMeshDS()->GetIDNodes();
 334       int idLast= lidl.First();
 335 //       SCRUTE(idLast);
 336
 337       // --- edge internal IDNodes (relies on good order storage, not checked)
 338
 339       int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
 340       //SCRUTE(nbPoints);
 341
 342       Standard_Real f,l;
 343       Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E,F,f,l);
 344
 345       const TColStd_ListOfInteger& indElt
 346         = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetIDNodes();
 347       TColStd_ListIteratorOfListOfInteger ite(indElt);
 348       //SCRUTE(nbPoints);
 349       //SCRUTE(indElt.Extent());
 350       ASSERT(nbPoints == indElt.Extent());
 351       bool isForward = (E.Orientation() == TopAbs_FORWARD);
 352       map<double,int> params;
 353       for (; ite.More(); ite.Next())
 354         {
 355           int nodeId = ite.Value();
 356           Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
 357           Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
 358           Handle (SMDS_EdgePosition) epos
 359             = Handle (SMDS_EdgePosition)::DownCast(node->GetPosition());
 360           double param = epos->GetUParameter();
 361           params[param] = nodeId;
 362 //        MESSAGE(" " << param << " " << params[param]);
 363         }
 364
 365       // --- load 2D values into MEFISTO structure,
 366       //     add IDNodes in mefistoToDS map
 367
 368       if (E.Orientation() == TopAbs_FORWARD)
 369         {
 370           gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
 371           uvslf [m].x = scalex * p.X();
 372           uvslf [m].y = scaley * p.Y();
 373           mefistoToDS[m+1] = idFirst;
 374           //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
 375           //MESSAGE("__ f "<<f<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
 376           m++;
 377           map<double,int>::iterator itp = params.begin();
 378           for (Standard_Integer i = 1; i<=nbPoints; i++) // nbPoints internal
 379             {
 380               double param = (*itp).first;
 381               gp_Pnt2d p = C2d->Value(param);
 382               uvslf [m].x = scalex * p.X();
 383               uvslf [m].y = scaley * p.Y();
 384               mefistoToDS[m+1] = (*itp).second;
 385 //            MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
 386 //            MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
 387               m++;
 388               itp++;
 389             }
 390         }
 391       else
 392         {
 393           gp_Pnt2d p = C2d->Value(l); // last point = Vertex Reversed
 394           uvslf [m].x = scalex * p.X();
 395           uvslf [m].y = scaley * p.Y();
 396           mefistoToDS[m+1] = idLast;
 397 //        MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
 398 //        MESSAGE("__ l "<<l<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
 399           m++;
 400           map<double,int>::reverse_iterator itp = params.rbegin();
 401           for (Standard_Integer i = nbPoints ; i >= 1; i--)
 402             {
 403               double param = (*itp).first;
 404               gp_Pnt2d p = C2d->Value(param);
 405               uvslf [m].x = scalex * p.X();
 406               uvslf [m].y = scaley * p.Y();
 407               mefistoToDS[m+1] = (*itp).second;
 408 //            MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
 409 //            MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
 410               m++;
 411               itp++;
 412             }
 413         }
 414     }
 415 }
 416
 417 //=============================================================================
 418 /*!
 419  *
 420  */
 421 //=============================================================================
 422
 423 // **** a mettre dans SMESH_Algo ou SMESH_2D_Algo
 424
 425 void SMESH_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh& aMesh,
 426                                           const TopoDS_Face& aFace,
 427                                           double& scalex,
 428                                           double& scaley)
 429 {
 430   //MESSAGE("SMESH_MEFISTO_2D::ComputeScaleOnFace");
 431   TopoDS_Face F = TopoDS::Face(aFace.Oriented(TopAbs_FORWARD));
 432   TopoDS_Wire W = BRepTools::OuterWire(F);
 433
 434   BRepTools_WireExplorer wexp(W,F);
 435
 436   double xmin =  1.e300; // min & max of face 2D parametric coord.
 437   double xmax = -1.e300;
 438   double ymin =  1.e300;
 439   double ymax = -1.e300;
 440   int nbp = 50;
 441   scalex = 1;
 442   scaley = 1;
 443   for (wexp.Init(W,F);wexp.More(); wexp.Next())
 444     {
 445       const TopoDS_Edge& E = wexp.Current();
 446       double f,l;
 447       Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E,F,f,l);
 448       for (int i = 0; i<= nbp; i++)
 449         {
 450           double param = f + (double(i)/double(nbp))*(l-f);
 451           gp_Pnt2d p = C2d->Value(param);
 452           if (p.X() < xmin) xmin = p.X();
 453           if (p.X() > xmax) xmax = p.X();
 454           if (p.Y() < ymin) ymin = p.Y();
 455           if (p.Y() > ymax) ymax = p.Y();
 456 //        MESSAGE(" "<< f<<" "<<l<<" "<<param<<" "<<xmin<<" "<<xmax<<" "<<ymin<<" "<<ymax);
 457         }
 458     }
 459 //   SCRUTE(xmin);
 460 //   SCRUTE(xmax);
 461 //   SCRUTE(ymin);
 462 //   SCRUTE(ymax);
 463   double xmoy = (xmax + xmin)/2.;
 464   double ymoy = (ymax + ymin)/2.;
 465
 466   Handle(Geom_Surface) S = BRep_Tool::Surface(F); // 3D surface
 467
 468   double length_x = 0;
 469   double length_y = 0;
 470   gp_Pnt PX0 = S->Value(xmin, ymoy);
 471   gp_Pnt PY0 = S->Value(xmoy, ymin);
 472   for (Standard_Integer i = 1; i<= nbp; i++)
 473     {
 474       double x = xmin + (double(i)/double(nbp))*(xmax-xmin);
 475       gp_Pnt PX = S->Value(x,ymoy);
 476       double y = ymin + (double(i)/double(nbp))*(ymax-ymin);
 477       gp_Pnt PY = S->Value(xmoy,y);
 478       length_x += PX.Distance(PX0);
 479       length_y += PY.Distance(PY0);
 480       PX0.SetCoord(PX.X(),PX.Y(),PX.Z());
 481       PY0.SetCoord(PY.X(),PY.Y(),PY.Z());
 482     }
 483 //   SCRUTE(length_x);
 484 //   SCRUTE(length_y);
 485   scalex = length_x/(xmax - xmin);
 486   scaley = length_y/(ymax - ymin);
 487 //   SCRUTE(scalex);
 488 //   SCRUTE(scaley);
 489   ASSERT(scalex);
 490   ASSERT(scaley);
 491 }
 492
 493 //=============================================================================
 494 /*!
 495  *
 496  */
 497 //=============================================================================
 498
 499 void SMESH_MEFISTO_2D::StoreResult (SMESH_Mesh& aMesh,
 500                                     Z nbst, R2* uvst, Z nbt, Z* nust,
 501                                     const TopoDS_Face& F, bool faceIsForward,
 502                                     map<int,int>& mefistoToDS)
 503 {
 504   double scalex;
 505   double scaley;
 506   ComputeScaleOnFace(aMesh, F, scalex, scaley);
 507
 508   Handle (SMESHDS_Mesh) meshDS = aMesh.GetMeshDS();
 509
 510   Z n,m;
 511   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 512
 513   for ( n=0; n<nbst; n++ )
 514     {
 515       double u = uvst[n][0]/scalex;
 516       double v = uvst[n][1]/scaley;
 517       gp_Pnt P = S->Value(u,v);
 518
 519       if (mefistoToDS.find(n+1) == mefistoToDS.end())
 520         {
 521           int nodeId = meshDS->AddNode(P.X(), P.Y(), P.Z());
 522           Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
 523           Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
 524           meshDS->SetNodeOnFace(node, F);
 525
 526           //MESSAGE(nodeId<<" "<<P.X()<<" "<<P.Y()<<" "<<P.Z());
 527           mefistoToDS[n+1] = nodeId;
 528           //MESSAGE(" "<<n<<" "<<mefistoToDS[n+1]);
 529           Handle (SMDS_FacePosition) fpos
 530             = Handle (SMDS_FacePosition)::DownCast(node->GetPosition());
 531           fpos->SetUParameter(u);
 532           fpos->SetVParameter(v);
 533         }
 534     }
 535
 536   m=0;
 537   int mt=0;
 538
 539   //SCRUTE(faceIsForward);
 540   for ( n=1; n<=nbt; n++ )
 541     {
 542       int inode1 = nust[m++];
 543       int inode2 = nust[m++];
 544       int inode3 = nust[m++];
 545
 546       int nodeId1 = mefistoToDS[inode1];
 547       int nodeId2 = mefistoToDS[inode2];
 548       int nodeId3 = mefistoToDS[inode3];
 549       //MESSAGE("-- "<<inode1<<" "<<inode2<<" "<<inode3<<" ++ "<<nodeId1<<" "<<nodeId2<<" "<<nodeId3);
 550
 551       // triangle points must be in trigonometric order if face is Forward
 552       // else they must be put clockwise
 553
 554       bool triangleIsWellOriented = faceIsForward;
 555       int faceId;
 556       if (triangleIsWellOriented)
 557         {
 558           faceId = meshDS->AddFace(nodeId1, nodeId2, nodeId3);
 559         }
 560       else
 561         {
 562           faceId = meshDS->AddFace(nodeId1, nodeId3, nodeId2);
 563         }
 564       Handle (SMDS_MeshElement) elt = meshDS->FindElement(faceId);
 565       meshDS->SetMeshElementOnShape(elt, F);
 566       m++;
 567     }
 568 }
 569
 570 //=============================================================================
 571 /*!
 572  *
 573  */
 574 //=============================================================================
 575
 576 double SMESH_MEFISTO_2D::ComputeEdgeElementLength(SMESH_Mesh& aMesh,
 577                                                   const TopoDS_Shape& aShape)
 578 {
 579   MESSAGE("SMESH_MEFISTO_2D::ComputeEdgeElementLength");
 580   // **** a mettre dans SMESH_2D_Algo ?
 581
 582   const TopoDS_Face& FF = TopoDS::Face(aShape);
 583   bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
 584   TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
 585
 586   double meanElementLength = 100;
 587   double wireLength =0;
 588   int wireElementsNumber =0;
 589   for (TopExp_Explorer exp (F, TopAbs_WIRE); exp.More(); exp.Next())
 590     {
 591       const TopoDS_Wire& W = TopoDS::Wire(exp.Current());
 592       for (TopExp_Explorer expe(W,TopAbs_EDGE); expe.More(); expe.Next())
 593         {
 594           const TopoDS_Edge& E = TopoDS::Edge(expe.Current());
 595           int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
 596           double length = EdgeLength(E);
 597           wireLength += length;
 598           wireElementsNumber += nb;
 599         }
 600     }
 601   if (wireElementsNumber)
 602     meanElementLength = wireLength/wireElementsNumber;
 603   //SCRUTE(meanElementLength);
 604   return meanElementLength;
 605 }
 606
 607 //=============================================================================
 608 /*!
 609  *
 610  */
 611 //=============================================================================
 612
 613 ostream & SMESH_MEFISTO_2D::SaveTo(ostream & save)
 614 {
 615   return save << this;
 616 }
 617
 618 //=============================================================================
 619 /*!
 620  *
 621  */
 622 //=============================================================================
 623
 624 istream & SMESH_MEFISTO_2D::LoadFrom(istream & load)
 625 {
 626   return load >> (*this);
 627 }
 628
 629 //=============================================================================
 630 /*!
 631  *
 632  */
 633 //=============================================================================
 634
 635 ostream & operator << (ostream & save, SMESH_MEFISTO_2D & hyp)
 636 {
 637   return save;
 638 }
 639
 640 //=============================================================================
 641 /*!
 642  *
 643  */
 644 //=============================================================================
 645
 646 istream & operator >> (istream & load, SMESH_MEFISTO_2D & hyp)
 647 {
 648   return load;
 649 }
 650