src/SMESH/SMESH_MesherHelper.cxx

   1 //  Copyright (C) 2007-2010  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  Copyright (C) 2003-2007  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.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 // File:      SMESH_MesherHelper.cxx
  24 // Created:   15.02.06 15:22:41
  25 // Author:    Sergey KUUL
  26 //
  27 #include "SMESH_MesherHelper.hxx"
  28
  29 #include "SMDS_FacePosition.hxx"
  30 #include "SMDS_EdgePosition.hxx"
  31 #include "SMDS_VolumeTool.hxx"
  32 #include "SMESH_subMesh.hxx"
  33
  34 #include <BRepAdaptor_Surface.hxx>
  35 #include <BRepTools.hxx>
  36 #include <BRepTools_WireExplorer.hxx>
  37 #include <BRep_Tool.hxx>
  38 #include <Geom2d_Curve.hxx>
  39 #include <GeomAPI_ProjectPointOnCurve.hxx>
  40 #include <GeomAPI_ProjectPointOnSurf.hxx>
  41 #include <Geom_Curve.hxx>
  42 #include <Geom_Surface.hxx>
  43 #include <ShapeAnalysis.hxx>
  44 #include <TopExp.hxx>
  45 #include <TopExp_Explorer.hxx>
  46 #include <TopTools_ListIteratorOfListOfShape.hxx>
  47 #include <TopTools_MapIteratorOfMapOfShape.hxx>
  48 #include <TopTools_MapOfShape.hxx>
  49 #include <TopoDS.hxx>
  50 #include <gp_Ax3.hxx>
  51 #include <gp_Pnt2d.hxx>
  52 #include <gp_Trsf.hxx>
  53
  54 #include <Standard_Failure.hxx>
  55 #include <Standard_ErrorHandler.hxx>
  56
  57 #include <utilities.h>
  58
  59 #include <limits>
  60
  61 using namespace std;
  62
  63 #define RETURN_BAD_RESULT(msg) { MESSAGE(msg); return false; }
  64
  65 namespace {
  66
  67   gp_XYZ XYZ(const SMDS_MeshNode* n) { return gp_XYZ(n->X(), n->Y(), n->Z()); }
  68
  69   enum { U_periodic = 1, V_periodic = 2 };
  70 }
  71
  72 //================================================================================
  73 /*!
  74  * \brief Constructor
  75  */
  76 //================================================================================
  77
  78 SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh)
  79   : myParIndex(0), myMesh(&theMesh), myShapeID(0), myCreateQuadratic(false)
  80 {
  81   myPar1[0] = myPar2[0] = myPar1[1] = myPar2[1] = 0;
  82   mySetElemOnShape = ( ! myMesh->HasShapeToMesh() );
  83 }
  84
  85 //=======================================================================
  86 //function : ~SMESH_MesherHelper
  87 //purpose  :
  88 //=======================================================================
  89
  90 SMESH_MesherHelper::~SMESH_MesherHelper()
  91 {
  92   {
  93     TID2ProjectorOnSurf::iterator i_proj = myFace2Projector.begin();
  94     for ( ; i_proj != myFace2Projector.end(); ++i_proj )
  95       delete i_proj->second;
  96   }
  97   {
  98     TID2ProjectorOnCurve::iterator i_proj = myEdge2Projector.begin();
  99     for ( ; i_proj != myEdge2Projector.end(); ++i_proj )
 100       delete i_proj->second;
 101   }
 102 }
 103
 104 //=======================================================================
 105 //function : IsQuadraticSubMesh
 106 //purpose  : Check submesh for given shape: if all elements on this shape
 107 //           are quadratic, quadratic elements will be created.
 108 //           Also fill myTLinkNodeMap
 109 //=======================================================================
 110
 111 bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
 112 {
 113   SMESHDS_Mesh* meshDS = GetMeshDS();
 114   // we can create quadratic elements only if all elements
 115   // created on subshapes of given shape are quadratic
 116   // also we have to fill myTLinkNodeMap
 117   myCreateQuadratic = true;
 118   mySeamShapeIds.clear();
 119   myDegenShapeIds.clear();
 120   TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
 121   SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge );
 122
 123   int nbOldLinks = myTLinkNodeMap.size();
 124
 125   TopExp_Explorer exp( aSh, subType );
 126   for (; exp.More() && myCreateQuadratic; exp.Next()) {
 127     if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) {
 128       if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) {
 129         while(it->more()) {
 130           const SMDS_MeshElement* e = it->next();
 131           if ( e->GetType() != elemType || !e->IsQuadratic() ) {
 132             myCreateQuadratic = false;
 133             break;
 134           }
 135           else {
 136             // fill TLinkNodeMap
 137             switch ( e->NbNodes() ) {
 138             case 3:
 139               AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break;
 140             case 6:
 141               AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3));
 142               AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4));
 143               AddTLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break;
 144             case 8:
 145               AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4));
 146               AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5));
 147               AddTLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6));
 148               AddTLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7));
 149               break;
 150             default:
 151               myCreateQuadratic = false;
 152               break;
 153             }
 154           }
 155         }
 156       }
 157     }
 158   }
 159
 160   if ( nbOldLinks == myTLinkNodeMap.size() )
 161     myCreateQuadratic = false;
 162
 163   if(!myCreateQuadratic) {
 164     myTLinkNodeMap.clear();
 165   }
 166   SetSubShape( aSh );
 167
 168   return myCreateQuadratic;
 169 }
 170
 171 //=======================================================================
 172 //function : SetSubShape
 173 //purpose  : Set geomerty to make elements on
 174 //=======================================================================
 175
 176 void SMESH_MesherHelper::SetSubShape(const int aShID)
 177 {
 178   if ( aShID == myShapeID )
 179     return;
 180   if ( aShID > 1 )
 181     SetSubShape( GetMeshDS()->IndexToShape( aShID ));
 182   else
 183     SetSubShape( TopoDS_Shape() );
 184 }
 185
 186 //=======================================================================
 187 //function : SetSubShape
 188 //purpose  : Set geomerty to create elements on
 189 //=======================================================================
 190
 191 void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
 192 {
 193   if ( myShape.IsSame( aSh ))
 194     return;
 195
 196   myShape = aSh;
 197   mySeamShapeIds.clear();
 198   myDegenShapeIds.clear();
 199
 200   if ( myShape.IsNull() ) {
 201     myShapeID  = 0;
 202     return;
 203   }
 204   SMESHDS_Mesh* meshDS = GetMeshDS();
 205   myShapeID = meshDS->ShapeToIndex(aSh);
 206   myParIndex = 0;
 207
 208   // treatment of periodic faces
 209   for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
 210   {
 211     const TopoDS_Face& face = TopoDS::Face( eF.Current() );
 212     BRepAdaptor_Surface surface( face );
 213     if ( surface.IsUPeriodic() || surface.IsVPeriodic() )
 214     {
 215       for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next())
 216       {
 217         // look for a seam edge
 218         const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
 219         if ( BRep_Tool::IsClosed( edge, face )) {
 220           // initialize myPar1, myPar2 and myParIndex
 221           gp_Pnt2d uv1, uv2;
 222           BRep_Tool::UVPoints( edge, face, uv1, uv2 );
 223           if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) ))
 224           {
 225             myParIndex |= U_periodic;
 226             myPar1[0] = surface.FirstUParameter();
 227             myPar2[0] = surface.LastUParameter();
 228           }
 229           else {
 230             myParIndex |= V_periodic;
 231             myPar1[1] = surface.FirstVParameter();
 232             myPar2[1] = surface.LastVParameter();
 233           }
 234           // store seam shape indices, negative if shape encounters twice
 235           int edgeID = meshDS->ShapeToIndex( edge );
 236           mySeamShapeIds.insert( IsSeamShape( edgeID ) ? -edgeID : edgeID );
 237           for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() ) {
 238             int vertexID = meshDS->ShapeToIndex( v.Current() );
 239             mySeamShapeIds.insert( IsSeamShape( vertexID ) ? -vertexID : vertexID );
 240           }
 241         }
 242
 243         // look for a degenerated edge
 244         if ( BRep_Tool::Degenerated( edge )) {
 245           myDegenShapeIds.insert( meshDS->ShapeToIndex( edge ));
 246           for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() )
 247             myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
 248         }
 249       }
 250     }
 251   }
 252 }
 253
 254 //=======================================================================
 255 //function : GetNodeUVneedInFaceNode
 256 //purpose  : Check if inFaceNode argument is necessary for call GetNodeUV(F,..)
 257 //           Return true if the face is periodic.
 258 //           If F is Null, answer about subshape set through IsQuadraticSubMesh() or
 259 //           * SetSubShape()
 260 //=======================================================================
 261
 262 bool SMESH_MesherHelper::GetNodeUVneedInFaceNode(const TopoDS_Face& F) const
 263 {
 264   if ( F.IsNull() ) return !mySeamShapeIds.empty();
 265
 266   if ( !F.IsNull() && !myShape.IsNull() && myShape.IsSame( F ))
 267     return !mySeamShapeIds.empty();
 268
 269   TopLoc_Location loc;
 270   Handle(Geom_Surface) aSurface = BRep_Tool::Surface( F,loc );
 271   if ( !aSurface.IsNull() )
 272     return ( aSurface->IsUPeriodic() || aSurface->IsVPeriodic() );
 273
 274   return false;
 275 }
 276
 277 //=======================================================================
 278 //function : IsMedium
 279 //purpose  :
 280 //=======================================================================
 281
 282 bool SMESH_MesherHelper::IsMedium(const SMDS_MeshNode*      node,
 283                                   const SMDSAbs_ElementType typeToCheck)
 284 {
 285   return SMESH_MeshEditor::IsMedium( node, typeToCheck );
 286 }
 287
 288 //=======================================================================
 289 //function : GetSubShapeByNode
 290 //purpose  : Return support shape of a node
 291 //=======================================================================
 292
 293 TopoDS_Shape SMESH_MesherHelper::GetSubShapeByNode(const SMDS_MeshNode* node,
 294                                                    SMESHDS_Mesh*        meshDS)
 295 {
 296   int shapeID = node->getshapeId();
 297   if ( 0 < shapeID && shapeID <= meshDS->MaxShapeIndex() )
 298     return meshDS->IndexToShape( shapeID );
 299   else
 300     return TopoDS_Shape();
 301 }
 302
 303
 304 //=======================================================================
 305 //function : AddTLinkNode
 306 //purpose  : add a link in my data structure
 307 //=======================================================================
 308
 309 void SMESH_MesherHelper::AddTLinkNode(const SMDS_MeshNode* n1,
 310                                       const SMDS_MeshNode* n2,
 311                                       const SMDS_MeshNode* n12)
 312 {
 313   // add new record to map
 314   SMESH_TLink link( n1, n2 );
 315   myTLinkNodeMap.insert( make_pair(link,n12));
 316 }
 317
 318 //=======================================================================
 319 //function : GetUVOnSeam
 320 //purpose  : Select UV on either of 2 pcurves of a seam edge, closest to the given UV
 321 //=======================================================================
 322
 323 gp_Pnt2d SMESH_MesherHelper::GetUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& uv2 ) const
 324 {
 325   gp_Pnt2d result = uv1;
 326   for ( int i = U_periodic; i <= V_periodic ; ++i )
 327   {
 328     if ( myParIndex & i )
 329     {
 330       double p1 = uv1.Coord( i );
 331       double dp1 = Abs( p1-myPar1[i-1]), dp2 = Abs( p1-myPar2[i-1]);
 332       if ( myParIndex == i ||
 333            dp1 < ( myPar2[i-1] - myPar2[i-1] ) / 100. ||
 334            dp2 < ( myPar2[i-1] - myPar2[i-1] ) / 100. )
 335       {
 336         double p2 = uv2.Coord( i );
 337         double p1Alt = ( dp1 < dp2 ) ? myPar2[i-1] : myPar1[i-1];
 338         if ( Abs( p2 - p1 ) > Abs( p2 - p1Alt ))
 339           result.SetCoord( i, p1Alt );
 340       }
 341     }
 342   }
 343   return result;
 344 }
 345
 346 //=======================================================================
 347 //function : GetNodeUV
 348 //purpose  : Return node UV on face
 349 //=======================================================================
 350
 351 gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
 352                                     const SMDS_MeshNode* n,
 353                                     const SMDS_MeshNode* n2,
 354                                     bool*                check) const
 355 {
 356   gp_Pnt2d uv( Precision::Infinite(), Precision::Infinite() );
 357   const SMDS_PositionPtr Pos = n->GetPosition();
 358   bool uvOK = false;
 359   if(Pos->GetTypeOfPosition()==SMDS_TOP_FACE)
 360   {
 361     // node has position on face
 362     const SMDS_FacePosition* fpos =
 363       static_cast<const SMDS_FacePosition*>(n->GetPosition());
 364     uv.SetCoord(fpos->GetUParameter(),fpos->GetVParameter());
 365     if ( check )
 366       uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ));
 367   }
 368   else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE)
 369   {
 370     // node has position on edge => it is needed to find
 371     // corresponding edge from face, get pcurve for this
 372     // edge and retrieve value from this pcurve
 373     const SMDS_EdgePosition* epos =
 374       static_cast<const SMDS_EdgePosition*>(n->GetPosition());
 375     int edgeID = n->getshapeId();
 376     TopoDS_Edge E = TopoDS::Edge(GetMeshDS()->IndexToShape(edgeID));
 377     double f, l, u = epos->GetUParameter();
 378     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
 379     bool validU = ( f < u && u < l );
 380     if ( validU )
 381       uv = C2d->Value( u );
 382     else
 383       uv.SetCoord(0.,0.);
 384     if ( check || !validU )
 385       uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ),/*force=*/ !validU );
 386
 387     // for a node on a seam edge select one of UVs on 2 pcurves
 388     if ( n2 && IsSeamShape( edgeID ) )
 389     {
 390       uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0, check ));
 391     }
 392     else
 393     { // adjust uv to period
 394       TopLoc_Location loc;
 395       Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
 396       Standard_Boolean isUPeriodic = S->IsUPeriodic();
 397       Standard_Boolean isVPeriodic = S->IsVPeriodic();
 398       if ( isUPeriodic || isVPeriodic ) {
 399         Standard_Real UF,UL,VF,VL;
 400         S->Bounds(UF,UL,VF,VL);
 401         if(isUPeriodic)
 402           uv.SetX( uv.X() + ShapeAnalysis::AdjustToPeriod(uv.X(),UF,UL));
 403         if(isVPeriodic)
 404           uv.SetY( uv.Y() + ShapeAnalysis::AdjustToPeriod(uv.Y(),VF,VL));
 405       }
 406     }
 407   }
 408   else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX)
 409   {
 410     if ( int vertexID = n->getshapeId() ) {
 411       const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
 412       try {
 413         uv = BRep_Tool::Parameters( V, F );
 414         uvOK = true;
 415       }
 416       catch (Standard_Failure& exc) {
 417       }
 418       if ( !uvOK ) {
 419         for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !uvOK && vert.More(); vert.Next() )
 420           uvOK = ( V == vert.Current() );
 421         if ( !uvOK ) {
 422 #ifdef _DEBUG_
 423           MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID
 424                << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
 425 #endif
 426           // get UV of a vertex closest to the node
 427           double dist = 1e100;
 428           gp_Pnt pn = XYZ( n );
 429           for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !uvOK && vert.More(); vert.Next() ) {
 430             TopoDS_Vertex curV = TopoDS::Vertex( vert.Current() );
 431             gp_Pnt p = BRep_Tool::Pnt( curV );
 432             double curDist = p.SquareDistance( pn );
 433             if ( curDist < dist ) {
 434               dist = curDist;
 435               uv = BRep_Tool::Parameters( curV, F );
 436               uvOK = ( dist < DBL_MIN );
 437             }
 438           }
 439         }
 440         else {
 441           uvOK = false;
 442           TopTools_ListIteratorOfListOfShape it( myMesh->GetAncestors( V ));
 443           for ( ; it.More(); it.Next() ) {
 444             if ( it.Value().ShapeType() == TopAbs_EDGE ) {
 445               const TopoDS_Edge & edge = TopoDS::Edge( it.Value() );
 446               double f,l;
 447               Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(edge, F, f, l);
 448               if ( !C2d.IsNull() ) {
 449                 double u = ( V == TopExp::FirstVertex( edge ) ) ?  f : l;
 450                 uv = C2d->Value( u );
 451                 uvOK = true;
 452                 break;
 453               }
 454             }
 455           }
 456         }
 457       }
 458       if ( n2 && IsSeamShape( vertexID ) )
 459         uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
 460     }
 461   }
 462
 463   if ( check )
 464     *check = uvOK;
 465
 466   return uv.XY();
 467 }
 468
 469 //=======================================================================
 470 //function : CheckNodeUV
 471 //purpose  : Check and fix node UV on a face
 472 //=======================================================================
 473
 474 bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
 475                                      const SMDS_MeshNode* n,
 476                                      gp_XY&               uv,
 477                                      const double         tol,
 478                                      const bool           force) const
 479 {
 480   if ( force || !myOkNodePosShapes.count( n->getshapeId() ))
 481   {
 482     double toldis = tol;
 483     double tolmin = 1.e-7*myMesh->GetMeshDS()->getMaxDim(); // nodes coordinates are stored in float format
 484     if (toldis < tolmin) toldis = tolmin;
 485     // check that uv is correct
 486     TopLoc_Location loc;
 487     Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc );
 488     gp_Pnt nodePnt = XYZ( n );
 489     if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
 490     if ( Precision::IsInfinite( uv.X() ) ||
 491          Precision::IsInfinite( uv.Y() ) ||
 492          nodePnt.Distance( surface->Value( uv.X(), uv.Y() )) > toldis )
 493     {
 494       // uv incorrect, project the node to surface
 495       GeomAPI_ProjectPointOnSurf& projector = GetProjector( F, loc, tol );
 496       projector.Perform( nodePnt );
 497       if ( !projector.IsDone() || projector.NbPoints() < 1 )
 498       {
 499         MESSAGE( "SMESH_MesherHelper::CheckNodeUV() failed to project" );
 500         return false;
 501       }
 502       Quantity_Parameter U,V;
 503       projector.LowerDistanceParameters(U,V);
 504       uv.SetCoord( U,V );
 505       if ( nodePnt.Distance( surface->Value( U, V )) > toldis )
 506       {
 507         MESSAGE( "SMESH_MesherHelper::CheckNodeUV(), invalid projection" );
 508         return false;
 509       }
 510     }
 511     else if ( uv.Modulus() > numeric_limits<double>::min() )
 512     {
 513       ((SMESH_MesherHelper*) this)->myOkNodePosShapes.insert( n->getshapeId() );
 514     }
 515   }
 516   return true;
 517 }
 518
 519 //=======================================================================
 520 //function : GetProjector
 521 //purpose  : Return projector intitialized by given face without location, which is returned
 522 //=======================================================================
 523
 524 GeomAPI_ProjectPointOnSurf& SMESH_MesherHelper::GetProjector(const TopoDS_Face& F,
 525                                                              TopLoc_Location&   loc,
 526                                                              double             tol ) const
 527 {
 528   Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc );
 529   int faceID = GetMeshDS()->ShapeToIndex( F );
 530   TID2ProjectorOnSurf& i2proj = const_cast< TID2ProjectorOnSurf&>( myFace2Projector );
 531   TID2ProjectorOnSurf::iterator i_proj = i2proj.find( faceID );
 532   if ( i_proj == i2proj.end() )
 533   {
 534     if ( tol == 0 ) tol = BRep_Tool::Tolerance( F );
 535     double U1, U2, V1, V2;
 536     surface->Bounds(U1, U2, V1, V2);
 537     GeomAPI_ProjectPointOnSurf* proj = new GeomAPI_ProjectPointOnSurf();
 538     proj->Init( surface, U1, U2, V1, V2, tol );
 539     i_proj = i2proj.insert( make_pair( faceID, proj )).first;
 540   }
 541   return *( i_proj->second );
 542 }
 543
 544 namespace
 545 {
 546   gp_XY AverageUV(const gp_XY& uv1, const gp_XY& uv2) { return ( uv1 + uv2 ) / 2.; }
 547   gp_XY_FunPtr(Added); // define gp_XY_Added pointer to function calling gp_XY::Added(gp_XY)
 548   gp_XY_FunPtr(Subtracted);
 549 }
 550
 551 //=======================================================================
 552 //function : applyIn2D
 553 //purpose  : Perform given operation on two 2d points in parameric space of given surface.
 554 //           It takes into account period of the surface. Use gp_XY_FunPtr macro
 555 //           to easily define pointer to function of gp_XY class.
 556 //=======================================================================
 557
 558 gp_XY SMESH_MesherHelper::applyIn2D(const Handle(Geom_Surface)& surface,
 559                                     const gp_XY&                uv1,
 560                                     const gp_XY&                uv2,
 561                                     xyFunPtr                    fun,
 562                                     const bool                  resultInPeriod)
 563 {
 564   Standard_Boolean isUPeriodic = surface.IsNull() ? false : surface->IsUPeriodic();
 565   Standard_Boolean isVPeriodic = surface.IsNull() ? false : surface->IsVPeriodic();
 566   if ( !isUPeriodic && !isVPeriodic )
 567     return fun(uv1,uv2);
 568
 569   // move uv2 not far than half-period from uv1
 570   double u2 =
 571     uv2.X()+(isUPeriodic ? ShapeAnalysis::AdjustByPeriod(uv2.X(),uv1.X(),surface->UPeriod()) :0);
 572   double v2 =
 573     uv2.Y()+(isVPeriodic ? ShapeAnalysis::AdjustByPeriod(uv2.Y(),uv1.Y(),surface->VPeriod()) :0);
 574
 575   // execute operation
 576   gp_XY res = fun( uv1, gp_XY(u2,v2) );
 577
 578   // move result within period
 579   if ( resultInPeriod )
 580   {
 581     Standard_Real UF,UL,VF,VL;
 582     surface->Bounds(UF,UL,VF,VL);
 583     if ( isUPeriodic )
 584       res.SetX( res.X() + ShapeAnalysis::AdjustToPeriod(res.X(),UF,UL));
 585     if ( isVPeriodic )
 586       res.SetY( res.Y() + ShapeAnalysis::AdjustToPeriod(res.Y(),VF,VL));
 587   }
 588
 589   return res;
 590 }
 591 //=======================================================================
 592 //function : GetMiddleUV
 593 //purpose  : Return middle UV taking in account surface period
 594 //=======================================================================
 595
 596 gp_XY SMESH_MesherHelper::GetMiddleUV(const Handle(Geom_Surface)& surface,
 597                                       const gp_XY&                p1,
 598                                       const gp_XY&                p2)
 599 {
 600   return applyIn2D( surface, p1, p2, & AverageUV );
 601 }
 602
 603 //=======================================================================
 604 //function : GetNodeU
 605 //purpose  : Return node U on edge
 606 //=======================================================================
 607
 608 double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
 609                                     const SMDS_MeshNode* n,
 610                                     const SMDS_MeshNode* inEdgeNode,
 611                                     bool*                check)
 612 {
 613   double param = 0;
 614   const SMDS_PositionPtr pos = n->GetPosition();
 615   if ( pos->GetTypeOfPosition()==SMDS_TOP_EDGE )
 616   {
 617     const SMDS_EdgePosition* epos = static_cast<const SMDS_EdgePosition*>( pos );
 618     param =  epos->GetUParameter();
 619   }
 620   else if( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
 621   {
 622     if ( inEdgeNode && TopExp::FirstVertex( E ).IsSame( TopExp::LastVertex( E ))) // issue 0020128
 623     {
 624       Standard_Real f,l;
 625       BRep_Tool::Range( E, f,l );
 626       double uInEdge = GetNodeU( E, inEdgeNode );
 627       param = ( fabs( uInEdge - f ) < fabs( l - uInEdge )) ? f : l;
 628     }
 629     else
 630     {
 631       SMESHDS_Mesh * meshDS = GetMeshDS();
 632       int vertexID = n->getshapeId();
 633       const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
 634       param =  BRep_Tool::Parameter( V, E );
 635     }
 636   }
 637   if ( check )
 638   {
 639     double tol = BRep_Tool::Tolerance( E );
 640     double f,l;  BRep_Tool::Range( E, f,l );
 641     bool force = ( param < f-tol || param > l+tol );
 642     if ( !force && pos->GetTypeOfPosition()==SMDS_TOP_EDGE )
 643       force = ( GetMeshDS()->ShapeToIndex( E ) != n->getshapeId() );
 644
 645     *check = CheckNodeU( E, n, param, 2*tol, force );
 646   }
 647   return param;
 648 }
 649
 650 //=======================================================================
 651 //function : CheckNodeU
 652 //purpose  : Check and fix node U on an edge
 653 //           Return false if U is bad and could not be fixed
 654 //=======================================================================
 655
 656 bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
 657                                     const SMDS_MeshNode* n,
 658                                     double&              u,
 659                                     const double         tol,
 660                                     const bool           force,
 661                                     double*              distance) const
 662 {
 663   if ( force || !myOkNodePosShapes.count( n->getshapeId() ))
 664   {
 665     double toldis = tol;
 666     double tolmin = 1.e-7*myMesh->GetMeshDS()->getMaxDim(); // nodes coordinates are stored in float format
 667     if (toldis < tolmin) toldis = tolmin;
 668     // check that u is correct
 669     TopLoc_Location loc; double f,l;
 670     Handle(Geom_Curve) curve = BRep_Tool::Curve( E,loc,f,l );
 671     if ( curve.IsNull() ) // degenerated edge
 672     {
 673       if ( u+tol < f || u-tol > l )
 674       {
 675         double r = Max( 0.5, 1 - tol*n->GetID()); // to get a unique u on edge
 676         u =  f*r + l*(1-r);
 677       }
 678     }
 679     else
 680     {
 681       gp_Pnt nodePnt = SMESH_MeshEditor::TNodeXYZ( n );
 682       if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
 683       double dist = nodePnt.Distance( curve->Value( u ));
 684       if ( distance ) *distance = dist;
 685       if ( dist > toldis )
 686       {
 687         // u incorrect, project the node to the curve
 688         int edgeID = GetMeshDS()->ShapeToIndex( E );
 689         TID2ProjectorOnCurve& i2proj = const_cast< TID2ProjectorOnCurve&>( myEdge2Projector );
 690         TID2ProjectorOnCurve::iterator i_proj =
 691           i2proj.insert( make_pair( edgeID, (GeomAPI_ProjectPointOnCurve*) 0 )).first;
 692         if ( !i_proj->second  )
 693         {
 694           i_proj->second = new GeomAPI_ProjectPointOnCurve();
 695           i_proj->second->Init( curve, f, l );
 696         }
 697         GeomAPI_ProjectPointOnCurve* projector = i_proj->second;
 698         projector->Perform( nodePnt );
 699         if ( projector->NbPoints() < 1 )
 700         {
 701           MESSAGE( "SMESH_MesherHelper::CheckNodeU() failed to project" );
 702           return false;
 703         }
 704         Quantity_Parameter U = projector->LowerDistanceParameter();
 705         u = double( U );
 706         dist = nodePnt.Distance( curve->Value( U ));
 707         if ( distance ) *distance = dist;
 708         if ( dist > toldis )
 709         {
 710           MESSAGE( "SMESH_MesherHelper::CheckNodeU(), invalid projection" );
 711           MESSAGE("distance " << nodePnt.Distance(curve->Value( U )) << " " << toldis);
 712           return false;
 713         }
 714         //u = double( U );
 715       }
 716       else if ( fabs( u ) > numeric_limits<double>::min() )
 717       {
 718         ((SMESH_MesherHelper*) this)->myOkNodePosShapes.insert( n->getshapeId() );
 719       }
 720       if (( u < f-tol || u > l+tol ) && force )
 721       {
 722         // node is on vertex but is set on periodic but trimmed edge (issue 0020890)
 723         try
 724         {
 725           // do not use IsPeriodic() as Geom_TrimmedCurve::IsPeriodic () returns false
 726           double period = curve->Period();
 727           u = ( u < f ) ? u + period : u - period;
 728         }
 729         catch (Standard_Failure& exc)
 730         {
 731           return false;
 732         }
 733       }
 734     }
 735   }
 736   return true;
 737 }
 738
 739 //=======================================================================
 740 //function : GetMediumNode
 741 //purpose  : Return existing or create new medium nodes between given ones
 742 //=======================================================================
 743
 744 const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
 745                                                        const SMDS_MeshNode* n2,
 746                                                        bool                 force3d)
 747 {
 748   // Find existing node
 749
 750   SMESH_TLink link(n1,n2);
 751   ItTLinkNode itLN = myTLinkNodeMap.find( link );
 752   if ( itLN != myTLinkNodeMap.end() ) {
 753     return (*itLN).second;
 754   }
 755
 756   // Create medium node
 757
 758   SMDS_MeshNode* n12;
 759   SMESHDS_Mesh* meshDS = GetMeshDS();
 760
 761   // get type of shape for the new medium node
 762   int faceID = -1, edgeID = -1;
 763   const SMDS_PositionPtr Pos1 = n1->GetPosition();
 764   const SMDS_PositionPtr Pos2 = n2->GetPosition();
 765
 766   TopoDS_Edge E; double u [2];
 767   TopoDS_Face F; gp_XY  uv[2];
 768   bool uvOK[2] = { false, false };
 769
 770   if( myShape.IsNull() )
 771   {
 772     if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) {
 773       faceID = n1->getshapeId();
 774     }
 775     else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) {
 776       faceID = n2->getshapeId();
 777     }
 778
 779     if( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
 780       edgeID = n1->getshapeId();
 781     }
 782     if( Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
 783       edgeID = n2->getshapeId();
 784     }
 785   }
 786   // get positions of the given nodes on shapes
 787   TopAbs_ShapeEnum shapeType = myShape.IsNull() ? TopAbs_SHAPE : myShape.ShapeType();
 788   if ( faceID>0 || shapeType == TopAbs_FACE)
 789   {
 790     if( myShape.IsNull() )
 791       F = TopoDS::Face(meshDS->IndexToShape(faceID));
 792     else {
 793       F = TopoDS::Face(myShape);
 794       faceID = myShapeID;
 795     }
 796     uv[0] = GetNodeUV(F,n1,n2, force3d ? 0 : &uvOK[0]);
 797     uv[1] = GetNodeUV(F,n2,n1, force3d ? 0 : &uvOK[1]);
 798   }
 799   else if (edgeID>0 || shapeType == TopAbs_EDGE)
 800   {
 801     if ( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE &&
 802          Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE &&
 803          n1->getshapeId() != n2->getshapeId() ) // issue 0021006
 804     return getMediumNodeOnComposedWire(n1,n2,force3d);
 805
 806     if( myShape.IsNull() )
 807       E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
 808     else {
 809       E = TopoDS::Edge(myShape);
 810       edgeID = myShapeID;
 811     }
 812     u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]);
 813     u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]);
 814   }
 815   if(!force3d)
 816   {
 817     // we try to create medium node using UV parameters of
 818     // nodes, else - medium between corresponding 3d points
 819     if( ! F.IsNull() )
 820     {
 821       if ( uvOK[0] && uvOK[1] )
 822       {
 823         if ( IsDegenShape( n1->getshapeId() ))
 824           if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 ));
 825           else                           uv[0].SetCoord( 2, uv[1].Coord( 2 ));
 826         else if ( IsDegenShape( n2->getshapeId() ))
 827           if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 ));
 828           else                           uv[1].SetCoord( 2, uv[0].Coord( 2 ));
 829
 830         TopLoc_Location loc;
 831         Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
 832         gp_XY UV = GetMiddleUV( S, uv[0], uv[1] );
 833         gp_Pnt P = S->Value( UV.X(), UV.Y() ).Transformed(loc);
 834         n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
 835         meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y());
 836         myTLinkNodeMap.insert(make_pair(link,n12));
 837         return n12;
 838       }
 839     }
 840     else if ( !E.IsNull() )
 841     {
 842       double f,l;
 843       Handle(Geom_Curve) C = BRep_Tool::Curve(E, f, l);
 844       if(!C.IsNull())
 845       {
 846         Standard_Boolean isPeriodic = C->IsPeriodic();
 847         double U;
 848         if(isPeriodic) {
 849           Standard_Real Period = C->Period();
 850           Standard_Real p = u[1]+ShapeAnalysis::AdjustByPeriod(u[1],u[0],Period);
 851           Standard_Real pmid = (u[0]+p)/2.;
 852           U = pmid+ShapeAnalysis::AdjustToPeriod(pmid,C->FirstParameter(),C->LastParameter());
 853         }
 854         else
 855           U = (u[0]+u[1])/2.;
 856
 857         gp_Pnt P = C->Value( U );
 858         n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
 859         meshDS->SetNodeOnEdge(n12, edgeID, U);
 860         myTLinkNodeMap.insert(make_pair(link,n12));
 861         return n12;
 862       }
 863     }
 864   }
 865
 866   // 3d variant
 867   double x = ( n1->X() + n2->X() )/2.;
 868   double y = ( n1->Y() + n2->Y() )/2.;
 869   double z = ( n1->Z() + n2->Z() )/2.;
 870   n12 = meshDS->AddNode(x,y,z);
 871
 872   if ( !F.IsNull() )
 873   {
 874     gp_XY UV = ( uv[0] + uv[1] ) / 2.;
 875     CheckNodeUV( F, n12, UV, 2*BRep_Tool::Tolerance( F ), /*force=*/true);
 876     meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y() );
 877   }
 878   else if ( !E.IsNull() )
 879   {
 880     double U = ( u[0] + u[1] ) / 2.;
 881     CheckNodeU( E, n12, U, 2*BRep_Tool::Tolerance( E ), /*force=*/true);
 882     meshDS->SetNodeOnEdge(n12, edgeID, U);
 883   }
 884   else if ( myShapeID > 0 )
 885   {
 886     meshDS->SetNodeInVolume(n12, myShapeID);
 887   }
 888
 889   myTLinkNodeMap.insert( make_pair( link, n12 ));
 890   return n12;
 891 }
 892
 893 //================================================================================
 894 /*!
 895  * \brief Makes a medium node if nodes reside different edges
 896  */
 897 //================================================================================
 898
 899 const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_MeshNode* n1,
 900                                                                      const SMDS_MeshNode* n2,
 901                                                                      bool                 force3d)
 902 {
 903   gp_Pnt middle = 0.5 * XYZ(n1) + 0.5 * XYZ(n2);
 904   SMDS_MeshNode* n12 = AddNode( middle.X(), middle.Y(), middle.Z() );
 905
 906   // To find position on edge and 3D position for n12,
 907   // project <middle> to 2 edges and select projection most close to <middle>
 908
 909   double u = 0, distMiddleProj = Precision::Infinite();
 910   int iOkEdge = 0;
 911   TopoDS_Edge edges[2];
 912   for ( int is2nd = 0; is2nd < 2; ++is2nd )
 913   {
 914     // get an edge
 915     const SMDS_MeshNode* n = is2nd ? n2 : n1;
 916     TopoDS_Shape shape = GetSubShapeByNode( n, GetMeshDS() );
 917     if ( shape.IsNull() || shape.ShapeType() != TopAbs_EDGE )
 918       continue;
 919
 920     // project to get U of projection and distance from middle to projection
 921     TopoDS_Edge edge = edges[ is2nd ] = TopoDS::Edge( shape );
 922     double node2MiddleDist = middle.Distance( XYZ(n) );
 923     double foundU = GetNodeU( edge, n ), foundDist = node2MiddleDist;
 924     CheckNodeU( edge, n12, foundU, 2*BRep_Tool::Tolerance(edge), /*force=*/true, &foundDist );
 925     if ( foundDist < node2MiddleDist )
 926     {
 927       distMiddleProj = foundDist;
 928       u = foundU;
 929       iOkEdge = is2nd;
 930     }
 931   }
 932   if ( Precision::IsInfinite( distMiddleProj ))
 933   {
 934     // both projections failed; set n12 on the edge of n1 with U of a common vertex
 935     TopoDS_Vertex vCommon;
 936     if ( TopExp::CommonVertex( edges[0], edges[1], vCommon ))
 937       u = BRep_Tool::Parameter( vCommon, edges[0] );
 938     else
 939     {
 940       double f,l, u0 = GetNodeU( edges[0], n1 );
 941       BRep_Tool::Range( edges[0],f,l );
 942       u = ( fabs(u0-f) < fabs(u0-l) ) ? f : l;
 943     }
 944     iOkEdge = 0;
 945     distMiddleProj = 0;
 946   }
 947
 948   // move n12 to position of a successfull projection
 949   double tol = BRep_Tool::Tolerance(edges[ iOkEdge ]);
 950   if ( !force3d && distMiddleProj > 2*tol )
 951   {
 952     TopLoc_Location loc; double f,l;
 953     Handle(Geom_Curve) curve = BRep_Tool::Curve( edges[iOkEdge],loc,f,l );
 954     gp_Pnt p = curve->Value( u );
 955     GetMeshDS()->MoveNode( n12, p.X(), p.Y(), p.Z() );
 956   }
 957
 958   GetMeshDS()->SetNodeOnEdge(n12, edges[iOkEdge], u);
 959
 960   myTLinkNodeMap.insert( make_pair( SMESH_TLink(n1,n2), n12 ));
 961
 962   return n12;
 963 }
 964
 965 //=======================================================================
 966 //function : AddNode
 967 //purpose  : Creates a node
 968 //=======================================================================
 969
 970 SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID)
 971 {
 972   SMESHDS_Mesh * meshDS = GetMeshDS();
 973   SMDS_MeshNode* node = 0;
 974   if ( ID )
 975     node = meshDS->AddNodeWithID( x, y, z, ID );
 976   else
 977     node = meshDS->AddNode( x, y, z );
 978   if ( mySetElemOnShape && myShapeID > 0 ) {
 979     switch ( myShape.ShapeType() ) {
 980     case TopAbs_SOLID:  meshDS->SetNodeInVolume( node, myShapeID); break;
 981     case TopAbs_SHELL:  meshDS->SetNodeInVolume( node, myShapeID); break;
 982     case TopAbs_FACE:   meshDS->SetNodeOnFace(   node, myShapeID); break;
 983     case TopAbs_EDGE:   meshDS->SetNodeOnEdge(   node, myShapeID); break;
 984     case TopAbs_VERTEX: meshDS->SetNodeOnVertex( node, myShapeID); break;
 985     default: ;
 986     }
 987   }
 988   return node;
 989 }
 990
 991 //=======================================================================
 992 //function : AddEdge
 993 //purpose  : Creates quadratic or linear edge
 994 //=======================================================================
 995
 996 SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1,
 997                                            const SMDS_MeshNode* n2,
 998                                            const int            id,
 999                                            const bool           force3d)
1000 {
1001   SMESHDS_Mesh * meshDS = GetMeshDS();
1002
1003   SMDS_MeshEdge* edge = 0;
1004   if (myCreateQuadratic) {
1005     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
1006     if(id)
1007       edge = meshDS->AddEdgeWithID(n1, n2, n12, id);
1008     else
1009       edge = meshDS->AddEdge(n1, n2, n12);
1010   }
1011   else {
1012     if(id)
1013       edge = meshDS->AddEdgeWithID(n1, n2, id);
1014     else
1015       edge = meshDS->AddEdge(n1, n2);
1016   }
1017
1018   if ( mySetElemOnShape && myShapeID > 0 )
1019     meshDS->SetMeshElementOnShape( edge, myShapeID );
1020
1021   return edge;
1022 }
1023
1024 //=======================================================================
1025 //function : AddFace
1026 //purpose  : Creates quadratic or linear triangle
1027 //=======================================================================
1028
1029 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
1030                                            const SMDS_MeshNode* n2,
1031                                            const SMDS_MeshNode* n3,
1032                                            const int id,
1033                                            const bool force3d)
1034 {
1035   SMESHDS_Mesh * meshDS = GetMeshDS();
1036   SMDS_MeshFace* elem = 0;
1037
1038   if( n1==n2 || n2==n3 || n3==n1 )
1039     return elem;
1040
1041   if(!myCreateQuadratic) {
1042     if(id)
1043       elem = meshDS->AddFaceWithID(n1, n2, n3, id);
1044     else
1045       elem = meshDS->AddFace(n1, n2, n3);
1046   }
1047   else {
1048     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
1049     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
1050     const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
1051
1052     if(id)
1053       elem = meshDS->AddFaceWithID(n1, n2, n3, n12, n23, n31, id);
1054     else
1055       elem = meshDS->AddFace(n1, n2, n3, n12, n23, n31);
1056   }
1057   if ( mySetElemOnShape && myShapeID > 0 )
1058     meshDS->SetMeshElementOnShape( elem, myShapeID );
1059
1060   return elem;
1061 }
1062
1063 //=======================================================================
1064 //function : AddFace
1065 //purpose  : Creates quadratic or linear quadrangle
1066 //=======================================================================
1067
1068 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
1069                                            const SMDS_MeshNode* n2,
1070                                            const SMDS_MeshNode* n3,
1071                                            const SMDS_MeshNode* n4,
1072                                            const int            id,
1073                                            const bool           force3d)
1074 {
1075   SMESHDS_Mesh * meshDS = GetMeshDS();
1076   SMDS_MeshFace* elem = 0;
1077
1078   if( n1==n2 ) {
1079     return AddFace(n1,n3,n4,id,force3d);
1080   }
1081   if( n1==n3 ) {
1082     return AddFace(n1,n2,n4,id,force3d);
1083   }
1084   if( n1==n4 ) {
1085     return AddFace(n1,n2,n3,id,force3d);
1086   }
1087   if( n2==n3 ) {
1088     return AddFace(n1,n2,n4,id,force3d);
1089   }
1090   if( n2==n4 ) {
1091     return AddFace(n1,n2,n3,id,force3d);
1092   }
1093   if( n3==n4 ) {
1094     return AddFace(n1,n2,n3,id,force3d);
1095   }
1096
1097   if(!myCreateQuadratic) {
1098     if(id)
1099       elem = meshDS->AddFaceWithID(n1, n2, n3, n4, id);
1100     else
1101       elem = meshDS->AddFace(n1, n2, n3, n4);
1102   }
1103   else {
1104     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
1105     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
1106     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
1107     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
1108
1109     if(id)
1110       elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, id);
1111     else
1112       elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
1113   }
1114   if ( mySetElemOnShape && myShapeID > 0 )
1115     meshDS->SetMeshElementOnShape( elem, myShapeID );
1116
1117   return elem;
1118 }
1119
1120 //=======================================================================
1121 //function : AddPolygonalFace
1122 //purpose  : Creates polygon, with additional nodes in quadratic mesh
1123 //=======================================================================
1124
1125 SMDS_MeshFace* SMESH_MesherHelper::AddPolygonalFace (const vector<const SMDS_MeshNode*>& nodes,
1126                                                      const int                           id,
1127                                                      const bool                          force3d)
1128 {
1129   SMESHDS_Mesh * meshDS = GetMeshDS();
1130   SMDS_MeshFace* elem = 0;
1131
1132   if(!myCreateQuadratic) {
1133     if(id)
1134       elem = meshDS->AddPolygonalFaceWithID(nodes, id);
1135     else
1136       elem = meshDS->AddPolygonalFace(nodes);
1137   }
1138   else {
1139     vector<const SMDS_MeshNode*> newNodes;
1140     for ( int i = 0; i < nodes.size(); ++i )
1141     {
1142       const SMDS_MeshNode* n1 = nodes[i];
1143       const SMDS_MeshNode* n2 = nodes[(i+1)/nodes.size()];
1144       const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
1145       newNodes.push_back( n1 );
1146       newNodes.push_back( n12 );
1147     }
1148     if(id)
1149       elem = meshDS->AddPolygonalFaceWithID(newNodes, id);
1150     else
1151       elem = meshDS->AddPolygonalFace(newNodes);
1152   }
1153   if ( mySetElemOnShape && myShapeID > 0 )
1154     meshDS->SetMeshElementOnShape( elem, myShapeID );
1155
1156   return elem;
1157 }
1158
1159 //=======================================================================
1160 //function : AddVolume
1161 //purpose  : Creates quadratic or linear prism
1162 //=======================================================================
1163
1164 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
1165                                                const SMDS_MeshNode* n2,
1166                                                const SMDS_MeshNode* n3,
1167                                                const SMDS_MeshNode* n4,
1168                                                const SMDS_MeshNode* n5,
1169                                                const SMDS_MeshNode* n6,
1170                                                const int id,
1171                                                const bool force3d)
1172 {
1173   SMESHDS_Mesh * meshDS = GetMeshDS();
1174   SMDS_MeshVolume* elem = 0;
1175   if(!myCreateQuadratic) {
1176     if(id)
1177       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, id);
1178     else
1179       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6);
1180   }
1181   else {
1182     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
1183     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
1184     const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
1185
1186     const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
1187     const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
1188     const SMDS_MeshNode* n64 = GetMediumNode(n6,n4,force3d);
1189
1190     const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
1191     const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
1192     const SMDS_MeshNode* n36 = GetMediumNode(n3,n6,force3d);
1193
1194     if(id)
1195       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
1196                                      n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
1197     else
1198       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
1199                                n12, n23, n31, n45, n56, n64, n14, n25, n36);
1200   }
1201   if ( mySetElemOnShape && myShapeID > 0 )
1202     meshDS->SetMeshElementOnShape( elem, myShapeID );
1203
1204   return elem;
1205 }
1206
1207 //=======================================================================
1208 //function : AddVolume
1209 //purpose  : Creates quadratic or linear tetrahedron
1210 //=======================================================================
1211
1212 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
1213                                                const SMDS_MeshNode* n2,
1214                                                const SMDS_MeshNode* n3,
1215                                                const SMDS_MeshNode* n4,
1216                                                const int id,
1217                                                const bool force3d)
1218 {
1219   SMESHDS_Mesh * meshDS = GetMeshDS();
1220   SMDS_MeshVolume* elem = 0;
1221   if(!myCreateQuadratic) {
1222     if(id)
1223       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, id);
1224     else
1225       elem = meshDS->AddVolume(n1, n2, n3, n4);
1226   }
1227   else {
1228     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
1229     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
1230     const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
1231
1232     const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
1233     const SMDS_MeshNode* n24 = GetMediumNode(n2,n4,force3d);
1234     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
1235
1236     if(id)
1237       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34, id);
1238     else
1239       elem = meshDS->AddVolume(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34);
1240   }
1241   if ( mySetElemOnShape && myShapeID > 0 )
1242     meshDS->SetMeshElementOnShape( elem, myShapeID );
1243
1244   return elem;
1245 }
1246
1247 //=======================================================================
1248 //function : AddVolume
1249 //purpose  : Creates quadratic or linear pyramid
1250 //=======================================================================
1251
1252 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
1253                                                const SMDS_MeshNode* n2,
1254                                                const SMDS_MeshNode* n3,
1255                                                const SMDS_MeshNode* n4,
1256                                                const SMDS_MeshNode* n5,
1257                                                const int id,
1258                                                const bool force3d)
1259 {
1260   SMDS_MeshVolume* elem = 0;
1261   if(!myCreateQuadratic) {
1262     if(id)
1263       elem = GetMeshDS()->AddVolumeWithID(n1, n2, n3, n4, n5, id);
1264     else
1265       elem = GetMeshDS()->AddVolume(n1, n2, n3, n4, n5);
1266   }
1267   else {
1268     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
1269     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
1270     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
1271     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
1272
1273     const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
1274     const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
1275     const SMDS_MeshNode* n35 = GetMediumNode(n3,n5,force3d);
1276     const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
1277
1278     if(id)
1279       elem = GetMeshDS()->AddVolumeWithID ( n1,  n2,  n3,  n4,  n5,
1280                                             n12, n23, n34, n41,
1281                                             n15, n25, n35, n45,
1282                                             id);
1283     else
1284       elem = GetMeshDS()->AddVolume( n1,  n2,  n3,  n4,  n5,
1285                                      n12, n23, n34, n41,
1286                                      n15, n25, n35, n45);
1287   }
1288   if ( mySetElemOnShape && myShapeID > 0 )
1289     GetMeshDS()->SetMeshElementOnShape( elem, myShapeID );
1290
1291   return elem;
1292 }
1293
1294 //=======================================================================
1295 //function : AddVolume
1296 //purpose  : Creates quadratic or linear hexahedron
1297 //=======================================================================
1298
1299 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
1300                                                const SMDS_MeshNode* n2,
1301                                                const SMDS_MeshNode* n3,
1302                                                const SMDS_MeshNode* n4,
1303                                                const SMDS_MeshNode* n5,
1304                                                const SMDS_MeshNode* n6,
1305                                                const SMDS_MeshNode* n7,
1306                                                const SMDS_MeshNode* n8,
1307                                                const int id,
1308                                                const bool force3d)
1309 {
1310   SMESHDS_Mesh * meshDS = GetMeshDS();
1311   SMDS_MeshVolume* elem = 0;
1312   if(!myCreateQuadratic) {
1313     if(id)
1314       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8, id);
1315     else
1316       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8);
1317   }
1318   else {
1319     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
1320     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
1321     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
1322     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
1323
1324     const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
1325     const SMDS_MeshNode* n67 = GetMediumNode(n6,n7,force3d);
1326     const SMDS_MeshNode* n78 = GetMediumNode(n7,n8,force3d);
1327     const SMDS_MeshNode* n85 = GetMediumNode(n8,n5,force3d);
1328
1329     const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
1330     const SMDS_MeshNode* n26 = GetMediumNode(n2,n6,force3d);
1331     const SMDS_MeshNode* n37 = GetMediumNode(n3,n7,force3d);
1332     const SMDS_MeshNode* n48 = GetMediumNode(n4,n8,force3d);
1333
1334     if(id)
1335       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
1336                                      n12, n23, n34, n41, n56, n67,
1337                                      n78, n85, n15, n26, n37, n48, id);
1338     else
1339       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
1340                                n12, n23, n34, n41, n56, n67,
1341                                n78, n85, n15, n26, n37, n48);
1342   }
1343   if ( mySetElemOnShape && myShapeID > 0 )
1344     meshDS->SetMeshElementOnShape( elem, myShapeID );
1345
1346   return elem;
1347 }
1348
1349 //=======================================================================
1350 //function : AddPolyhedralVolume
1351 //purpose  : Creates polyhedron. In quadratic mesh, adds medium nodes
1352 //=======================================================================
1353
1354 SMDS_MeshVolume*
1355 SMESH_MesherHelper::AddPolyhedralVolume (const std::vector<const SMDS_MeshNode*>& nodes,
1356                                          const std::vector<int>&                  quantities,
1357                                          const int                                id,
1358                                          const bool                               force3d)
1359 {
1360   SMESHDS_Mesh * meshDS = GetMeshDS();
1361   SMDS_MeshVolume* elem = 0;
1362   if(!myCreateQuadratic)
1363   {
1364     if(id)
1365       elem = meshDS->AddPolyhedralVolumeWithID(nodes, quantities, id);
1366     else
1367       elem = meshDS->AddPolyhedralVolume(nodes, quantities);
1368   }
1369   else
1370   {
1371     vector<const SMDS_MeshNode*> newNodes;
1372     vector<int> newQuantities;
1373     for ( int iFace=0, iN=0; iFace < quantities.size(); ++iFace)
1374     {
1375       int nbNodesInFace = quantities[iFace];
1376       newQuantities.push_back(0);
1377       for ( int i = 0; i < nbNodesInFace; ++i )
1378       {
1379         const SMDS_MeshNode* n1 = nodes[ iN + i ];
1380         newNodes.push_back( n1 );
1381         newQuantities.back()++;
1382
1383         const SMDS_MeshNode* n2 = nodes[ iN + ( i+1==nbNodesInFace ? 0 : i+1 )];
1384 //         if ( n1->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE &&
1385 //              n2->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE )
1386         {
1387           const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
1388           newNodes.push_back( n12 );
1389           newQuantities.back()++;
1390         }
1391       }
1392       iN += nbNodesInFace;
1393     }
1394     if(id)
1395       elem = meshDS->AddPolyhedralVolumeWithID( newNodes, newQuantities, id );
1396     else
1397       elem = meshDS->AddPolyhedralVolume( newNodes, newQuantities );
1398   }
1399   if ( mySetElemOnShape && myShapeID > 0 )
1400     meshDS->SetMeshElementOnShape( elem, myShapeID );
1401
1402   return elem;
1403 }
1404
1405 //=======================================================================
1406 //function : LoadNodeColumns
1407 //purpose  : Load nodes bound to face into a map of node columns
1408 //=======================================================================
1409
1410 bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
1411                                          const TopoDS_Face& theFace,
1412                                          const TopoDS_Edge& theBaseEdge,
1413                                          SMESHDS_Mesh*      theMesh)
1414 {
1415   SMESHDS_SubMesh* faceSubMesh = theMesh->MeshElements( theFace );
1416   if ( !faceSubMesh || faceSubMesh->NbElements() == 0 )
1417     return false;
1418
1419   // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them
1420
1421   map< double, const SMDS_MeshNode*> sortedBaseNodes;
1422   if ( !SMESH_Algo::GetSortedNodesOnEdge( theMesh, theBaseEdge,/*noMedium=*/true, sortedBaseNodes)
1423        || sortedBaseNodes.size() < 2 )
1424     return false;
1425
1426   int nbRows = faceSubMesh->NbElements() / ( sortedBaseNodes.size()-1 ) + 1;
1427   map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNodes.begin();
1428   double f = u_n->first, range = sortedBaseNodes.rbegin()->first - f;
1429   for ( ; u_n != sortedBaseNodes.end(); u_n++ )
1430   {
1431     double par = ( u_n->first - f ) / range;
1432     vector<const SMDS_MeshNode*>& nCol = theParam2ColumnMap[ par ];
1433     nCol.resize( nbRows );
1434     nCol[0] = u_n->second;
1435   }
1436
1437   // fill theParam2ColumnMap column by column by passing from nodes on
1438   // theBaseEdge up via mesh faces on theFace
1439
1440   TParam2ColumnMap::iterator par_nVec_2 = theParam2ColumnMap.begin();
1441   TParam2ColumnMap::iterator par_nVec_1 = par_nVec_2++;
1442   TIDSortedElemSet emptySet, avoidSet;
1443   for ( ; par_nVec_2 != theParam2ColumnMap.end(); ++par_nVec_1, ++par_nVec_2 )
1444   {
1445     vector<const SMDS_MeshNode*>& nCol1 = par_nVec_1->second;
1446     vector<const SMDS_MeshNode*>& nCol2 = par_nVec_2->second;
1447
1448     int i1, i2, iRow = 0;
1449     const SMDS_MeshNode *n1 = nCol1[0], *n2 = nCol2[0];
1450     // find face sharing node n1 and n2 and belonging to faceSubMesh
1451     while ( const SMDS_MeshElement* face =
1452             SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2))
1453     {
1454       if ( faceSubMesh->Contains( face ))
1455       {
1456         int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes();
1457         if ( nbNodes != 4 )
1458           return false;
1459         n1 = face->GetNode( (i2+2) % 4 ); // opposite corner of quadrangle face
1460         n2 = face->GetNode( (i1+2) % 4 );
1461         if ( ++iRow >= nbRows )
1462           return false;
1463         nCol1[ iRow ] = n1;
1464         nCol2[ iRow ] = n2;
1465         avoidSet.clear();
1466       }
1467       avoidSet.insert( face );
1468     }
1469     if ( iRow + 1 < nbRows ) // compact if necessary
1470       nCol1.resize( iRow + 1 ), nCol2.resize( iRow + 1 );
1471   }
1472   return true;
1473 }
1474
1475 //=======================================================================
1476 //function : NbAncestors
1477 //purpose  : Return number of unique ancestors of the shape
1478 //=======================================================================
1479
1480 int SMESH_MesherHelper::NbAncestors(const TopoDS_Shape& shape,
1481                                     const SMESH_Mesh&   mesh,
1482                                     TopAbs_ShapeEnum    ancestorType/*=TopAbs_SHAPE*/)
1483 {
1484   TopTools_MapOfShape ancestors;
1485   TopTools_ListIteratorOfListOfShape ansIt( mesh.GetAncestors(shape) );
1486   for ( ; ansIt.More(); ansIt.Next() ) {
1487     if ( ancestorType == TopAbs_SHAPE || ansIt.Value().ShapeType() == ancestorType )
1488       ancestors.Add( ansIt.Value() );
1489   }
1490   return ancestors.Extent();
1491 }
1492
1493 //=======================================================================
1494 //function : GetSubShapeOri
1495 //purpose  : Return orientation of sub-shape in the main shape
1496 //=======================================================================
1497
1498 TopAbs_Orientation SMESH_MesherHelper::GetSubShapeOri(const TopoDS_Shape& shape,
1499                                                       const TopoDS_Shape& subShape)
1500 {
1501   TopAbs_Orientation ori = TopAbs_Orientation(-1);
1502   if ( !shape.IsNull() && !subShape.IsNull() )
1503   {
1504     TopExp_Explorer e( shape, subShape.ShapeType() );
1505     if ( shape.Orientation() >= TopAbs_INTERNAL ) // TopAbs_INTERNAL or TopAbs_EXTERNAL
1506       e.Init( shape.Oriented(TopAbs_FORWARD), subShape.ShapeType() );
1507     for ( ; e.More(); e.Next())
1508       if ( subShape.IsSame( e.Current() ))
1509         break;
1510     if ( e.More() )
1511       ori = e.Current().Orientation();
1512   }
1513   return ori;
1514 }
1515
1516 //=======================================================================
1517 //function : IsSubShape
1518 //purpose  :
1519 //=======================================================================
1520
1521 bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape,
1522                                      const TopoDS_Shape& mainShape )
1523 {
1524   if ( !shape.IsNull() && !mainShape.IsNull() )
1525   {
1526     for ( TopExp_Explorer exp( mainShape, shape.ShapeType());
1527           exp.More();
1528           exp.Next() )
1529       if ( shape.IsSame( exp.Current() ))
1530         return true;
1531   }
1532   SCRUTE((shape.IsNull()));
1533   SCRUTE((mainShape.IsNull()));
1534   return false;
1535 }
1536
1537 //=======================================================================
1538 //function : IsSubShape
1539 //purpose  :
1540 //=======================================================================
1541
1542 bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, SMESH_Mesh* aMesh )
1543 {
1544   if ( shape.IsNull() || !aMesh )
1545     return false;
1546   return
1547     aMesh->GetMeshDS()->ShapeToIndex( shape ) ||
1548     // PAL16202
1549     shape.ShapeType() == TopAbs_COMPOUND && aMesh->GetMeshDS()->IsGroupOfSubShapes( shape );
1550 }
1551
1552 //================================================================================
1553 /*!
1554  * \brief Return maximal tolerance of shape
1555  */
1556 //================================================================================
1557
1558 double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
1559 {
1560   double tol = Precision::Confusion();
1561   TopExp_Explorer exp;
1562   for ( exp.Init( shape, TopAbs_FACE ); exp.More(); exp.Next() )
1563     tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Face( exp.Current())));
1564   for ( exp.Init( shape, TopAbs_EDGE ); exp.More(); exp.Next() )
1565     tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Edge( exp.Current())));
1566   for ( exp.Init( shape, TopAbs_VERTEX ); exp.More(); exp.Next() )
1567     tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Vertex( exp.Current())));
1568
1569   return tol;
1570 }
1571
1572 //=======================================================================
1573 //function : IsQuadraticMesh
1574 //purpose  : Check mesh without geometry for: if all elements on this shape are quadratic,
1575 //           quadratic elements will be created.
1576 //           Used then generated 3D mesh without geometry.
1577 //=======================================================================
1578
1579 SMESH_MesherHelper:: MType SMESH_MesherHelper::IsQuadraticMesh()
1580 {
1581   int NbAllEdgsAndFaces=0;
1582   int NbQuadFacesAndEdgs=0;
1583   int NbFacesAndEdges=0;
1584   //All faces and edges
1585   NbAllEdgsAndFaces = myMesh->NbEdges() + myMesh->NbFaces();
1586
1587   //Quadratic faces and edges
1588   NbQuadFacesAndEdgs = myMesh->NbEdges(ORDER_QUADRATIC) + myMesh->NbFaces(ORDER_QUADRATIC);
1589
1590   //Linear faces and edges
1591   NbFacesAndEdges = myMesh->NbEdges(ORDER_LINEAR) + myMesh->NbFaces(ORDER_LINEAR);
1592
1593   if (NbAllEdgsAndFaces == NbQuadFacesAndEdgs) {
1594     //Quadratic mesh
1595     return SMESH_MesherHelper::QUADRATIC;
1596   }
1597   else if (NbAllEdgsAndFaces == NbFacesAndEdges) {
1598     //Linear mesh
1599     return SMESH_MesherHelper::LINEAR;
1600   }
1601   else
1602     //Mesh with both type of elements
1603     return SMESH_MesherHelper::COMP;
1604 }
1605
1606 //=======================================================================
1607 //function : GetOtherParam
1608 //purpose  : Return an alternative parameter for a node on seam
1609 //=======================================================================
1610
1611 double SMESH_MesherHelper::GetOtherParam(const double param) const
1612 {
1613   int i = myParIndex & U_periodic ? 0 : 1;
1614   return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i];
1615 }
1616
1617 //#include <Perf_Meter.hxx>
1618
1619 //=======================================================================
1620 namespace { // Structures used by FixQuadraticElements()
1621 //=======================================================================
1622
1623 #define __DMP__(txt) \
1624   //cout << txt
1625 #define MSG(txt) __DMP__(txt<<endl)
1626 #define MSGBEG(txt) __DMP__(txt)
1627
1628   //const double straightTol2 = 1e-33; // to detect straing links
1629   bool isStraightLink(double linkLen2, double middleNodeMove2)
1630   {
1631     // straight if <node move> < 1/15 * <link length>
1632     return middleNodeMove2 < 1/15./15. * linkLen2;
1633   }
1634
1635   struct QFace;
1636   // ---------------------------------------
1637   /*!
1638    * \brief Quadratic link knowing its faces
1639    */
1640   struct QLink: public SMESH_TLink
1641   {
1642     const SMDS_MeshNode*          _mediumNode;
1643     mutable vector<const QFace* > _faces;
1644     mutable gp_Vec                _nodeMove;
1645     mutable int                   _nbMoves;
1646
1647     QLink(const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshNode* nm):
1648       SMESH_TLink( n1,n2 ), _mediumNode(nm), _nodeMove(0,0,0), _nbMoves(0) {
1649       _faces.reserve(4);
1650       //if ( MediumPos() != SMDS_TOP_3DSPACE )
1651         _nodeMove = MediumPnt() - MiddlePnt();
1652     }
1653     void SetContinuesFaces() const;
1654     const QFace* GetContinuesFace( const QFace* face ) const;
1655     bool OnBoundary() const;
1656     gp_XYZ MiddlePnt() const { return ( XYZ( node1() ) + XYZ( node2() )) / 2.; }
1657     gp_XYZ MediumPnt() const { return XYZ( _mediumNode ); }
1658
1659     SMDS_TypeOfPosition MediumPos() const
1660     { return _mediumNode->GetPosition()->GetTypeOfPosition(); }
1661     SMDS_TypeOfPosition EndPos(bool isSecond) const
1662     { return (isSecond ? node2() : node1())->GetPosition()->GetTypeOfPosition(); }
1663     const SMDS_MeshNode* EndPosNode(SMDS_TypeOfPosition pos) const
1664     { return EndPos(0) == pos ? node1() : EndPos(1) == pos ? node2() : 0; }
1665
1666     void Move(const gp_Vec& move, bool sum=false) const
1667     { _nodeMove += move; _nbMoves += sum ? (_nbMoves==0) : 1; }
1668     gp_XYZ Move() const { return _nodeMove.XYZ() / _nbMoves; }
1669     bool IsMoved() const { return (_nbMoves > 0 && !IsStraight()); }
1670     bool IsStraight() const
1671     { return isStraightLink( (XYZ(node1())-XYZ(node2())).SquareModulus(),
1672                              _nodeMove.SquareMagnitude());
1673     }
1674     bool operator<(const QLink& other) const {
1675       return (node1()->GetID() == other.node1()->GetID() ?
1676               node2()->GetID() < other.node2()->GetID() :
1677               node1()->GetID() < other.node1()->GetID());
1678     }
1679     struct PtrComparator {
1680       bool operator() (const QLink* l1, const QLink* l2 ) const { return *l1 < *l2; }
1681     };
1682   };
1683   // ---------------------------------------------------------
1684   /*!
1685    * \brief Link in the chain of links; it connects two faces
1686    */
1687   struct TChainLink
1688   {
1689     const QLink*         _qlink;
1690     mutable const QFace* _qfaces[2];
1691
1692     TChainLink(const QLink* qlink=0):_qlink(qlink) {
1693       _qfaces[0] = _qfaces[1] = 0;
1694     }
1695     void SetFace(const QFace* face) const { int iF = _qfaces[0] ? 1 : 0; _qfaces[iF]=face; }
1696
1697     bool IsBoundary() const { return !_qfaces[1]; }
1698
1699     void RemoveFace( const QFace* face ) const
1700     { _qfaces[(face == _qfaces[1])] = 0; if (!_qfaces[0]) std::swap(_qfaces[0],_qfaces[1]); }
1701
1702     const QFace* NextFace( const QFace* f ) const
1703     { return _qfaces[0]==f ? _qfaces[1] : _qfaces[0]; }
1704
1705     const SMDS_MeshNode* NextNode( const SMDS_MeshNode* n ) const
1706     { return n == _qlink->node1() ? _qlink->node2() : _qlink->node1(); }
1707
1708     bool operator<(const TChainLink& other) const { return *_qlink < *other._qlink; }
1709
1710     operator bool() const { return (_qlink); }
1711
1712     const QLink* operator->() const { return _qlink; }
1713
1714     gp_Vec Normal() const;
1715   };
1716   // --------------------------------------------------------------------
1717   typedef list< TChainLink > TChain;
1718   typedef set < TChainLink > TLinkSet;
1719   typedef TLinkSet::const_iterator TLinkInSet;
1720
1721   const int theFirstStep = 5;
1722
1723   enum { ERR_OK, ERR_TRI, ERR_PRISM, ERR_UNKNOWN }; // errors of QFace::GetLinkChain()
1724   // --------------------------------------------------------------------
1725   /*!
1726    * \brief Face shared by two volumes and bound by QLinks
1727    */
1728   struct QFace: public TIDSortedElemSet
1729   {
1730     mutable const SMDS_MeshElement* _volumes[2];
1731     mutable vector< const QLink* >  _sides;
1732     mutable bool                    _sideIsAdded[4]; // added in chain of links
1733     gp_Vec                          _normal;
1734 #ifdef _DEBUG_
1735     mutable const SMDS_MeshElement* _face;
1736 #endif
1737
1738     QFace( const vector< const QLink*>& links, const SMDS_MeshElement* face=0 );
1739
1740     void SetVolume(const SMDS_MeshElement* v) const { _volumes[ _volumes[0] ? 1 : 0 ] = v; }
1741
1742     int NbVolumes() const { return !_volumes[0] ? 0 : !_volumes[1] ? 1 : 2; }
1743
1744     void AddSelfToLinks() const {
1745       for ( int i = 0; i < _sides.size(); ++i )
1746         _sides[i]->_faces.push_back( this );
1747     }
1748     int LinkIndex( const QLink* side ) const {
1749       for (int i=0; i<_sides.size(); ++i ) if ( _sides[i] == side ) return i;
1750       return -1;
1751     }
1752     bool GetLinkChain( int iSide, TChain& chain, SMDS_TypeOfPosition pos, int& err) const;
1753
1754     bool GetLinkChain( TChainLink& link, TChain& chain, SMDS_TypeOfPosition pos, int& err) const
1755     {
1756       int i = LinkIndex( link._qlink );
1757       if ( i < 0 ) return true;
1758       _sideIsAdded[i] = true;
1759       link.SetFace( this );
1760       // continue from opposite link
1761       return GetLinkChain( (i+2)%_sides.size(), chain, pos, err );
1762     }
1763     bool IsBoundary() const { return !_volumes[1]; }
1764
1765     bool Contains( const SMDS_MeshNode* node ) const { return count(node); }
1766
1767     TLinkInSet GetBoundaryLink( const TLinkSet&      links,
1768                                 const TChainLink&    avoidLink,
1769                                 TLinkInSet *         notBoundaryLink = 0,
1770                                 const SMDS_MeshNode* nodeToContain = 0,
1771                                 bool *               isAdjacentUsed = 0,
1772                                 int                  nbRecursionsLeft = -1) const;
1773
1774     TLinkInSet GetLinkByNode( const TLinkSet&      links,
1775                               const TChainLink&    avoidLink,
1776                               const SMDS_MeshNode* nodeToContain) const;
1777
1778     const SMDS_MeshNode* GetNodeInFace() const {
1779       for ( int iL = 0; iL < _sides.size(); ++iL )
1780         if ( _sides[iL]->MediumPos() == SMDS_TOP_FACE ) return _sides[iL]->_mediumNode;
1781       return 0;
1782     }
1783
1784     gp_Vec LinkNorm(const int i, SMESH_MesherHelper* theFaceHelper=0) const;
1785
1786     double MoveByBoundary( const TChainLink&   theLink,
1787                            const gp_Vec&       theRefVec,
1788                            const TLinkSet&     theLinks,
1789                            SMESH_MesherHelper* theFaceHelper=0,
1790                            const double        thePrevLen=0,
1791                            const int           theStep=theFirstStep,
1792                            gp_Vec*             theLinkNorm=0,
1793                            double              theSign=1.0) const;
1794   };
1795
1796   //================================================================================
1797   /*!
1798    * \brief Dump QLink and QFace
1799    */
1800   ostream& operator << (ostream& out, const QLink& l)
1801   {
1802     out <<"QLink nodes: "
1803         << l.node1()->GetID() << " - "
1804         << l._mediumNode->GetID() << " - "
1805         << l.node2()->GetID() << endl;
1806     return out;
1807   }
1808   ostream& operator << (ostream& out, const QFace& f)
1809   {
1810     out <<"QFace nodes: "/*<< &f << "  "*/;
1811     for ( TIDSortedElemSet::const_iterator n = f.begin(); n != f.end(); ++n )
1812       out << (*n)->GetID() << " ";
1813     out << " \tvolumes: "
1814         << (f._volumes[0] ? f._volumes[0]->GetID() : 0) << " "
1815         << (f._volumes[1] ? f._volumes[1]->GetID() : 0);
1816     out << "  \tNormal: "<< f._normal.X() <<", "<<f._normal.Y() <<", "<<f._normal.Z() << endl;
1817     return out;
1818   }
1819
1820   //================================================================================
1821   /*!
1822    * \brief Construct QFace from QLinks
1823    */
1824   //================================================================================
1825
1826   QFace::QFace( const vector< const QLink*>& links, const SMDS_MeshElement* face )
1827   {
1828     _volumes[0] = _volumes[1] = 0;
1829     _sides = links;
1830     _sideIsAdded[0]=_sideIsAdded[1]=_sideIsAdded[2]=_sideIsAdded[3]=false;
1831     _normal.SetCoord(0,0,0);
1832     for ( int i = 1; i < _sides.size(); ++i ) {
1833       const QLink *l1 = _sides[i-1], *l2 = _sides[i];
1834       insert( l1->node1() ); insert( l1->node2() );
1835       // compute normal
1836       gp_Vec v1( XYZ( l1->node2()), XYZ( l1->node1()));
1837       gp_Vec v2( XYZ( l2->node1()), XYZ( l2->node2()));
1838       if ( l1->node1() != l2->node1() && l1->node2() != l2->node2() )
1839         v1.Reverse();
1840       _normal += v1 ^ v2;
1841     }
1842     double normSqSize = _normal.SquareMagnitude();
1843     if ( normSqSize > numeric_limits<double>::min() )
1844       _normal /= sqrt( normSqSize );
1845     else
1846       _normal.SetCoord(1e-33,0,0);
1847
1848 #ifdef _DEBUG_
1849     _face = face;
1850 #endif
1851   }
1852   //================================================================================
1853   /*!
1854    * \brief Make up a chain of links
1855    *  \param iSide - link to add first
1856    *  \param chain - chain to fill in
1857    *  \param pos   - postion of medium nodes the links should have
1858    *  \param error - out, specifies what is wrong
1859    *  \retval bool - false if valid chain can't be built; "valid" means that links
1860    *                 of the chain belongs to rectangles bounding hexahedrons
1861    */
1862   //================================================================================
1863
1864   bool QFace::GetLinkChain( int iSide, TChain& chain, SMDS_TypeOfPosition pos, int& error) const
1865   {
1866     if ( iSide >= _sides.size() ) // wrong argument iSide
1867       return false;
1868     if ( _sideIsAdded[ iSide ]) // already in chain
1869       return true;
1870
1871     if ( _sides.size() != 4 ) { // triangle - visit all my continous faces
1872       MSGBEG( *this );
1873       TLinkSet links;
1874       list< const QFace* > faces( 1, this );
1875       while ( !faces.empty() ) {
1876         const QFace* face = faces.front();
1877         for ( int i = 0; i < face->_sides.size(); ++i ) {
1878           if ( !face->_sideIsAdded[i] && face->_sides[i] ) {
1879             face->_sideIsAdded[i] = true;
1880             // find a face side in the chain
1881             TLinkInSet chLink = links.insert( TChainLink(face->_sides[i])).first;
1882 //             TChain::iterator chLink = chain.begin();
1883 //             for ( ; chLink != chain.end(); ++chLink )
1884 //               if ( chLink->_qlink == face->_sides[i] )
1885 //                 break;
1886 //             if ( chLink == chain.end() )
1887 //               chLink = chain.insert( chain.begin(), TChainLink(face->_sides[i]));
1888             // add a face to a chained link and put a continues face in the queue
1889             chLink->SetFace( face );
1890             if ( face->_sides[i]->MediumPos() >= pos )
1891               if ( const QFace* contFace = face->_sides[i]->GetContinuesFace( face ))
1892                 faces.push_back( contFace );
1893           }
1894         }
1895         faces.pop_front();
1896       }
1897       if ( error < ERR_TRI )
1898         error = ERR_TRI;
1899       chain.insert( chain.end(), links.begin(),links.end() );
1900       return false;
1901     }
1902     _sideIsAdded[iSide] = true; // not to add this link to chain again
1903     const QLink* link = _sides[iSide];
1904     if ( !link)
1905       return true;
1906
1907     // add link into chain
1908     TChain::iterator chLink = chain.insert( chain.begin(), TChainLink(link));
1909     chLink->SetFace( this );
1910     MSGBEG( *this );
1911
1912     // propagate from quadrangle to neighbour faces
1913     if ( link->MediumPos() >= pos ) {
1914       int nbLinkFaces = link->_faces.size();
1915       if ( nbLinkFaces == 4 || nbLinkFaces < 4 && link->OnBoundary()) {
1916         // hexahedral mesh or boundary quadrangles - goto a continous face
1917         if ( const QFace* f = link->GetContinuesFace( this ))
1918           return f->GetLinkChain( *chLink, chain, pos, error );
1919       }
1920       else {
1921         TChainLink chLink(link); // side face of prismatic mesh - visit all faces of iSide
1922         for ( int i = 0; i < nbLinkFaces; ++i )
1923           if ( link->_faces[i] )
1924             link->_faces[i]->GetLinkChain( chLink, chain, pos, error );
1925         if ( error < ERR_PRISM )
1926           error = ERR_PRISM;
1927         return false;
1928       }
1929     }
1930     return true;
1931   }
1932
1933   //================================================================================
1934   /*!
1935    * \brief Return a boundary link of the triangle face
1936    *  \param links - set of all links
1937    *  \param avoidLink - link not to return
1938    *  \param notBoundaryLink - out, neither the returned link nor avoidLink
1939    *  \param nodeToContain - node the returned link must contain; if provided, search
1940    *                         also performed on adjacent faces
1941    *  \param isAdjacentUsed - returns true if link is found in adjacent faces
1942    *  \param nbRecursionsLeft - to limit recursion
1943    */
1944   //================================================================================
1945
1946   TLinkInSet QFace::GetBoundaryLink( const TLinkSet&      links,
1947                                      const TChainLink&    avoidLink,
1948                                      TLinkInSet *         notBoundaryLink,
1949                                      const SMDS_MeshNode* nodeToContain,
1950                                      bool *               isAdjacentUsed,
1951                                      int                  nbRecursionsLeft) const
1952   {
1953     TLinkInSet linksEnd = links.end(), boundaryLink = linksEnd;
1954
1955     typedef list< pair< const QFace*, TLinkInSet > > TFaceLinkList;
1956     TFaceLinkList adjacentFaces;
1957
1958     for ( int iL = 0; iL < _sides.size(); ++iL )
1959     {
1960       if ( avoidLink._qlink == _sides[iL] )
1961         continue;
1962       TLinkInSet link = links.find( _sides[iL] );
1963       if ( link == linksEnd ) continue;
1964       if ( (*link)->MediumPos() > SMDS_TOP_FACE )
1965         continue; // We work on faces here, don't go inside a solid
1966
1967       // check link
1968       if ( link->IsBoundary() ) {
1969         if ( !nodeToContain ||
1970              (*link)->node1() == nodeToContain ||
1971              (*link)->node2() == nodeToContain )
1972         {
1973           boundaryLink = link;
1974           if ( !notBoundaryLink ) break;
1975         }
1976       }
1977       else if ( notBoundaryLink ) {
1978         *notBoundaryLink = link;
1979         if ( boundaryLink != linksEnd ) break;
1980       }
1981
1982       if ( boundaryLink == linksEnd && nodeToContain ) // collect adjacent faces
1983         if ( const QFace* adj = link->NextFace( this ))
1984           if ( adj->Contains( nodeToContain ))
1985             adjacentFaces.push_back( make_pair( adj, link ));
1986     }
1987
1988     if ( isAdjacentUsed ) *isAdjacentUsed = false;
1989     if ( boundaryLink == linksEnd && nodeToContain && nbRecursionsLeft) // check adjacent faces
1990     {
1991       if ( nbRecursionsLeft < 0 )
1992         nbRecursionsLeft = nodeToContain->NbInverseElements();
1993       TFaceLinkList::iterator adj = adjacentFaces.begin();
1994       for ( ; boundaryLink == linksEnd && adj != adjacentFaces.end(); ++adj )
1995         boundaryLink = adj->first->GetBoundaryLink( links, *(adj->second), 0, nodeToContain,
1996                                                     isAdjacentUsed, nbRecursionsLeft-1);
1997       if ( isAdjacentUsed ) *isAdjacentUsed = true;
1998     }
1999     return boundaryLink;
2000   }
2001   //================================================================================
2002   /*!
2003    * \brief Return a link ending at the given node but not avoidLink
2004    */
2005   //================================================================================
2006
2007   TLinkInSet QFace::GetLinkByNode( const TLinkSet&      links,
2008                                    const TChainLink&    avoidLink,
2009                                    const SMDS_MeshNode* nodeToContain) const
2010   {
2011     for ( int i = 0; i < _sides.size(); ++i )
2012       if ( avoidLink._qlink != _sides[i] &&
2013            (_sides[i]->node1() == nodeToContain || _sides[i]->node2() == nodeToContain ))
2014         return links.find( _sides[ i ]);
2015     return links.end();
2016   }
2017
2018   //================================================================================
2019   /*!
2020    * \brief Return normal to the i-th side pointing outside the face
2021    */
2022   //================================================================================
2023
2024   gp_Vec QFace::LinkNorm(const int i, SMESH_MesherHelper* /*uvHelper*/) const
2025   {
2026     gp_Vec norm, vecOut;
2027 //     if ( uvHelper ) {
2028 //       TopoDS_Face face = TopoDS::Face( uvHelper->GetSubShape());
2029 //       const SMDS_MeshNode* inFaceNode = uvHelper->GetNodeUVneedInFaceNode() ? GetNodeInFace() : 0;
2030 //       gp_XY uv1 = uvHelper->GetNodeUV( face, _sides[i]->node1(), inFaceNode );
2031 //       gp_XY uv2 = uvHelper->GetNodeUV( face, _sides[i]->node2(), inFaceNode );
2032 //       norm.SetCoord( uv1.Y() - uv2.Y(), uv2.X() - uv1.X(), 0 );
2033
2034 //       const QLink* otherLink = _sides[(i + 1) % _sides.size()];
2035 //       const SMDS_MeshNode* otherNode =
2036 //         otherLink->node1() == _sides[i]->node1() ? otherLink->node2() : otherLink->node1();
2037 //       gp_XY pIn = uvHelper->GetNodeUV( face, otherNode, inFaceNode );
2038 //       vecOut.SetCoord( uv1.X() - pIn.X(), uv1.Y() - pIn.Y(), 0 );
2039 //     }
2040 //     else {
2041       norm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2()));
2042       gp_XYZ pIn = ( XYZ( _sides[0]->node1() ) +
2043                      XYZ( _sides[0]->node2() ) +
2044                      XYZ( _sides[1]->node1() )) / 3.;
2045       vecOut.SetXYZ( _sides[i]->MiddlePnt() - pIn );
2046       //}
2047     if ( norm * vecOut < 0 )
2048       norm.Reverse();
2049     double mag2 = norm.SquareMagnitude();
2050     if ( mag2 > numeric_limits<double>::min() )
2051       norm /= sqrt( mag2 );
2052     return norm;
2053   }
2054   //================================================================================
2055   /*!
2056    * \brief Move medium node of theLink according to its distance from boundary
2057    *  \param theLink - link to fix
2058    *  \param theRefVec - movement of boundary
2059    *  \param theLinks - all adjacent links of continous triangles
2060    *  \param theFaceHelper - helper is not used so far
2061    *  \param thePrevLen - distance from the boundary
2062    *  \param theStep - number of steps till movement propagation limit
2063    *  \param theLinkNorm - out normal to theLink
2064    *  \param theSign - 1 or -1 depending on movement of boundary
2065    *  \retval double - distance from boundary to propagation limit or other boundary
2066    */
2067   //================================================================================
2068
2069   double QFace::MoveByBoundary( const TChainLink&   theLink,
2070                                 const gp_Vec&       theRefVec,
2071                                 const TLinkSet&     theLinks,
2072                                 SMESH_MesherHelper* theFaceHelper,
2073                                 const double        thePrevLen,
2074                                 const int           theStep,
2075                                 gp_Vec*             theLinkNorm,
2076                                 double              theSign) const
2077   {
2078     if ( !theStep )
2079       return thePrevLen; // propagation limit reached
2080
2081     int iL; // index of theLink
2082     for ( iL = 0; iL < _sides.size(); ++iL )
2083       if ( theLink._qlink == _sides[ iL ])
2084         break;
2085
2086     MSG(string(theStep,'.')<<" Ref( "<<theRefVec.X()<<","<<theRefVec.Y()<<","<<theRefVec.Z()<<" )"
2087         <<" thePrevLen " << thePrevLen);
2088     MSG(string(theStep,'.')<<" "<<*theLink._qlink);
2089
2090     gp_Vec linkNorm = -LinkNorm( iL/*, theFaceHelper*/ ); // normal to theLink
2091     double refProj = theRefVec * linkNorm; // project movement vector to normal of theLink
2092     if ( theStep == theFirstStep )
2093       theSign = refProj < 0. ? -1. : 1.;
2094     else if ( theSign * refProj < 0.4 * theRefVec.Magnitude())
2095       return thePrevLen; // to propagate movement forward only, not in side dir or backward
2096
2097     int iL1 = (iL + 1) % 3, iL2 = (iL + 2) % 3; // indices of the two other links of triangle
2098     TLinkInSet link1 = theLinks.find( _sides[iL1] );
2099     TLinkInSet link2 = theLinks.find( _sides[iL2] );
2100     if ( link1 == theLinks.end() || link2 == theLinks.end() )
2101       return thePrevLen;
2102     const QFace* f1 = link1->NextFace( this ); // adjacent faces
2103     const QFace* f2 = link2->NextFace( this );
2104
2105     // propagate to adjacent faces till limit step or boundary
2106     double len1 = thePrevLen + (theLink->MiddlePnt() - _sides[iL1]->MiddlePnt()).Modulus();
2107     double len2 = thePrevLen + (theLink->MiddlePnt() - _sides[iL2]->MiddlePnt()).Modulus();
2108     gp_Vec linkDir1(0,0,0); // initialize to avoid valgrind error ("Conditional jump...")
2109     gp_Vec linkDir2(0,0,0);
2110     try {
2111       OCC_CATCH_SIGNALS;
2112       if ( f1 )
2113         len1 = f1->MoveByBoundary
2114           ( *link1, theRefVec, theLinks, theFaceHelper, len1, theStep-1, &linkDir1, theSign);
2115       else
2116         linkDir1 = LinkNorm( iL1/*, theFaceHelper*/ );
2117     } catch (...) {
2118       MSG( " --------------- EXCEPTION");
2119       return thePrevLen;
2120     }
2121     try {
2122       OCC_CATCH_SIGNALS;
2123       if ( f2 )
2124         len2 = f2->MoveByBoundary
2125           ( *link2, theRefVec, theLinks, theFaceHelper, len2, theStep-1, &linkDir2, theSign);
2126       else
2127         linkDir2 = LinkNorm( iL2/*, theFaceHelper*/ );
2128     } catch (...) {
2129       MSG( " --------------- EXCEPTION");
2130       return thePrevLen;
2131     }
2132
2133     double fullLen = 0;
2134     if ( theStep != theFirstStep )
2135     {
2136       // choose chain length by direction of propagation most codirected with theRefVec
2137       bool choose1 = ( theRefVec * linkDir1 * theSign > theRefVec * linkDir2 * theSign );
2138       fullLen = choose1 ? len1 : len2;
2139       double r = thePrevLen / fullLen;
2140
2141       gp_Vec move = linkNorm * refProj * ( 1 - r );
2142       theLink->Move( move, true );
2143
2144       MSG(string(theStep,'.')<<" Move "<< theLink->_mediumNode->GetID()<<
2145           " by " << refProj * ( 1 - r ) << " following " <<
2146           (choose1 ? *link1->_qlink : *link2->_qlink));
2147
2148       if ( theLinkNorm ) *theLinkNorm = linkNorm;
2149     }
2150     return fullLen;
2151   }
2152
2153   //================================================================================
2154   /*!
2155    * \brief Find pairs of continues faces
2156    */
2157   //================================================================================
2158
2159   void QLink::SetContinuesFaces() const
2160   {
2161     //       x0         x - QLink, [-|] - QFace, v - volume
2162     //   v0  |   v1
2163     //       |          Between _faces of link x2 two vertical faces are continues
2164     // x1----x2-----x3  and two horizontal faces are continues. We set vertical faces
2165     //       |          to _faces[0] and _faces[1] and horizontal faces to
2166     //   v2  |   v3     _faces[2] and _faces[3] (or vise versa).
2167     //       x4
2168
2169     if ( _faces.empty() )
2170       return;
2171     int iFaceCont = -1;
2172     for ( int iF = 1; iFaceCont < 0 && iF < _faces.size(); ++iF )
2173     {
2174       // look for a face bounding none of volumes bound by _faces[0]
2175       bool sameVol = false;
2176       int nbVol = _faces[iF]->NbVolumes();
2177       for ( int iV = 0; !sameVol && iV < nbVol; ++iV )
2178         sameVol = ( _faces[iF]->_volumes[iV] == _faces[0]->_volumes[0] ||
2179                     _faces[iF]->_volumes[iV] == _faces[0]->_volumes[1]);
2180       if ( !sameVol )
2181         iFaceCont = iF;
2182     }
2183     if ( iFaceCont > 0 ) // continues faces found, set one by the other
2184     {
2185       if ( iFaceCont != 1 )
2186         std::swap( _faces[1], _faces[iFaceCont] );
2187     }
2188     else if ( _faces.size() > 1 ) // not found, set NULL by the first face
2189     {
2190       _faces.insert( ++_faces.begin(), 0 );
2191     }
2192   }
2193   //================================================================================
2194   /*!
2195    * \brief Return a face continues to the given one
2196    */
2197   //================================================================================
2198
2199   const QFace* QLink::GetContinuesFace( const QFace* face ) const
2200   {
2201     for ( int i = 0; i < _faces.size(); ++i ) {
2202       if ( _faces[i] == face ) {
2203         int iF = i < 2 ? 1-i : 5-i;
2204         return iF < _faces.size() ? _faces[iF] : 0;
2205       }
2206     }
2207     return 0;
2208   }
2209   //================================================================================
2210   /*!
2211    * \brief True if link is on mesh boundary
2212    */
2213   //================================================================================
2214
2215   bool QLink::OnBoundary() const
2216   {
2217     for ( int i = 0; i < _faces.size(); ++i )
2218       if (_faces[i] && _faces[i]->IsBoundary()) return true;
2219     return false;
2220   }
2221   //================================================================================
2222   /*!
2223    * \brief Return normal of link of the chain
2224    */
2225   //================================================================================
2226
2227   gp_Vec TChainLink::Normal() const {
2228     gp_Vec norm;
2229     if (_qfaces[0]) norm  = _qfaces[0]->_normal;
2230     if (_qfaces[1]) norm += _qfaces[1]->_normal;
2231     return norm;
2232   }
2233   //================================================================================
2234   /*!
2235    * \brief Move medium nodes of vertical links of pentahedrons adjacent by side faces
2236    */
2237   //================================================================================
2238
2239   void fixPrism( TChain& allLinks )
2240   {
2241     // separate boundary links from internal ones
2242     typedef set<const QLink*/*, QLink::PtrComparator*/> QLinkSet;
2243     QLinkSet interLinks, bndLinks1, bndLink2;
2244
2245     bool isCurved = false;
2246     for ( TChain::iterator lnk = allLinks.begin(); lnk != allLinks.end(); ++lnk ) {
2247       if ( (*lnk)->OnBoundary() )
2248         bndLinks1.insert( lnk->_qlink );
2249       else
2250         interLinks.insert( lnk->_qlink );
2251       isCurved = isCurved || !(*lnk)->IsStraight();
2252     }
2253     if ( !isCurved )
2254       return; // no need to move
2255
2256     QLinkSet *curBndLinks = &bndLinks1, *newBndLinks = &bndLink2;
2257
2258     while ( !interLinks.empty() && !curBndLinks->empty() )
2259     {
2260       // propagate movement from boundary links to connected internal links
2261       QLinkSet::iterator bnd = curBndLinks->begin(), bndEnd = curBndLinks->end();
2262       for ( ; bnd != bndEnd; ++bnd )
2263       {
2264         const QLink* bndLink = *bnd;
2265         for ( int i = 0; i < bndLink->_faces.size(); ++i ) // loop on faces of bndLink
2266         {
2267           const QFace* face = bndLink->_faces[i]; // quadrange lateral face of a prism
2268           if ( !face ) continue;
2269           // find and move internal link opposite to bndLink within the face
2270           int interInd = ( face->LinkIndex( bndLink ) + 2 ) % face->_sides.size();
2271           const QLink* interLink = face->_sides[ interInd ];
2272           QLinkSet::iterator pInterLink = interLinks.find( interLink );
2273           if ( pInterLink == interLinks.end() ) continue; // not internal link
2274           interLink->Move( bndLink->_nodeMove );
2275           // treated internal links become new boundary ones
2276           interLinks. erase( pInterLink );
2277           newBndLinks->insert( interLink );
2278         }
2279       }
2280       curBndLinks->clear();
2281       std::swap( curBndLinks, newBndLinks );
2282     }
2283   }
2284
2285   //================================================================================
2286   /*!
2287    * \brief Fix links of continues triangles near curved boundary
2288    */
2289   //================================================================================
2290
2291   void fixTriaNearBoundary( TChain & allLinks, SMESH_MesherHelper& /*helper*/)
2292   {
2293     if ( allLinks.empty() ) return;
2294
2295     TLinkSet linkSet( allLinks.begin(), allLinks.end());
2296     TLinkInSet linkIt = linkSet.begin(), linksEnd = linkSet.end();
2297
2298     // move in 2d if we are on geom face
2299 //     TopoDS_Face face;
2300 //     TopLoc_Location loc;
2301 //     SMESH_MesherHelper faceHelper( *helper.GetMesh());
2302 //     while ( linkIt->IsBoundary()) ++linkIt;
2303 //     if ( linkIt == linksEnd ) return;
2304 //     if ( (*linkIt)->MediumPos() == SMDS_TOP_FACE ) {
2305 //       bool checkPos = true;
2306 //       TopoDS_Shape f = helper.GetSubShapeByNode( (*linkIt)->_mediumNode, helper.GetMeshDS() );
2307 //       if ( !f.IsNull() && f.ShapeType() == TopAbs_FACE ) {
2308 //         face = TopoDS::Face( f );
2309 //         helper.GetNodeUV( face, (*linkIt)->_mediumNode, 0, &checkPos);
2310 //         if (checkPos)
2311 //           face.Nullify();
2312 //         else
2313 //           faceHelper.SetSubShape( face );
2314 //       }
2315 //     }
2316     for ( linkIt = linkSet.begin(); linkIt != linksEnd; ++linkIt)
2317     {
2318       if ( linkIt->IsBoundary() && !(*linkIt)->IsStraight() && linkIt->_qfaces[0])
2319       {
2320 //         if ( !face.IsNull() ) {
2321 //           const SMDS_MeshNode* inFaceNode =
2322 //             faceHelper.GetNodeUVneedInFaceNode() ? linkIt->_qfaces[0]->GetNodeInFace() : 0;
2323 //           gp_XY uvm = helper.GetNodeUV( face, (*linkIt)->_mediumNode, inFaceNode );
2324 //           gp_XY uv1 = helper.GetNodeUV( face, (*linkIt)->node1(), inFaceNode);
2325 //           gp_XY uv2 = helper.GetNodeUV( face, (*linkIt)->node2(), inFaceNode);
2326 //           gp_XY uvMove = uvm - helper.GetMiddleUV( BRep_Tool::Surface(face,loc), uv1, uv2);
2327 //           gp_Vec move( uvMove.X(), uvMove.Y(), 0 );
2328 //           linkIt->_qfaces[0]->MoveByBoundary( *linkIt, move, linkSet, &faceHelper );
2329 //         }
2330 //         else {
2331           linkIt->_qfaces[0]->MoveByBoundary( *linkIt, (*linkIt)->_nodeMove, linkSet );
2332           //}
2333       }
2334     }
2335   }
2336
2337   //================================================================================
2338   /*!
2339    * \brief Detect rectangular structure of links and build chains from them
2340    */
2341   //================================================================================
2342
2343   enum TSplitTriaResult {
2344     _OK, _NO_CORNERS, _FEW_ROWS, _MANY_ROWS, _NO_SIDELINK, _BAD_MIDQUAD, _NOT_RECT,
2345     _NO_MIDQUAD, _NO_UPTRIA, _BAD_SET_SIZE, _BAD_CORNER, _BAD_START, _NO_BOTLINK, _TWISTED_CHAIN };
2346
2347   TSplitTriaResult splitTrianglesIntoChains( TChain &            allLinks,
2348                                              vector< TChain> &   resultChains,
2349                                              SMDS_TypeOfPosition pos )
2350   {
2351     // put links in the set and evalute number of result chains by number of boundary links
2352     TLinkSet linkSet;
2353     int nbBndLinks = 0;
2354     for ( TChain::iterator lnk = allLinks.begin(); lnk != allLinks.end(); ++lnk ) {
2355       linkSet.insert( *lnk );
2356       nbBndLinks += lnk->IsBoundary();
2357     }
2358     resultChains.clear();
2359     resultChains.reserve( nbBndLinks / 2 );
2360
2361     TLinkInSet linkIt, linksEnd = linkSet.end();
2362
2363     // find a boundary link with corner node; corner node has position pos-2
2364     // i.e. SMDS_TOP_VERTEX for links on faces and SMDS_TOP_EDGE for
2365     // links in volume
2366     SMDS_TypeOfPosition cornerPos = SMDS_TypeOfPosition(pos-2);
2367     const SMDS_MeshNode* corner = 0;
2368     for ( linkIt = linkSet.begin(); linkIt != linksEnd; ++linkIt )
2369       if ( linkIt->IsBoundary() && (corner = (*linkIt)->EndPosNode(cornerPos)))
2370         break;
2371     if ( !corner)
2372       return _NO_CORNERS;
2373
2374     TLinkInSet           startLink = linkIt;
2375     const SMDS_MeshNode* startCorner = corner;
2376     vector< TChain* >    rowChains;
2377     int iCol = 0;
2378
2379     while ( startLink != linksEnd) // loop on columns
2380     {
2381       // We suppose we have a rectangular structure like shown here. We have found a
2382       //               corner of the rectangle (startCorner) and a boundary link sharing
2383       //    |/  |/  |  the startCorner (startLink). We are going to loop on rows of the
2384       //  --o---o---o  structure making several chains at once. One chain (columnChain)
2385       //    |\  |  /|  starts at startLink and continues upward (we look at the structure
2386       //  \ | \ | / |  from such point that startLink is on the bottom of the structure).
2387       //   \|  \|/  |  While going upward we also fill horizontal chains (rowChains) we
2388       //  --o---o---o  encounter.
2389       //   /|\  |\  |
2390       //  / | \ | \ |  startCorner
2391       //    |  \|  \|,'
2392       //  --o---o---o
2393       //          `.startLink
2394
2395       if ( resultChains.size() == nbBndLinks / 2 )
2396         return _NOT_RECT;
2397       resultChains.push_back( TChain() );
2398       TChain& columnChain = resultChains.back();
2399
2400       TLinkInSet botLink = startLink; // current horizontal link to go up from
2401       corner = startCorner; // current corner the botLink ends at
2402       int iRow = 0;
2403       while ( botLink != linksEnd ) // loop on rows
2404       {
2405         // add botLink to the columnChain
2406         columnChain.push_back( *botLink );
2407
2408         const QFace* botTria = botLink->_qfaces[0]; // bottom triangle bound by botLink
2409         if ( !botTria )
2410         { // the column ends
2411           if ( botLink == startLink )
2412             return _TWISTED_CHAIN; // issue 0020951
2413           linkSet.erase( botLink );
2414           if ( iRow != rowChains.size() )
2415             return _FEW_ROWS; // different nb of rows in columns
2416           break;
2417         }
2418         // find the link dividing the quadrangle (midQuadLink) and vertical boundary
2419         // link ending at <corner> (sideLink); there are two cases:
2420         // 1) midQuadLink does not end at <corner>, then we easily find it by botTria,
2421         //   since midQuadLink is not at boundary while sideLink is.
2422         // 2) midQuadLink ends at <corner>
2423         bool isCase2;
2424         TLinkInSet midQuadLink = linksEnd;
2425         TLinkInSet sideLink = botTria->GetBoundaryLink( linkSet, *botLink, &midQuadLink,
2426                                                         corner, &isCase2 );
2427         if ( isCase2 ) { // find midQuadLink among links of botTria
2428           midQuadLink = botTria->GetLinkByNode( linkSet, *botLink, corner );
2429           if ( midQuadLink->IsBoundary() )
2430             return _BAD_MIDQUAD;
2431         }
2432         if ( sideLink == linksEnd || midQuadLink == linksEnd || sideLink == midQuadLink )
2433           return sideLink == linksEnd ? _NO_SIDELINK : _NO_MIDQUAD;
2434
2435         // fill chains
2436         columnChain.push_back( *midQuadLink );
2437         if ( iRow >= rowChains.size() ) {
2438           if ( iCol > 0 )
2439             return _MANY_ROWS; // different nb of rows in columns
2440           if ( resultChains.size() == nbBndLinks / 2 )
2441             return _NOT_RECT;
2442           resultChains.push_back( TChain() );
2443           rowChains.push_back( & resultChains.back() );
2444         }
2445         rowChains[iRow]->push_back( *sideLink );
2446         rowChains[iRow]->push_back( *midQuadLink );
2447
2448         const QFace* upTria = midQuadLink->NextFace( botTria ); // upper tria of the rectangle
2449         if ( !upTria)
2450           return _NO_UPTRIA;
2451         if ( iRow == 0 ) {
2452           // prepare startCorner and startLink for the next column
2453           startCorner = startLink->NextNode( startCorner );
2454           if (isCase2)
2455             startLink = botTria->GetBoundaryLink( linkSet, *botLink, 0, startCorner );
2456           else
2457             startLink = upTria->GetBoundaryLink( linkSet, *midQuadLink, 0, startCorner );
2458           // check if no more columns remains
2459           if ( startLink != linksEnd ) {
2460             const SMDS_MeshNode* botNode = startLink->NextNode( startCorner );
2461             if ( (isCase2 ? botTria : upTria)->Contains( botNode ))
2462               startLink = linksEnd; // startLink bounds upTria or botTria
2463             else if ( startLink == botLink || startLink == midQuadLink || startLink == sideLink )
2464               return _BAD_START;
2465           }
2466         }
2467         // find bottom link and corner for the next row
2468         corner = sideLink->NextNode( corner );
2469         // next bottom link ends at the new corner
2470         linkSet.erase( botLink );
2471         botLink = upTria->GetLinkByNode( linkSet, (isCase2 ? *sideLink : *midQuadLink), corner );
2472         if ( botLink == linksEnd || botLink == midQuadLink || botLink == sideLink)
2473           return _NO_BOTLINK;
2474         if ( midQuadLink == startLink || sideLink == startLink )
2475           return _TWISTED_CHAIN; // issue 0020951
2476         linkSet.erase( midQuadLink );
2477         linkSet.erase( sideLink );
2478
2479         // make faces neighboring the found ones be boundary
2480         if ( startLink != linksEnd ) {
2481           const QFace* tria = isCase2 ? botTria : upTria;
2482           for ( int iL = 0; iL < 3; ++iL ) {
2483             linkIt = linkSet.find( tria->_sides[iL] );
2484             if ( linkIt != linksEnd )
2485               linkIt->RemoveFace( tria );
2486           }
2487         }
2488         if ( botLink->_qfaces[0] == upTria || botLink->_qfaces[1] == upTria )
2489           botLink->RemoveFace( upTria ); // make next botTria first in vector
2490
2491         iRow++;
2492       } // loop on rows
2493
2494       iCol++;
2495     }
2496     // In the linkSet, there must remain the last links of rowChains; add them
2497     if ( linkSet.size() != rowChains.size() )
2498       return _BAD_SET_SIZE;
2499     for ( int iRow = 0; iRow < rowChains.size(); ++iRow ) {
2500       // find the link (startLink) ending at startCorner
2501       corner = 0;
2502       for ( startLink = linkSet.begin(); startLink != linksEnd; ++startLink ) {
2503         if ( (*startLink)->node1() == startCorner ) {
2504           corner = (*startLink)->node2(); break;
2505         }
2506         else if ( (*startLink)->node2() == startCorner) {
2507           corner = (*startLink)->node1(); break;
2508         }
2509       }
2510       if ( startLink == linksEnd )
2511         return _BAD_CORNER;
2512       rowChains[ iRow ]->push_back( *startLink );
2513       linkSet.erase( startLink );
2514       startCorner = corner;
2515     }
2516
2517     return _OK;
2518   }
2519 } //namespace
2520
2521 //=======================================================================
2522 /*!
2523  * \brief Move medium nodes of faces and volumes to fix distorted elements
2524  * \param volumeOnly - to fix nodes on faces or not, if the shape is solid
2525  *
2526  * Issue 0020307: EDF 992 SMESH : Linea/Quadratic with Medium Node on Geometry
2527  */
2528 //=======================================================================
2529
2530 void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
2531 {
2532   // 0. Apply algorithm to solids or geom faces
2533   // ----------------------------------------------
2534   if ( myShape.IsNull() ) {
2535     if ( !myMesh->HasShapeToMesh() ) return;
2536     SetSubShape( myMesh->GetShapeToMesh() );
2537
2538 #ifdef _DEBUG_
2539     int nbSolids = 0;
2540     TopTools_IndexedMapOfShape solids;
2541     TopExp::MapShapes(myShape,TopAbs_SOLID,solids);
2542     nbSolids = solids.Extent();
2543 #endif
2544     TopTools_MapOfShape faces; // faces not in solid or in not meshed solid
2545     for ( TopExp_Explorer f(myShape,TopAbs_FACE,TopAbs_SOLID); f.More(); f.Next() ) {
2546       faces.Add( f.Current() ); // not in solid
2547     }
2548     for ( TopExp_Explorer s(myShape,TopAbs_SOLID); s.More(); s.Next() ) {
2549       if ( myMesh->GetSubMesh( s.Current() )->IsEmpty() ) { // get faces of solid
2550         for ( TopExp_Explorer f( s.Current(), TopAbs_FACE); f.More(); f.Next() )
2551           faces.Add( f.Current() ); // in not meshed solid
2552       }
2553       else { // fix nodes in the solid and its faces
2554         MSG("FIX SOLID " << nbSolids-- << " #" << GetMeshDS()->ShapeToIndex(s.Current()));
2555         SMESH_MesherHelper h(*myMesh);
2556         h.SetSubShape( s.Current() );
2557         h.FixQuadraticElements(false);
2558       }
2559     }
2560     // fix nodes on geom faces
2561 #ifdef _DEBUG_
2562     int nbfaces = faces.Extent();
2563 #endif
2564     for ( TopTools_MapIteratorOfMapOfShape fIt( faces ); fIt.More(); fIt.Next() ) {
2565       MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key()));
2566       SMESH_MesherHelper h(*myMesh);
2567       h.SetSubShape( fIt.Key() );
2568       h.FixQuadraticElements(true);
2569     }
2570     //perf_print_all_meters(1);
2571     return;
2572   }
2573
2574   // 1. Find out type of elements and get iterator on them
2575   // ---------------------------------------------------
2576
2577   SMDS_ElemIteratorPtr elemIt;
2578   SMDSAbs_ElementType elemType = SMDSAbs_All;
2579
2580   SMESH_subMesh* submesh = myMesh->GetSubMeshContaining( myShapeID );
2581   if ( !submesh )
2582     return;
2583   if ( SMESHDS_SubMesh* smDS = submesh->GetSubMeshDS() ) {
2584     elemIt = smDS->GetElements();
2585     if ( elemIt->more() ) {
2586       elemType = elemIt->next()->GetType();
2587       elemIt = smDS->GetElements();
2588     }
2589   }
2590   if ( !elemIt || !elemIt->more() || elemType < SMDSAbs_Face )
2591     return;
2592
2593   // 2. Fill in auxiliary data structures
2594   // ----------------------------------
2595
2596   set< QLink > links;
2597   set< QFace > faces;
2598   set< QLink >::iterator pLink;
2599   set< QFace >::iterator pFace;
2600
2601   bool isCurved = false;
2602   //bool hasRectFaces = false;
2603   set<int> nbElemNodeSet;
2604
2605   if ( elemType == SMDSAbs_Volume )
2606   {
2607     SMDS_VolumeTool volTool;
2608     while ( elemIt->more() ) // loop on volumes
2609     {
2610       const SMDS_MeshElement* vol = elemIt->next();
2611       if ( !vol->IsQuadratic() || !volTool.Set( vol ))
2612         return; //continue;
2613       for ( int iF = 0; iF < volTool.NbFaces(); ++iF ) // loop on faces of volume
2614       {
2615         int nbN = volTool.NbFaceNodes( iF );
2616         nbElemNodeSet.insert( nbN );
2617         const SMDS_MeshNode** faceNodes = volTool.GetFaceNodes( iF );
2618         vector< const QLink* > faceLinks( nbN/2 );
2619         for ( int iN = 0; iN < nbN; iN += 2 ) // loop on links of a face
2620         {
2621           // store QLink
2622           QLink link( faceNodes[iN], faceNodes[iN+2], faceNodes[iN+1] );
2623           pLink = links.insert( link ).first;
2624           faceLinks[ iN/2 ] = & *pLink;
2625           if ( !isCurved )
2626             isCurved = !link.IsStraight();
2627           if ( link.MediumPos() == SMDS_TOP_3DSPACE && !link.IsStraight() )
2628             return; // already fixed
2629         }
2630         // store QFace
2631         pFace = faces.insert( QFace( faceLinks )).first;
2632         if ( pFace->NbVolumes() == 0 )
2633           pFace->AddSelfToLinks();
2634         pFace->SetVolume( vol );
2635 //         hasRectFaces = hasRectFaces ||
2636 //           ( volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_HEXA ||
2637 //             volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_PENTA );
2638 #ifdef _DEBUG_
2639         if ( nbN == 6 )
2640           pFace->_face = GetMeshDS()->FindFace(faceNodes[0],faceNodes[2],faceNodes[4]);
2641         else
2642           pFace->_face = GetMeshDS()->FindFace(faceNodes[0],faceNodes[2],
2643                                                faceNodes[4],faceNodes[6] );
2644 #endif
2645       }
2646     }
2647     set< QLink >::iterator pLink = links.begin();
2648     for ( ; pLink != links.end(); ++pLink )
2649       pLink->SetContinuesFaces();
2650   }
2651   else
2652   {
2653     while ( elemIt->more() ) // loop on faces
2654     {
2655       const SMDS_MeshElement* face = elemIt->next();
2656       if ( !face->IsQuadratic() )
2657         continue;
2658       nbElemNodeSet.insert( face->NbNodes() );
2659       int nbN = face->NbNodes()/2;
2660       vector< const QLink* > faceLinks( nbN );
2661       for ( int iN = 0; iN < nbN; ++iN ) // loop on links of a face
2662       {
2663         // store QLink
2664         QLink link( face->GetNode(iN), face->GetNode((iN+1)%nbN), face->GetNode(iN+nbN) );
2665         pLink = links.insert( link ).first;
2666         faceLinks[ iN ] = & *pLink;
2667         if ( !isCurved )
2668           isCurved = !link.IsStraight();
2669       }
2670       // store QFace
2671       pFace = faces.insert( QFace( faceLinks )).first;
2672       pFace->AddSelfToLinks();
2673       //hasRectFaces = ( hasRectFaces || nbN == 4 );
2674     }
2675   }
2676   if ( !isCurved )
2677     return; // no curved edges of faces
2678
2679   // 3. Compute displacement of medium nodes
2680   // -------------------------------------
2681
2682   // two loops on faces: the first is to treat boundary links, the second is for internal ones
2683   TopLoc_Location loc;
2684   // not treat boundary of volumic submesh
2685   int isInside = ( elemType == SMDSAbs_Volume && volumeOnly ) ? 1 : 0;
2686   for ( ; isInside < 2; ++isInside ) {
2687     MSG( "--------------- LOOP (inside=" << isInside << ") ------------------");
2688     SMDS_TypeOfPosition pos = isInside ? SMDS_TOP_3DSPACE : SMDS_TOP_FACE;
2689     SMDS_TypeOfPosition bndPos = isInside ? SMDS_TOP_FACE : SMDS_TOP_EDGE;
2690
2691     for ( pFace = faces.begin(); pFace != faces.end(); ++pFace ) {
2692       if ( bool(isInside) == pFace->IsBoundary() )
2693         continue;
2694       for ( int dir = 0; dir < 2; ++dir ) // 2 directions of propagation from quadrangle
2695       {
2696         MSG( "CHAIN");
2697         // make chain of links connected via continues faces
2698         int error = ERR_OK;
2699         TChain rawChain;
2700         if ( !pFace->GetLinkChain( dir, rawChain, pos, error) && error ==ERR_UNKNOWN ) continue;
2701         rawChain.reverse();
2702         if ( !pFace->GetLinkChain( dir+2, rawChain, pos, error ) && error ==ERR_UNKNOWN ) continue;
2703
2704         vector< TChain > chains;
2705         if ( error == ERR_OK ) { // chain contains continues rectangles
2706           chains.resize(1);
2707           chains[0].splice( chains[0].begin(), rawChain );
2708         }
2709         else if ( error == ERR_TRI ) {  // chain contains continues triangles
2710           TSplitTriaResult res = splitTrianglesIntoChains( rawChain, chains, pos );
2711           if ( res != _OK ) { // not quadrangles split into triangles
2712             fixTriaNearBoundary( rawChain, *this );
2713             break;
2714           }
2715         }
2716         else if ( error == ERR_PRISM ) { // quadrangle side faces of prisms
2717           fixPrism( rawChain );
2718           break;
2719         }
2720         else {
2721           continue;
2722         }
2723         for ( int iC = 0; iC < chains.size(); ++iC )
2724         {
2725           TChain& chain = chains[iC];
2726           if ( chain.empty() ) continue;
2727           if ( chain.front()->IsStraight() && chain.back()->IsStraight() ) {
2728             MSG("3D straight - ignore");
2729             continue;
2730           }
2731           if ( chain.front()->MediumPos() > bndPos ||
2732                chain.back()->MediumPos() > bndPos ) {
2733             MSG("Internal chain - ignore");
2734             continue;
2735           }
2736           // mesure chain length and compute link position along the chain
2737           double chainLen = 0;
2738           vector< double > linkPos;
2739           MSGBEG( "Link medium nodes: ");
2740           TChain::iterator link0 = chain.begin(), link1 = chain.begin(), link2;
2741           for ( ++link1; link1 != chain.end(); ++link1, ++link0 ) {
2742             MSGBEG( (*link0)->_mediumNode->GetID() << "-" <<(*link1)->_mediumNode->GetID()<<" ");
2743             double len = ((*link0)->MiddlePnt() - (*link1)->MiddlePnt()).Modulus();
2744             while ( len < numeric_limits<double>::min() ) { // remove degenerated link
2745               link1 = chain.erase( link1 );
2746               if ( link1 == chain.end() )
2747                 break;
2748               len = ((*link0)->MiddlePnt() - (*link1)->MiddlePnt()).Modulus();
2749             }
2750             chainLen += len;
2751             linkPos.push_back( chainLen );
2752           }
2753           MSG("");
2754           if ( linkPos.size() < 2 )
2755             continue;
2756
2757           gp_Vec move0 = chain.front()->_nodeMove;
2758           gp_Vec move1 = chain.back ()->_nodeMove;
2759
2760           TopoDS_Face face;
2761           bool checkUV = true;
2762           if ( !isInside ) {
2763             // compute node displacement of end links in parametric space of face
2764             const SMDS_MeshNode* nodeOnFace = (*(++chain.begin()))->_mediumNode;
2765             TopoDS_Shape f = GetSubShapeByNode( nodeOnFace, GetMeshDS() );
2766             if ( !f.IsNull() && f.ShapeType() == TopAbs_FACE )
2767             {
2768               face = TopoDS::Face( f );
2769               Handle(Geom_Surface) surf = BRep_Tool::Surface(face,loc);
2770               bool isStraight[2];
2771               for ( int is1 = 0; is1 < 2; ++is1 ) // move0 or move1
2772               {
2773                 TChainLink& link = is1 ? chain.back() : chain.front();
2774                 gp_XY uvm = GetNodeUV( face, link->_mediumNode, nodeOnFace, &checkUV);
2775                 gp_XY uv1 = GetNodeUV( face, link->node1(), nodeOnFace, &checkUV);
2776                 gp_XY uv2 = GetNodeUV( face, link->node2(), nodeOnFace, &checkUV);
2777                 gp_XY uv12 = GetMiddleUV( surf, uv1, uv2);
2778                 // uvMove = uvm - uv12
2779                 gp_XY uvMove = applyIn2D(surf, uvm, uv12, gp_XY_Subtracted, /*inPeriod=*/false);
2780                 ( is1 ? move1 : move0 ).SetCoord( uvMove.X(), uvMove.Y(), 0 );
2781                 if ( !is1 ) // correct nodeOnFace for move1 (issue 0020919)
2782                   nodeOnFace = (*(++chain.rbegin()))->_mediumNode;
2783                 isStraight[is1] = isStraightLink( (uv2-uv1).SquareModulus(),uvMove.SquareModulus());
2784               }
2785 //               if ( move0.SquareMagnitude() < straightTol2 &&
2786 //                    move1.SquareMagnitude() < straightTol2 ) {
2787               if ( isStraight[0] && isStraight[1] ) {
2788                 MSG("2D straight - ignore");
2789                 continue; // straight - no need to move nodes of internal links
2790               }
2791             }
2792           }
2793           gp_Trsf trsf;
2794           if ( isInside || face.IsNull() )
2795           {
2796             // compute node displacement of end links in their local coord systems
2797             {
2798               TChainLink& ln0 = chain.front(), ln1 = *(++chain.begin());
2799               trsf.SetTransformation( gp_Ax3( gp::Origin(), ln0.Normal(),
2800                                               gp_Vec( ln0->MiddlePnt(), ln1->MiddlePnt() )));
2801               move0.Transform(trsf);
2802             }
2803             {
2804               TChainLink& ln0 = *(++chain.rbegin()), ln1 = chain.back();
2805               trsf.SetTransformation( gp_Ax3( gp::Origin(), ln1.Normal(),
2806                                               gp_Vec( ln0->MiddlePnt(), ln1->MiddlePnt() )));
2807               move1.Transform(trsf);
2808             }
2809           }
2810           // compute displacement of medium nodes
2811           link2 = chain.begin();
2812           link0 = link2++;
2813           link1 = link2++;
2814           for ( int i = 0; link2 != chain.end(); ++link0, ++link1, ++link2, ++i )
2815           {
2816             double r = linkPos[i] / chainLen;
2817             // displacement in local coord system
2818             gp_Vec move = (1. - r) * move0 + r * move1;
2819             if ( isInside || face.IsNull()) {
2820               // transform to global
2821               gp_Vec x01( (*link0)->MiddlePnt(), (*link1)->MiddlePnt() );
2822               gp_Vec x12( (*link1)->MiddlePnt(), (*link2)->MiddlePnt() );
2823               gp_Vec x = x01.Normalized() + x12.Normalized();
2824               trsf.SetTransformation( gp_Ax3( gp::Origin(), link1->Normal(), x), gp_Ax3() );
2825               move.Transform(trsf);
2826             }
2827             else {
2828               // compute 3D displacement by 2D one
2829               Handle(Geom_Surface) s = BRep_Tool::Surface(face,loc);
2830               gp_XY oldUV   = GetNodeUV( face, (*link1)->_mediumNode, 0, &checkUV);
2831               gp_XY newUV   = applyIn2D( s, oldUV, gp_XY( move.X(),move.Y()), gp_XY_Added);
2832               gp_Pnt newPnt = s->Value( newUV.X(), newUV.Y());
2833               move = gp_Vec( XYZ((*link1)->_mediumNode), newPnt.Transformed(loc) );
2834 #ifdef _DEBUG_
2835               if ( (XYZ((*link1)->node1()) - XYZ((*link1)->node2())).SquareModulus() <
2836                    move.SquareMagnitude())
2837               {
2838                 gp_XY uv0 = GetNodeUV( face, (*link0)->_mediumNode, 0, &checkUV);
2839                 gp_XY uv2 = GetNodeUV( face, (*link2)->_mediumNode, 0, &checkUV);
2840                 MSG( "TOO LONG MOVE \t" <<
2841                      "uv0: "<<uv0.X()<<", "<<uv0.Y()<<" \t" <<
2842                      "uv2: "<<uv2.X()<<", "<<uv2.Y()<<" \t" <<
2843                      "uvOld: "<<oldUV.X()<<", "<<oldUV.Y()<<" \t" <<
2844                      "newUV: "<<newUV.X()<<", "<<newUV.Y()<<" \t");
2845               }
2846 #endif
2847             }
2848             (*link1)->Move( move );
2849             MSG( "Move " << (*link1)->_mediumNode->GetID() << " following "
2850                  << chain.front()->_mediumNode->GetID() <<"-"
2851                  << chain.back ()->_mediumNode->GetID() <<
2852                  " by " << move.Magnitude());
2853           }
2854         } // loop on chains of links
2855       } // loop on 2 directions of propagation from quadrangle
2856     } // loop on faces
2857   }
2858
2859   // 4. Move nodes
2860   // -----------
2861
2862   for ( pLink = links.begin(); pLink != links.end(); ++pLink ) {
2863     if ( pLink->IsMoved() ) {
2864       //gp_Pnt p = pLink->MediumPnt() + pLink->Move();
2865       gp_Pnt p = pLink->MiddlePnt() + pLink->Move();
2866       GetMeshDS()->MoveNode( pLink->_mediumNode, p.X(), p.Y(), p.Z());
2867     }
2868   }
2869 }
2870
2871 //=======================================================================
2872 /*!
2873  * \brief Iterator on ancestors of the given type
2874  */
2875 //=======================================================================
2876
2877 struct TAncestorsIterator : public SMDS_Iterator<const TopoDS_Shape*>
2878 {
2879   TopTools_ListIteratorOfListOfShape _ancIter;
2880   TopAbs_ShapeEnum                   _type;
2881   TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type)
2882     : _ancIter( ancestors ), _type( type )
2883   {
2884     if ( _ancIter.More() && _ancIter.Value().ShapeType() != _type ) next();
2885   }
2886   virtual bool more()
2887   {
2888     return _ancIter.More();
2889   }
2890   virtual const TopoDS_Shape* next()
2891   {
2892     const TopoDS_Shape* s = _ancIter.More() ? & _ancIter.Value() : 0;
2893     if ( _ancIter.More() )
2894       for ( _ancIter.Next();  _ancIter.More(); _ancIter.Next())
2895         if ( _ancIter.Value().ShapeType() == _type )
2896           break;
2897     return s;
2898   }
2899 };
2900
2901 //=======================================================================
2902 /*!
2903  * \brief Return iterator on ancestors of the given type
2904  */
2905 //=======================================================================
2906
2907 PShapeIteratorPtr SMESH_MesherHelper::GetAncestors(const TopoDS_Shape& shape,
2908                                                    const SMESH_Mesh&   mesh,
2909                                                    TopAbs_ShapeEnum    ancestorType)
2910 {
2911   return PShapeIteratorPtr( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType));
2912 }