Salome HOME
0022523: [CEA 1096] Add a colomn "biquadratic" in "Mesh Information"
[modules/smesh.git] / src / StdMeshers / StdMeshers_MEFISTO_2D.cxx
index 3cdff0adf84774e14769ff1906b242c5320428ce..75db563136dd7a0ab154bcb6785bd8771646fc61 100644 (file)
@@ -1,66 +1,76 @@
-//  SMESH SMESH : implementaion of SMESH idl descriptions
+// Copyright (C) 2007-2014  CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 //
-//  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
-//  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS 
-// 
-//  This library is free software; you can redistribute it and/or 
-//  modify it under the terms of the GNU Lesser General Public 
-//  License as published by the Free Software Foundation; either 
-//  version 2.1 of the License. 
-// 
-//  This library is distributed in the hope that it will be useful, 
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of 
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
-//  Lesser General Public License for more details. 
-// 
-//  You should have received a copy of the GNU Lesser General Public 
-//  License along with this library; if not, write to the Free Software 
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA 
-// 
-//  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org 
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
 //
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
 //
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+//  SMESH SMESH : implementaion of SMESH idl descriptions
 //  File   : StdMeshers_MEFISTO_2D.cxx
 //           Moved here from SMESH_MEFISTO_2D.cxx
 //  Author : Paul RASCLE, EDF
 //  Module : SMESH
-//  $Header$
-
-using namespace std;
+//
 #include "StdMeshers_MEFISTO_2D.hxx"
+
+#include "SMDS_EdgePosition.hxx"
+#include "SMDS_MeshElement.hxx"
+#include "SMDS_MeshNode.hxx"
+#include "SMESH_Comment.hxx"
 #include "SMESH_Gen.hxx"
 #include "SMESH_Mesh.hxx"
+#include "SMESH_MesherHelper.hxx"
 #include "SMESH_subMesh.hxx"
-
-#include "StdMeshers_MaxElementArea.hxx"
+#include "StdMeshers_FaceSide.hxx"
 #include "StdMeshers_LengthFromEdges.hxx"
+#include "StdMeshers_MaxElementArea.hxx"
+#include "StdMeshers_ViscousLayers2D.hxx"
+
+#include "utilities.h"
 
 #include "Rn.h"
 #include "aptrte.h"
 
-#include "SMDS_MeshElement.hxx"
-#include "SMDS_MeshNode.hxx"
-#include "SMDS_EdgePosition.hxx"
-#include "SMDS_FacePosition.hxx"
-
-#include "utilities.h"
-
-#include <TopoDS_Face.hxx>
-#include <TopoDS_Edge.hxx>
-#include <TopoDS_Shape.hxx>
-#include <Geom_Surface.hxx>
-#include <GeomAdaptor_Curve.hxx>
-#include <Geom2d_Curve.hxx>
-#include <gp_Pnt2d.hxx>
-#include <BRep_Tool.hxx>
+#include <BRepGProp.hxx>
 #include <BRepTools.hxx>
-#include <BRepTools_WireExplorer.hxx>
-#include <GCPnts_AbscissaPoint.hxx>
-#include <GCPnts_UniformAbscissa.hxx>
-#include <TColStd_ListIteratorOfListOfInteger.hxx>
+#include <BRep_Tool.hxx>
+#include <GProp_GProps.hxx>
+#include <Geom2d_Curve.hxx>
+#include <Geom_Curve.hxx>
+#include <Geom_Surface.hxx>
+#include <Precision.hxx>
+#include <TopExp.hxx>
+#include <TopExp_Explorer.hxx>
 #include <TopTools_ListIteratorOfListOfShape.hxx>
 #include <TopTools_ListOfShape.hxx>
+#include <TopTools_MapOfShape.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopoDS_Iterator.hxx>
+#include <TopoDS_Wire.hxx>
+#include <gp_Pnt2d.hxx>
+
+using namespace std;
+
+#ifdef _DEBUG_
+//#define DUMP_POINTS // to print coordinates of MEFISTO input
+#endif
 
 //=============================================================================
 /*!
@@ -68,20 +78,21 @@ using namespace std;
  */
 //=============================================================================
 
-StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D(int hypId, int studyId,
-                                             SMESH_Gen * gen):SMESH_2D_Algo(hypId, studyId, gen)
+StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D(int hypId, int studyId, SMESH_Gen * gen):
+  SMESH_2D_Algo(hypId, studyId, gen)
 {
   MESSAGE("StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D");
   _name = "MEFISTO_2D";
   _shapeType = (1 << TopAbs_FACE);
   _compatibleHypothesis.push_back("MaxElementArea");
   _compatibleHypothesis.push_back("LengthFromEdges");
+  _compatibleHypothesis.push_back("ViscousLayers2D");
 
   _edgeLength = 0;
   _maxElementArea = 0;
   _hypMaxElementArea = NULL;
   _hypLengthFromEdges = NULL;
-  myTool = 0;
+  _helper = 0;
 }
 
 //=============================================================================
@@ -102,75 +113,68 @@ StdMeshers_MEFISTO_2D::~StdMeshers_MEFISTO_2D()
 //=============================================================================
 
 bool StdMeshers_MEFISTO_2D::CheckHypothesis
-                         (SMESH_Mesh& aMesh,
-                          const TopoDS_Shape& aShape,
+                         (SMESH_Mesh&                          aMesh,
+                          const TopoDS_Shape&                  aShape,
                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 {
-       //MESSAGE("StdMeshers_MEFISTO_2D::CheckHypothesis");
+  _hypMaxElementArea = NULL;
+  _hypLengthFromEdges = NULL;
+  _edgeLength = 0;
+  _maxElementArea = 0;
 
-       _hypMaxElementArea = NULL;
-       _hypLengthFromEdges = NULL;
+  list <const SMESHDS_Hypothesis * >::const_iterator itl;
+  const SMESHDS_Hypothesis *theHyp;
 
-       list <const SMESHDS_Hypothesis * >::const_iterator itl;
-       const SMESHDS_Hypothesis *theHyp;
+  const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
+  int nbHyp = hyps.size();
+  if (!nbHyp)
+  {
+    aStatus = SMESH_Hypothesis::HYP_OK; //SMESH_Hypothesis::HYP_MISSING;
+    return true;  // (PAL13464) can work with no hypothesis, LengthFromEdges is default one
+  }
 
-       const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
-       int nbHyp = hyps.size();
-        if (!nbHyp)
-        {
-          aStatus = SMESH_Hypothesis::HYP_MISSING;
-          return false;  // can't work with no hypothesis
-        }
+  itl = hyps.begin();
+  theHyp = (*itl); // use only the first hypothesis
+
+  string hypName = theHyp->GetName();
+
+  bool isOk = false;
 
-       itl = hyps.begin();
-       theHyp = (*itl); // use only the first hypothesis
-
-       string hypName = theHyp->GetName();
-       //int hypId = theHyp->GetID();
-       //SCRUTE(hypName);
-
-       bool isOk = false;
-
-       if (hypName == "MaxElementArea")
-       {
-               _hypMaxElementArea = static_cast<const StdMeshers_MaxElementArea *>(theHyp);
-               ASSERT(_hypMaxElementArea);
-               _maxElementArea = _hypMaxElementArea->GetMaxArea();
-               _edgeLength = 0;
-               isOk = true;
-                aStatus = SMESH_Hypothesis::HYP_OK;
-       }
-
-       else if (hypName == "LengthFromEdges")
-       {
-               _hypLengthFromEdges = static_cast<const StdMeshers_LengthFromEdges *>(theHyp);
-               ASSERT(_hypLengthFromEdges);
-               _edgeLength = 0;
-               _maxElementArea = 0;
-               isOk = true;
-                aStatus = SMESH_Hypothesis::HYP_OK;
-       }
-        else
-          aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
-
-       if (isOk)
-       {
-               isOk = false;
-               if (_maxElementArea > 0)
-               {
-//                     _edgeLength = 2 * sqrt(_maxElementArea);        // triangles : minorant
-                       _edgeLength = 2 * sqrt(_maxElementArea/sqrt(3.0));
-                       isOk = true;
-               }
-               else
-                       isOk = (_hypLengthFromEdges != NULL);   // **** check mode
-                if (!isOk)
-                  aStatus = SMESH_Hypothesis::HYP_BAD_PARAMETER;
-       }
-
-       //SCRUTE(_edgeLength);
-       //SCRUTE(_maxElementArea);
-       return isOk;
+  if (hypName == "MaxElementArea")
+  {
+    _hypMaxElementArea = static_cast<const StdMeshers_MaxElementArea *>(theHyp);
+    ASSERT(_hypMaxElementArea);
+    _maxElementArea = _hypMaxElementArea->GetMaxArea();
+    isOk = true;
+    aStatus = SMESH_Hypothesis::HYP_OK;
+  }
+
+  else if (hypName == "LengthFromEdges")
+  {
+    _hypLengthFromEdges = static_cast<const StdMeshers_LengthFromEdges *>(theHyp);
+    ASSERT(_hypLengthFromEdges);
+    isOk = true;
+    aStatus = SMESH_Hypothesis::HYP_OK;
+  }
+  else
+    aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
+
+  if (isOk)
+  {
+    isOk = false;
+    if (_maxElementArea > 0)
+    {
+      //_edgeLength = 2 * sqrt(_maxElementArea);        // triangles : minorant
+      _edgeLength = sqrt(2. * _maxElementArea/sqrt(3.0));
+      isOk = true;
+    }
+    else
+      isOk = (_hypLengthFromEdges != NULL);     // **** check mode
+    if (!isOk)
+      aStatus = SMESH_Hypothesis::HYP_BAD_PARAMETER;
+  }
+
+  return isOk;
 }
 
 //=============================================================================
@@ -183,21 +187,52 @@ bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh
 {
   MESSAGE("StdMeshers_MEFISTO_2D::Compute");
 
-  if (_hypLengthFromEdges)
-    _edgeLength = ComputeEdgeElementLength(aMesh, aShape);
-  
-  bool isOk = false;
-  //const SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
-  //SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
+  TopoDS_Face F = TopoDS::Face(aShape.Oriented(TopAbs_FORWARD));
+
+  // helper builds quadratic mesh if necessary
+  SMESH_MesherHelper helper(aMesh);
+  _helper = &helper;
+  _quadraticMesh = _helper->IsQuadraticSubMesh(aShape);
+  const bool skipMediumNodes = _quadraticMesh;
+
+  // build viscous layers if required
+  SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F );
+  if ( !proxyMesh )
+    return false;
+
+  // get all edges of a face
+  TError problem;
+  TWireVector wires =
+    StdMeshers_FaceSide::GetFaceWires( F, aMesh, skipMediumNodes, problem, proxyMesh );
+  int nbWires = wires.size();
+  if ( problem && !problem->IsOK() ) return error( problem );
+  if ( nbWires == 0 ) return error( "Problem in StdMeshers_FaceSide::GetFaceWires()");
+  if ( wires[0]->NbSegments() < 3 ) // ex: a circle with 2 segments
+    return error(COMPERR_BAD_INPUT_MESH,
+                 SMESH_Comment("Too few segments: ")<<wires[0]->NbSegments());
+
+  // compute average edge length
+  if (!_hypMaxElementArea)
+  {
+    _edgeLength = 0;
+    int nbSegments = 0;
+    for ( int iW = 0; iW < nbWires; ++iW )
+    {
+      StdMeshers_FaceSidePtr wire = wires[ iW ];
+      _edgeLength += wire->Length();
+      nbSegments  += wire->NbSegments();
+    }
+    if ( nbSegments )
+      _edgeLength /= nbSegments;
+  }
 
-  const TopoDS_Face & FF = TopoDS::Face(aShape);
-  bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
-  TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
+  if (/*_hypLengthFromEdges &&*/ _edgeLength < DBL_MIN )
+    _edgeLength = 100;
 
-  Z nblf;                                              //nombre de lignes fermees (enveloppe en tete)
-  Z *nudslf = NULL;                    //numero du dernier sommet de chaque ligne fermee
-  R2 *uvslf = NULL;
-  Z nbpti = 0;                         //nombre points internes futurs sommets de la triangulation
+  Z nblf;                 //nombre de lignes fermees (enveloppe en tete)
+  Z *nudslf = NULL;       //numero du dernier sommet de chaque ligne fermee
+  R2 *uvslf = NULL;       
+  Z nbpti = 0;            //nombre points internes futurs sommets de la triangulation
   R2 *uvpti = NULL;
   
   Z nbst;
@@ -206,101 +241,147 @@ bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh
   Z *nust = NULL;
   Z ierr = 0;
 
-  Z nutysu = 1;                                // 1: il existe un fonction areteideale_()
-  // Z  nutysu=0;              // 0: on utilise aretmx
-  R aretmx = _edgeLength;              // longueur max aretes future triangulation
+  Z nutysu = 1;           // 1: il existe un fonction areteideale_()
+  // Z  nutysu=0;         // 0: on utilise aretmx
+  R aretmx = _edgeLength; // longueur max aretes future triangulation
+  if ( _hypMaxElementArea )
+    aretmx *= 1.5;
   
-  nblf = NumberOfWires(F);
+  nblf = nbWires;
   
   nudslf = new Z[1 + nblf];
   nudslf[0] = 0;
   int iw = 1;
   int nbpnt = 0;
 
-  myTool = new SMESH_MesherHelper(aMesh);
-  _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
-
-  if ( _quadraticMesh && _hypLengthFromEdges )
-    aretmx *= 2.;
-
-  myOuterWire = BRepTools::OuterWire(F);
-  nbpnt += NumberOfPoints(aMesh, myOuterWire);
-  if ( nbpnt < 3 ) { // ex: a circle with 2 segments
-    delete myTool; myTool = 0;
-    return false;
-  }
-  nudslf[iw++] = nbpnt;
-
-  for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next()) {
-    const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
-    if (!myOuterWire.IsSame(W)) {
-      nbpnt += NumberOfPoints(aMesh, W);
-      nudslf[iw++] = nbpnt;
-    }
+  // count nb of input points
+  for ( int iW = 0; iW < nbWires; ++iW )
+  {
+    nbpnt += wires[iW]->NbPoints() - 1;
+    nudslf[iw++] = nbpnt;
   }
 
-  // avoid passing same uv points for a vertex common to 2 wires
-  TopTools_IndexedDataMapOfShapeListOfShape VWMap;
-  if ( iw - 1 > 1 ) // nbofWires > 1
-    TopExp::MapShapesAndAncestors( F , TopAbs_VERTEX, TopAbs_WIRE, VWMap );
-
   uvslf = new R2[nudslf[nblf]];
-  int m = 0;
 
   double scalex, scaley;
   ComputeScaleOnFace(aMesh, F, scalex, scaley);
 
-  map<int, const SMDS_MeshNode*> mefistoToDS;  // correspondence mefisto index--> points IDNodes
-  if ( !LoadPoints(aMesh, F, myOuterWire, uvslf, m,
-                   mefistoToDS, scalex, scaley, VWMap) ) {
-    delete myTool; myTool = 0;
-    return false;
-  }
+  // correspondence mefisto index --> Nodes
+  vector< const SMDS_MeshNode*> mefistoToDS(nbpnt, (const SMDS_MeshNode*)0);
 
-  for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
+  bool isOk = false;
+
+  // fill input points UV
+  if ( LoadPoints(wires, uvslf, mefistoToDS, scalex, scaley) )
   {
-    const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
-    if (!myOuterWire.IsSame(W))
+    // Compute
+    aptrte(nutysu, aretmx,
+           nblf, nudslf, uvslf, nbpti, uvpti, nbst, uvst, nbt, nust, ierr);
+
+    if (ierr == 0)
     {
-      if (! LoadPoints(aMesh, F, W, uvslf, m,
-                       mefistoToDS, scalex, scaley, VWMap )) {
-        delete myTool; myTool = 0;
-        return false;
-      }
+      MESSAGE("... End Triangulation Generated Triangle Number " << nbt);
+      MESSAGE("                                    Node Number " << nbst);
+      StoreResult(nbst, uvst, nbt, nust, mefistoToDS, scalex, scaley);
+      isOk = true;
+    }
+    else
+    {
+      error(ierr,"Error in Triangulation (aptrte())");
     }
   }
+  if (nudslf != NULL) delete[]nudslf;
+  if (uvslf != NULL)  delete[]uvslf;
+  if (uvst != NULL)   delete[]uvst;
+  if (nust != NULL)   delete[]nust;
 
-  uvst = NULL;
-  nust = NULL;
-  aptrte(nutysu, aretmx,
-         nblf, nudslf, uvslf, nbpti, uvpti, nbst, uvst, nbt, nust, ierr);
+  return isOk;
+}
 
-  if (ierr == 0)
-  {
-    MESSAGE("... End Triangulation Generated Triangle Number " << nbt);
-    MESSAGE("                                    Node Number " << nbst);
-    StoreResult(aMesh, nbst, uvst, nbt, nust, F,
-                faceIsForward, mefistoToDS, scalex, scaley);
-    isOk = true;
+
+//=============================================================================
+/*!
+ *  
+ */
+//=============================================================================
+
+bool StdMeshers_MEFISTO_2D::Evaluate(SMESH_Mesh & aMesh,
+                                     const TopoDS_Shape & aShape,
+                                     MapShapeNbElems& aResMap)
+{
+  MESSAGE("StdMeshers_MEFISTO_2D::Evaluate");
+
+  TopoDS_Face F = TopoDS::Face(aShape.Oriented(TopAbs_FORWARD));
+
+  double aLen = 0.0;
+  int NbSeg = 0;
+  bool IsQuadratic = false;
+  bool IsFirst = true;
+  TopExp_Explorer exp(F,TopAbs_EDGE);
+  for(; exp.More(); exp.Next()) {
+    TopoDS_Edge E = TopoDS::Edge(exp.Current());
+    MapShapeNbElemsItr anIt = aResMap.find( aMesh.GetSubMesh(E) );
+    if( anIt == aResMap.end() ) continue;
+    std::vector<int> aVec = (*anIt).second;
+    int nbe = Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
+    NbSeg += nbe;
+    if(IsFirst) {
+      IsQuadratic = ( aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge] );
+      IsFirst = false;
+    }
+    double a,b;
+    TopLoc_Location L;
+    Handle(Geom_Curve) C = BRep_Tool::Curve(E,L,a,b);
+    gp_Pnt P1;
+    C->D0(a,P1);
+    double dp = (b-a)/nbe;
+    for(int i=1; i<=nbe; i++) {
+      gp_Pnt P2;
+      C->D0(a+i*dp,P2);
+      aLen += P1.Distance(P2);
+      P1 = P2;
+    }
   }
-  else
-  {
-    MESSAGE("Error in Triangulation");
-    isOk = false;
+  if(NbSeg<1) {
+    std::vector<int> aResVec(SMDSEntity_Last);
+    for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
+    SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+    aResMap.insert(std::make_pair(sm,aResVec));
+    SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+    smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
+                                          "Submesh can not be evaluated",this));
+    return false;
   }
-  if (nudslf != NULL)
-    delete[]nudslf;
-  if (uvslf != NULL)
-    delete[]uvslf;
-  if (uvst != NULL)
-    delete[]uvst;
-  if (nust != NULL)
-    delete[]nust;
-  delete myTool; myTool = 0;
+  aLen = aLen/NbSeg; // middle length
 
-  return isOk;
+  _edgeLength = Precision::Infinite();
+  double tmpLength = Min( _edgeLength, aLen );
+
+  GProp_GProps G;
+  BRepGProp::SurfaceProperties(aShape,G);
+  double anArea = G.Mass();
+
+  int nbFaces = Precision::IsInfinite( tmpLength ) ? 0 :
+    (int)( anArea/(tmpLength*tmpLength*sqrt(3.)/4) );
+  int nbNodes = (int) ( nbFaces*3 - (NbSeg-1)*2 ) / 6;
+
+  std::vector<int> aVec(SMDSEntity_Last);
+  for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
+  if(IsQuadratic) {
+    aVec[SMDSEntity_Quad_Triangle] = nbFaces;
+    aVec[SMDSEntity_Node] = (int)( nbNodes + nbFaces*3 - (NbSeg-1) );
+  }
+  else {
+    aVec[SMDSEntity_Node] = nbNodes;
+    aVec[SMDSEntity_Triangle] = nbFaces;
+  }
+  SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+  aResMap.insert(std::make_pair(sm,aVec));
+
+  return true;
 }
 
+
 //=======================================================================
 //function : fixOverlappedLinkUV
 //purpose  : prevent failure due to overlapped adjacent links
@@ -357,25 +438,30 @@ static bool fixOverlappedLinkUV( R2& uv0, const R2& uv1, const R2& uv2 )
 //purpose  : 
 //=======================================================================
 
-static bool fixCommonVertexUV (gp_Pnt2d &           theUV,
+static bool fixCommonVertexUV (R2 &                 theUV,
                                const TopoDS_Vertex& theV,
-                               const TopoDS_Wire&   theW,
-                               const TopoDS_Wire&   theOW,
                                const TopoDS_Face&   theF,
                                const TopTools_IndexedDataMapOfShapeListOfShape & theVWMap,
                                SMESH_Mesh &         theMesh,
-                               bool CreateQuadratic)
+                               const double         theScaleX,
+                               const double         theScaleY,
+                               const bool           theCreateQuadratic)
 {
-  if( theW.IsSame( theOW ) ||
-      !theVWMap.Contains( theV )) return false;
+  if( !theVWMap.Contains( theV )) return false;
 
   // check if there is another wire sharing theV
   const TopTools_ListOfShape& WList = theVWMap.FindFromKey( theV );
   TopTools_ListIteratorOfListOfShape aWIt;
+  TopTools_MapOfShape aWires;
   for ( aWIt.Initialize( WList ); aWIt.More(); aWIt.Next() )
-    if ( !theW.IsSame( aWIt.Value() ))
-      break;
-  if ( !aWIt.More() ) return false;
+    aWires.Add( aWIt.Value() );
+  if ( aWires.Extent() < 2 ) return false;
+
+  TopoDS_Shape anOuterWire = BRepTools::OuterWire(theF);
+  TopoDS_Shape anInnerWire;
+  for ( aWIt.Initialize( WList ); aWIt.More() && anInnerWire.IsNull(); aWIt.Next() )
+    if ( !anOuterWire.IsSame( aWIt.Value() ))
+      anInnerWire = aWIt.Value();
 
   TopTools_ListOfShape EList;
   list< double >       UList;
@@ -383,7 +469,7 @@ static bool fixCommonVertexUV (gp_Pnt2d &           theUV,
   // find edges of theW sharing theV
   // and find 2d normal to them at theV
   gp_Vec2d N(0.,0.);
-  TopoDS_Iterator itE( theW );
+  TopoDS_Iterator itE( anInnerWire );
   for (  ; itE.More(); itE.Next() )
   {
     const TopoDS_Edge& E = TopoDS::Edge( itE.Value() );
@@ -401,7 +487,7 @@ static bool fixCommonVertexUV (gp_Pnt2d &           theUV,
       gp_Vec2d d1;
       gp_Pnt2d p;
       C2d->D1( u, p, d1 );
-      gp_Vec2d n( d1.Y(), -d1.X() );
+      gp_Vec2d n( d1.Y() * theScaleX, -d1.X() * theScaleY);
       if ( E.Orientation() == TopAbs_REVERSED )
         n.Reverse();
       N += n.Normalized();
@@ -411,6 +497,7 @@ static bool fixCommonVertexUV (gp_Pnt2d &           theUV,
   // define step size by which to move theUV
 
   gp_Pnt2d nextUV; // uv of next node on edge, most distant of the four
+  gp_Pnt2d thisUV( theUV.x, theUV.y );
   double maxDist = -DBL_MAX;
   TopTools_ListIteratorOfListOfShape aEIt (EList);
   list< double >::iterator aUIt = UList.begin();
@@ -431,10 +518,10 @@ static bool fixCommonVertexUV (gp_Pnt2d &           theUV,
       while ( nIt->more() ) {
         const SMDS_MeshNode* node = nIt->next();
         // check if node is medium
-        if ( CreateQuadratic && SMESH_MesherHelper::IsMedium( node, SMDSAbs_Edge ))
+        if ( theCreateQuadratic && SMESH_MesherHelper::IsMedium( node, SMDSAbs_Edge ))
           continue;
         const SMDS_EdgePosition* epos =
-          static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
+          static_cast<const SMDS_EdgePosition*>(node->GetPosition());
         double u = epos->GetUParameter();
         if ( u < umin )
           umin = u;
@@ -444,29 +531,33 @@ static bool fixCommonVertexUV (gp_Pnt2d &           theUV,
     }
     bool isFirstCommon = ( *aUIt == f );
     gp_Pnt2d uv = C2d->Value( isFirstCommon ? umin : umax );
-    double dist = theUV.SquareDistance( uv );
+    double dist = thisUV.SquareDistance( uv );
     if ( dist > maxDist ) {
       maxDist = dist;
       nextUV  = uv;
     }
   }
   R2 uv0, uv1, uv2;
-  uv0.x = theUV.X();    uv0.y = theUV.Y(); 
+  uv0.x = thisUV.X();   uv0.y = thisUV.Y();
   uv1.x = nextUV.X();   uv1.y = nextUV.Y(); 
-  uv2.x = uv0.x;        uv2.y = uv0.y;
+  uv2.x = thisUV.X();   uv2.y = thisUV.Y();
+
+  uv1.x *= theScaleX;   uv1.y *= theScaleY; 
+
   if ( fixOverlappedLinkUV( uv0, uv1, uv2 ))
   {
-    double step = theUV.Distance( gp_Pnt2d( uv0.x, uv0.y ));
+    double step = thisUV.Distance( gp_Pnt2d( uv0.x, uv0.y ));
 
     // move theUV along the normal by the step
 
     N *= step;
 
-    MESSAGE("--fixCommonVertexUV move(" << theUV.X() << " " << theUV.Y()
+    MESSAGE("--fixCommonVertexUV move(" << theUV.x << " " << theUV.x
             << ") by (" << N.X() << " " << N.Y() << ")" 
             << endl << "--- MAX DIST " << maxDist);
 
-    theUV.SetXY( theUV.XY() + N.XY() );
+    theUV.x += N.X();
+    theUV.y += N.Y();
 
     return true;
   }
@@ -479,113 +570,109 @@ static bool fixCommonVertexUV (gp_Pnt2d &           theUV,
  */
 //=============================================================================
 
-bool StdMeshers_MEFISTO_2D::LoadPoints(SMESH_Mesh &        aMesh,
-                                       const TopoDS_Face & FF,
-                                       const TopoDS_Wire & WW,
-                                       R2 *                uvslf,
-                                       int &               m,
-                                       map<int, const SMDS_MeshNode*>&mefistoToDS,
-                                       double              scalex,
-                                       double              scaley,
-                                       const TopTools_IndexedDataMapOfShapeListOfShape& VWMap)
+bool StdMeshers_MEFISTO_2D::LoadPoints(TWireVector &                 wires,
+                                       R2 *                          uvslf,
+                                       vector<const SMDS_MeshNode*>& mefistoToDS,
+                                       double                        scalex,
+                                       double                        scaley)
 {
-  //  MESSAGE("StdMeshers_MEFISTO_2D::LoadPoints");
-
-  TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
-
-  list< int > mOnVertex;
-
-  gp_XY scale( scalex, scaley );
-
-  TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
-  BRepTools_WireExplorer wexp(W, F);
-  for ( wexp.Init(W, F); wexp.More(); wexp.Next() )
+  // to avoid passing same uv points for a vertex common to 2 wires
+  TopoDS_Face F;
+  TopTools_IndexedDataMapOfShapeListOfShape VWMap;
+  if ( wires.size() > 1 )
   {
-    const TopoDS_Edge & E = wexp.Current();
-    bool isForward = (E.Orientation() == TopAbs_FORWARD);
-
-    // --- IDNodes of first and last Vertex
+    F = TopoDS::Face( _helper->GetSubShape() );
+    TopExp::MapShapesAndAncestors( F, TopAbs_VERTEX, TopAbs_WIRE, VWMap );
+    int nbVertices = 0;
+    for ( int iW = 0; iW < wires.size(); ++iW )
+      nbVertices += wires[ iW ]->NbEdges();
+    if ( nbVertices == VWMap.Extent() )
+      VWMap.Clear(); // wires have no common vertices
+  }
 
-    TopoDS_Vertex VFirst, VLast, V;
-    TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
-    V = isForward ? VFirst : VLast;
+  int m = 0;
 
-    ASSERT(!V.IsNull());
-    SMDS_NodeIteratorPtr lid= aMesh.GetSubMesh(V)->GetSubMeshDS()->GetNodes();
-    if ( !lid->more() ) {
-      MESSAGE (" NO NODE BUILT ON VERTEX ");
-      return false;
+  for ( int iW = 0; iW < wires.size(); ++iW )
+  {
+    const vector<UVPtStruct>& uvPtVec = wires[ iW ]->GetUVPtStruct();
+    if ( uvPtVec.size() != wires[ iW ]->NbPoints() ) {
+      return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Unexpected nb of points on wire ")
+                   << iW << ": " << uvPtVec.size()<<" != "<<wires[ iW ]->NbPoints()
+                   << ", probably because of invalid node parameters on geom edges");
+    }
+    if ( m + uvPtVec.size()-1 > mefistoToDS.size() ) {
+      MESSAGE("Wrong mefistoToDS.size: "<<mefistoToDS.size()<<" < "<<m + uvPtVec.size()-1);
+      return error("Internal error");
     }
-    const SMDS_MeshNode* idFirst = lid->next();
-
-    // --- edge internal IDNodes (relies on good order storage, not checked)
-
-    map<double, const SMDS_MeshNode*> params;
-    const SMDS_MeshNode * node;
 
-    SMDS_NodeIteratorPtr nodeIt= aMesh.GetSubMesh(E)->GetSubMeshDS()->GetNodes();
-    while ( nodeIt->more() )
+    list< int > mOnVertex;
+    vector<UVPtStruct>::const_iterator uvPt = uvPtVec.begin();
+    for ( ++uvPt; uvPt != uvPtVec.end(); ++uvPt )
     {
-      node = nodeIt->next();
-      if ( _quadraticMesh && SMESH_MesherHelper::IsMedium( node, SMDSAbs_Edge ))
-        continue;
-      const SMDS_EdgePosition* epos =
-        static_cast<const SMDS_EdgePosition*>(node->GetPosition().get());
-      double param = epos->GetUParameter();
-      if ( !isForward ) param = -param;
-      if ( !params.insert( make_pair( param, node )).second )
+      // bind mefisto ID to node
+      mefistoToDS[m] = uvPt->node;
+      // set UV
+      uvslf[m].x = uvPt->u * scalex;
+      uvslf[m].y = uvPt->v * scaley;
+      switch ( uvPt->node->GetPosition()->GetTypeOfPosition())
       {
-        MESSAGE( "BAD NODE ON EDGE POSITIONS" );
-        return false;
+      case SMDS_TOP_VERTEX:
+        mOnVertex.push_back( m );
+        break;
+      case SMDS_TOP_EDGE:
+        // In order to detect degenerated faces easily, we replace
+        // nodes on a degenerated edge by node on the vertex of that edge
+        if ( _helper->IsDegenShape( uvPt->node->getshapeId() ))
+        {
+          int edgeID = uvPt->node->getshapeId();
+          SMESH_subMesh* edgeSM = _helper->GetMesh()->GetSubMeshContaining( edgeID );
+          SMESH_subMeshIteratorPtr smIt = edgeSM->getDependsOnIterator( /*includeSelf=*/0,
+                                                                        /*complexShapeFirst=*/0);
+          if ( smIt->more() )
+          {
+            SMESH_subMesh* vertexSM = smIt->next();
+            SMDS_NodeIteratorPtr nIt = vertexSM->GetSubMeshDS()->GetNodes();
+            if ( nIt->more() )
+              mefistoToDS[m] = nIt->next();
+          }
+        }
+        break;
+      default:;
       }
+      m++;
     }
 
-    // --- load 2D values into MEFISTO structure,
-    //     add IDNodes in mefistoToDS map
-
-    double f, l, uFirst, u;
-    Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
-    uFirst = isForward ? f : l;
-
-    // vertex node
-    gp_Pnt2d p = C2d->Value( uFirst ).XY().Multiplied( scale );
-    if ( fixCommonVertexUV( p, V, W, myOuterWire, F, VWMap, aMesh, _quadraticMesh ))
-      myNodesOnCommonV.push_back( idFirst );
-    mOnVertex.push_back( m );
-    uvslf[m].x = p.X();
-    uvslf[m].y = p.Y();
-    mefistoToDS[m + 1] = idFirst;
-    //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
-    //MESSAGE("__ f "<<uFirst<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
-    m++;
-
-    // internal nodes
-    map<double, const SMDS_MeshNode*>::iterator u_n = params.begin();
-    for ( int i = 0; u_n != params.end(); ++u_n, ++i )
+    int mFirst = mOnVertex.front(), mLast = m - 1;
+    list< int >::iterator mIt = mOnVertex.begin();
+    for ( ; mIt != mOnVertex.end(); ++mIt)
     {
-      u = isForward ? u_n->first : - u_n->first;
-      gp_Pnt2d p = C2d->Value( u ).XY().Multiplied( scale );
-      uvslf[m].x = p.X();
-      uvslf[m].y = p.Y();
-      mefistoToDS[m + 1] = u_n->second;
-      //MESSAGE(" "<<m<<" "<<mefistoToDS[m+1]);
-      //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[m].x <<" "<<uvslf[m].y);
-      m++;
+      int m = *mIt;
+      if ( iW && !VWMap.IsEmpty()) { // except outer wire
+        // avoid passing same uv point for a vertex common to 2 wires
+        int vID = mefistoToDS[m]->getshapeId();
+        TopoDS_Vertex V = TopoDS::Vertex( _helper->GetMeshDS()->IndexToShape( vID ));
+        if ( fixCommonVertexUV( uvslf[m], V, F, VWMap, *_helper->GetMesh(),
+                                scalex, scaley, _quadraticMesh )) {
+          myNodesOnCommonV.push_back( mefistoToDS[m] );
+          continue;
+        }
+      }
+      // prevent failure on overlapped adjacent links,
+      // check only links ending in vertex nodes
+      int mB = m - 1, mA = m + 1; // indices Before and After
+      if ( mB < mFirst ) mB = mLast;
+      if ( mA > mLast )  mA = mFirst;
+      fixOverlappedLinkUV (uvslf[ mB ], uvslf[ m ], uvslf[ mA ]);
     }
-  } // for  wexp
-
-  // prevent failure on overlapped adjacent links,
-  // check only links ending in vertex nodes
-  int mFirst = mOnVertex.front(), mLast = m - 1;
-  list< int >::iterator mIt = mOnVertex.begin();
-  for ( ; mIt != mOnVertex.end(); ++mIt ) {
-    int i = *mIt;
-    int iB = i - 1, iA = i + 1; // indices Before and After
-    if ( iB < mFirst ) iB = mLast;
-    if ( iA > mLast )  iA = mFirst;
-    fixOverlappedLinkUV (uvslf[ iB ], uvslf[ i ], uvslf[ iA ]);
   }
 
+#ifdef DUMP_POINTS
+  cout << "MEFISTO INPUT************" << endl;
+  for ( int i =0; i < m; ++i )
+    cout << i << ": \t" << uvslf[i].x << ", " << uvslf[i].y
+         << " Node " << mefistoToDS[i]->GetID()<< endl;
+#endif
+
   return true;
 }
 
@@ -595,12 +682,12 @@ bool StdMeshers_MEFISTO_2D::LoadPoints(SMESH_Mesh &        aMesh,
  */
 //=============================================================================
 
-void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh,
-       const TopoDS_Face & aFace, double &scalex, double &scaley)
+void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh &        aMesh,
+                                               const TopoDS_Face & aFace,
+                                               double &            scalex,
+                                               double &            scaley)
 {
-  //MESSAGE("StdMeshers_MEFISTO_2D::ComputeScaleOnFace");
-  TopoDS_Face F = TopoDS::Face(aFace.Oriented(TopAbs_FORWARD));
-  TopoDS_Wire W = BRepTools::OuterWire(F);
+  TopoDS_Wire W = BRepTools::OuterWire(aFace);
 
   double xmin = 1.e300;         // min & max of face 2D parametric coord.
   double xmax = -1.e300;
@@ -615,7 +702,7 @@ void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh,
   {
     const TopoDS_Edge & E = TopoDS::Edge( wexp.Current() );
     double f, l;
-    Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
+    Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, aFace, f, l);
     if ( C2d.IsNull() ) continue;
     double du = (l - f) / double (nbp);
     for (int i = 0; i <= nbp; i++)
@@ -642,7 +729,8 @@ void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh,
   double xsize = xmax - xmin;
   double ysize = ymax - ymin;
 
-  Handle(Geom_Surface) S = BRep_Tool::Surface(F);       // 3D surface
+  TopLoc_Location L;
+  Handle(Geom_Surface) S = BRep_Tool::Surface(aFace,L);       // 3D surface
 
   double length_x = 0;
   double length_y = 0;
@@ -682,64 +770,85 @@ void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh,
   ASSERT(scaley);
 }
 
+// namespace
+// {
+//   bool isDegenTria( const SMDS_MeshNode * nn[3] )
+//   {
+//     SMESH_TNodeXYZ p1( nn[0] );
+//     SMESH_TNodeXYZ p2( nn[1] );
+//     SMESH_TNodeXYZ p3( nn[2] );
+//     gp_XYZ vec1 = p2 - p1;
+//     gp_XYZ vec2 = p3 - p1;
+//     gp_XYZ cross = vec1 ^ vec2;
+//     const double eps = 1e-100;
+//     return ( fabs( cross.X() ) < eps &&
+//              fabs( cross.Y() ) < eps &&
+//              fabs( cross.Z() ) < eps );
+//   }
+// }
+
 //=============================================================================
 /*!
  *  
  */
 //=============================================================================
 
-void StdMeshers_MEFISTO_2D::StoreResult(SMESH_Mesh & aMesh,
-                                        Z nbst, R2 * uvst, Z nbt, Z * nust,
-                                        const TopoDS_Face & F, bool faceIsForward,
-                                        map<int, const SMDS_MeshNode*>&mefistoToDS,
+void StdMeshers_MEFISTO_2D::StoreResult(Z nbst, R2 * uvst, Z nbt, Z * nust,
+                                        vector< const SMDS_MeshNode*>&mefistoToDS,
                                         double scalex, double scaley)
 {
-  SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
-  int faceID = meshDS->ShapeToIndex( F );
+  _helper->SetElementsOnShape( true );
 
-  Z n, m;
-  Handle(Geom_Surface) S = BRep_Tool::Surface(F);
+  TopoDS_Face F = TopoDS::Face( _helper->GetSubShape() );
+  Handle(Geom_Surface) S = BRep_Tool::Surface( F );
 
-  for (n = 0; n < nbst; n++)
+  //const size_t nbInputNodes = mefistoToDS.size();
+
+  Z n = mefistoToDS.size(); // nb input points
+  mefistoToDS.resize( nbst );
+  for ( ; n < nbst; n++)
   {
-    if (mefistoToDS.find(n + 1) == mefistoToDS.end())
+    if (!mefistoToDS[n])
     {
       double u = uvst[n][0] / scalex;
       double v = uvst[n][1] / scaley;
       gp_Pnt P = S->Value(u, v);
 
-      SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
-      meshDS->SetNodeOnFace(node, faceID, u, v);
-
-      //MESSAGE(P.X()<<" "<<P.Y()<<" "<<P.Z());
-      mefistoToDS[n + 1] = node;
-      //MESSAGE("NEW: "<<n<<" "<<mefistoToDS[n+1]);
+      mefistoToDS[n] = _helper->AddNode( P.X(), P.Y(), P.Z(), 0, u, v );
     }
   }
 
-  m = 0;
+  m = 0;
 
   // triangle points must be in trigonometric order if face is Forward
   // else they must be put clockwise
 
-  bool triangleIsWellOriented = faceIsForward;
+  int i1 = 1, i2 = 2;
+  if ( F.Orientation() != TopAbs_FORWARD )
+    std::swap( i1, i2 );
 
+  const SMDS_MeshNode * nn[3];
   for (n = 1; n <= nbt; n++)
   {
-    const SMDS_MeshNode * n1 = mefistoToDS[ nust[m++] ];
-    const SMDS_MeshNode * n2 = mefistoToDS[ nust[m++] ];
-    const SMDS_MeshNode * n3 = mefistoToDS[ nust[m++] ];
-
-    SMDS_MeshElement * elt;
-    if (triangleIsWellOriented)
-      //elt = meshDS->AddFace(n1, n2, n3);
-      elt = myTool->AddFace(n1, n2, n3);
-    else
-      //elt = meshDS->AddFace(n1, n3, n2);
-      elt = myTool->AddFace(n1, n3, n2);
-
-    meshDS->SetMeshElementOnShape(elt, faceID);
+    // const bool allNodesAreOld = ( nust[m + 0] <= nbInputNodes &&
+    //                               nust[m + 1] <= nbInputNodes &&
+    //                               nust[m + 2] <= nbInputNodes );
+    nn[ 0 ] = mefistoToDS[ nust[m++] - 1 ];
+    nn[ 1 ] = mefistoToDS[ nust[m++] - 1 ];
+    nn[ 2 ] = mefistoToDS[ nust[m++] - 1 ];
     m++;
+
+    // avoid creating degenetrated faces
+    bool isDegen = ( _helper->HasDegeneratedEdges() &&
+                     ( nn[0] == nn[1] || nn[1] == nn[2] || nn[2] == nn[0] ));
+
+    // It was an attemp to fix a problem of a zero area face whose all nodes
+    // are on one staight EDGE. But omitting this face makes a hole in the mesh :(
+    // if ( !isDegen && allNodesAreOld )
+    //   isDegen = isDegenTria( nn );
+
+    if ( !isDegen )
+      _helper->AddFace( nn[0], nn[i1], nn[i2] );
   }
 
   // remove bad elements built on vertices shared by wires
@@ -759,86 +868,8 @@ void StdMeshers_MEFISTO_2D::StoreResult(SMESH_Mesh & aMesh,
           nbSame++;
       if (nbSame > 1) {
         MESSAGE( "RM bad element " << elem->GetID());
-        meshDS->RemoveElement( elem );
+        _helper->GetMeshDS()->RemoveElement( elem );
       }
     }
   }
-
-}
-
-//=============================================================================
-/*!
- *  
- */
-//=============================================================================
-
-double StdMeshers_MEFISTO_2D::ComputeEdgeElementLength(SMESH_Mesh & aMesh,
-       const TopoDS_Shape & aShape)
-{
-       //MESSAGE("StdMeshers_MEFISTO_2D::ComputeEdgeElementLength");
-       // **** a mettre dans SMESH_2D_Algo ?
-
-       //const TopoDS_Face & FF = TopoDS::Face(aShape);
-       //bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
-       //TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
-
-       double meanElementLength = 100;
-       double wireLength = 0;
-       int wireElementsNumber = 0;
-               for (TopExp_Explorer expe(aShape, TopAbs_EDGE); expe.More(); expe.Next())
-               {
-                       const TopoDS_Edge & E = TopoDS::Edge(expe.Current());
-                       int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
-                       double length = EdgeLength(E);
-                       wireLength += length;
-                       wireElementsNumber += nb;
-               }
-       if (wireElementsNumber)
-               meanElementLength = wireLength / wireElementsNumber;
-       //SCRUTE(meanElementLength);
-       return meanElementLength;
-}
-
-//=============================================================================
-/*!
- *  
- */
-//=============================================================================
-
-ostream & StdMeshers_MEFISTO_2D::SaveTo(ostream & save)
-{
-  return save;
-}
-
-//=============================================================================
-/*!
- *  
- */
-//=============================================================================
-
-istream & StdMeshers_MEFISTO_2D::LoadFrom(istream & load)
-{
-  return load;
-}
-
-//=============================================================================
-/*!
- *  
- */
-//=============================================================================
-
-ostream & operator <<(ostream & save, StdMeshers_MEFISTO_2D & hyp)
-{
-  return hyp.SaveTo( save );
-}
-
-//=============================================================================
-/*!
- *  
- */
-//=============================================================================
-
-istream & operator >>(istream & load, StdMeshers_MEFISTO_2D & hyp)
-{
-  return hyp.LoadFrom( load );
 }