Back to index

salome-smesh  6.5.0
SMESH_MesherHelper.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 // File:      SMESH_MesherHelper.cxx
00024 // Created:   15.02.06 15:22:41
00025 // Author:    Sergey KUUL
00026 //
00027 #include "SMESH_MesherHelper.hxx"
00028 
00029 #include "SMDS_FacePosition.hxx" 
00030 #include "SMDS_EdgePosition.hxx"
00031 #include "SMDS_VolumeTool.hxx"
00032 #include "SMESH_subMesh.hxx"
00033 #include "SMESH_ProxyMesh.hxx"
00034 
00035 #include <BRepAdaptor_Surface.hxx>
00036 #include <BRepTools.hxx>
00037 #include <BRepTools_WireExplorer.hxx>
00038 #include <BRep_Tool.hxx>
00039 #include <Geom2d_Curve.hxx>
00040 #include <GeomAPI_ProjectPointOnCurve.hxx>
00041 #include <GeomAPI_ProjectPointOnSurf.hxx>
00042 #include <Geom_Curve.hxx>
00043 #include <Geom_RectangularTrimmedSurface.hxx>
00044 #include <Geom_Surface.hxx>
00045 #include <ShapeAnalysis.hxx>
00046 #include <TopExp.hxx>
00047 #include <TopExp_Explorer.hxx>
00048 #include <TopTools_ListIteratorOfListOfShape.hxx>
00049 #include <TopTools_MapIteratorOfMapOfShape.hxx>
00050 #include <TopTools_MapOfShape.hxx>
00051 #include <TopoDS.hxx>
00052 #include <gp_Ax3.hxx>
00053 #include <gp_Pnt2d.hxx>
00054 #include <gp_Trsf.hxx>
00055 
00056 #include <Standard_Failure.hxx>
00057 #include <Standard_ErrorHandler.hxx>
00058 
00059 #include <utilities.h>
00060 
00061 #include <limits>
00062 
00063 using namespace std;
00064 
00065 #define RETURN_BAD_RESULT(msg) { MESSAGE(msg); return false; }
00066 
00067 namespace {
00068 
00069   gp_XYZ XYZ(const SMDS_MeshNode* n) { return gp_XYZ(n->X(), n->Y(), n->Z()); }
00070 
00071   enum { U_periodic = 1, V_periodic = 2 };
00072 }
00073 
00074 //================================================================================
00078 //================================================================================
00079 
00080 SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh)
00081   : myParIndex(0), myMesh(&theMesh), myShapeID(0), myCreateQuadratic(false)
00082 {
00083   myPar1[0] = myPar2[0] = myPar1[1] = myPar2[1] = 0;
00084   mySetElemOnShape = ( ! myMesh->HasShapeToMesh() );
00085 }
00086 
00087 //=======================================================================
00088 //function : ~SMESH_MesherHelper
00089 //purpose  : 
00090 //=======================================================================
00091 
00092 SMESH_MesherHelper::~SMESH_MesherHelper()
00093 {
00094   {
00095     TID2ProjectorOnSurf::iterator i_proj = myFace2Projector.begin();
00096     for ( ; i_proj != myFace2Projector.end(); ++i_proj )
00097       delete i_proj->second;
00098   }
00099   {
00100     TID2ProjectorOnCurve::iterator i_proj = myEdge2Projector.begin();
00101     for ( ; i_proj != myEdge2Projector.end(); ++i_proj )
00102       delete i_proj->second;
00103   }
00104 }
00105 
00106 //=======================================================================
00107 //function : IsQuadraticSubMesh
00108 //purpose  : Check submesh for given shape: if all elements on this shape 
00109 //           are quadratic, quadratic elements will be created.
00110 //           Also fill myTLinkNodeMap
00111 //=======================================================================
00112 
00113 bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
00114 {
00115   SMESHDS_Mesh* meshDS = GetMeshDS();
00116   // we can create quadratic elements only if all elements
00117   // created on sub-shapes of given shape are quadratic
00118   // also we have to fill myTLinkNodeMap
00119   myCreateQuadratic = true;
00120   mySeamShapeIds.clear();
00121   myDegenShapeIds.clear();
00122   TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
00123   SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge );
00124 
00125   int nbOldLinks = myTLinkNodeMap.size();
00126 
00127   if ( !myMesh->HasShapeToMesh() )
00128   {
00129     if (( myCreateQuadratic = myMesh->NbFaces( ORDER_QUADRATIC )))
00130     {
00131       SMDS_FaceIteratorPtr fIt = meshDS->facesIterator();
00132       while ( fIt->more() )
00133         AddTLinks( static_cast< const SMDS_MeshFace* >( fIt->next() ));
00134     }
00135   }
00136   else
00137   {
00138     TopExp_Explorer exp( aSh, subType );
00139     TopTools_MapOfShape checkedSubShapes;
00140     for (; exp.More() && myCreateQuadratic; exp.Next()) {
00141       if ( !checkedSubShapes.Add( exp.Current() ))
00142         continue; // needed if aSh is compound of solids
00143       if ( SMESHDS_SubMesh * subMesh = meshDS->MeshElements( exp.Current() )) {
00144         if ( SMDS_ElemIteratorPtr it = subMesh->GetElements() ) {
00145           while(it->more()) {
00146             const SMDS_MeshElement* e = it->next();
00147             if ( e->GetType() != elemType || !e->IsQuadratic() ) {
00148               myCreateQuadratic = false;
00149               break;
00150             }
00151             else {
00152               // fill TLinkNodeMap
00153               switch ( e->NbNodes() ) {
00154               case 3:
00155                 AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(2)); break;
00156               case 6:
00157                 AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(3));
00158                 AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(4));
00159                 AddTLinkNode(e->GetNode(2),e->GetNode(0),e->GetNode(5)); break;
00160               case 8:
00161                 AddTLinkNode(e->GetNode(0),e->GetNode(1),e->GetNode(4));
00162                 AddTLinkNode(e->GetNode(1),e->GetNode(2),e->GetNode(5));
00163                 AddTLinkNode(e->GetNode(2),e->GetNode(3),e->GetNode(6));
00164                 AddTLinkNode(e->GetNode(3),e->GetNode(0),e->GetNode(7));
00165                 break;
00166               default:
00167                 myCreateQuadratic = false;
00168                 break;
00169               }
00170             }
00171           }
00172         }
00173       }
00174     }
00175   }
00176 
00177   if ( nbOldLinks == myTLinkNodeMap.size() )
00178     myCreateQuadratic = false;
00179 
00180   if(!myCreateQuadratic) {
00181     myTLinkNodeMap.clear();
00182   }
00183   SetSubShape( aSh );
00184 
00185   return myCreateQuadratic;
00186 }
00187 
00188 //=======================================================================
00189 //function : SetSubShape
00190 //purpose  : Set geomerty to make elements on
00191 //=======================================================================
00192 
00193 void SMESH_MesherHelper::SetSubShape(const int aShID)
00194 {
00195   if ( aShID == myShapeID )
00196     return;
00197   if ( aShID > 0 )
00198     SetSubShape( GetMeshDS()->IndexToShape( aShID ));
00199   else
00200     SetSubShape( TopoDS_Shape() );
00201 }
00202 
00203 //=======================================================================
00204 //function : SetSubShape
00205 //purpose  : Set geomerty to create elements on
00206 //=======================================================================
00207 
00208 void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
00209 {
00210   if ( myShape.IsSame( aSh ))
00211     return;
00212 
00213   myShape = aSh;
00214   mySeamShapeIds.clear();
00215   myDegenShapeIds.clear();
00216 
00217   if ( myShape.IsNull() ) {
00218     myShapeID  = 0;
00219     return;
00220   }
00221   SMESHDS_Mesh* meshDS = GetMeshDS();
00222   myShapeID = meshDS->ShapeToIndex(aSh);
00223   myParIndex = 0;
00224 
00225   // treatment of periodic faces
00226   for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
00227   {
00228     const TopoDS_Face& face = TopoDS::Face( eF.Current() );
00229     TopLoc_Location loc;
00230     Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
00231 
00232     if ( surface->IsUPeriodic() || surface->IsVPeriodic() ||
00233          surface->IsUClosed()   || surface->IsVClosed() )
00234     {
00235       //while ( surface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface )))
00236       //surface = Handle(Geom_RectangularTrimmedSurface)::DownCast( surface )->BasisSurface();
00237       GeomAdaptor_Surface surf( surface );
00238 
00239       for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next())
00240       {
00241         // look for a seam edge
00242         const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
00243         if ( BRep_Tool::IsClosed( edge, face )) {
00244           // initialize myPar1, myPar2 and myParIndex
00245           gp_Pnt2d uv1, uv2;
00246           BRep_Tool::UVPoints( edge, face, uv1, uv2 );
00247           if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) ))
00248           {
00249             myParIndex |= U_periodic;
00250             myPar1[0] = surf.FirstUParameter();
00251             myPar2[0] = surf.LastUParameter();
00252           }
00253           else {
00254             myParIndex |= V_periodic;
00255             myPar1[1] = surf.FirstVParameter();
00256             myPar2[1] = surf.LastVParameter();
00257           }
00258           // store seam shape indices, negative if shape encounters twice
00259           int edgeID = meshDS->ShapeToIndex( edge );
00260           mySeamShapeIds.insert( IsSeamShape( edgeID ) ? -edgeID : edgeID );
00261           for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() ) {
00262             int vertexID = meshDS->ShapeToIndex( v.Current() );
00263             mySeamShapeIds.insert( IsSeamShape( vertexID ) ? -vertexID : vertexID );
00264           }
00265         }
00266 
00267         // look for a degenerated edge
00268         if ( BRep_Tool::Degenerated( edge )) {
00269           myDegenShapeIds.insert( meshDS->ShapeToIndex( edge ));
00270           for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() )
00271             myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
00272         }
00273       }
00274     }
00275   }
00276 }
00277 
00278 //=======================================================================
00279 //function : GetNodeUVneedInFaceNode
00280 //purpose  : Check if inFaceNode argument is necessary for call GetNodeUV(F,..)
00281 //           Return true if the face is periodic.
00282 //           If F is Null, answer about sub-shape set through IsQuadraticSubMesh() or
00283 //           * SetSubShape()
00284 //=======================================================================
00285 
00286 bool SMESH_MesherHelper::GetNodeUVneedInFaceNode(const TopoDS_Face& F) const
00287 {
00288   if ( F.IsNull() ) return !mySeamShapeIds.empty();
00289 
00290   if ( !F.IsNull() && !myShape.IsNull() && myShape.IsSame( F ))
00291     return !mySeamShapeIds.empty();
00292 
00293   TopLoc_Location loc;
00294   Handle(Geom_Surface) aSurface = BRep_Tool::Surface( F,loc );
00295   if ( !aSurface.IsNull() )
00296     return ( aSurface->IsUPeriodic() || aSurface->IsVPeriodic() );
00297 
00298   return false;
00299 }
00300 
00301 //=======================================================================
00302 //function : IsMedium
00303 //purpose  : 
00304 //=======================================================================
00305 
00306 bool SMESH_MesherHelper::IsMedium(const SMDS_MeshNode*      node,
00307                                   const SMDSAbs_ElementType typeToCheck)
00308 {
00309   return SMESH_MeshEditor::IsMedium( node, typeToCheck );
00310 }
00311 
00312 //=======================================================================
00313 //function : GetSubShapeByNode
00314 //purpose  : Return support shape of a node
00315 //=======================================================================
00316 
00317 TopoDS_Shape SMESH_MesherHelper::GetSubShapeByNode(const SMDS_MeshNode* node,
00318                                                    const SMESHDS_Mesh*  meshDS)
00319 {
00320   int shapeID = node->getshapeId();
00321   if ( 0 < shapeID && shapeID <= meshDS->MaxShapeIndex() )
00322     return meshDS->IndexToShape( shapeID );
00323   else
00324     return TopoDS_Shape();
00325 }
00326 
00327 
00328 //=======================================================================
00329 //function : AddTLinkNode
00330 //purpose  : add a link in my data structure
00331 //=======================================================================
00332 
00333 void SMESH_MesherHelper::AddTLinkNode(const SMDS_MeshNode* n1,
00334                                       const SMDS_MeshNode* n2,
00335                                       const SMDS_MeshNode* n12)
00336 {
00337   // add new record to map
00338   SMESH_TLink link( n1, n2 );
00339   myTLinkNodeMap.insert( make_pair(link,n12));
00340 }
00341 
00342 //================================================================================
00346 //================================================================================
00347 
00348 void SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge)
00349 {
00350   if ( edge->IsQuadratic() )
00351     AddTLinkNode(edge->GetNode(0), edge->GetNode(1), edge->GetNode(2));
00352 }
00353 
00354 //================================================================================
00358 //================================================================================
00359 
00360 void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
00361 {
00362   if ( !f->IsPoly() )
00363     switch ( f->NbNodes() ) {
00364     case 6:
00365       AddTLinkNode(f->GetNode(0),f->GetNode(1),f->GetNode(3));
00366       AddTLinkNode(f->GetNode(1),f->GetNode(2),f->GetNode(4));
00367       AddTLinkNode(f->GetNode(2),f->GetNode(0),f->GetNode(5)); break;
00368     case 8:
00369       AddTLinkNode(f->GetNode(0),f->GetNode(1),f->GetNode(4));
00370       AddTLinkNode(f->GetNode(1),f->GetNode(2),f->GetNode(5));
00371       AddTLinkNode(f->GetNode(2),f->GetNode(3),f->GetNode(6));
00372       AddTLinkNode(f->GetNode(3),f->GetNode(0),f->GetNode(7));
00373     default:;
00374     }
00375 }
00376 
00377 //================================================================================
00381 //================================================================================
00382 
00383 void SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume)
00384 {
00385   if ( volume->IsQuadratic() )
00386   {
00387     SMDS_VolumeTool vTool( volume );
00388     const SMDS_MeshNode** nodes = vTool.GetNodes();
00389     set<int> addedLinks;
00390     for ( int iF = 1; iF < vTool.NbFaces(); ++iF )
00391     {
00392       const int nbN = vTool.NbFaceNodes( iF );
00393       const int* iNodes = vTool.GetFaceNodesIndices( iF );
00394       for ( int i = 0; i < nbN; )
00395       {
00396         int iN1  = iNodes[i++];
00397         int iN12 = iNodes[i++];
00398         int iN2  = iNodes[i++];
00399         if ( iN1 > iN2 ) std::swap( iN1, iN2 );
00400         int linkID = iN1 * vTool.NbNodes() + iN2;
00401         pair< set<int>::iterator, bool > it_isNew = addedLinks.insert( linkID );
00402         if ( it_isNew.second )
00403           AddTLinkNode( nodes[iN1], nodes[iN2], nodes[iN12] );
00404         else
00405           addedLinks.erase( it_isNew.first ); // each link encounters only twice
00406       }
00407     }
00408   }
00409 }
00410 
00411 //================================================================================
00416 //================================================================================
00417 
00418 bool SMESH_MesherHelper::toCheckPosOnShape(int shapeID ) const
00419 {
00420   map< int,bool >::const_iterator id_ok = myNodePosShapesValidity.find( shapeID );
00421   return ( id_ok == myNodePosShapesValidity.end() || !id_ok->second );
00422 }
00423 
00424 //================================================================================
00429 //================================================================================
00430 
00431 void SMESH_MesherHelper::setPosOnShapeValidity(int shapeID, bool ok ) const
00432 {
00433   ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok));
00434 }
00435 
00436 //=======================================================================
00437 //function : GetUVOnSeam
00438 //purpose  : Select UV on either of 2 pcurves of a seam edge, closest to the given UV
00439 //=======================================================================
00440 
00441 gp_Pnt2d SMESH_MesherHelper::GetUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& uv2 ) const
00442 {
00443   gp_Pnt2d result = uv1;
00444   for ( int i = U_periodic; i <= V_periodic ; ++i )
00445   {
00446     if ( myParIndex & i )
00447     {
00448       double p1 = uv1.Coord( i );
00449       double dp1 = Abs( p1-myPar1[i-1]), dp2 = Abs( p1-myPar2[i-1]);
00450       if ( myParIndex == i ||
00451            dp1 < ( myPar2[i-1] - myPar2[i-1] ) / 100. ||
00452            dp2 < ( myPar2[i-1] - myPar2[i-1] ) / 100. )
00453       {
00454         double p2 = uv2.Coord( i );
00455         double p1Alt = ( dp1 < dp2 ) ? myPar2[i-1] : myPar1[i-1];
00456         if ( Abs( p2 - p1 ) > Abs( p2 - p1Alt ))
00457           result.SetCoord( i, p1Alt );
00458       }
00459     }
00460   }
00461   return result;
00462 }
00463 
00464 //=======================================================================
00465 //function : GetNodeUV
00466 //purpose  : Return node UV on face
00467 //=======================================================================
00468 
00469 gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
00470                                     const SMDS_MeshNode* n,
00471                                     const SMDS_MeshNode* n2,
00472                                     bool*                check) const
00473 {
00474   gp_Pnt2d uv( Precision::Infinite(), Precision::Infinite() );
00475   const SMDS_PositionPtr Pos = n->GetPosition();
00476   bool uvOK = false;
00477   if(Pos->GetTypeOfPosition()==SMDS_TOP_FACE)
00478   {
00479     // node has position on face
00480     const SMDS_FacePosition* fpos =
00481       static_cast<const SMDS_FacePosition*>(n->GetPosition());
00482     uv.SetCoord(fpos->GetUParameter(),fpos->GetVParameter());
00483     if ( check )
00484       uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ));
00485   }
00486   else if(Pos->GetTypeOfPosition()==SMDS_TOP_EDGE)
00487   {
00488     // node has position on edge => it is needed to find
00489     // corresponding edge from face, get pcurve for this
00490     // edge and retrieve value from this pcurve
00491     const SMDS_EdgePosition* epos =
00492       static_cast<const SMDS_EdgePosition*>(n->GetPosition());
00493     int edgeID = n->getshapeId();
00494     TopoDS_Edge E = TopoDS::Edge(GetMeshDS()->IndexToShape(edgeID));
00495     double f, l, u = epos->GetUParameter();
00496     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
00497     bool validU = ( f < u && u < l );
00498     if ( validU )
00499       uv = C2d->Value( u );
00500     else
00501       uv.SetCoord( Precision::Infinite(),0.);
00502     if ( check || !validU )
00503       uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ),/*force=*/ !validU );
00504 
00505     // for a node on a seam edge select one of UVs on 2 pcurves
00506     if ( n2 && IsSeamShape( edgeID ) )
00507     {
00508       uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0, check ));
00509     }
00510     else
00511     { // adjust uv to period
00512       TopLoc_Location loc;
00513       Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
00514       Standard_Boolean isUPeriodic = S->IsUPeriodic();
00515       Standard_Boolean isVPeriodic = S->IsVPeriodic();
00516       if ( isUPeriodic || isVPeriodic ) {
00517         Standard_Real UF,UL,VF,VL;
00518         S->Bounds(UF,UL,VF,VL);
00519         if(isUPeriodic)
00520           uv.SetX( uv.X() + ShapeAnalysis::AdjustToPeriod(uv.X(),UF,UL));
00521         if(isVPeriodic)
00522           uv.SetY( uv.Y() + ShapeAnalysis::AdjustToPeriod(uv.Y(),VF,VL));
00523       }
00524     }
00525   }
00526   else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX)
00527   {
00528     if ( int vertexID = n->getshapeId() ) {
00529       const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
00530       try {
00531         uv = BRep_Tool::Parameters( V, F );
00532         uvOK = true;
00533       }
00534       catch (Standard_Failure& exc) {
00535       }
00536       if ( !uvOK ) {
00537         for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !uvOK && vert.More(); vert.Next() )
00538           uvOK = ( V == vert.Current() );
00539         if ( !uvOK ) {
00540 #ifdef _DEBUG_
00541           MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID
00542                << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
00543 #endif
00544           // get UV of a vertex closest to the node
00545           double dist = 1e100;
00546           gp_Pnt pn = XYZ( n );
00547           for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !uvOK && vert.More(); vert.Next() ) {
00548             TopoDS_Vertex curV = TopoDS::Vertex( vert.Current() );
00549             gp_Pnt p = BRep_Tool::Pnt( curV );
00550             double curDist = p.SquareDistance( pn );
00551             if ( curDist < dist ) {
00552               dist = curDist;
00553               uv = BRep_Tool::Parameters( curV, F );
00554               uvOK = ( dist < DBL_MIN );
00555             }
00556           }
00557         }
00558         else {
00559           uvOK = false;
00560           TopTools_ListIteratorOfListOfShape it( myMesh->GetAncestors( V ));
00561           for ( ; it.More(); it.Next() ) {
00562             if ( it.Value().ShapeType() == TopAbs_EDGE ) {
00563               const TopoDS_Edge & edge = TopoDS::Edge( it.Value() );
00564               double f,l;
00565               Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(edge, F, f, l);
00566               if ( !C2d.IsNull() ) {
00567                 double u = ( V == TopExp::FirstVertex( edge ) ) ?  f : l;
00568                 uv = C2d->Value( u );
00569                 uvOK = true;
00570                 break;
00571               }
00572             }
00573           }
00574         }
00575       }
00576       if ( n2 && IsSeamShape( vertexID ) )
00577         uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
00578     }
00579   }
00580   else
00581   {
00582     uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ));
00583   }
00584 
00585   if ( check )
00586     *check = uvOK;
00587 
00588   return uv.XY();
00589 }
00590 
00591 //=======================================================================
00592 //function : CheckNodeUV
00593 //purpose  : Check and fix node UV on a face
00594 //=======================================================================
00595 
00596 bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
00597                                      const SMDS_MeshNode* n,
00598                                      gp_XY&               uv,
00599                                      const double         tol,
00600                                      const bool           force,
00601                                      double               distXYZ[4]) const
00602 {
00603   int shapeID = n->getshapeId();
00604   bool infinit = ( Precision::IsInfinite( uv.X() ) || Precision::IsInfinite( uv.Y() ));
00605   if ( force || toCheckPosOnShape( shapeID ) || infinit )
00606   {
00607     // check that uv is correct
00608     TopLoc_Location loc;
00609     Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc );
00610     gp_Pnt nodePnt = XYZ( n ), surfPnt(0,0,0);
00611     double dist = 0;
00612     if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
00613     if ( infinit ||
00614          (dist = nodePnt.Distance( surfPnt = surface->Value( uv.X(), uv.Y() ))) > tol )
00615     {
00616       setPosOnShapeValidity( shapeID, false );
00617       if ( !infinit && distXYZ ) {
00618         surfPnt.Transform( loc );
00619         distXYZ[0] = dist;
00620         distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z();
00621       }
00622       // uv incorrect, project the node to surface
00623       GeomAPI_ProjectPointOnSurf& projector = GetProjector( F, loc, tol );
00624       projector.Perform( nodePnt );
00625       if ( !projector.IsDone() || projector.NbPoints() < 1 )
00626       {
00627         MESSAGE( "SMESH_MesherHelper::CheckNodeUV() failed to project" );
00628         return false;
00629       }
00630       Quantity_Parameter U,V;
00631       projector.LowerDistanceParameters(U,V);
00632       uv.SetCoord( U,V );
00633       surfPnt = surface->Value( U, V );
00634       dist = nodePnt.Distance( surfPnt );
00635       if ( distXYZ ) {
00636         surfPnt.Transform( loc );
00637         distXYZ[0] = dist;
00638         distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z();
00639       }
00640       if ( dist > tol )
00641       {
00642         MESSAGE( "SMESH_MesherHelper::CheckNodeUV(), invalid projection" );
00643         return false;
00644       }
00645       // store the fixed UV on the face
00646       if ( myShape.IsSame(F) && shapeID == myShapeID )
00647         const_cast<SMDS_MeshNode*>(n)->SetPosition
00648           ( SMDS_PositionPtr( new SMDS_FacePosition( U, V )));
00649     }
00650     else if ( uv.Modulus() > numeric_limits<double>::min() )
00651     {
00652       setPosOnShapeValidity( shapeID, true );
00653     }
00654   }
00655   return true;
00656 }
00657 
00658 //=======================================================================
00659 //function : GetProjector
00660 //purpose  : Return projector intitialized by given face without location, which is returned
00661 //=======================================================================
00662 
00663 GeomAPI_ProjectPointOnSurf& SMESH_MesherHelper::GetProjector(const TopoDS_Face& F,
00664                                                              TopLoc_Location&   loc,
00665                                                              double             tol ) const
00666 {
00667   Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc );
00668   int faceID = GetMeshDS()->ShapeToIndex( F );
00669   TID2ProjectorOnSurf& i2proj = const_cast< TID2ProjectorOnSurf&>( myFace2Projector );
00670   TID2ProjectorOnSurf::iterator i_proj = i2proj.find( faceID );
00671   if ( i_proj == i2proj.end() )
00672   {
00673     if ( tol == 0 ) tol = BRep_Tool::Tolerance( F );
00674     double U1, U2, V1, V2;
00675     surface->Bounds(U1, U2, V1, V2);
00676     GeomAPI_ProjectPointOnSurf* proj = new GeomAPI_ProjectPointOnSurf();
00677     proj->Init( surface, U1, U2, V1, V2, tol );
00678     i_proj = i2proj.insert( make_pair( faceID, proj )).first;
00679   }
00680   return *( i_proj->second );
00681 }
00682 
00683 namespace
00684 {
00685   gp_XY AverageUV(const gp_XY& uv1, const gp_XY& uv2) { return ( uv1 + uv2 ) / 2.; }
00686   gp_XY_FunPtr(Added); // define gp_XY_Added pointer to function calling gp_XY::Added(gp_XY)
00687   gp_XY_FunPtr(Subtracted); 
00688 }
00689 
00690 //=======================================================================
00691 //function : applyIn2D
00692 //purpose  : Perform given operation on two 2d points in parameric space of given surface.
00693 //           It takes into account period of the surface. Use gp_XY_FunPtr macro
00694 //           to easily define pointer to function of gp_XY class.
00695 //=======================================================================
00696 
00697 gp_XY SMESH_MesherHelper::applyIn2D(const Handle(Geom_Surface)& surface,
00698                                     const gp_XY&                uv1,
00699                                     const gp_XY&                uv2,
00700                                     xyFunPtr                    fun,
00701                                     const bool                  resultInPeriod)
00702 {
00703   Standard_Boolean isUPeriodic = surface.IsNull() ? false : surface->IsUPeriodic();
00704   Standard_Boolean isVPeriodic = surface.IsNull() ? false : surface->IsVPeriodic();
00705   if ( !isUPeriodic && !isVPeriodic )
00706     return fun(uv1,uv2);
00707 
00708   // move uv2 not far than half-period from uv1
00709   double u2 = 
00710     uv2.X()+(isUPeriodic ? ShapeAnalysis::AdjustByPeriod(uv2.X(),uv1.X(),surface->UPeriod()) :0);
00711   double v2 = 
00712     uv2.Y()+(isVPeriodic ? ShapeAnalysis::AdjustByPeriod(uv2.Y(),uv1.Y(),surface->VPeriod()) :0);
00713 
00714   // execute operation
00715   gp_XY res = fun( uv1, gp_XY(u2,v2) );
00716 
00717   // move result within period
00718   if ( resultInPeriod )
00719   {
00720     Standard_Real UF,UL,VF,VL;
00721     surface->Bounds(UF,UL,VF,VL);
00722     if ( isUPeriodic )
00723       res.SetX( res.X() + ShapeAnalysis::AdjustToPeriod(res.X(),UF,UL));
00724     if ( isVPeriodic )
00725       res.SetY( res.Y() + ShapeAnalysis::AdjustToPeriod(res.Y(),VF,VL));
00726   }
00727 
00728   return res;
00729 }
00730 //=======================================================================
00731 //function : GetMiddleUV
00732 //purpose  : Return middle UV taking in account surface period
00733 //=======================================================================
00734 
00735 gp_XY SMESH_MesherHelper::GetMiddleUV(const Handle(Geom_Surface)& surface,
00736                                       const gp_XY&                p1,
00737                                       const gp_XY&                p2)
00738 {
00739   // NOTE:
00740   // the proper place of getting basic surface seems to be in applyIn2D()
00741   // but we put it here to decrease a risk of regressions just before releasing a version
00742   Handle(Geom_Surface) surf = surface;
00743   while ( !surf.IsNull() && surf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface )))
00744     surf = Handle(Geom_RectangularTrimmedSurface)::DownCast( surf )->BasisSurface();
00745 
00746   return applyIn2D( surf, p1, p2, & AverageUV );
00747 }
00748 
00749 //=======================================================================
00750 //function : GetNodeU
00751 //purpose  : Return node U on edge
00752 //=======================================================================
00753 
00754 double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
00755                                     const SMDS_MeshNode* n,
00756                                     const SMDS_MeshNode* inEdgeNode,
00757                                     bool*                check)
00758 {
00759   double param = 0;
00760   const SMDS_PositionPtr pos = n->GetPosition();
00761   if ( pos->GetTypeOfPosition()==SMDS_TOP_EDGE )
00762   {
00763     const SMDS_EdgePosition* epos = static_cast<const SMDS_EdgePosition*>( pos );
00764     param =  epos->GetUParameter();
00765   }
00766   else if( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
00767   {
00768     if ( inEdgeNode && TopExp::FirstVertex( E ).IsSame( TopExp::LastVertex( E ))) // issue 0020128
00769     {
00770       Standard_Real f,l;
00771       BRep_Tool::Range( E, f,l );
00772       double uInEdge = GetNodeU( E, inEdgeNode );
00773       param = ( fabs( uInEdge - f ) < fabs( l - uInEdge )) ? f : l;
00774     }
00775     else
00776     {
00777       SMESHDS_Mesh * meshDS = GetMeshDS();
00778       int vertexID = n->getshapeId();
00779       const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
00780       param =  BRep_Tool::Parameter( V, E );
00781     }
00782   }
00783   if ( check )
00784   {
00785     double tol = BRep_Tool::Tolerance( E );
00786     double f,l;  BRep_Tool::Range( E, f,l );
00787     bool force = ( param < f-tol || param > l+tol );
00788     if ( !force && pos->GetTypeOfPosition()==SMDS_TOP_EDGE )
00789       force = ( GetMeshDS()->ShapeToIndex( E ) != n->getshapeId() );
00790 
00791     *check = CheckNodeU( E, n, param, 2*tol, force );
00792   }
00793   return param;
00794 }
00795 
00796 //=======================================================================
00797 //function : CheckNodeU
00798 //purpose  : Check and fix node U on an edge
00799 //           Return false if U is bad and could not be fixed
00800 //=======================================================================
00801 
00802 bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
00803                                     const SMDS_MeshNode* n,
00804                                     double&              u,
00805                                     const double         tol,
00806                                     const bool           force,
00807                                     double               distXYZ[4]) const
00808 {
00809   int shapeID = n->getshapeId();
00810   if ( force || toCheckPosOnShape( shapeID ))
00811   {
00812     TopLoc_Location loc; double f,l;
00813     Handle(Geom_Curve) curve = BRep_Tool::Curve( E,loc,f,l );
00814     if ( curve.IsNull() ) // degenerated edge
00815     {
00816       if ( u+tol < f || u-tol > l )
00817       {
00818         double r = Max( 0.5, 1 - tol*n->GetID()); // to get a unique u on edge
00819         u =  f*r + l*(1-r);
00820       }
00821     }
00822     else
00823     {
00824       gp_Pnt nodePnt = SMESH_TNodeXYZ( n );
00825       if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
00826       gp_Pnt curvPnt = curve->Value( u );
00827       double dist = nodePnt.Distance( curvPnt );
00828       if ( distXYZ ) {
00829         curvPnt.Transform( loc );
00830         distXYZ[0] = dist;
00831         distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
00832       }
00833       if ( dist > tol )
00834       {
00835         setPosOnShapeValidity( shapeID, false );
00836         // u incorrect, project the node to the curve
00837         int edgeID = GetMeshDS()->ShapeToIndex( E );
00838         TID2ProjectorOnCurve& i2proj = const_cast< TID2ProjectorOnCurve&>( myEdge2Projector );
00839         TID2ProjectorOnCurve::iterator i_proj =
00840           i2proj.insert( make_pair( edgeID, (GeomAPI_ProjectPointOnCurve*) 0 )).first;
00841         if ( !i_proj->second  )
00842         {
00843           i_proj->second = new GeomAPI_ProjectPointOnCurve();
00844           i_proj->second->Init( curve, f, l );
00845         }
00846         GeomAPI_ProjectPointOnCurve* projector = i_proj->second;
00847         projector->Perform( nodePnt );
00848         if ( projector->NbPoints() < 1 )
00849         {
00850           MESSAGE( "SMESH_MesherHelper::CheckNodeU() failed to project" );
00851           return false;
00852         }
00853         Quantity_Parameter U = projector->LowerDistanceParameter();
00854         u = double( U );
00855         curvPnt = curve->Value( u );
00856         dist = nodePnt.Distance( curvPnt );
00857         if ( distXYZ ) {
00858           curvPnt.Transform( loc );
00859           distXYZ[0] = dist;
00860           distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
00861         }
00862         if ( dist > tol )
00863         {
00864           MESSAGE( "SMESH_MesherHelper::CheckNodeU(), invalid projection" );
00865           MESSAGE("distance " << dist << " " << tol );
00866           return false;
00867         }
00868         // store the fixed U on the edge
00869         if ( myShape.IsSame(E) && shapeID == myShapeID )
00870           const_cast<SMDS_MeshNode*>(n)->SetPosition
00871             ( SMDS_PositionPtr( new SMDS_EdgePosition( U )));
00872       }
00873       else if ( fabs( u ) > numeric_limits<double>::min() )
00874       {
00875         setPosOnShapeValidity( shapeID, true );
00876       }
00877       if (( u < f-tol || u > l+tol ) && force )
00878       {
00879         // node is on vertex but is set on periodic but trimmed edge (issue 0020890)
00880         try
00881         {
00882           // do not use IsPeriodic() as Geom_TrimmedCurve::IsPeriodic () returns false
00883           double period = curve->Period();
00884           u = ( u < f ) ? u + period : u - period;
00885         }
00886         catch (Standard_Failure& exc)
00887         {
00888           return false;
00889         }
00890       }
00891     }
00892   }
00893   return true;
00894 }
00895 
00896 //=======================================================================
00897 //function : GetMediumPos
00898 //purpose  : Return index and type of the shape  (EDGE or FACE only) to
00899 //          set a medium node on
00900 //=======================================================================
00901 
00902 std::pair<int, TopAbs_ShapeEnum> SMESH_MesherHelper::GetMediumPos(const SMDS_MeshNode* n1,
00903                                                                   const SMDS_MeshNode* n2)
00904 {
00905   TopAbs_ShapeEnum shapeType = TopAbs_SHAPE;
00906   int              shapeID = -1;
00907   TopoDS_Shape     shape;
00908 
00909   if (( myShapeID == n1->getshapeId() || myShapeID == n2->getshapeId() ) && myShapeID > 0 )
00910   {
00911     shapeType = myShape.ShapeType();
00912     shapeID   = myShapeID;
00913   }
00914   else if ( n1->getshapeId() == n2->getshapeId() )
00915   {
00916     shapeID = n2->getshapeId();
00917     shape = GetSubShapeByNode( n1, GetMeshDS() );
00918   }
00919   else
00920   {
00921     const SMDS_TypeOfPosition Pos1 = n1->GetPosition()->GetTypeOfPosition();
00922     const SMDS_TypeOfPosition Pos2 = n2->GetPosition()->GetTypeOfPosition();
00923 
00924     if ( Pos1 == SMDS_TOP_3DSPACE || Pos2 == SMDS_TOP_3DSPACE )
00925     {
00926     }
00927     else if ( Pos1 == SMDS_TOP_FACE || Pos2 == SMDS_TOP_FACE )
00928     {
00929       if ( Pos1 != SMDS_TOP_FACE || Pos2 != SMDS_TOP_FACE )
00930       {
00931         if ( Pos1 != SMDS_TOP_FACE ) std::swap( n1,n2 );
00932         TopoDS_Shape F = GetSubShapeByNode( n1, GetMeshDS() );
00933         TopoDS_Shape S = GetSubShapeByNode( n2, GetMeshDS() );
00934         if ( IsSubShape( S, F ))
00935         {
00936           shapeType = TopAbs_FACE;
00937           shapeID   = n1->getshapeId();
00938         }
00939       }
00940     }
00941     else if ( Pos1 == SMDS_TOP_EDGE && Pos2 == SMDS_TOP_EDGE )
00942     {
00943       TopoDS_Shape E1 = GetSubShapeByNode( n1, GetMeshDS() );
00944       TopoDS_Shape E2 = GetSubShapeByNode( n2, GetMeshDS() );
00945       shape = GetCommonAncestor( E1, E2, *myMesh, TopAbs_FACE );
00946     }
00947     else if ( Pos1 == SMDS_TOP_VERTEX && Pos2 == SMDS_TOP_VERTEX )
00948     {
00949       TopoDS_Shape V1 = GetSubShapeByNode( n1, GetMeshDS() );
00950       TopoDS_Shape V2 = GetSubShapeByNode( n2, GetMeshDS() );
00951       shape = GetCommonAncestor( V1, V2, *myMesh, TopAbs_EDGE );
00952       if ( shape.IsNull() ) shape = GetCommonAncestor( V1, V2, *myMesh, TopAbs_FACE );
00953     }
00954     else // VERTEX and EDGE
00955     {
00956       if ( Pos1 != SMDS_TOP_VERTEX ) std::swap( n1,n2 );
00957       TopoDS_Shape V = GetSubShapeByNode( n1, GetMeshDS() );
00958       TopoDS_Shape E = GetSubShapeByNode( n2, GetMeshDS() );
00959       if ( IsSubShape( V, E ))
00960         shape = E;
00961       else
00962         shape = GetCommonAncestor( V, E, *myMesh, TopAbs_FACE );
00963     }
00964   }
00965 
00966   if ( !shape.IsNull() )
00967   {
00968     if ( shapeID < 1 )
00969       shapeID = GetMeshDS()->ShapeToIndex( shape );
00970     shapeType = shape.ShapeType();
00971   }
00972   return make_pair( shapeID, shapeType );
00973 }
00974 
00975 //=======================================================================
00976 //function : GetMediumNode
00977 //purpose  : Return existing or create new medium nodes between given ones
00978 //=======================================================================
00979 
00980 const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
00981                                                        const SMDS_MeshNode* n2,
00982                                                        bool                 force3d)
00983 {
00984   // Find existing node
00985 
00986   SMESH_TLink link(n1,n2);
00987   ItTLinkNode itLN = myTLinkNodeMap.find( link );
00988   if ( itLN != myTLinkNodeMap.end() ) {
00989     return (*itLN).second;
00990   }
00991 
00992   // Create medium node
00993 
00994   SMDS_MeshNode* n12;
00995   SMESHDS_Mesh* meshDS = GetMeshDS();
00996 
00997   if ( IsSeamShape( n1->getshapeId() ))
00998     // to get a correct UV of a node on seam, the second node must have checked UV
00999     std::swap( n1, n2 );
01000 
01001   // get type of shape for the new medium node
01002   int faceID = -1, edgeID = -1;
01003   TopoDS_Edge E; double u [2];
01004   TopoDS_Face F; gp_XY  uv[2];
01005   bool uvOK[2] = { false, false };
01006 
01007   pair<int, TopAbs_ShapeEnum> pos = GetMediumPos( n1, n2 );
01008 
01009   // get positions of the given nodes on shapes
01010   if ( pos.second == TopAbs_FACE )
01011   {
01012     F = TopoDS::Face(meshDS->IndexToShape( faceID = pos.first ));
01013     uv[0] = GetNodeUV(F,n1,n2, force3d ? 0 : &uvOK[0]);
01014     uv[1] = GetNodeUV(F,n2,n1, force3d ? 0 : &uvOK[1]);
01015   }
01016   else if ( pos.second == TopAbs_EDGE )
01017   {
01018     const SMDS_PositionPtr Pos1 = n1->GetPosition();
01019     const SMDS_PositionPtr Pos2 = n2->GetPosition();
01020     if ( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE &&
01021          Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE &&
01022          n1->getshapeId() != n2->getshapeId() )
01023     {
01024       // issue 0021006
01025       return getMediumNodeOnComposedWire(n1,n2,force3d);
01026     }
01027     E = TopoDS::Edge(meshDS->IndexToShape( edgeID = pos.first ));
01028     u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]);
01029     u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]);
01030   }
01031 
01032   if ( !force3d & uvOK[0] && uvOK[1] )
01033   {
01034     // we try to create medium node using UV parameters of
01035     // nodes, else - medium between corresponding 3d points
01036     if( ! F.IsNull() )
01037     {
01038       //if ( uvOK[0] && uvOK[1] )
01039       {
01040         if ( IsDegenShape( n1->getshapeId() )) {
01041           if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 ));
01042           else                           uv[0].SetCoord( 2, uv[1].Coord( 2 ));
01043         }
01044         else if ( IsDegenShape( n2->getshapeId() )) {
01045           if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 ));
01046           else                           uv[1].SetCoord( 2, uv[0].Coord( 2 ));
01047         }
01048 
01049         TopLoc_Location loc;
01050         Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
01051         gp_XY UV = GetMiddleUV( S, uv[0], uv[1] );
01052         gp_Pnt P = S->Value( UV.X(), UV.Y() ).Transformed(loc);
01053         n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
01054         meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y());
01055         myTLinkNodeMap.insert(make_pair(link,n12));
01056         return n12;
01057       }
01058     }
01059     else if ( !E.IsNull() )
01060     {
01061       double f,l;
01062       Handle(Geom_Curve) C = BRep_Tool::Curve(E, f, l);
01063       if(!C.IsNull())
01064       {
01065         Standard_Boolean isPeriodic = C->IsPeriodic();
01066         double U;
01067         if(isPeriodic) {
01068           Standard_Real Period = C->Period();
01069           Standard_Real p = u[1]+ShapeAnalysis::AdjustByPeriod(u[1],u[0],Period);
01070           Standard_Real pmid = (u[0]+p)/2.;
01071           U = pmid+ShapeAnalysis::AdjustToPeriod(pmid,C->FirstParameter(),C->LastParameter());
01072         }
01073         else
01074           U = (u[0]+u[1])/2.;
01075 
01076         gp_Pnt P = C->Value( U );
01077         n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
01078         meshDS->SetNodeOnEdge(n12, edgeID, U);
01079         myTLinkNodeMap.insert(make_pair(link,n12));
01080         return n12;
01081       }
01082     }
01083   }
01084 
01085   // 3d variant
01086   double x = ( n1->X() + n2->X() )/2.;
01087   double y = ( n1->Y() + n2->Y() )/2.;
01088   double z = ( n1->Z() + n2->Z() )/2.;
01089   n12 = meshDS->AddNode(x,y,z);
01090 
01091   if ( !F.IsNull() )
01092   {
01093     gp_XY UV = ( uv[0] + uv[1] ) / 2.;
01094     CheckNodeUV( F, n12, UV, 2*BRep_Tool::Tolerance( F ), /*force=*/true);
01095     meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y() );
01096   }
01097   else if ( !E.IsNull() )
01098   {
01099     double U = ( u[0] + u[1] ) / 2.;
01100     CheckNodeU( E, n12, U, 2*BRep_Tool::Tolerance( E ), /*force=*/true);
01101     meshDS->SetNodeOnEdge(n12, edgeID, U);
01102   }
01103   else if ( myShapeID > 0 )
01104   {
01105     meshDS->SetNodeInVolume(n12, myShapeID);
01106   }
01107 
01108   myTLinkNodeMap.insert( make_pair( link, n12 ));
01109   return n12;
01110 }
01111 
01112 //================================================================================
01116 //================================================================================
01117 
01118 const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_MeshNode* n1,
01119                                                                      const SMDS_MeshNode* n2,
01120                                                                      bool                 force3d)
01121 {
01122   gp_Pnt middle = 0.5 * XYZ(n1) + 0.5 * XYZ(n2);
01123   SMDS_MeshNode* n12 = AddNode( middle.X(), middle.Y(), middle.Z() );
01124 
01125   // To find position on edge and 3D position for n12,
01126   // project <middle> to 2 edges and select projection most close to <middle>
01127 
01128   double u = 0, distMiddleProj = Precision::Infinite(), distXYZ[4];
01129   int iOkEdge = 0;
01130   TopoDS_Edge edges[2];
01131   for ( int is2nd = 0; is2nd < 2; ++is2nd )
01132   {
01133     // get an edge
01134     const SMDS_MeshNode* n = is2nd ? n2 : n1;
01135     TopoDS_Shape shape = GetSubShapeByNode( n, GetMeshDS() );
01136     if ( shape.IsNull() || shape.ShapeType() != TopAbs_EDGE )
01137       continue;
01138 
01139     // project to get U of projection and distance from middle to projection
01140     TopoDS_Edge edge = edges[ is2nd ] = TopoDS::Edge( shape );
01141     double node2MiddleDist = middle.Distance( XYZ(n) );
01142     double foundU = GetNodeU( edge, n );
01143     CheckNodeU( edge, n12, foundU, 2*BRep_Tool::Tolerance(edge), /*force=*/true, distXYZ );
01144     if ( distXYZ[0] < node2MiddleDist )
01145     {
01146       distMiddleProj = distXYZ[0];
01147       u = foundU;
01148       iOkEdge = is2nd;
01149     }
01150   }
01151   if ( Precision::IsInfinite( distMiddleProj ))
01152   {
01153     // both projections failed; set n12 on the edge of n1 with U of a common vertex
01154     TopoDS_Vertex vCommon;
01155     if ( TopExp::CommonVertex( edges[0], edges[1], vCommon ))
01156       u = BRep_Tool::Parameter( vCommon, edges[0] );
01157     else
01158     {
01159       double f,l, u0 = GetNodeU( edges[0], n1 );
01160       BRep_Tool::Range( edges[0],f,l );
01161       u = ( fabs(u0-f) < fabs(u0-l) ) ? f : l;
01162     }
01163     iOkEdge = 0;
01164     distMiddleProj = 0;
01165   }
01166 
01167   // move n12 to position of a successfull projection
01168   double tol = BRep_Tool::Tolerance(edges[ iOkEdge ]);
01169   if ( !force3d && distMiddleProj > 2*tol )
01170   {
01171     TopLoc_Location loc; double f,l;
01172     Handle(Geom_Curve) curve = BRep_Tool::Curve( edges[iOkEdge],loc,f,l );
01173     gp_Pnt p = curve->Value( u );
01174     GetMeshDS()->MoveNode( n12, p.X(), p.Y(), p.Z() );
01175   }
01176 
01177   GetMeshDS()->SetNodeOnEdge(n12, edges[iOkEdge], u);
01178 
01179   myTLinkNodeMap.insert( make_pair( SMESH_TLink(n1,n2), n12 ));
01180 
01181   return n12;
01182 }
01183 
01184 //=======================================================================
01185 //function : AddNode
01186 //purpose  : Creates a node
01187 //=======================================================================
01188 
01189 SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID)
01190 {
01191   SMESHDS_Mesh * meshDS = GetMeshDS();
01192   SMDS_MeshNode* node = 0;
01193   if ( ID )
01194     node = meshDS->AddNodeWithID( x, y, z, ID );
01195   else
01196     node = meshDS->AddNode( x, y, z );
01197   if ( mySetElemOnShape && myShapeID > 0 ) {
01198     switch ( myShape.ShapeType() ) {
01199     case TopAbs_SOLID:  meshDS->SetNodeInVolume( node, myShapeID); break;
01200     case TopAbs_SHELL:  meshDS->SetNodeInVolume( node, myShapeID); break;
01201     case TopAbs_FACE:   meshDS->SetNodeOnFace(   node, myShapeID); break;
01202     case TopAbs_EDGE:   meshDS->SetNodeOnEdge(   node, myShapeID); break;
01203     case TopAbs_VERTEX: meshDS->SetNodeOnVertex( node, myShapeID); break;
01204     default: ;
01205     }
01206   }
01207   return node;
01208 }
01209 
01210 //=======================================================================
01211 //function : AddEdge
01212 //purpose  : Creates quadratic or linear edge
01213 //=======================================================================
01214 
01215 SMDS_MeshEdge* SMESH_MesherHelper::AddEdge(const SMDS_MeshNode* n1,
01216                                            const SMDS_MeshNode* n2,
01217                                            const int            id,
01218                                            const bool           force3d)
01219 {
01220   SMESHDS_Mesh * meshDS = GetMeshDS();
01221   
01222   SMDS_MeshEdge* edge = 0;
01223   if (myCreateQuadratic) {
01224     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01225     if(id)
01226       edge = meshDS->AddEdgeWithID(n1, n2, n12, id);
01227     else
01228       edge = meshDS->AddEdge(n1, n2, n12);
01229   }
01230   else {
01231     if(id)
01232       edge = meshDS->AddEdgeWithID(n1, n2, id);
01233     else
01234       edge = meshDS->AddEdge(n1, n2);
01235   }
01236 
01237   if ( mySetElemOnShape && myShapeID > 0 )
01238     meshDS->SetMeshElementOnShape( edge, myShapeID );
01239 
01240   return edge;
01241 }
01242 
01243 //=======================================================================
01244 //function : AddFace
01245 //purpose  : Creates quadratic or linear triangle
01246 //=======================================================================
01247 
01248 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
01249                                            const SMDS_MeshNode* n2,
01250                                            const SMDS_MeshNode* n3,
01251                                            const int id,
01252                                            const bool force3d)
01253 {
01254   SMESHDS_Mesh * meshDS = GetMeshDS();
01255   SMDS_MeshFace* elem = 0;
01256 
01257   if( n1==n2 || n2==n3 || n3==n1 )
01258     return elem;
01259 
01260   if(!myCreateQuadratic) {
01261     if(id)
01262       elem = meshDS->AddFaceWithID(n1, n2, n3, id);
01263     else
01264       elem = meshDS->AddFace(n1, n2, n3);
01265   }
01266   else {
01267     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01268     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
01269     const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
01270 
01271     if(id)
01272       elem = meshDS->AddFaceWithID(n1, n2, n3, n12, n23, n31, id);
01273     else
01274       elem = meshDS->AddFace(n1, n2, n3, n12, n23, n31);
01275   }
01276   if ( mySetElemOnShape && myShapeID > 0 )
01277     meshDS->SetMeshElementOnShape( elem, myShapeID );
01278 
01279   return elem;
01280 }
01281 
01282 //=======================================================================
01283 //function : AddFace
01284 //purpose  : Creates quadratic or linear quadrangle
01285 //=======================================================================
01286 
01287 SMDS_MeshFace* SMESH_MesherHelper::AddFace(const SMDS_MeshNode* n1,
01288                                            const SMDS_MeshNode* n2,
01289                                            const SMDS_MeshNode* n3,
01290                                            const SMDS_MeshNode* n4,
01291                                            const int            id,
01292                                            const bool           force3d)
01293 {
01294   SMESHDS_Mesh * meshDS = GetMeshDS();
01295   SMDS_MeshFace* elem = 0;
01296 
01297   if( n1==n2 ) {
01298     return AddFace(n1,n3,n4,id,force3d);
01299   }
01300   if( n1==n3 ) {
01301     return AddFace(n1,n2,n4,id,force3d);
01302   }
01303   if( n1==n4 ) {
01304     return AddFace(n1,n2,n3,id,force3d);
01305   }
01306   if( n2==n3 ) {
01307     return AddFace(n1,n2,n4,id,force3d);
01308   }
01309   if( n2==n4 ) {
01310     return AddFace(n1,n2,n3,id,force3d);
01311   }
01312   if( n3==n4 ) {
01313     return AddFace(n1,n2,n3,id,force3d);
01314   }
01315 
01316   if(!myCreateQuadratic) {
01317     if(id)
01318       elem = meshDS->AddFaceWithID(n1, n2, n3, n4, id);
01319     else
01320       elem = meshDS->AddFace(n1, n2, n3, n4);
01321   }
01322   else {
01323     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01324     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
01325     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
01326     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
01327 
01328     if(id)
01329       elem = meshDS->AddFaceWithID(n1, n2, n3, n4, n12, n23, n34, n41, id);
01330     else
01331       elem = meshDS->AddFace(n1, n2, n3, n4, n12, n23, n34, n41);
01332   }
01333   if ( mySetElemOnShape && myShapeID > 0 )
01334     meshDS->SetMeshElementOnShape( elem, myShapeID );
01335 
01336   return elem;
01337 }
01338 
01339 //=======================================================================
01340 //function : AddPolygonalFace
01341 //purpose  : Creates polygon, with additional nodes in quadratic mesh
01342 //=======================================================================
01343 
01344 SMDS_MeshFace* SMESH_MesherHelper::AddPolygonalFace (const vector<const SMDS_MeshNode*>& nodes,
01345                                                      const int                           id,
01346                                                      const bool                          force3d)
01347 {
01348   SMESHDS_Mesh * meshDS = GetMeshDS();
01349   SMDS_MeshFace* elem = 0;
01350 
01351   if(!myCreateQuadratic) {
01352     if(id)
01353       elem = meshDS->AddPolygonalFaceWithID(nodes, id);
01354     else
01355       elem = meshDS->AddPolygonalFace(nodes);
01356   }
01357   else {
01358     vector<const SMDS_MeshNode*> newNodes;
01359     for ( int i = 0; i < nodes.size(); ++i )
01360     {
01361       const SMDS_MeshNode* n1 = nodes[i];
01362       const SMDS_MeshNode* n2 = nodes[(i+1)%nodes.size()];
01363       const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01364       newNodes.push_back( n1 );
01365       newNodes.push_back( n12 );
01366     }
01367     if(id)
01368       elem = meshDS->AddPolygonalFaceWithID(newNodes, id);
01369     else
01370       elem = meshDS->AddPolygonalFace(newNodes);
01371   }
01372   if ( mySetElemOnShape && myShapeID > 0 )
01373     meshDS->SetMeshElementOnShape( elem, myShapeID );
01374 
01375   return elem;
01376 }
01377 
01378 //=======================================================================
01379 //function : AddVolume
01380 //purpose  : Creates quadratic or linear prism
01381 //=======================================================================
01382 
01383 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
01384                                                const SMDS_MeshNode* n2,
01385                                                const SMDS_MeshNode* n3,
01386                                                const SMDS_MeshNode* n4,
01387                                                const SMDS_MeshNode* n5,
01388                                                const SMDS_MeshNode* n6,
01389                                                const int id,
01390                                                const bool force3d)
01391 {
01392   SMESHDS_Mesh * meshDS = GetMeshDS();
01393   SMDS_MeshVolume* elem = 0;
01394   if(!myCreateQuadratic) {
01395     if(id)
01396       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, id);
01397     else
01398       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6);
01399   }
01400   else {
01401     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01402     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
01403     const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
01404 
01405     const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
01406     const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
01407     const SMDS_MeshNode* n64 = GetMediumNode(n6,n4,force3d);
01408 
01409     const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
01410     const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
01411     const SMDS_MeshNode* n36 = GetMediumNode(n3,n6,force3d);
01412 
01413     if(id)
01414       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, 
01415                                      n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
01416     else
01417       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
01418                                n12, n23, n31, n45, n56, n64, n14, n25, n36);
01419   }
01420   if ( mySetElemOnShape && myShapeID > 0 )
01421     meshDS->SetMeshElementOnShape( elem, myShapeID );
01422 
01423   return elem;
01424 }
01425 
01426 //=======================================================================
01427 //function : AddVolume
01428 //purpose  : Creates quadratic or linear tetrahedron
01429 //=======================================================================
01430 
01431 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
01432                                                const SMDS_MeshNode* n2,
01433                                                const SMDS_MeshNode* n3,
01434                                                const SMDS_MeshNode* n4,
01435                                                const int id, 
01436                                                const bool force3d)
01437 {
01438   SMESHDS_Mesh * meshDS = GetMeshDS();
01439   SMDS_MeshVolume* elem = 0;
01440   if(!myCreateQuadratic) {
01441     if(id)
01442       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, id);
01443     else
01444       elem = meshDS->AddVolume(n1, n2, n3, n4);
01445   }
01446   else {
01447     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01448     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
01449     const SMDS_MeshNode* n31 = GetMediumNode(n3,n1,force3d);
01450 
01451     const SMDS_MeshNode* n14 = GetMediumNode(n1,n4,force3d);
01452     const SMDS_MeshNode* n24 = GetMediumNode(n2,n4,force3d);
01453     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
01454 
01455     if(id)
01456       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34, id);
01457     else
01458       elem = meshDS->AddVolume(n1, n2, n3, n4, n12, n23, n31, n14, n24, n34);
01459   }
01460   if ( mySetElemOnShape && myShapeID > 0 )
01461     meshDS->SetMeshElementOnShape( elem, myShapeID );
01462 
01463   return elem;
01464 }
01465 
01466 //=======================================================================
01467 //function : AddVolume
01468 //purpose  : Creates quadratic or linear pyramid
01469 //=======================================================================
01470 
01471 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
01472                                                const SMDS_MeshNode* n2,
01473                                                const SMDS_MeshNode* n3,
01474                                                const SMDS_MeshNode* n4,
01475                                                const SMDS_MeshNode* n5,
01476                                                const int id, 
01477                                                const bool force3d)
01478 {
01479   SMDS_MeshVolume* elem = 0;
01480   if(!myCreateQuadratic) {
01481     if(id)
01482       elem = GetMeshDS()->AddVolumeWithID(n1, n2, n3, n4, n5, id);
01483     else
01484       elem = GetMeshDS()->AddVolume(n1, n2, n3, n4, n5);
01485   }
01486   else {
01487     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01488     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
01489     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
01490     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
01491 
01492     const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
01493     const SMDS_MeshNode* n25 = GetMediumNode(n2,n5,force3d);
01494     const SMDS_MeshNode* n35 = GetMediumNode(n3,n5,force3d);
01495     const SMDS_MeshNode* n45 = GetMediumNode(n4,n5,force3d);
01496 
01497     if(id)
01498       elem = GetMeshDS()->AddVolumeWithID ( n1,  n2,  n3,  n4,  n5,
01499                                             n12, n23, n34, n41,
01500                                             n15, n25, n35, n45,
01501                                             id);
01502     else
01503       elem = GetMeshDS()->AddVolume( n1,  n2,  n3,  n4,  n5,
01504                                      n12, n23, n34, n41,
01505                                      n15, n25, n35, n45);
01506   }
01507   if ( mySetElemOnShape && myShapeID > 0 )
01508     GetMeshDS()->SetMeshElementOnShape( elem, myShapeID );
01509 
01510   return elem;
01511 }
01512 
01513 //=======================================================================
01514 //function : AddVolume
01515 //purpose  : Creates quadratic or linear hexahedron
01516 //=======================================================================
01517 
01518 SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
01519                                                const SMDS_MeshNode* n2,
01520                                                const SMDS_MeshNode* n3,
01521                                                const SMDS_MeshNode* n4,
01522                                                const SMDS_MeshNode* n5,
01523                                                const SMDS_MeshNode* n6,
01524                                                const SMDS_MeshNode* n7,
01525                                                const SMDS_MeshNode* n8,
01526                                                const int id,
01527                                                const bool force3d)
01528 {
01529   SMESHDS_Mesh * meshDS = GetMeshDS();
01530   SMDS_MeshVolume* elem = 0;
01531   if(!myCreateQuadratic) {
01532     if(id)
01533       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8, id);
01534     else
01535       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8);
01536   }
01537   else {
01538     const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01539     const SMDS_MeshNode* n23 = GetMediumNode(n2,n3,force3d);
01540     const SMDS_MeshNode* n34 = GetMediumNode(n3,n4,force3d);
01541     const SMDS_MeshNode* n41 = GetMediumNode(n4,n1,force3d);
01542 
01543     const SMDS_MeshNode* n56 = GetMediumNode(n5,n6,force3d);
01544     const SMDS_MeshNode* n67 = GetMediumNode(n6,n7,force3d);
01545     const SMDS_MeshNode* n78 = GetMediumNode(n7,n8,force3d);
01546     const SMDS_MeshNode* n85 = GetMediumNode(n8,n5,force3d);
01547 
01548     const SMDS_MeshNode* n15 = GetMediumNode(n1,n5,force3d);
01549     const SMDS_MeshNode* n26 = GetMediumNode(n2,n6,force3d);
01550     const SMDS_MeshNode* n37 = GetMediumNode(n3,n7,force3d);
01551     const SMDS_MeshNode* n48 = GetMediumNode(n4,n8,force3d);
01552 
01553     if(id)
01554       elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8,
01555                                      n12, n23, n34, n41, n56, n67,
01556                                      n78, n85, n15, n26, n37, n48, id);
01557     else
01558       elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, n7, n8,
01559                                n12, n23, n34, n41, n56, n67,
01560                                n78, n85, n15, n26, n37, n48);
01561   }
01562   if ( mySetElemOnShape && myShapeID > 0 )
01563     meshDS->SetMeshElementOnShape( elem, myShapeID );
01564 
01565   return elem;
01566 }
01567 
01568 //=======================================================================
01569 //function : AddPolyhedralVolume
01570 //purpose  : Creates polyhedron. In quadratic mesh, adds medium nodes
01571 //=======================================================================
01572 
01573 SMDS_MeshVolume*
01574 SMESH_MesherHelper::AddPolyhedralVolume (const std::vector<const SMDS_MeshNode*>& nodes,
01575                                          const std::vector<int>&                  quantities,
01576                                          const int                                id,
01577                                          const bool                               force3d)
01578 {
01579   SMESHDS_Mesh * meshDS = GetMeshDS();
01580   SMDS_MeshVolume* elem = 0;
01581   if(!myCreateQuadratic)
01582   {
01583     if(id)
01584       elem = meshDS->AddPolyhedralVolumeWithID(nodes, quantities, id);
01585     else
01586       elem = meshDS->AddPolyhedralVolume(nodes, quantities);
01587   }
01588   else
01589   {
01590     vector<const SMDS_MeshNode*> newNodes;
01591     vector<int> newQuantities;
01592     for ( int iFace=0, iN=0; iFace < quantities.size(); ++iFace)
01593     {
01594       int nbNodesInFace = quantities[iFace];
01595       newQuantities.push_back(0);
01596       for ( int i = 0; i < nbNodesInFace; ++i )
01597       {
01598         const SMDS_MeshNode* n1 = nodes[ iN + i ];
01599         newNodes.push_back( n1 );
01600         newQuantities.back()++;
01601         
01602         const SMDS_MeshNode* n2 = nodes[ iN + ( i+1==nbNodesInFace ? 0 : i+1 )];
01603 //         if ( n1->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE &&
01604 //              n2->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE )
01605         {
01606           const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
01607           newNodes.push_back( n12 );
01608           newQuantities.back()++;
01609         }
01610       }
01611       iN += nbNodesInFace;
01612     }
01613     if(id)
01614       elem = meshDS->AddPolyhedralVolumeWithID( newNodes, newQuantities, id );
01615     else
01616       elem = meshDS->AddPolyhedralVolume( newNodes, newQuantities );
01617   }
01618   if ( mySetElemOnShape && myShapeID > 0 )
01619     meshDS->SetMeshElementOnShape( elem, myShapeID );
01620 
01621   return elem;
01622 }
01623 
01624 //=======================================================================
01625 //function : LoadNodeColumns
01626 //purpose  : Load nodes bound to face into a map of node columns
01627 //=======================================================================
01628 
01629 bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
01630                                          const TopoDS_Face& theFace,
01631                                          const TopoDS_Edge& theBaseEdge,
01632                                          SMESHDS_Mesh*      theMesh,
01633                                          SMESH_ProxyMesh*   theProxyMesh)
01634 {
01635   return LoadNodeColumns(theParam2ColumnMap,
01636                          theFace,
01637                          std::list<TopoDS_Edge>(1,theBaseEdge),
01638                          theMesh,
01639                          theProxyMesh);
01640 }
01641 
01642 //=======================================================================
01643 //function : LoadNodeColumns
01644 //purpose  : Load nodes bound to face into a map of node columns
01645 //=======================================================================
01646 
01647 bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2ColumnMap,
01648                                          const TopoDS_Face&            theFace,
01649                                          const std::list<TopoDS_Edge>& theBaseSide,
01650                                          SMESHDS_Mesh*                 theMesh,
01651                                          SMESH_ProxyMesh*              theProxyMesh)
01652 {
01653   // get a right submesh of theFace
01654 
01655   const SMESHDS_SubMesh* faceSubMesh = 0;
01656   if ( theProxyMesh )
01657   {
01658     faceSubMesh = theProxyMesh->GetSubMesh( theFace );
01659     if ( !faceSubMesh ||
01660          faceSubMesh->NbElements() == 0 ||
01661          theProxyMesh->IsTemporary( faceSubMesh->GetElements()->next() ))
01662     {
01663       // can use a proxy sub-mesh with not temporary elements only
01664       faceSubMesh = 0;
01665       theProxyMesh = 0;
01666     }
01667   }
01668   if ( !faceSubMesh )
01669     faceSubMesh = theMesh->MeshElements( theFace );
01670   if ( !faceSubMesh || faceSubMesh->NbElements() == 0 )
01671     return false;
01672 
01673   // get data of edges for normalization of params
01674 
01675   vector< double > length;
01676   double fullLen = 0;
01677   list<TopoDS_Edge>::const_iterator edge;
01678   {
01679     for ( edge = theBaseSide.begin(); edge != theBaseSide.end(); ++edge )
01680     {
01681       double len = std::max( 1e-10, SMESH_Algo::EdgeLength( *edge ));
01682       fullLen += len;
01683       length.push_back( len );
01684     }
01685   }
01686 
01687   // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them
01688   edge = theBaseSide.begin();
01689   for ( int iE = 0; edge != theBaseSide.end(); ++edge, ++iE )
01690   {
01691     map< double, const SMDS_MeshNode*> sortedBaseNodes;
01692     SMESH_Algo::GetSortedNodesOnEdge( theMesh, *edge,/*noMedium=*/true, sortedBaseNodes);
01693     if ( sortedBaseNodes.empty() ) continue;
01694 
01695     double f, l;
01696     BRep_Tool::Range( *edge, f, l );
01697     if ( edge->Orientation() == TopAbs_REVERSED ) std::swap( f, l );
01698     const double coeff = 1. / ( l - f ) * length[iE] / fullLen;
01699     const double prevPar = theParam2ColumnMap.empty() ? 0 : theParam2ColumnMap.rbegin()->first;
01700     map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNodes.begin();
01701     for ( ; u_n != sortedBaseNodes.end(); u_n++ )
01702     {
01703       double par = prevPar + coeff * ( u_n->first - f );
01704       TParam2ColumnMap::iterator u2nn =
01705         theParam2ColumnMap.insert( theParam2ColumnMap.end(), make_pair( par, TNodeColumn()));
01706       u2nn->second.push_back( u_n->second );
01707     }
01708   }
01709   TParam2ColumnMap::iterator par_nVec_2, par_nVec_1 = theParam2ColumnMap.begin();
01710   if ( theProxyMesh )
01711   {
01712     for ( ; par_nVec_1 != theParam2ColumnMap.end(); ++par_nVec_1 )
01713     {
01714       const SMDS_MeshNode* & n = par_nVec_1->second[0];
01715       n = theProxyMesh->GetProxyNode( n );
01716     }
01717   }
01718 
01719   int nbRows = 1 + faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 );
01720 
01721   // fill theParam2ColumnMap column by column by passing from nodes on
01722   // theBaseEdge up via mesh faces on theFace
01723 
01724   par_nVec_2 = theParam2ColumnMap.begin();
01725   par_nVec_1 = par_nVec_2++;
01726   TIDSortedElemSet emptySet, avoidSet;
01727   for ( ; par_nVec_2 != theParam2ColumnMap.end(); ++par_nVec_1, ++par_nVec_2 )
01728   {
01729     vector<const SMDS_MeshNode*>& nCol1 = par_nVec_1->second;
01730     vector<const SMDS_MeshNode*>& nCol2 = par_nVec_2->second;
01731     nCol1.resize( nbRows );
01732     nCol2.resize( nbRows );
01733 
01734     int i1, i2, iRow = 0;
01735     const SMDS_MeshNode *n1 = nCol1[0], *n2 = nCol2[0];
01736     // find face sharing node n1 and n2 and belonging to faceSubMesh
01737     while ( const SMDS_MeshElement* face =
01738             SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2))
01739     {
01740       if ( faceSubMesh->Contains( face ))
01741       {
01742         int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes();
01743         if ( nbNodes != 4 )
01744           return false;
01745         n1 = face->GetNode( (i2+2) % 4 ); // opposite corner of quadrangle face
01746         n2 = face->GetNode( (i1+2) % 4 );
01747         if ( ++iRow >= nbRows )
01748           return false;
01749         nCol1[ iRow ] = n1;
01750         nCol2[ iRow ] = n2;
01751         avoidSet.clear();
01752       }
01753       avoidSet.insert( face );
01754     }
01755     // set a real height
01756     nCol1.resize( iRow + 1 );
01757     nCol2.resize( iRow + 1 );
01758   }
01759   return theParam2ColumnMap.size() > 1 && theParam2ColumnMap.begin()->second.size() > 1;
01760 }
01761 
01762 //=======================================================================
01763 //function : NbAncestors
01764 //purpose  : Return number of unique ancestors of the shape
01765 //=======================================================================
01766 
01767 int SMESH_MesherHelper::NbAncestors(const TopoDS_Shape& shape,
01768                                     const SMESH_Mesh&   mesh,
01769                                     TopAbs_ShapeEnum    ancestorType/*=TopAbs_SHAPE*/)
01770 {
01771   TopTools_MapOfShape ancestors;
01772   TopTools_ListIteratorOfListOfShape ansIt( mesh.GetAncestors(shape) );
01773   for ( ; ansIt.More(); ansIt.Next() ) {
01774     if ( ancestorType == TopAbs_SHAPE || ansIt.Value().ShapeType() == ancestorType )
01775       ancestors.Add( ansIt.Value() );
01776   }
01777   return ancestors.Extent();
01778 }
01779 
01780 //=======================================================================
01781 //function : GetSubShapeOri
01782 //purpose  : Return orientation of sub-shape in the main shape
01783 //=======================================================================
01784 
01785 TopAbs_Orientation SMESH_MesherHelper::GetSubShapeOri(const TopoDS_Shape& shape,
01786                                                       const TopoDS_Shape& subShape)
01787 {
01788   TopAbs_Orientation ori = TopAbs_Orientation(-1);
01789   if ( !shape.IsNull() && !subShape.IsNull() )
01790   {
01791     TopExp_Explorer e( shape, subShape.ShapeType() );
01792     if ( shape.Orientation() >= TopAbs_INTERNAL ) // TopAbs_INTERNAL or TopAbs_EXTERNAL
01793       e.Init( shape.Oriented(TopAbs_FORWARD), subShape.ShapeType() );
01794     for ( ; e.More(); e.Next())
01795       if ( subShape.IsSame( e.Current() ))
01796         break;
01797     if ( e.More() )
01798       ori = e.Current().Orientation();
01799   }
01800   return ori;
01801 }
01802 
01803 //=======================================================================
01804 //function : IsSubShape
01805 //purpose  : 
01806 //=======================================================================
01807 
01808 bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape,
01809                                      const TopoDS_Shape& mainShape )
01810 {
01811   if ( !shape.IsNull() && !mainShape.IsNull() )
01812   {
01813     for ( TopExp_Explorer exp( mainShape, shape.ShapeType());
01814           exp.More();
01815           exp.Next() )
01816       if ( shape.IsSame( exp.Current() ))
01817         return true;
01818   }
01819   SCRUTE((shape.IsNull()));
01820   SCRUTE((mainShape.IsNull()));
01821   return false;
01822 }
01823 
01824 //=======================================================================
01825 //function : IsSubShape
01826 //purpose  : 
01827 //=======================================================================
01828 
01829 bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, SMESH_Mesh* aMesh )
01830 {
01831   if ( shape.IsNull() || !aMesh )
01832     return false;
01833   return
01834     aMesh->GetMeshDS()->ShapeToIndex( shape ) ||
01835     // PAL16202
01836     (shape.ShapeType() == TopAbs_COMPOUND && aMesh->GetMeshDS()->IsGroupOfSubShapes( shape ));
01837 }
01838 
01839 //================================================================================
01843 //================================================================================
01844 
01845 double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
01846 {
01847   double tol = Precision::Confusion();
01848   TopExp_Explorer exp;
01849   for ( exp.Init( shape, TopAbs_FACE ); exp.More(); exp.Next() )
01850     tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Face( exp.Current())));
01851   for ( exp.Init( shape, TopAbs_EDGE ); exp.More(); exp.Next() )
01852     tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Edge( exp.Current())));
01853   for ( exp.Init( shape, TopAbs_VERTEX ); exp.More(); exp.Next() )
01854     tol = Max( tol, BRep_Tool::Tolerance( TopoDS::Vertex( exp.Current())));
01855 
01856   return tol;
01857 }
01858 
01859 //================================================================================
01865 //================================================================================
01866 
01867 bool SMESH_MesherHelper::IsClosedEdge( const TopoDS_Edge& anEdge )
01868 {
01869   if ( anEdge.Orientation() >= TopAbs_INTERNAL )
01870     return IsClosedEdge( TopoDS::Edge( anEdge.Oriented( TopAbs_FORWARD )));
01871   return TopExp::FirstVertex( anEdge ).IsSame( TopExp::LastVertex( anEdge ));
01872 }
01873 
01874 //================================================================================
01879 //================================================================================
01880 
01881 TopoDS_Vertex SMESH_MesherHelper::IthVertex( const bool  is2nd,
01882                                              TopoDS_Edge anEdge,
01883                                              const bool  CumOri )
01884 {
01885   if ( anEdge.Orientation() >= TopAbs_INTERNAL )
01886     anEdge.Orientation( TopAbs_FORWARD );
01887 
01888   const TopAbs_Orientation tgtOri = is2nd ? TopAbs_REVERSED : TopAbs_FORWARD;
01889   TopoDS_Iterator vIt( anEdge, CumOri );
01890   while ( vIt.More() && vIt.Value().Orientation() != tgtOri )
01891     vIt.Next();
01892 
01893   return ( vIt.More() ? TopoDS::Vertex(vIt.Value()) : TopoDS_Vertex() );
01894 }
01895 
01896 //=======================================================================
01897 //function : IsQuadraticMesh
01898 //purpose  : Check mesh without geometry for: if all elements on this shape are quadratic,
01899 //           quadratic elements will be created.
01900 //           Used then generated 3D mesh without geometry.
01901 //=======================================================================
01902 
01903 SMESH_MesherHelper:: MType SMESH_MesherHelper::IsQuadraticMesh()
01904 {
01905   int NbAllEdgsAndFaces=0;
01906   int NbQuadFacesAndEdgs=0;
01907   int NbFacesAndEdges=0;
01908   //All faces and edges
01909   NbAllEdgsAndFaces = myMesh->NbEdges() + myMesh->NbFaces();
01910   
01911   //Quadratic faces and edges
01912   NbQuadFacesAndEdgs = myMesh->NbEdges(ORDER_QUADRATIC) + myMesh->NbFaces(ORDER_QUADRATIC);
01913 
01914   //Linear faces and edges
01915   NbFacesAndEdges = myMesh->NbEdges(ORDER_LINEAR) + myMesh->NbFaces(ORDER_LINEAR);
01916   
01917   if (NbAllEdgsAndFaces == NbQuadFacesAndEdgs) {
01918     //Quadratic mesh
01919     return SMESH_MesherHelper::QUADRATIC;
01920   }
01921   else if (NbAllEdgsAndFaces == NbFacesAndEdges) {
01922     //Linear mesh
01923     return SMESH_MesherHelper::LINEAR;
01924   }
01925   else
01926     //Mesh with both type of elements
01927     return SMESH_MesherHelper::COMP;
01928 }
01929 
01930 //=======================================================================
01931 //function : GetOtherParam
01932 //purpose  : Return an alternative parameter for a node on seam
01933 //=======================================================================
01934 
01935 double SMESH_MesherHelper::GetOtherParam(const double param) const
01936 {
01937   int i = myParIndex & U_periodic ? 0 : 1;
01938   return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i];
01939 }
01940 
01941 namespace {
01942 
01943   //=======================================================================
01947   //=======================================================================
01948 
01949   struct TAncestorsIterator : public SMDS_Iterator<const TopoDS_Shape*>
01950   {
01951     TopTools_ListIteratorOfListOfShape _ancIter;
01952     TopAbs_ShapeEnum                   _type;
01953     TopTools_MapOfShape                _encountered;
01954     TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type)
01955       : _ancIter( ancestors ), _type( type )
01956     {
01957       if ( _ancIter.More() ) {
01958         if ( _ancIter.Value().ShapeType() != _type ) next();
01959         else _encountered.Add( _ancIter.Value() );
01960       }
01961     }
01962     virtual bool more()
01963     {
01964       return _ancIter.More();
01965     }
01966     virtual const TopoDS_Shape* next()
01967     {
01968       const TopoDS_Shape* s = _ancIter.More() ? & _ancIter.Value() : 0;
01969       if ( _ancIter.More() )
01970         for ( _ancIter.Next();  _ancIter.More(); _ancIter.Next())
01971           if ( _ancIter.Value().ShapeType() == _type && _encountered.Add( _ancIter.Value() ))
01972             break;
01973       return s;
01974     }
01975   };
01976 
01977 } // namespace
01978 
01979 //=======================================================================
01983 //=======================================================================
01984 
01985 PShapeIteratorPtr SMESH_MesherHelper::GetAncestors(const TopoDS_Shape& shape,
01986                                                    const SMESH_Mesh&   mesh,
01987                                                    TopAbs_ShapeEnum    ancestorType)
01988 {
01989   return PShapeIteratorPtr( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType));
01990 }
01991 
01992 //=======================================================================
01993 //function : GetCommonAncestor
01994 //purpose  : Find a common ancestors of two shapes of the given type
01995 //=======================================================================
01996 
01997 TopoDS_Shape SMESH_MesherHelper::GetCommonAncestor(const TopoDS_Shape& shape1,
01998                                                    const TopoDS_Shape& shape2,
01999                                                    const SMESH_Mesh&   mesh,
02000                                                    TopAbs_ShapeEnum    ancestorType)
02001 {
02002   TopoDS_Shape commonAnc;
02003   if ( !shape1.IsNull() && !shape2.IsNull() )
02004   {
02005     PShapeIteratorPtr ancIt = GetAncestors( shape1, mesh, ancestorType );
02006     while ( const TopoDS_Shape* anc = ancIt->next() )
02007       if ( IsSubShape( shape2, *anc ))
02008       {
02009         commonAnc = *anc;
02010         break;
02011       }
02012   }
02013   return commonAnc;
02014 }
02015 
02016 //#include <Perf_Meter.hxx>
02017 
02018 //=======================================================================
02019 namespace { // Structures used by FixQuadraticElements()
02020 //=======================================================================
02021 
02022 #define __DMP__(txt) \
02023   //cout << txt
02024 #define MSG(txt) __DMP__(txt<<endl)
02025 #define MSGBEG(txt) __DMP__(txt)
02026 
02027   //const double straightTol2 = 1e-33; // to detect straing links
02028   bool isStraightLink(double linkLen2, double middleNodeMove2)
02029   {
02030     // straight if <node move> < 1/15 * <link length>
02031     return middleNodeMove2 < 1/15./15. * linkLen2;
02032   }
02033 
02034   struct QFace;
02035   // ---------------------------------------
02039   struct QLink: public SMESH_TLink
02040   {
02041     const SMDS_MeshNode*          _mediumNode;
02042     mutable vector<const QFace* > _faces;
02043     mutable gp_Vec                _nodeMove;
02044     mutable int                   _nbMoves;
02045 
02046     QLink(const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshNode* nm):
02047       SMESH_TLink( n1,n2 ), _mediumNode(nm), _nodeMove(0,0,0), _nbMoves(0) {
02048       _faces.reserve(4);
02049       //if ( MediumPos() != SMDS_TOP_3DSPACE )
02050         _nodeMove = MediumPnt() - MiddlePnt();
02051     }
02052     void SetContinuesFaces() const;
02053     const QFace* GetContinuesFace( const QFace* face ) const;
02054     bool OnBoundary() const;
02055     gp_XYZ MiddlePnt() const { return ( XYZ( node1() ) + XYZ( node2() )) / 2.; }
02056     gp_XYZ MediumPnt() const { return XYZ( _mediumNode ); }
02057 
02058     SMDS_TypeOfPosition MediumPos() const
02059     { return _mediumNode->GetPosition()->GetTypeOfPosition(); }
02060     SMDS_TypeOfPosition EndPos(bool isSecond) const
02061     { return (isSecond ? node2() : node1())->GetPosition()->GetTypeOfPosition(); }
02062     const SMDS_MeshNode* EndPosNode(SMDS_TypeOfPosition pos) const
02063     { return EndPos(0) == pos ? node1() : EndPos(1) == pos ? node2() : 0; }
02064 
02065     void Move(const gp_Vec& move, bool sum=false) const
02066     { _nodeMove += move; _nbMoves += sum ? (_nbMoves==0) : 1; }
02067     gp_XYZ Move() const { return _nodeMove.XYZ() / _nbMoves; }
02068     bool IsMoved() const { return (_nbMoves > 0 /*&& !IsStraight()*/); }
02069     bool IsStraight() const
02070     { return isStraightLink( (XYZ(node1())-XYZ(node2())).SquareModulus(),
02071                              _nodeMove.SquareMagnitude());
02072     }
02073     bool operator<(const QLink& other) const {
02074       return (node1()->GetID() == other.node1()->GetID() ?
02075               node2()->GetID() < other.node2()->GetID() :
02076               node1()->GetID() < other.node1()->GetID());
02077     }
02078 //     struct PtrComparator {
02079 //       bool operator() (const QLink* l1, const QLink* l2 ) const { return *l1 < *l2; }
02080 //     };
02081   };
02082   // ---------------------------------------------------------
02086   struct TChainLink
02087   {
02088     const QLink*         _qlink;
02089     mutable const QFace* _qfaces[2];
02090 
02091     TChainLink(const QLink* qlink=0):_qlink(qlink) {
02092       _qfaces[0] = _qfaces[1] = 0;
02093     }
02094     void SetFace(const QFace* face) const { int iF = _qfaces[0] ? 1 : 0; _qfaces[iF]=face; }
02095 
02096     bool IsBoundary() const { return !_qfaces[1]; }
02097 
02098     void RemoveFace( const QFace* face ) const
02099     { _qfaces[(face == _qfaces[1])] = 0; if (!_qfaces[0]) std::swap(_qfaces[0],_qfaces[1]); }
02100 
02101     const QFace* NextFace( const QFace* f ) const
02102     { return _qfaces[0]==f ? _qfaces[1] : _qfaces[0]; }
02103 
02104     const SMDS_MeshNode* NextNode( const SMDS_MeshNode* n ) const
02105     { return n == _qlink->node1() ? _qlink->node2() : _qlink->node1(); }
02106 
02107     bool operator<(const TChainLink& other) const { return *_qlink < *other._qlink; }
02108 
02109     operator bool() const { return (_qlink); }
02110 
02111     const QLink* operator->() const { return _qlink; }
02112 
02113     gp_Vec Normal() const;
02114 
02115     bool IsStraight() const;
02116   };
02117   // --------------------------------------------------------------------
02118   typedef list< TChainLink > TChain;
02119   typedef set < TChainLink > TLinkSet;
02120   typedef TLinkSet::const_iterator TLinkInSet;
02121 
02122   const int theFirstStep = 5;
02123 
02124   enum { ERR_OK, ERR_TRI, ERR_PRISM, ERR_UNKNOWN }; // errors of QFace::GetLinkChain()
02125   // --------------------------------------------------------------------
02129   struct QFace: public TIDSortedNodeSet
02130   {
02131     mutable const SMDS_MeshElement* _volumes[2];
02132     mutable vector< const QLink* >  _sides;
02133     mutable bool                    _sideIsAdded[4]; // added in chain of links
02134     gp_Vec                          _normal;
02135 #ifdef _DEBUG_
02136     mutable const SMDS_MeshElement* _face;
02137 #endif
02138 
02139     QFace( const vector< const QLink*>& links, const SMDS_MeshElement* face=0 );
02140 
02141     void SetVolume(const SMDS_MeshElement* v) const { _volumes[ _volumes[0] ? 1 : 0 ] = v; }
02142 
02143     int NbVolumes() const { return !_volumes[0] ? 0 : !_volumes[1] ? 1 : 2; }
02144 
02145     void AddSelfToLinks() const {
02146       for ( int i = 0; i < _sides.size(); ++i )
02147         _sides[i]->_faces.push_back( this );
02148     }
02149     int LinkIndex( const QLink* side ) const {
02150       for (int i=0; i<_sides.size(); ++i ) if ( _sides[i] == side ) return i;
02151       return -1;
02152     }
02153     bool GetLinkChain( int iSide, TChain& chain, SMDS_TypeOfPosition pos, int& err) const;
02154 
02155     bool GetLinkChain( TChainLink& link, TChain& chain, SMDS_TypeOfPosition pos, int& err) const
02156     {
02157       int i = LinkIndex( link._qlink );
02158       if ( i < 0 ) return true;
02159       _sideIsAdded[i] = true;
02160       link.SetFace( this );
02161       // continue from opposite link
02162       return GetLinkChain( (i+2)%_sides.size(), chain, pos, err );
02163     }
02164     bool IsBoundary() const { return !_volumes[1]; }
02165 
02166     bool Contains( const SMDS_MeshNode* node ) const { return count(node); }
02167 
02168     bool IsSpoiled(const QLink* bentLink ) const;
02169 
02170     TLinkInSet GetBoundaryLink( const TLinkSet&      links,
02171                                 const TChainLink&    avoidLink,
02172                                 TLinkInSet *         notBoundaryLink = 0,
02173                                 const SMDS_MeshNode* nodeToContain = 0,
02174                                 bool *               isAdjacentUsed = 0,
02175                                 int                  nbRecursionsLeft = -1) const;
02176 
02177     TLinkInSet GetLinkByNode( const TLinkSet&      links,
02178                               const TChainLink&    avoidLink,
02179                               const SMDS_MeshNode* nodeToContain) const;
02180 
02181     const SMDS_MeshNode* GetNodeInFace() const {
02182       for ( int iL = 0; iL < _sides.size(); ++iL )
02183         if ( _sides[iL]->MediumPos() == SMDS_TOP_FACE ) return _sides[iL]->_mediumNode;
02184       return 0;
02185     }
02186 
02187     gp_Vec LinkNorm(const int i, SMESH_MesherHelper* theFaceHelper=0) const;
02188 
02189     double MoveByBoundary( const TChainLink&   theLink,
02190                            const gp_Vec&       theRefVec,
02191                            const TLinkSet&     theLinks,
02192                            SMESH_MesherHelper* theFaceHelper=0,
02193                            const double        thePrevLen=0,
02194                            const int           theStep=theFirstStep,
02195                            gp_Vec*             theLinkNorm=0,
02196                            double              theSign=1.0) const;
02197   };
02198 
02199   //================================================================================
02203   ostream& operator << (ostream& out, const QLink& l)
02204   {
02205     out <<"QLink nodes: "
02206         << l.node1()->GetID() << " - "
02207         << l._mediumNode->GetID() << " - "
02208         << l.node2()->GetID() << endl;
02209     return out;
02210   }
02211   ostream& operator << (ostream& out, const QFace& f)
02212   {
02213     out <<"QFace nodes: "/*<< &f << "  "*/;
02214     for ( TIDSortedNodeSet::const_iterator n = f.begin(); n != f.end(); ++n )
02215       out << (*n)->GetID() << " ";
02216     out << " \tvolumes: "
02217         << (f._volumes[0] ? f._volumes[0]->GetID() : 0) << " "
02218         << (f._volumes[1] ? f._volumes[1]->GetID() : 0);
02219     out << "  \tNormal: "<< f._normal.X() <<", "<<f._normal.Y() <<", "<<f._normal.Z() << endl;
02220     return out;
02221   }
02222 
02223   //================================================================================
02227   //================================================================================
02228 
02229   QFace::QFace( const vector< const QLink*>& links, const SMDS_MeshElement* face )
02230   {
02231     _volumes[0] = _volumes[1] = 0;
02232     _sides = links;
02233     _sideIsAdded[0]=_sideIsAdded[1]=_sideIsAdded[2]=_sideIsAdded[3]=false;
02234     _normal.SetCoord(0,0,0);
02235     for ( int i = 1; i < _sides.size(); ++i ) {
02236       const QLink *l1 = _sides[i-1], *l2 = _sides[i];
02237       insert( l1->node1() ); insert( l1->node2() );
02238       // compute normal
02239       gp_Vec v1( XYZ( l1->node2()), XYZ( l1->node1()));
02240       gp_Vec v2( XYZ( l2->node1()), XYZ( l2->node2()));
02241       if ( l1->node1() != l2->node1() && l1->node2() != l2->node2() )
02242         v1.Reverse(); 
02243       _normal += v1 ^ v2;
02244     }
02245     double normSqSize = _normal.SquareMagnitude();
02246     if ( normSqSize > numeric_limits<double>::min() )
02247       _normal /= sqrt( normSqSize );
02248     else
02249       _normal.SetCoord(1e-33,0,0);
02250 
02251 #ifdef _DEBUG_
02252     _face = face;
02253 #endif
02254   }
02255   //================================================================================
02265   //================================================================================
02266 
02267   bool QFace::GetLinkChain( int iSide, TChain& chain, SMDS_TypeOfPosition pos, int& error) const
02268   {
02269     if ( iSide >= _sides.size() ) // wrong argument iSide
02270       return false;
02271     if ( _sideIsAdded[ iSide ]) // already in chain
02272       return true;
02273 
02274     if ( _sides.size() != 4 ) { // triangle - visit all my continous faces
02275       MSGBEG( *this );
02276       TLinkSet links;
02277       list< const QFace* > faces( 1, this );
02278       while ( !faces.empty() ) {
02279         const QFace* face = faces.front();
02280         for ( int i = 0; i < face->_sides.size(); ++i ) {
02281           if ( !face->_sideIsAdded[i] && face->_sides[i] ) {
02282             face->_sideIsAdded[i] = true;
02283             // find a face side in the chain
02284             TLinkInSet chLink = links.insert( TChainLink(face->_sides[i])).first;
02285 //             TChain::iterator chLink = chain.begin();
02286 //             for ( ; chLink != chain.end(); ++chLink )
02287 //               if ( chLink->_qlink == face->_sides[i] )
02288 //                 break;
02289 //             if ( chLink == chain.end() )
02290 //               chLink = chain.insert( chain.begin(), TChainLink(face->_sides[i]));
02291             // add a face to a chained link and put a continues face in the queue
02292             chLink->SetFace( face );
02293             if ( face->_sides[i]->MediumPos() == pos )
02294               if ( const QFace* contFace = face->_sides[i]->GetContinuesFace( face ))
02295                 if ( contFace->_sides.size() == 3 )
02296                   faces.push_back( contFace );
02297           }
02298         }
02299         faces.pop_front();
02300       }
02301       if ( error < ERR_TRI )
02302         error = ERR_TRI;
02303       chain.insert( chain.end(), links.begin(),links.end() );
02304       return false;
02305     }
02306     _sideIsAdded[iSide] = true; // not to add this link to chain again
02307     const QLink* link = _sides[iSide];
02308     if ( !link)
02309       return true;
02310 
02311     // add link into chain
02312     TChain::iterator chLink = chain.insert( chain.begin(), TChainLink(link));
02313     chLink->SetFace( this );
02314     MSGBEG( *this );
02315 
02316     // propagate from quadrangle to neighbour faces
02317     if ( link->MediumPos() >= pos ) {
02318       int nbLinkFaces = link->_faces.size();
02319       if ( nbLinkFaces == 4 || (/*nbLinkFaces < 4 && */link->OnBoundary())) {
02320         // hexahedral mesh or boundary quadrangles - goto a continous face
02321         if ( const QFace* f = link->GetContinuesFace( this ))
02322           if ( f->_sides.size() == 4 )
02323             return f->GetLinkChain( *chLink, chain, pos, error );
02324       }
02325       else {
02326         TChainLink chLink(link); // side face of prismatic mesh - visit all faces of iSide
02327         for ( int i = 0; i < nbLinkFaces; ++i )
02328           if ( link->_faces[i] )
02329             link->_faces[i]->GetLinkChain( chLink, chain, pos, error );
02330         if ( error < ERR_PRISM )
02331           error = ERR_PRISM;
02332         return false;
02333       }
02334     }
02335     return true;
02336   }
02337 
02338   //================================================================================
02349   //================================================================================
02350 
02351   TLinkInSet QFace::GetBoundaryLink( const TLinkSet&      links,
02352                                      const TChainLink&    avoidLink,
02353                                      TLinkInSet *         notBoundaryLink,
02354                                      const SMDS_MeshNode* nodeToContain,
02355                                      bool *               isAdjacentUsed,
02356                                      int                  nbRecursionsLeft) const
02357   {
02358     TLinkInSet linksEnd = links.end(), boundaryLink = linksEnd;
02359 
02360     typedef list< pair< const QFace*, TLinkInSet > > TFaceLinkList;
02361     TFaceLinkList adjacentFaces;
02362 
02363     for ( int iL = 0; iL < _sides.size(); ++iL )
02364     {
02365       if ( avoidLink._qlink == _sides[iL] )
02366         continue;
02367       TLinkInSet link = links.find( _sides[iL] );
02368       if ( link == linksEnd ) continue;
02369       if ( (*link)->MediumPos() > SMDS_TOP_FACE )
02370         continue; // We work on faces here, don't go inside a solid
02371 
02372       // check link
02373       if ( link->IsBoundary() ) {
02374         if ( !nodeToContain ||
02375              (*link)->node1() == nodeToContain ||
02376              (*link)->node2() == nodeToContain )
02377         {
02378           boundaryLink = link;
02379           if ( !notBoundaryLink ) break;
02380         }
02381       }
02382       else if ( notBoundaryLink ) {
02383         *notBoundaryLink = link;
02384         if ( boundaryLink != linksEnd ) break;
02385       }
02386 
02387       if ( boundaryLink == linksEnd && nodeToContain ) // collect adjacent faces
02388         if ( const QFace* adj = link->NextFace( this ))
02389           if ( adj->Contains( nodeToContain ))
02390             adjacentFaces.push_back( make_pair( adj, link ));
02391     }
02392 
02393     if ( isAdjacentUsed ) *isAdjacentUsed = false;
02394     if ( boundaryLink == linksEnd && nodeToContain && nbRecursionsLeft) // check adjacent faces
02395     {
02396       if ( nbRecursionsLeft < 0 )
02397         nbRecursionsLeft = nodeToContain->NbInverseElements();
02398       TFaceLinkList::iterator adj = adjacentFaces.begin();
02399       for ( ; boundaryLink == linksEnd && adj != adjacentFaces.end(); ++adj )
02400         boundaryLink = adj->first->GetBoundaryLink( links, *(adj->second), 0, nodeToContain,
02401                                                     isAdjacentUsed, nbRecursionsLeft-1);
02402       if ( isAdjacentUsed ) *isAdjacentUsed = true;
02403     }
02404     return boundaryLink;
02405   }
02406   //================================================================================
02410   //================================================================================
02411 
02412   TLinkInSet QFace::GetLinkByNode( const TLinkSet&      links,
02413                                    const TChainLink&    avoidLink,
02414                                    const SMDS_MeshNode* nodeToContain) const
02415   {
02416     for ( int i = 0; i < _sides.size(); ++i )
02417       if ( avoidLink._qlink != _sides[i] &&
02418            (_sides[i]->node1() == nodeToContain || _sides[i]->node2() == nodeToContain ))
02419         return links.find( _sides[ i ]);
02420     return links.end();
02421   }
02422 
02423   //================================================================================
02427   //================================================================================
02428 
02429   gp_Vec QFace::LinkNorm(const int i, SMESH_MesherHelper* /*uvHelper*/) const
02430   {
02431     gp_Vec norm, vecOut;
02432 //     if ( uvHelper ) {
02433 //       TopoDS_Face face = TopoDS::Face( uvHelper->GetSubShape());
02434 //       const SMDS_MeshNode* inFaceNode = uvHelper->GetNodeUVneedInFaceNode() ? GetNodeInFace() : 0;
02435 //       gp_XY uv1 = uvHelper->GetNodeUV( face, _sides[i]->node1(), inFaceNode );
02436 //       gp_XY uv2 = uvHelper->GetNodeUV( face, _sides[i]->node2(), inFaceNode );
02437 //       norm.SetCoord( uv1.Y() - uv2.Y(), uv2.X() - uv1.X(), 0 );
02438 
02439 //       const QLink* otherLink = _sides[(i + 1) % _sides.size()];
02440 //       const SMDS_MeshNode* otherNode =
02441 //         otherLink->node1() == _sides[i]->node1() ? otherLink->node2() : otherLink->node1();
02442 //       gp_XY pIn = uvHelper->GetNodeUV( face, otherNode, inFaceNode );
02443 //       vecOut.SetCoord( uv1.X() - pIn.X(), uv1.Y() - pIn.Y(), 0 );
02444 //     }
02445 //     else {
02446       norm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2()));
02447       gp_XYZ pIn = ( XYZ( _sides[0]->node1() ) +
02448                      XYZ( _sides[0]->node2() ) +
02449                      XYZ( _sides[1]->node1() )) / 3.;
02450       vecOut.SetXYZ( _sides[i]->MiddlePnt() - pIn );
02451       //}
02452     if ( norm * vecOut < 0 )
02453       norm.Reverse();
02454     double mag2 = norm.SquareMagnitude();
02455     if ( mag2 > numeric_limits<double>::min() )
02456       norm /= sqrt( mag2 );
02457     return norm;
02458   }
02459   //================================================================================
02472   //================================================================================
02473 
02474   double QFace::MoveByBoundary( const TChainLink&   theLink,
02475                                 const gp_Vec&       theRefVec,
02476                                 const TLinkSet&     theLinks,
02477                                 SMESH_MesherHelper* theFaceHelper,
02478                                 const double        thePrevLen,
02479                                 const int           theStep,
02480                                 gp_Vec*             theLinkNorm,
02481                                 double              theSign) const
02482   {
02483     if ( !theStep )
02484       return thePrevLen; // propagation limit reached
02485 
02486     int iL; // index of theLink
02487     for ( iL = 0; iL < _sides.size(); ++iL )
02488       if ( theLink._qlink == _sides[ iL ])
02489         break;
02490 
02491     MSG(string(theStep,'.')<<" Ref( "<<theRefVec.X()<<","<<theRefVec.Y()<<","<<theRefVec.Z()<<" )"
02492         <<" thePrevLen " << thePrevLen);
02493     MSG(string(theStep,'.')<<" "<<*theLink._qlink);
02494 
02495     gp_Vec linkNorm = -LinkNorm( iL/*, theFaceHelper*/ ); // normal to theLink
02496     double refProj = theRefVec * linkNorm; // project movement vector to normal of theLink
02497     if ( theStep == theFirstStep )
02498       theSign = refProj < 0. ? -1. : 1.;
02499     else if ( theSign * refProj < 0.4 * theRefVec.Magnitude())
02500       return thePrevLen; // to propagate movement forward only, not in side dir or backward
02501 
02502     int iL1 = (iL + 1) % 3, iL2 = (iL + 2) % 3; // indices of the two other links of triangle
02503     TLinkInSet link1 = theLinks.find( _sides[iL1] );
02504     TLinkInSet link2 = theLinks.find( _sides[iL2] );
02505     if ( link1 == theLinks.end() || link2 == theLinks.end() )
02506       return thePrevLen;
02507     const QFace* f1 = link1->NextFace( this ); // adjacent faces
02508     const QFace* f2 = link2->NextFace( this );
02509 
02510     // propagate to adjacent faces till limit step or boundary
02511     double len1 = thePrevLen + (theLink->MiddlePnt() - _sides[iL1]->MiddlePnt()).Modulus();
02512     double len2 = thePrevLen + (theLink->MiddlePnt() - _sides[iL2]->MiddlePnt()).Modulus();
02513     gp_Vec linkDir1(0,0,0); // initialize to avoid valgrind error ("Conditional jump...")
02514     gp_Vec linkDir2(0,0,0);
02515     try {
02516       OCC_CATCH_SIGNALS;
02517       if ( f1 && theLink->MediumPos() <= (*link1)->MediumPos() )
02518         len1 = f1->MoveByBoundary
02519           ( *link1, theRefVec, theLinks, theFaceHelper, len1, theStep-1, &linkDir1, theSign);
02520       else
02521         linkDir1 = LinkNorm( iL1/*, theFaceHelper*/ );
02522     } catch (...) {
02523       MSG( " --------------- EXCEPTION");
02524       return thePrevLen;
02525     }
02526     try {
02527       OCC_CATCH_SIGNALS;
02528       if ( f2 && theLink->MediumPos() <= (*link2)->MediumPos() )
02529         len2 = f2->MoveByBoundary
02530           ( *link2, theRefVec, theLinks, theFaceHelper, len2, theStep-1, &linkDir2, theSign);
02531       else
02532         linkDir2 = LinkNorm( iL2/*, theFaceHelper*/ );
02533     } catch (...) {
02534       MSG( " --------------- EXCEPTION");
02535       return thePrevLen;
02536     }
02537 
02538     double fullLen = 0;
02539     if ( theStep != theFirstStep )
02540     {
02541       // choose chain length by direction of propagation most codirected with theRefVec
02542       bool choose1 = ( theRefVec * linkDir1 * theSign > theRefVec * linkDir2 * theSign );
02543       fullLen = choose1 ? len1 : len2;
02544       double r = thePrevLen / fullLen;
02545 
02546       gp_Vec move = linkNorm * refProj * ( 1 - r );
02547       theLink->Move( move, true );
02548 
02549       MSG(string(theStep,'.')<<" Move "<< theLink->_mediumNode->GetID()<<
02550           " by " << refProj * ( 1 - r ) << " following " <<
02551           (choose1 ? *link1->_qlink : *link2->_qlink));
02552 
02553       if ( theLinkNorm ) *theLinkNorm = linkNorm;
02554     }
02555     return fullLen;
02556   }
02557 
02558   //================================================================================
02562   //================================================================================
02563 
02564   bool QFace::IsSpoiled(const QLink* bentLink ) const
02565   {
02566     // code is valid for convex faces only
02567     gp_XYZ gc(0,0,0);
02568     for ( TIDSortedNodeSet::const_iterator n = begin(); n!=end(); ++n)
02569       gc += XYZ( *n ) / size();
02570     for (unsigned i = 0; i < _sides.size(); ++i )
02571     {
02572       if ( _sides[i] == bentLink ) continue;
02573       gp_Vec linkNorm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2()));
02574       gp_Vec vecOut( gc, _sides[i]->MiddlePnt() );
02575       if ( linkNorm * vecOut < 0 )
02576         linkNorm.Reverse();
02577       double mag2 = linkNorm.SquareMagnitude();
02578       if ( mag2 > numeric_limits<double>::min() )
02579         linkNorm /= sqrt( mag2 );
02580       gp_Vec vecBent    ( _sides[i]->MiddlePnt(), bentLink->MediumPnt());
02581       gp_Vec vecStraight( _sides[i]->MiddlePnt(), bentLink->MiddlePnt());
02582       if ( vecBent * linkNorm > -0.1*vecStraight.Magnitude() )
02583         return true;
02584     }
02585     return false;
02586     
02587   }
02588 
02589   //================================================================================
02593   //================================================================================
02594 
02595   void QLink::SetContinuesFaces() const
02596   {
02597     //       x0         x - QLink, [-|] - QFace, v - volume
02598     //   v0  |   v1   
02599     //       |          Between _faces of link x2 two vertical faces are continues
02600     // x1----x2-----x3  and two horizontal faces are continues. We set vertical faces
02601     //       |          to _faces[0] and _faces[1] and horizontal faces to
02602     //   v2  |   v3     _faces[2] and _faces[3] (or vise versa).
02603     //       x4
02604 
02605     if ( _faces.empty() )
02606       return;
02607     int iFaceCont = -1, nbBoundary = 0, iBoundary[2]={-1,-1};
02608     if ( _faces[0]->IsBoundary() )
02609       iBoundary[ nbBoundary++ ] = 0;
02610     for ( int iF = 1; iFaceCont < 0 && iF < _faces.size(); ++iF )
02611     {
02612       // look for a face bounding none of volumes bound by _faces[0]
02613       bool sameVol = false;
02614       int nbVol = _faces[iF]->NbVolumes();
02615       for ( int iV = 0; !sameVol && iV < nbVol; ++iV )
02616         sameVol = ( _faces[iF]->_volumes[iV] == _faces[0]->_volumes[0] ||
02617                     _faces[iF]->_volumes[iV] == _faces[0]->_volumes[1]);
02618       if ( !sameVol )
02619         iFaceCont = iF;
02620       if ( _faces[iF]->IsBoundary() )
02621         iBoundary[ nbBoundary++ ] = iF;
02622     }
02623     // Set continues faces: arrange _faces to have
02624     // _faces[0] continues to _faces[1]
02625     // _faces[2] continues to _faces[3]
02626     if ( nbBoundary == 2 ) // bnd faces are continues
02627     {
02628       if (( iBoundary[0] < 2 ) != ( iBoundary[1] < 2 ))
02629       {
02630         int iNear0 = iBoundary[0] < 2 ? 1-iBoundary[0] : 5-iBoundary[0];
02631         std::swap( _faces[ iBoundary[1] ], _faces[iNear0] );
02632       }
02633     }
02634     else if ( iFaceCont > 0 ) // continues faces found
02635     {
02636       if ( iFaceCont != 1 )
02637         std::swap( _faces[1], _faces[iFaceCont] );
02638     }
02639     else if ( _faces.size() > 1 ) // not found, set NULL by the first face
02640     {
02641       _faces.insert( ++_faces.begin(), 0 );
02642     }
02643   }
02644   //================================================================================
02648   //================================================================================
02649 
02650   const QFace* QLink::GetContinuesFace( const QFace* face ) const
02651   {
02652     for ( int i = 0; i < _faces.size(); ++i ) {
02653       if ( _faces[i] == face ) {
02654         int iF = i < 2 ? 1-i : 5-i;
02655         return iF < _faces.size() ? _faces[iF] : 0;
02656       }
02657     }
02658     return 0;
02659   }
02660   //================================================================================
02664   //================================================================================
02665 
02666   bool QLink::OnBoundary() const
02667   {
02668     for ( int i = 0; i < _faces.size(); ++i )
02669       if (_faces[i] && _faces[i]->IsBoundary()) return true;
02670     return false;
02671   }
02672   //================================================================================
02676   //================================================================================
02677 
02678   gp_Vec TChainLink::Normal() const {
02679     gp_Vec norm;
02680     if (_qfaces[0]) norm  = _qfaces[0]->_normal;
02681     if (_qfaces[1]) norm += _qfaces[1]->_normal;
02682     return norm;
02683   }
02684   //================================================================================
02688   //================================================================================
02689 
02690   bool TChainLink::IsStraight() const
02691   {
02692     bool isStraight = _qlink->IsStraight();
02693     if ( isStraight && _qfaces[0] && !_qfaces[1] )
02694     {
02695       int i = _qfaces[0]->LinkIndex( _qlink );
02696       int iOpp = ( i + 2 ) % _qfaces[0]->_sides.size();
02697       gp_XYZ mid1 = _qlink->MiddlePnt();
02698       gp_XYZ mid2 = _qfaces[0]->_sides[ iOpp ]->MiddlePnt();
02699       double faceSize2 = (mid1-mid2).SquareModulus();
02700       isStraight = _qlink->_nodeMove.SquareMagnitude() < 1/10./10. * faceSize2;
02701     }
02702     return isStraight;
02703   }
02704   
02705   //================================================================================
02709   //================================================================================
02710 
02711   void fixPrism( TChain& allLinks )
02712   {
02713     // separate boundary links from internal ones
02714     typedef set<const QLink*/*, QLink::PtrComparator*/> QLinkSet;
02715     QLinkSet interLinks, bndLinks1, bndLink2;
02716 
02717     bool isCurved = false;
02718     for ( TChain::iterator lnk = allLinks.begin(); lnk != allLinks.end(); ++lnk ) {
02719       if ( (*lnk)->OnBoundary() )
02720         bndLinks1.insert( lnk->_qlink );
02721       else
02722         interLinks.insert( lnk->_qlink );
02723       isCurved = isCurved || !lnk->IsStraight();
02724     }
02725     if ( !isCurved )
02726       return; // no need to move
02727 
02728     QLinkSet *curBndLinks = &bndLinks1, *newBndLinks = &bndLink2;
02729 
02730     while ( !interLinks.empty() && !curBndLinks->empty() )
02731     {
02732       // propagate movement from boundary links to connected internal links
02733       QLinkSet::iterator bnd = curBndLinks->begin(), bndEnd = curBndLinks->end();
02734       for ( ; bnd != bndEnd; ++bnd )
02735       {
02736         const QLink* bndLink = *bnd;
02737         for ( int i = 0; i < bndLink->_faces.size(); ++i ) // loop on faces of bndLink
02738         {
02739           const QFace* face = bndLink->_faces[i]; // quadrange lateral face of a prism
02740           if ( !face ) continue;
02741           // find and move internal link opposite to bndLink within the face
02742           int interInd = ( face->LinkIndex( bndLink ) + 2 ) % face->_sides.size();
02743           const QLink* interLink = face->_sides[ interInd ];
02744           QLinkSet::iterator pInterLink = interLinks.find( interLink );
02745           if ( pInterLink == interLinks.end() ) continue; // not internal link
02746           interLink->Move( bndLink->_nodeMove );
02747           // treated internal links become new boundary ones
02748           interLinks. erase( pInterLink );
02749           newBndLinks->insert( interLink );
02750         }
02751       }
02752       curBndLinks->clear();
02753       std::swap( curBndLinks, newBndLinks );
02754     }
02755   }
02756 
02757   //================================================================================
02761   //================================================================================
02762 
02763   void fixTriaNearBoundary( TChain & allLinks, SMESH_MesherHelper& /*helper*/)
02764   {
02765     if ( allLinks.empty() ) return;
02766 
02767     TLinkSet linkSet( allLinks.begin(), allLinks.end());
02768     TLinkInSet linkIt = linkSet.begin(), linksEnd = linkSet.end();
02769 
02770     for ( linkIt = linkSet.begin(); linkIt != linksEnd; ++linkIt)
02771     {
02772       if ( linkIt->IsBoundary() && !linkIt->IsStraight() && linkIt->_qfaces[0])
02773       {
02774         // move iff a boundary link is bent towards inside of a face (issue 0021084)
02775         const QFace* face = linkIt->_qfaces[0];
02776         gp_XYZ pIn = ( face->_sides[0]->MiddlePnt() +
02777                        face->_sides[1]->MiddlePnt() +
02778                        face->_sides[2]->MiddlePnt() ) / 3.;
02779         gp_XYZ insideDir( pIn - (*linkIt)->MiddlePnt());
02780         bool linkBentInside = ((*linkIt)->_nodeMove.Dot( insideDir ) > 0 );
02781         //if ( face->IsSpoiled( linkIt->_qlink ))
02782         if ( linkBentInside )
02783           face->MoveByBoundary( *linkIt, (*linkIt)->_nodeMove, linkSet );
02784       }
02785     }
02786   }
02787 
02788   //================================================================================
02792   //================================================================================
02793 
02794   enum TSplitTriaResult {
02795     _OK, _NO_CORNERS, _FEW_ROWS, _MANY_ROWS, _NO_SIDELINK, _BAD_MIDQUAD, _NOT_RECT,
02796     _NO_MIDQUAD, _NO_UPTRIA, _BAD_SET_SIZE, _BAD_CORNER, _BAD_START, _NO_BOTLINK, _TWISTED_CHAIN };
02797 
02798   TSplitTriaResult splitTrianglesIntoChains( TChain &            allLinks,
02799                                              vector< TChain> &   resultChains,
02800                                              SMDS_TypeOfPosition pos )
02801   {
02802     // put links in the set and evalute number of result chains by number of boundary links
02803     TLinkSet linkSet;
02804     int nbBndLinks = 0;
02805     for ( TChain::iterator lnk = allLinks.begin(); lnk != allLinks.end(); ++lnk ) {
02806       linkSet.insert( *lnk );
02807       nbBndLinks += lnk->IsBoundary();
02808     }
02809     resultChains.clear();
02810     resultChains.reserve( nbBndLinks / 2 );
02811 
02812     TLinkInSet linkIt, linksEnd = linkSet.end();
02813 
02814     // find a boundary link with corner node; corner node has position pos-2
02815     // i.e. SMDS_TOP_VERTEX for links on faces and SMDS_TOP_EDGE for
02816     // links in volume
02817     SMDS_TypeOfPosition cornerPos = SMDS_TypeOfPosition(pos-2);
02818     const SMDS_MeshNode* corner = 0;
02819     for ( linkIt = linkSet.begin(); linkIt != linksEnd; ++linkIt )
02820       if ( linkIt->IsBoundary() && (corner = (*linkIt)->EndPosNode(cornerPos)))
02821         break;
02822     if ( !corner)
02823       return _NO_CORNERS;
02824 
02825     TLinkInSet           startLink = linkIt;
02826     const SMDS_MeshNode* startCorner = corner;
02827     vector< TChain* >    rowChains;
02828     int iCol = 0;
02829 
02830     while ( startLink != linksEnd) // loop on columns
02831     {
02832       // We suppose we have a rectangular structure like shown here. We have found a
02833       //               corner of the rectangle (startCorner) and a boundary link sharing  
02834       //    |/  |/  |  the startCorner (startLink). We are going to loop on rows of the   
02835       //  --o---o---o  structure making several chains at once. One chain (columnChain)   
02836       //    |\  |  /|  starts at startLink and continues upward (we look at the structure 
02837       //  \ | \ | / |  from such point that startLink is on the bottom of the structure). 
02838       //   \|  \|/  |  While going upward we also fill horizontal chains (rowChains) we   
02839       //  --o---o---o  encounter.                                                         
02840       //   /|\  |\  |
02841       //  / | \ | \ |  startCorner
02842       //    |  \|  \|,'
02843       //  --o---o---o
02844       //          `.startLink
02845 
02846       if ( resultChains.size() == nbBndLinks / 2 )
02847         return _NOT_RECT;
02848       resultChains.push_back( TChain() );
02849       TChain& columnChain = resultChains.back();
02850 
02851       TLinkInSet botLink = startLink; // current horizontal link to go up from
02852       corner = startCorner; // current corner the botLink ends at
02853       int iRow = 0;
02854       while ( botLink != linksEnd ) // loop on rows
02855       {
02856         // add botLink to the columnChain
02857         columnChain.push_back( *botLink );
02858 
02859         const QFace* botTria = botLink->_qfaces[0]; // bottom triangle bound by botLink
02860         if ( !botTria )
02861         { // the column ends
02862           if ( botLink == startLink )
02863             return _TWISTED_CHAIN; // issue 0020951
02864           linkSet.erase( botLink );
02865           if ( iRow != rowChains.size() )
02866             return _FEW_ROWS; // different nb of rows in columns
02867           break;
02868         }
02869         // find the link dividing the quadrangle (midQuadLink) and vertical boundary
02870         // link ending at <corner> (sideLink); there are two cases:
02871         // 1) midQuadLink does not end at <corner>, then we easily find it by botTria,
02872         //   since midQuadLink is not at boundary while sideLink is.
02873         // 2) midQuadLink ends at <corner>
02874         bool isCase2;
02875         TLinkInSet midQuadLink = linksEnd;
02876         TLinkInSet sideLink = botTria->GetBoundaryLink( linkSet, *botLink, &midQuadLink,
02877                                                         corner, &isCase2 );
02878         if ( isCase2 ) { // find midQuadLink among links of botTria
02879           midQuadLink = botTria->GetLinkByNode( linkSet, *botLink, corner );
02880           if ( midQuadLink->IsBoundary() )
02881             return _BAD_MIDQUAD;
02882         }
02883         if ( sideLink == linksEnd || midQuadLink == linksEnd || sideLink == midQuadLink )
02884           return sideLink == linksEnd ? _NO_SIDELINK : _NO_MIDQUAD;
02885 
02886         // fill chains
02887         columnChain.push_back( *midQuadLink );
02888         if ( iRow >= rowChains.size() ) {
02889           if ( iCol > 0 )
02890             return _MANY_ROWS; // different nb of rows in columns
02891           if ( resultChains.size() == nbBndLinks / 2 )
02892             return _NOT_RECT;
02893           resultChains.push_back( TChain() );
02894           rowChains.push_back( & resultChains.back() );
02895         }
02896         rowChains[iRow]->push_back( *sideLink );
02897         rowChains[iRow]->push_back( *midQuadLink );
02898 
02899         const QFace* upTria = midQuadLink->NextFace( botTria ); // upper tria of the rectangle
02900         if ( !upTria)
02901           return _NO_UPTRIA;
02902         if ( iRow == 0 ) {
02903           // prepare startCorner and startLink for the next column
02904           startCorner = startLink->NextNode( startCorner );
02905           if (isCase2)
02906             startLink = botTria->GetBoundaryLink( linkSet, *botLink, 0, startCorner );
02907           else
02908             startLink = upTria->GetBoundaryLink( linkSet, *midQuadLink, 0, startCorner );
02909           // check if no more columns remains
02910           if ( startLink != linksEnd ) {
02911             const SMDS_MeshNode* botNode = startLink->NextNode( startCorner );
02912             if ( (isCase2 ? botTria : upTria)->Contains( botNode ))
02913               startLink = linksEnd; // startLink bounds upTria or botTria
02914             else if ( startLink == botLink || startLink == midQuadLink || startLink == sideLink )
02915               return _BAD_START;
02916           }
02917         }
02918         // find bottom link and corner for the next row
02919         corner = sideLink->NextNode( corner );
02920         // next bottom link ends at the new corner
02921         linkSet.erase( botLink );
02922         botLink = upTria->GetLinkByNode( linkSet, (isCase2 ? *sideLink : *midQuadLink), corner );
02923         if ( botLink == linksEnd || botLink == midQuadLink || botLink == sideLink)
02924           return _NO_BOTLINK;
02925         if ( midQuadLink == startLink || sideLink == startLink )
02926           return _TWISTED_CHAIN; // issue 0020951
02927         linkSet.erase( midQuadLink );
02928         linkSet.erase( sideLink );
02929 
02930         // make faces neighboring the found ones be boundary
02931         if ( startLink != linksEnd ) {
02932           const QFace* tria = isCase2 ? botTria : upTria;
02933           for ( int iL = 0; iL < 3; ++iL ) {
02934             linkIt = linkSet.find( tria->_sides[iL] );
02935             if ( linkIt != linksEnd )
02936               linkIt->RemoveFace( tria );
02937           }
02938         }
02939         if ( botLink->_qfaces[0] == upTria || botLink->_qfaces[1] == upTria )
02940           botLink->RemoveFace( upTria ); // make next botTria first in vector
02941 
02942         iRow++;
02943       } // loop on rows
02944 
02945       iCol++;
02946     }
02947     // In the linkSet, there must remain the last links of rowChains; add them
02948     if ( linkSet.size() != rowChains.size() )
02949       return _BAD_SET_SIZE;
02950     for ( int iRow = 0; iRow < rowChains.size(); ++iRow ) {
02951       // find the link (startLink) ending at startCorner
02952       corner = 0;
02953       for ( startLink = linkSet.begin(); startLink != linksEnd; ++startLink ) {
02954         if ( (*startLink)->node1() == startCorner ) {
02955           corner = (*startLink)->node2(); break;
02956         }
02957         else if ( (*startLink)->node2() == startCorner) {
02958           corner = (*startLink)->node1(); break;
02959         }
02960       }
02961       if ( startLink == linksEnd )
02962         return _BAD_CORNER;
02963       rowChains[ iRow ]->push_back( *startLink );
02964       linkSet.erase( startLink );
02965       startCorner = corner;
02966     }
02967 
02968     return _OK;
02969   }
02970 } //namespace
02971 
02972 //=======================================================================
02979 //=======================================================================
02980 
02981 void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
02982 {
02983   // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
02984   if ( getenv("NO_FixQuadraticElements") )
02985     return;
02986 
02987   // 0. Apply algorithm to solids or geom faces
02988   // ----------------------------------------------
02989   if ( myShape.IsNull() ) {
02990     if ( !myMesh->HasShapeToMesh() ) return;
02991     SetSubShape( myMesh->GetShapeToMesh() );
02992 
02993 #ifdef _DEBUG_
02994     int nbSolids = 0;
02995     TopTools_IndexedMapOfShape solids;
02996     TopExp::MapShapes(myShape,TopAbs_SOLID,solids);
02997     nbSolids = solids.Extent();
02998 #endif
02999     TopTools_MapOfShape faces; // faces not in solid or in not meshed solid
03000     for ( TopExp_Explorer f(myShape,TopAbs_FACE,TopAbs_SOLID); f.More(); f.Next() ) {
03001       faces.Add( f.Current() ); // not in solid
03002     }
03003     for ( TopExp_Explorer s(myShape,TopAbs_SOLID); s.More(); s.Next() ) {
03004       if ( myMesh->GetSubMesh( s.Current() )->IsEmpty() ) { // get faces of solid
03005         for ( TopExp_Explorer f( s.Current(), TopAbs_FACE); f.More(); f.Next() )
03006           faces.Add( f.Current() ); // in not meshed solid
03007       }
03008       else { // fix nodes in the solid and its faces
03009 #ifdef _DEBUG_
03010         MSG("FIX SOLID " << nbSolids-- << " #" << GetMeshDS()->ShapeToIndex(s.Current()));
03011 #endif
03012         SMESH_MesherHelper h(*myMesh);
03013         h.SetSubShape( s.Current() );
03014         h.FixQuadraticElements(false);
03015       }
03016     }
03017     // fix nodes on geom faces
03018 #ifdef _DEBUG_
03019     int nbfaces = faces.Extent(); /*avoid "unused varianbles": */ nbfaces++, nbfaces--; 
03020 #endif
03021     for ( TopTools_MapIteratorOfMapOfShape fIt( faces ); fIt.More(); fIt.Next() ) {
03022       MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key()));
03023       SMESH_MesherHelper h(*myMesh);
03024       h.SetSubShape( fIt.Key() );
03025       h.FixQuadraticElements(true);
03026     }
03027     //perf_print_all_meters(1);
03028     return;
03029   }
03030 
03031   // 1. Find out type of elements and get iterator on them
03032   // ---------------------------------------------------
03033 
03034   SMDS_ElemIteratorPtr elemIt;
03035   SMDSAbs_ElementType elemType = SMDSAbs_All;
03036 
03037   SMESH_subMesh* submesh = myMesh->GetSubMeshContaining( myShapeID );
03038   if ( !submesh )
03039     return;
03040   if ( SMESHDS_SubMesh* smDS = submesh->GetSubMeshDS() ) {
03041     elemIt = smDS->GetElements();
03042     if ( elemIt->more() ) {
03043       elemType = elemIt->next()->GetType();
03044       elemIt = smDS->GetElements();
03045     }
03046   }
03047   if ( !elemIt || !elemIt->more() || elemType < SMDSAbs_Face )
03048     return;
03049 
03050   // 2. Fill in auxiliary data structures
03051   // ----------------------------------
03052 
03053   set< QLink > links;
03054   set< QFace > faces;
03055   set< QLink >::iterator pLink;
03056   set< QFace >::iterator pFace;
03057 
03058   bool isCurved = false;
03059   //bool hasRectFaces = false;
03060   //set<int> nbElemNodeSet;
03061   SMDS_VolumeTool volTool;
03062 
03063   TIDSortedNodeSet apexOfPyramid;
03064   const int apexIndex = 4;
03065 
03066   if ( elemType == SMDSAbs_Volume )
03067   {
03068     while ( elemIt->more() ) // loop on volumes
03069     {
03070       const SMDS_MeshElement* vol = elemIt->next();
03071       if ( !vol->IsQuadratic() || !volTool.Set( vol ))
03072         return;
03073       double volMinSize2 = -1.;
03074       for ( int iF = 0; iF < volTool.NbFaces(); ++iF ) // loop on faces of volume
03075       {
03076         int nbN = volTool.NbFaceNodes( iF );
03077         //nbElemNodeSet.insert( nbN );
03078         const SMDS_MeshNode** faceNodes = volTool.GetFaceNodes( iF );
03079         vector< const QLink* > faceLinks( nbN/2 );
03080         for ( int iN = 0; iN < nbN; iN += 2 ) // loop on links of a face
03081         {
03082           // store QLink
03083           QLink link( faceNodes[iN], faceNodes[iN+2], faceNodes[iN+1] );
03084           pLink = links.insert( link ).first;
03085           faceLinks[ iN/2 ] = & *pLink;
03086 
03087           if ( link.MediumPos() == SMDS_TOP_3DSPACE )
03088           {
03089             if ( !link.IsStraight() )
03090               return; // already fixed
03091           }
03092           else if ( !isCurved )
03093           {
03094             if ( volMinSize2 < 0 ) volMinSize2 = volTool.MinLinearSize2();
03095             isCurved = !isStraightLink( volMinSize2, link._nodeMove.SquareMagnitude() );
03096           }
03097         }
03098         // store QFace
03099         pFace = faces.insert( QFace( faceLinks )).first;
03100         if ( pFace->NbVolumes() == 0 )
03101           pFace->AddSelfToLinks();
03102         pFace->SetVolume( vol );
03103 //         hasRectFaces = hasRectFaces ||
03104 //           ( volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_HEXA ||
03105 //             volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_PENTA );
03106 #ifdef _DEBUG_
03107         if ( nbN == 6 )
03108           pFace->_face = GetMeshDS()->FindFace(faceNodes[0],faceNodes[2],faceNodes[4]);
03109         else
03110           pFace->_face = GetMeshDS()->FindFace(faceNodes[0],faceNodes[2],
03111                                                faceNodes[4],faceNodes[6] );
03112 #endif
03113       }
03114       // collect pyramid apexes for further correction
03115       if ( vol->NbCornerNodes() == 5 )
03116         apexOfPyramid.insert( vol->GetNode( apexIndex ));
03117     }
03118     set< QLink >::iterator pLink = links.begin();
03119     for ( ; pLink != links.end(); ++pLink )
03120       pLink->SetContinuesFaces();
03121   }
03122   else
03123   {
03124     while ( elemIt->more() ) // loop on faces
03125     {
03126       const SMDS_MeshElement* face = elemIt->next();
03127       if ( !face->IsQuadratic() )
03128         continue;
03129       //nbElemNodeSet.insert( face->NbNodes() );
03130       int nbN = face->NbNodes()/2;
03131       vector< const QLink* > faceLinks( nbN );
03132       for ( int iN = 0; iN < nbN; ++iN ) // loop on links of a face
03133       {
03134         // store QLink
03135         QLink link( face->GetNode(iN), face->GetNode((iN+1)%nbN), face->GetNode(iN+nbN) );
03136         pLink = links.insert( link ).first;
03137         faceLinks[ iN ] = & *pLink;
03138         if ( !isCurved )
03139           isCurved = !link.IsStraight();
03140       }
03141       // store QFace
03142       pFace = faces.insert( QFace( faceLinks )).first;
03143       pFace->AddSelfToLinks();
03144       //hasRectFaces = ( hasRectFaces || nbN == 4 );
03145     }
03146   }
03147   if ( !isCurved )
03148     return; // no curved edges of faces
03149 
03150   // 3. Compute displacement of medium nodes
03151   // ---------------------------------------
03152 
03153   // two loops on QFaces: the first is to treat boundary links, the second is for internal ones
03154   TopLoc_Location loc;
03155   // not treat boundary of volumic submesh
03156   int isInside = ( elemType == SMDSAbs_Volume && volumeOnly ) ? 1 : 0;
03157   for ( ; isInside < 2; ++isInside ) {
03158     MSG( "--------------- LOOP (inside=" << isInside << ") ------------------");
03159     SMDS_TypeOfPosition pos = isInside ? SMDS_TOP_3DSPACE : SMDS_TOP_FACE;
03160     SMDS_TypeOfPosition bndPos = isInside ? SMDS_TOP_FACE : SMDS_TOP_EDGE;
03161 
03162     for ( pFace = faces.begin(); pFace != faces.end(); ++pFace ) {
03163       if ( bool(isInside) == pFace->IsBoundary() )
03164         continue;
03165       for ( int dir = 0; dir < 2; ++dir ) // 2 directions of propagation from the quadrangle
03166       {
03167         MSG( "CHAIN");
03168         // make chain of links connected via continues faces
03169         int error = ERR_OK;
03170         TChain rawChain;
03171         if ( !pFace->GetLinkChain( dir, rawChain, pos, error) && error ==ERR_UNKNOWN ) continue;
03172         rawChain.reverse();
03173         if ( !pFace->GetLinkChain( dir+2, rawChain, pos, error ) && error ==ERR_UNKNOWN ) continue;
03174 
03175         vector< TChain > chains;
03176         if ( error == ERR_OK ) { // chain contains continues rectangles
03177           chains.resize(1);
03178           chains[0].splice( chains[0].begin(), rawChain );
03179         }
03180         else if ( error == ERR_TRI ) {  // chain contains continues triangles
03181           TSplitTriaResult res = splitTrianglesIntoChains( rawChain, chains, pos );
03182           if ( res != _OK ) { // not quadrangles split into triangles
03183             fixTriaNearBoundary( rawChain, *this );
03184             break;
03185           }
03186         }
03187         else if ( error == ERR_PRISM ) { // quadrangle side faces of prisms
03188           fixPrism( rawChain );
03189           break;
03190         }
03191         else {
03192           continue;
03193         }
03194         for ( int iC = 0; iC < chains.size(); ++iC )
03195         {
03196           TChain& chain = chains[iC];
03197           if ( chain.empty() ) continue;
03198           if ( chain.front().IsStraight() && chain.back().IsStraight() ) {
03199             MSG("3D straight - ignore");
03200             continue;
03201           }
03202           if ( chain.front()->MediumPos() > bndPos ||
03203                chain.back() ->MediumPos() > bndPos ) {
03204             MSG("Internal chain - ignore");
03205             continue;
03206           }
03207           // mesure chain length and compute link position along the chain
03208           double chainLen = 0;
03209           vector< double > linkPos;
03210           MSGBEG( "Link medium nodes: ");
03211           TChain::iterator link0 = chain.begin(), link1 = chain.begin(), link2;
03212           for ( ++link1; link1 != chain.end(); ++link1, ++link0 ) {
03213             MSGBEG( (*link0)->_mediumNode->GetID() << "-" <<(*link1)->_mediumNode->GetID()<<" ");
03214             double len = ((*link0)->MiddlePnt() - (*link1)->MiddlePnt()).Modulus();
03215             while ( len < numeric_limits<double>::min() ) { // remove degenerated link
03216               link1 = chain.erase( link1 );
03217               if ( link1 == chain.end() )
03218                 break;
03219               len = ((*link0)->MiddlePnt() - (*link1)->MiddlePnt()).Modulus();
03220             }
03221             chainLen += len;
03222             linkPos.push_back( chainLen );
03223           }
03224           MSG("");
03225           if ( linkPos.size() < 2 )
03226             continue;
03227 
03228           gp_Vec move0 = chain.front()->_nodeMove;
03229           gp_Vec move1 = chain.back ()->_nodeMove;
03230 
03231           TopoDS_Face face;
03232           bool checkUV = true;
03233           if ( !isInside )
03234           {
03235             // compute node displacement of end links of chain in parametric space of face
03236             TChainLink& linkOnFace = *(++chain.begin());
03237             const SMDS_MeshNode* nodeOnFace = linkOnFace->_mediumNode;
03238             TopoDS_Shape f = GetSubShapeByNode( nodeOnFace, GetMeshDS() );
03239             if ( !f.IsNull() && f.ShapeType() == TopAbs_FACE )
03240             {
03241               face = TopoDS::Face( f );
03242               Handle(Geom_Surface) surf = BRep_Tool::Surface(face,loc);
03243               bool isStraight[2];
03244               for ( int is1 = 0; is1 < 2; ++is1 ) // move0 or move1
03245               {
03246                 TChainLink& link = is1 ? chain.back() : chain.front();
03247                 gp_XY uvm = GetNodeUV( face, link->_mediumNode, nodeOnFace, &checkUV);
03248                 gp_XY uv1 = GetNodeUV( face, link->node1(), nodeOnFace, &checkUV);
03249                 gp_XY uv2 = GetNodeUV( face, link->node2(), nodeOnFace, &checkUV);
03250                 gp_XY uv12 = GetMiddleUV( surf, uv1, uv2);
03251                 // uvMove = uvm - uv12
03252                 gp_XY uvMove = applyIn2D(surf, uvm, uv12, gp_XY_Subtracted, /*inPeriod=*/false);
03253                 ( is1 ? move1 : move0 ).SetCoord( uvMove.X(), uvMove.Y(), 0 );
03254                 if ( !is1 ) // correct nodeOnFace for move1 (issue 0020919)
03255                   nodeOnFace = (*(++chain.rbegin()))->_mediumNode;
03256                 isStraight[is1] = isStraightLink( (uv2-uv1).SquareModulus(),
03257                                                   10 * uvMove.SquareModulus());
03258               }
03259               if ( isStraight[0] && isStraight[1] ) {
03260                 MSG("2D straight - ignore");
03261                 continue; // straight - no need to move nodes of internal links
03262               }
03263 
03264               // check if a chain is already fixed
03265               gp_XY uvm = GetNodeUV( face, linkOnFace->_mediumNode, 0, &checkUV);
03266               gp_XY uv1 = GetNodeUV( face, linkOnFace->node1(), nodeOnFace, &checkUV);
03267               gp_XY uv2 = GetNodeUV( face, linkOnFace->node2(), nodeOnFace, &checkUV);
03268               gp_XY uv12 = GetMiddleUV( surf, uv1, uv2);
03269               if (( uvm - uv12 ).SquareModulus() > 1e-10 )
03270               {
03271                 MSG("Already fixed - ignore");
03272                 continue;
03273               }
03274             }
03275           }
03276           gp_Trsf trsf;
03277           if ( isInside || face.IsNull() )
03278           {
03279             // compute node displacement of end links in their local coord systems
03280             {
03281               TChainLink& ln0 = chain.front(), ln1 = *(++chain.begin());
03282               trsf.SetTransformation( gp_Ax3( gp::Origin(), ln0.Normal(),
03283                                               gp_Vec( ln0->MiddlePnt(), ln1->MiddlePnt() )));
03284               move0.Transform(trsf);
03285             }
03286             {
03287               TChainLink& ln0 = *(++chain.rbegin()), ln1 = chain.back();
03288               trsf.SetTransformation( gp_Ax3( gp::Origin(), ln1.Normal(),
03289                                               gp_Vec( ln0->MiddlePnt(), ln1->MiddlePnt() )));
03290               move1.Transform(trsf);
03291             }
03292           }
03293           // compute displacement of medium nodes
03294           link2 = chain.begin();
03295           link0 = link2++;
03296           link1 = link2++;
03297           for ( int i = 0; link2 != chain.end(); ++link0, ++link1, ++link2, ++i )
03298           {
03299             double r = linkPos[i] / chainLen;
03300             // displacement in local coord system
03301             gp_Vec move = (1. - r) * move0 + r * move1;
03302             if ( isInside || face.IsNull()) {
03303               // transform to global
03304               gp_Vec x01( (*link0)->MiddlePnt(), (*link1)->MiddlePnt() );
03305               gp_Vec x12( (*link1)->MiddlePnt(), (*link2)->MiddlePnt() );
03306               gp_Vec x = x01.Normalized() + x12.Normalized();
03307               trsf.SetTransformation( gp_Ax3( gp::Origin(), link1->Normal(), x), gp_Ax3() );
03308               move.Transform(trsf);
03309             }
03310             else {
03311               // compute 3D displacement by 2D one
03312               Handle(Geom_Surface) s = BRep_Tool::Surface(face,loc);
03313               gp_XY oldUV   = GetNodeUV( face, (*link1)->_mediumNode, 0, &checkUV);
03314               gp_XY newUV   = applyIn2D( s, oldUV, gp_XY( move.X(),move.Y()), gp_XY_Added);
03315               gp_Pnt newPnt = s->Value( newUV.X(), newUV.Y());
03316               move = gp_Vec( XYZ((*link1)->_mediumNode), newPnt.Transformed(loc) );
03317 #ifdef _DEBUG_
03318               if ( (XYZ((*link1)->node1()) - XYZ((*link1)->node2())).SquareModulus() <
03319                    move.SquareMagnitude())
03320               {
03321                 gp_XY uv0 = GetNodeUV( face, (*link0)->_mediumNode, 0, &checkUV);
03322                 gp_XY uv2 = GetNodeUV( face, (*link2)->_mediumNode, 0, &checkUV);
03323                 MSG( "TOO LONG MOVE \t" <<
03324                      "uv0: "<<uv0.X()<<", "<<uv0.Y()<<" \t" <<
03325                      "uv2: "<<uv2.X()<<", "<<uv2.Y()<<" \t" <<
03326                      "uvOld: "<<oldUV.X()<<", "<<oldUV.Y()<<" \t" <<
03327                      "newUV: "<<newUV.X()<<", "<<newUV.Y()<<" \t");
03328               }
03329 #endif
03330             }
03331             (*link1)->Move( move );
03332             MSG( "Move " << (*link1)->_mediumNode->GetID() << " following "
03333                  << chain.front()->_mediumNode->GetID() <<"-"
03334                  << chain.back ()->_mediumNode->GetID() <<
03335                  " by " << move.Magnitude());
03336           }
03337         } // loop on chains of links
03338       } // loop on 2 directions of propagation from quadrangle
03339     } // loop on faces
03340   }
03341 
03342   // 4. Move nodes
03343   // -------------
03344 
03345 //   vector<const SMDS_MeshElement*> vols( 100 );
03346 //   vector<double>                  volSize( 100 );
03347 //   int nbVols;
03348 //   bool ok;
03349   for ( pLink = links.begin(); pLink != links.end(); ++pLink ) {
03350     if ( pLink->IsMoved() ) {
03351       gp_Pnt p = pLink->MiddlePnt() + pLink->Move();
03352       GetMeshDS()->MoveNode( pLink->_mediumNode, p.X(), p.Y(), p.Z());
03353       //
03354 //       gp_Pnt pNew = pLink->MiddlePnt() + pLink->Move();
03355 //       if ( pLink->MediumPos() != SMDS_TOP_3DSPACE )
03356 //       {
03357 //         // avoid making distorted volumes near boundary
03358 //         SMDS_ElemIteratorPtr volIt =
03359 //           (*pLink)._mediumNode->GetInverseElementIterator( SMDSAbs_Volume );
03360 //         for ( nbVols = 0; volIt->more() && volTool.Set( volIt->next() ); ++nbVols )
03361 //         {
03362 //           vols   [ nbVols ] = volTool.Element();
03363 //           volSize[ nbVols ] = volTool.GetSize();
03364 //         }
03365 //         gp_Pnt pOld = pLink->MediumPnt();
03366 //         const_cast<SMDS_MeshNode*>( pLink->_mediumNode )->setXYZ( pNew.X(), pNew.Y(), pNew.Z() );
03367 //         ok = true;
03368 //         while ( nbVols-- && ok )
03369 //         {
03370 //           volTool.Set( vols[ nbVols ]);
03371 //           ok = ( volSize[ nbVols ] * volTool.GetSize() > 1e-20 ); 
03372 //         }
03373 //         if ( !ok )
03374 //         {
03375 //           const_cast<SMDS_MeshNode*>( pLink->_mediumNode )->setXYZ( pOld.X(), pOld.Y(), pOld.Z() );
03376 //           MSG( "Do NOT move \t" << pLink->_mediumNode->GetID()
03377 //                << " because of distortion of volume " << vols[ nbVols+1 ]->GetID());
03378 //           continue;
03379 //         }
03380 //       }
03381 //       GetMeshDS()->MoveNode( pLink->_mediumNode, pNew.X(), pNew.Y(), pNew.Z() );
03382     }
03383   }
03384 
03385   //return;
03386 
03387   // issue 0020982
03388   // Move the apex of pyramid together with the most curved link
03389 
03390   TIDSortedNodeSet::iterator apexIt = apexOfPyramid.begin();
03391   for ( ; apexIt != apexOfPyramid.end(); ++apexIt )
03392   {
03393     SMESH_TNodeXYZ apex = *apexIt;
03394 
03395     gp_Vec maxMove( 0,0,0 );
03396     double maxMoveSize2 = 0;
03397 
03398     // shift of node index to get medium nodes between the base nodes
03399     const int base2MediumShift = 5;
03400 
03401     // find maximal movement of medium node
03402     SMDS_ElemIteratorPtr volIt = apex._node->GetInverseElementIterator( SMDSAbs_Volume );
03403     vector< const SMDS_MeshElement* > pyramids;
03404     while ( volIt->more() )
03405     {
03406       const SMDS_MeshElement* pyram = volIt->next();
03407       if ( pyram->GetEntityType() != SMDSEntity_Quad_Pyramid ) continue;
03408       pyramids.push_back( pyram );
03409 
03410       for ( int iBase = 0; iBase < apexIndex; ++iBase )
03411       {
03412         SMESH_TNodeXYZ medium = pyram->GetNode( iBase + base2MediumShift );
03413         if ( medium._node->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE )
03414         {
03415           SMESH_TNodeXYZ n1 = pyram->GetNode( iBase );
03416           SMESH_TNodeXYZ n2 = pyram->GetNode( ( iBase+1 ) % 4 );
03417           gp_Pnt middle = 0.5 * ( n1 + n2 );
03418           gp_Vec move( middle, medium );
03419           double moveSize2 = move.SquareMagnitude();
03420           if ( moveSize2 > maxMoveSize2 )
03421             maxMove = move, maxMoveSize2 = moveSize2;
03422         }
03423       }
03424     }
03425 
03426     // move the apex
03427     if ( maxMoveSize2 > 1e-20 )
03428     {
03429       apex += maxMove.XYZ();
03430       GetMeshDS()->MoveNode( apex._node, apex.X(), apex.Y(), apex.Z());
03431 
03432       // move medium nodes neighboring the apex to the middle
03433       const int base2MediumShift_2 = 9;
03434       for ( unsigned i = 0; i < pyramids.size(); ++i )
03435         for ( int iBase = 0; iBase < apexIndex; ++iBase )
03436         {
03437           SMESH_TNodeXYZ         base = pyramids[i]->GetNode( iBase );
03438           const SMDS_MeshNode* medium = pyramids[i]->GetNode( iBase + base2MediumShift_2 );
03439           gp_XYZ middle = 0.5 * ( apex + base );
03440           GetMeshDS()->MoveNode( medium, middle.X(), middle.Y(), middle.Z());
03441         }
03442     }
03443   }
03444 }
03445