Salome HOME
Regression of XSMESH_TEST/SMESHCOMMON/SMESH_TEST/Grids/smesh/bugs12/M6
[modules/smesh.git] / src / SMESH / SMESH_MeshEditor.cxx
1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
5 //
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
10 //
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14 // Lesser General Public License for more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
19 //
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 //
22
23 // File      : SMESH_MeshEditor.cxx
24 // Created   : Mon Apr 12 16:10:22 2004
25 // Author    : Edward AGAPOV (eap)
26
27 #include "SMESH_MeshEditor.hxx"
28
29 #include "SMDS_FaceOfNodes.hxx"
30 #include "SMDS_VolumeTool.hxx"
31 #include "SMDS_EdgePosition.hxx"
32 #include "SMDS_FacePosition.hxx"
33 #include "SMDS_SpacePosition.hxx"
34 #include "SMDS_MeshGroup.hxx"
35 #include "SMDS_LinearEdge.hxx"
36 #include "SMDS_Downward.hxx"
37 #include "SMDS_SetIterator.hxx"
38
39 #include "SMESHDS_Group.hxx"
40 #include "SMESHDS_Mesh.hxx"
41
42 #include "SMESH_Algo.hxx"
43 #include "SMESH_ControlsDef.hxx"
44 #include "SMESH_Group.hxx"
45 #include "SMESH_MesherHelper.hxx"
46 #include "SMESH_OctreeNode.hxx"
47 #include "SMESH_subMesh.hxx"
48
49 #include <Basics_OCCTVersion.hxx>
50
51 #include "utilities.h"
52
53 #include <BRepAdaptor_Surface.hxx>
54 #include <BRepBuilderAPI_MakeEdge.hxx>
55 #include <BRepClass3d_SolidClassifier.hxx>
56 #include <BRep_Tool.hxx>
57 #include <ElCLib.hxx>
58 #include <Extrema_GenExtPS.hxx>
59 #include <Extrema_POnCurv.hxx>
60 #include <Extrema_POnSurf.hxx>
61 #include <GC_MakeSegment.hxx>
62 #include <Geom2d_Curve.hxx>
63 #include <GeomAPI_ExtremaCurveCurve.hxx>
64 #include <GeomAdaptor_Surface.hxx>
65 #include <Geom_Curve.hxx>
66 #include <Geom_Line.hxx>
67 #include <Geom_Surface.hxx>
68 #include <IntAna_IntConicQuad.hxx>
69 #include <IntAna_Quadric.hxx>
70 #include <Precision.hxx>
71 #include <TColStd_ListOfInteger.hxx>
72 #include <TopAbs_State.hxx>
73 #include <TopExp.hxx>
74 #include <TopExp_Explorer.hxx>
75 #include <TopTools_ListIteratorOfListOfShape.hxx>
76 #include <TopTools_ListOfShape.hxx>
77 #include <TopTools_SequenceOfShape.hxx>
78 #include <TopoDS.hxx>
79 #include <TopoDS_Face.hxx>
80 #include <TopoDS_Solid.hxx>
81 #include <gp.hxx>
82 #include <gp_Ax1.hxx>
83 #include <gp_Dir.hxx>
84 #include <gp_Lin.hxx>
85 #include <gp_Pln.hxx>
86 #include <gp_Trsf.hxx>
87 #include <gp_Vec.hxx>
88 #include <gp_XY.hxx>
89 #include <gp_XYZ.hxx>
90
91 #include <cmath>
92
93 #include <map>
94 #include <set>
95 #include <numeric>
96 #include <limits>
97 #include <algorithm>
98 #include <sstream>
99
100 #include <boost/tuple/tuple.hpp>
101
102 #include <Standard_Failure.hxx>
103 #include <Standard_ErrorHandler.hxx>
104
105 #define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
106
107 using namespace std;
108 using namespace SMESH::Controls;
109
110 typedef SMDS_SetIterator< SMDS_pElement, TIDSortedElemSet::const_iterator> TSetIterator;
111
112 //=======================================================================
113 //function : SMESH_MeshEditor
114 //purpose  :
115 //=======================================================================
116
117 SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh )
118   :myMesh( theMesh ) // theMesh may be NULL
119 {
120 }
121
122 //================================================================================
123 /*!
124  * \brief Clears myLastCreatedNodes and myLastCreatedElems
125  */
126 //================================================================================
127
128 void SMESH_MeshEditor::CrearLastCreated()
129 {
130   myLastCreatedNodes.Clear();
131   myLastCreatedElems.Clear();
132 }
133
134
135 //=======================================================================
136 /*!
137  * \brief Add element
138  */
139 //=======================================================================
140
141 SMDS_MeshElement*
142 SMESH_MeshEditor::AddElement(const vector<const SMDS_MeshNode*> & node,
143                              const SMDSAbs_ElementType            type,
144                              const bool                           isPoly,
145                              const int                            ID,
146                              const double                         ballDiameter)
147 {
148   //MESSAGE("AddElement " <<node.size() << " " << type << " " << isPoly << " " << ID);
149   SMDS_MeshElement* e = 0;
150   int nbnode = node.size();
151   SMESHDS_Mesh* mesh = GetMeshDS();
152   switch ( type ) {
153   case SMDSAbs_Face:
154     if ( !isPoly ) {
155       if      (nbnode == 3) {
156         if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID);
157         else           e = mesh->AddFace      (node[0], node[1], node[2] );
158       }
159       else if (nbnode == 4) {
160         if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID);
161         else           e = mesh->AddFace      (node[0], node[1], node[2], node[3] );
162       }
163       else if (nbnode == 6) {
164         if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
165                                                node[4], node[5], ID);
166         else           e = mesh->AddFace      (node[0], node[1], node[2], node[3],
167                                                node[4], node[5] );
168       }
169       else if (nbnode == 8) {
170         if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
171                                                node[4], node[5], node[6], node[7], ID);
172         else           e = mesh->AddFace      (node[0], node[1], node[2], node[3],
173                                                node[4], node[5], node[6], node[7] );
174       }
175       else if (nbnode == 9) {
176         if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3],
177                                                node[4], node[5], node[6], node[7], node[8], ID);
178         else           e = mesh->AddFace      (node[0], node[1], node[2], node[3],
179                                                node[4], node[5], node[6], node[7], node[8] );
180       }
181     } else {
182       if ( ID >= 1 ) e = mesh->AddPolygonalFaceWithID(node, ID);
183       else           e = mesh->AddPolygonalFace      (node    );
184     }
185     break;
186
187   case SMDSAbs_Volume:
188     if ( !isPoly ) {
189       if      (nbnode == 4) {
190         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID);
191         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3] );
192       }
193       else if (nbnode == 5) {
194         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
195                                                  node[4], ID);
196         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
197                                                  node[4] );
198       }
199       else if (nbnode == 6) {
200         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
201                                                  node[4], node[5], ID);
202         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
203                                                  node[4], node[5] );
204       }
205       else if (nbnode == 8) {
206         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
207                                                  node[4], node[5], node[6], node[7], ID);
208         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
209                                                  node[4], node[5], node[6], node[7] );
210       }
211       else if (nbnode == 10) {
212         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
213                                                  node[4], node[5], node[6], node[7],
214                                                  node[8], node[9], ID);
215         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
216                                                  node[4], node[5], node[6], node[7],
217                                                  node[8], node[9] );
218       }
219       else if (nbnode == 12) {
220         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
221                                                  node[4], node[5], node[6], node[7],
222                                                  node[8], node[9], node[10], node[11], ID);
223         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
224                                                  node[4], node[5], node[6], node[7],
225                                                  node[8], node[9], node[10], node[11] );
226       }
227       else if (nbnode == 13) {
228         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
229                                                  node[4], node[5], node[6], node[7],
230                                                  node[8], node[9], node[10],node[11],
231                                                  node[12],ID);
232         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
233                                                  node[4], node[5], node[6], node[7],
234                                                  node[8], node[9], node[10],node[11],
235                                                  node[12] );
236       }
237       else if (nbnode == 15) {
238         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
239                                                  node[4], node[5], node[6], node[7],
240                                                  node[8], node[9], node[10],node[11],
241                                                  node[12],node[13],node[14],ID);
242         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
243                                                  node[4], node[5], node[6], node[7],
244                                                  node[8], node[9], node[10],node[11],
245                                                  node[12],node[13],node[14] );
246       }
247       else if (nbnode == 20) {
248         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
249                                                  node[4], node[5], node[6], node[7],
250                                                  node[8], node[9], node[10],node[11],
251                                                  node[12],node[13],node[14],node[15],
252                                                  node[16],node[17],node[18],node[19],ID);
253         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
254                                                  node[4], node[5], node[6], node[7],
255                                                  node[8], node[9], node[10],node[11],
256                                                  node[12],node[13],node[14],node[15],
257                                                  node[16],node[17],node[18],node[19] );
258       }
259       else if (nbnode == 27) {
260         if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3],
261                                                  node[4], node[5], node[6], node[7],
262                                                  node[8], node[9], node[10],node[11],
263                                                  node[12],node[13],node[14],node[15],
264                                                  node[16],node[17],node[18],node[19],
265                                                  node[20],node[21],node[22],node[23],
266                                                  node[24],node[25],node[26], ID);
267         else           e = mesh->AddVolume      (node[0], node[1], node[2], node[3],
268                                                  node[4], node[5], node[6], node[7],
269                                                  node[8], node[9], node[10],node[11],
270                                                  node[12],node[13],node[14],node[15],
271                                                  node[16],node[17],node[18],node[19],
272                                                  node[20],node[21],node[22],node[23],
273                                                  node[24],node[25],node[26] );
274       }
275     }
276     break;
277
278   case SMDSAbs_Edge:
279     if ( nbnode == 2 ) {
280       if ( ID >= 1 ) e = mesh->AddEdgeWithID(node[0], node[1], ID);
281       else           e = mesh->AddEdge      (node[0], node[1] );
282     }
283     else if ( nbnode == 3 ) {
284       if ( ID >= 1 ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID);
285       else           e = mesh->AddEdge      (node[0], node[1], node[2] );
286     }
287     break;
288
289   case SMDSAbs_0DElement:
290     if ( nbnode == 1 ) {
291       if ( ID >= 1 ) e = mesh->Add0DElementWithID(node[0], ID);
292       else           e = mesh->Add0DElement      (node[0] );
293     }
294     break;
295
296   case SMDSAbs_Node:
297     if ( ID >= 1 ) e = mesh->AddNodeWithID(node[0]->X(), node[0]->Y(), node[0]->Z(), ID);
298     else           e = mesh->AddNode      (node[0]->X(), node[0]->Y(), node[0]->Z());
299     break;
300
301   case SMDSAbs_Ball:
302     if ( ID >= 1 ) e = mesh->AddBallWithID(node[0], ballDiameter, ID);
303     else           e = mesh->AddBall      (node[0], ballDiameter);
304     break;
305
306   default:;
307   }
308   if ( e ) myLastCreatedElems.Append( e );
309   return e;
310 }
311
312 //=======================================================================
313 /*!
314  * \brief Add element
315  */
316 //=======================================================================
317
318 SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector<int> &       nodeIDs,
319                                                const SMDSAbs_ElementType type,
320                                                const bool                isPoly,
321                                                const int                 ID)
322 {
323   vector<const SMDS_MeshNode*> nodes;
324   nodes.reserve( nodeIDs.size() );
325   vector<int>::const_iterator id = nodeIDs.begin();
326   while ( id != nodeIDs.end() ) {
327     if ( const SMDS_MeshNode* node = GetMeshDS()->FindNode( *id++ ))
328       nodes.push_back( node );
329     else
330       return 0;
331   }
332   return AddElement( nodes, type, isPoly, ID );
333 }
334
335 //=======================================================================
336 //function : Remove
337 //purpose  : Remove a node or an element.
338 //           Modify a compute state of sub-meshes which become empty
339 //=======================================================================
340
341 int SMESH_MeshEditor::Remove (const list< int >& theIDs,
342                               const bool         isNodes )
343 {
344   myLastCreatedElems.Clear();
345   myLastCreatedNodes.Clear();
346
347   SMESHDS_Mesh* aMesh = GetMeshDS();
348   set< SMESH_subMesh *> smmap;
349
350   int removed = 0;
351   list<int>::const_iterator it = theIDs.begin();
352   for ( ; it != theIDs.end(); it++ ) {
353     const SMDS_MeshElement * elem;
354     if ( isNodes )
355       elem = aMesh->FindNode( *it );
356     else
357       elem = aMesh->FindElement( *it );
358     if ( !elem )
359       continue;
360
361     // Notify VERTEX sub-meshes about modification
362     if ( isNodes ) {
363       const SMDS_MeshNode* node = cast2Node( elem );
364       if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
365         if ( int aShapeID = node->getshapeId() )
366           if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
367             smmap.insert( sm );
368     }
369     // Find sub-meshes to notify about modification
370     //     SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
371     //     while ( nodeIt->more() ) {
372     //       const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
373     //       const SMDS_PositionPtr& aPosition = node->GetPosition();
374     //       if ( aPosition.get() ) {
375     //         if ( int aShapeID = aPosition->GetShapeId() ) {
376     //           if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( aShapeID ) )
377     //             smmap.insert( sm );
378     //         }
379     //       }
380     //     }
381
382     // Do remove
383     if ( isNodes )
384       aMesh->RemoveNode( static_cast< const SMDS_MeshNode* >( elem ));
385     else
386       aMesh->RemoveElement( elem );
387     removed++;
388   }
389
390   // Notify sub-meshes about modification
391   if ( !smmap.empty() ) {
392     set< SMESH_subMesh *>::iterator smIt;
393     for ( smIt = smmap.begin(); smIt != smmap.end(); smIt++ )
394       (*smIt)->ComputeStateEngine( SMESH_subMesh::MESH_ENTITY_REMOVED );
395   }
396
397   //   // Check if the whole mesh becomes empty
398   //   if ( SMESH_subMesh * sm = GetMesh()->GetSubMeshContaining( 1 ) )
399   //     sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
400
401   return removed;
402 }
403
404 //================================================================================
405 /*!
406  * \brief Create 0D elements on all nodes of the given object except those
407  *        nodes on which a 0D element already exists.
408  *  \param elements - Elements on whose nodes to create 0D elements; if empty, 
409  *                    the all mesh is treated
410  *  \param all0DElems - returns all 0D elements found or created on nodes of \a elements
411  */
412 //================================================================================
413
414 void SMESH_MeshEditor::Create0DElementsOnAllNodes( const TIDSortedElemSet& elements,
415                                                    TIDSortedElemSet&       all0DElems )
416 {
417   typedef SMDS_SetIterator<const SMDS_MeshElement*, TIDSortedElemSet::const_iterator> TSetIterator;
418   SMDS_ElemIteratorPtr elemIt;
419   if ( elements.empty() )
420     elemIt = GetMeshDS()->elementsIterator( SMDSAbs_Node );
421   else
422     elemIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() ));
423
424   while ( elemIt->more() )
425   {
426     const SMDS_MeshElement* e = elemIt->next();
427     SMDS_ElemIteratorPtr nodeIt = e->nodesIterator();
428     while ( nodeIt->more() )
429     {
430       const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
431       SMDS_ElemIteratorPtr it0D = n->GetInverseElementIterator( SMDSAbs_0DElement );
432       if ( it0D->more() )
433         all0DElems.insert( it0D->next() );
434       else {
435         myLastCreatedElems.Append( GetMeshDS()->Add0DElement( n ));
436         all0DElems.insert( myLastCreatedElems.Last() );
437       }
438     }
439   }
440 }
441
442 //=======================================================================
443 //function : FindShape
444 //purpose  : Return an index of the shape theElem is on
445 //           or zero if a shape not found
446 //=======================================================================
447
448 int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
449 {
450   myLastCreatedElems.Clear();
451   myLastCreatedNodes.Clear();
452
453   SMESHDS_Mesh * aMesh = GetMeshDS();
454   if ( aMesh->ShapeToMesh().IsNull() )
455     return 0;
456
457   int aShapeID = theElem->getshapeId();
458   if ( aShapeID < 1 )
459     return 0;
460
461   if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ))
462     if ( sm->Contains( theElem ))
463       return aShapeID;
464
465   if ( theElem->GetType() == SMDSAbs_Node ) {
466     MESSAGE( ":( Error: invalid myShapeId of node " << theElem->GetID() );
467   }
468   else {
469     MESSAGE( ":( Error: invalid myShapeId of element " << theElem->GetID() );
470   }
471
472   TopoDS_Shape aShape; // the shape a node of theElem is on
473   if ( theElem->GetType() != SMDSAbs_Node )
474   {
475     SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator();
476     while ( nodeIt->more() ) {
477       const SMDS_MeshNode* node = static_cast<const SMDS_MeshNode*>( nodeIt->next() );
478       if ((aShapeID = node->getshapeId()) > 0) {
479         if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ) ) {
480           if ( sm->Contains( theElem ))
481             return aShapeID;
482           if ( aShape.IsNull() )
483             aShape = aMesh->IndexToShape( aShapeID );
484         }
485       }
486     }
487   }
488
489   // None of nodes is on a proper shape,
490   // find the shape among ancestors of aShape on which a node is
491   if ( !aShape.IsNull() ) {
492     TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape ));
493     for ( ; ancIt.More(); ancIt.Next() ) {
494       SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() );
495       if ( sm && sm->Contains( theElem ))
496         return aMesh->ShapeToIndex( ancIt.Value() );
497     }
498   }
499   else
500   {
501     const map<int,SMESHDS_SubMesh*>& id2sm = GetMeshDS()->SubMeshes();
502     map<int,SMESHDS_SubMesh*>::const_iterator id_sm = id2sm.begin();
503     for ( ; id_sm != id2sm.end(); ++id_sm )
504       if ( id_sm->second->Contains( theElem ))
505         return id_sm->first;
506   }
507
508   //MESSAGE ("::FindShape() - SHAPE NOT FOUND")
509   return 0;
510 }
511
512 //=======================================================================
513 //function : IsMedium
514 //purpose  :
515 //=======================================================================
516
517 bool SMESH_MeshEditor::IsMedium(const SMDS_MeshNode*      node,
518                                 const SMDSAbs_ElementType typeToCheck)
519 {
520   bool isMedium = false;
521   SMDS_ElemIteratorPtr it = node->GetInverseElementIterator(typeToCheck);
522   while (it->more() && !isMedium ) {
523     const SMDS_MeshElement* elem = it->next();
524     isMedium = elem->IsMediumNode(node);
525   }
526   return isMedium;
527 }
528
529 //=======================================================================
530 //function : ShiftNodesQuadTria
531 //purpose  : auxilary
532 //           Shift nodes in the array corresponded to quadratic triangle
533 //           example: (0,1,2,3,4,5) -> (1,2,0,4,5,3)
534 //=======================================================================
535 static void ShiftNodesQuadTria(const SMDS_MeshNode* aNodes[])
536 {
537   const SMDS_MeshNode* nd1 = aNodes[0];
538   aNodes[0] = aNodes[1];
539   aNodes[1] = aNodes[2];
540   aNodes[2] = nd1;
541   const SMDS_MeshNode* nd2 = aNodes[3];
542   aNodes[3] = aNodes[4];
543   aNodes[4] = aNodes[5];
544   aNodes[5] = nd2;
545 }
546
547 //=======================================================================
548 //function : edgeConnectivity
549 //purpose  : auxilary
550 //           return number of the edges connected with the theNode.
551 //           if theEdges has connections with the other type of the
552 //           elements, return -1
553 //=======================================================================
554 static int nbEdgeConnectivity(const SMDS_MeshNode* theNode)
555 {
556   SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator();
557   int nb=0;
558   while(elemIt->more()) {
559     elemIt->next();
560     nb++;
561   }
562   return nb;
563 }
564
565
566 //=======================================================================
567 //function : GetNodesFromTwoTria
568 //purpose  : auxilary
569 //           Shift nodes in the array corresponded to quadratic triangle
570 //           example: (0,1,2,3,4,5) -> (1,2,0,4,5,3)
571 //=======================================================================
572 static bool GetNodesFromTwoTria(const SMDS_MeshElement * theTria1,
573                                 const SMDS_MeshElement * theTria2,
574                                 const SMDS_MeshNode* N1[],
575                                 const SMDS_MeshNode* N2[])
576 {
577   SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
578   int i=0;
579   while(i<6) {
580     N1[i] = static_cast<const SMDS_MeshNode*>( it->next() );
581     i++;
582   }
583   if(it->more()) return false;
584   it = theTria2->nodesIterator();
585   i=0;
586   while(i<6) {
587     N2[i] = static_cast<const SMDS_MeshNode*>( it->next() );
588     i++;
589   }
590   if(it->more()) return false;
591
592   int sames[3] = {-1,-1,-1};
593   int nbsames = 0;
594   int j;
595   for(i=0; i<3; i++) {
596     for(j=0; j<3; j++) {
597       if(N1[i]==N2[j]) {
598         sames[i] = j;
599         nbsames++;
600         break;
601       }
602     }
603   }
604   if(nbsames!=2) return false;
605   if(sames[0]>-1) {
606     ShiftNodesQuadTria(N1);
607     if(sames[1]>-1) {
608       ShiftNodesQuadTria(N1);
609     }
610   }
611   i = sames[0] + sames[1] + sames[2];
612   for(; i<2; i++) {
613     ShiftNodesQuadTria(N2);
614   }
615   // now we receive following N1 and N2 (using numeration as above image)
616   // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
617   // i.e. first nodes from both arrays determ new diagonal
618   return true;
619 }
620
621 //=======================================================================
622 //function : InverseDiag
623 //purpose  : Replace two neighbour triangles with ones built on the same 4 nodes
624 //           but having other common link.
625 //           Return False if args are improper
626 //=======================================================================
627
628 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1,
629                                     const SMDS_MeshElement * theTria2 )
630 {
631   MESSAGE("InverseDiag");
632   myLastCreatedElems.Clear();
633   myLastCreatedNodes.Clear();
634
635   if (!theTria1 || !theTria2)
636     return false;
637
638   const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( theTria1 );
639   if (!F1) return false;
640   const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( theTria2 );
641   if (!F2) return false;
642   if ((theTria1->GetEntityType() == SMDSEntity_Triangle) &&
643       (theTria2->GetEntityType() == SMDSEntity_Triangle)) {
644
645     //  1 +--+ A  theTria1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
646     //    | /|    theTria2: ( B A 2 ) B->1 ( 1 A 2 )   |\ |
647     //    |/ |                                         | \|
648     //  B +--+ 2                                     B +--+ 2
649
650     // put nodes in array and find out indices of the same ones
651     const SMDS_MeshNode* aNodes [6];
652     int sameInd [] = { 0, 0, 0, 0, 0, 0 };
653     int i = 0;
654     SMDS_ElemIteratorPtr it = theTria1->nodesIterator();
655     while ( it->more() ) {
656       aNodes[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
657
658       if ( i > 2 ) // theTria2
659         // find same node of theTria1
660         for ( int j = 0; j < 3; j++ )
661           if ( aNodes[ i ] == aNodes[ j ]) {
662             sameInd[ j ] = i;
663             sameInd[ i ] = j;
664             break;
665           }
666       // next
667       i++;
668       if ( i == 3 ) {
669         if ( it->more() )
670           return false; // theTria1 is not a triangle
671         it = theTria2->nodesIterator();
672       }
673       if ( i == 6 && it->more() )
674         return false; // theTria2 is not a triangle
675     }
676
677     // find indices of 1,2 and of A,B in theTria1
678     int iA = 0, iB = 0, i1 = 0, i2 = 0;
679     for ( i = 0; i < 6; i++ ) {
680       if ( sameInd [ i ] == 0 ) {
681         if ( i < 3 ) i1 = i;
682         else         i2 = i;
683       }
684       else if (i < 3) {
685         if ( iA ) iB = i;
686         else      iA = i;
687       }
688     }
689     // nodes 1 and 2 should not be the same
690     if ( aNodes[ i1 ] == aNodes[ i2 ] )
691       return false;
692
693     // theTria1: A->2
694     aNodes[ iA ] = aNodes[ i2 ];
695     // theTria2: B->1
696     aNodes[ sameInd[ iB ]] = aNodes[ i1 ];
697
698     GetMeshDS()->ChangeElementNodes( theTria1, aNodes, 3 );
699     GetMeshDS()->ChangeElementNodes( theTria2, &aNodes[ 3 ], 3 );
700
701     return true;
702
703   } // end if(F1 && F2)
704
705   // check case of quadratic faces
706   if (theTria1->GetEntityType() != SMDSEntity_Quad_Triangle)
707     return false;
708   if (theTria2->GetEntityType() != SMDSEntity_Quad_Triangle)
709     return false;
710
711   //       5
712   //  1 +--+--+ 2  theTria1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
713   //    |    /|    theTria2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
714   //    |   / |
715   //  7 +  +  + 6
716   //    | /9  |
717   //    |/    |
718   //  4 +--+--+ 3
719   //       8
720
721   const SMDS_MeshNode* N1 [6];
722   const SMDS_MeshNode* N2 [6];
723   if(!GetNodesFromTwoTria(theTria1,theTria2,N1,N2))
724     return false;
725   // now we receive following N1 and N2 (using numeration as above image)
726   // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
727   // i.e. first nodes from both arrays determ new diagonal
728
729   const SMDS_MeshNode* N1new [6];
730   const SMDS_MeshNode* N2new [6];
731   N1new[0] = N1[0];
732   N1new[1] = N2[0];
733   N1new[2] = N2[1];
734   N1new[3] = N1[4];
735   N1new[4] = N2[3];
736   N1new[5] = N1[5];
737   N2new[0] = N1[0];
738   N2new[1] = N1[1];
739   N2new[2] = N2[0];
740   N2new[3] = N1[3];
741   N2new[4] = N2[5];
742   N2new[5] = N1[4];
743   // replaces nodes in faces
744   GetMeshDS()->ChangeElementNodes( theTria1, N1new, 6 );
745   GetMeshDS()->ChangeElementNodes( theTria2, N2new, 6 );
746
747   return true;
748 }
749
750 //=======================================================================
751 //function : findTriangles
752 //purpose  : find triangles sharing theNode1-theNode2 link
753 //=======================================================================
754
755 static bool findTriangles(const SMDS_MeshNode *    theNode1,
756                           const SMDS_MeshNode *    theNode2,
757                           const SMDS_MeshElement*& theTria1,
758                           const SMDS_MeshElement*& theTria2)
759 {
760   if ( !theNode1 || !theNode2 ) return false;
761
762   theTria1 = theTria2 = 0;
763
764   set< const SMDS_MeshElement* > emap;
765   SMDS_ElemIteratorPtr it = theNode1->GetInverseElementIterator(SMDSAbs_Face);
766   while (it->more()) {
767     const SMDS_MeshElement* elem = it->next();
768     if ( elem->NbNodes() == 3 )
769       emap.insert( elem );
770   }
771   it = theNode2->GetInverseElementIterator(SMDSAbs_Face);
772   while (it->more()) {
773     const SMDS_MeshElement* elem = it->next();
774     if ( emap.find( elem ) != emap.end() ) {
775       if ( theTria1 ) {
776         // theTria1 must be element with minimum ID
777         if( theTria1->GetID() < elem->GetID() ) {
778           theTria2 = elem;
779         }
780         else {
781           theTria2 = theTria1;
782           theTria1 = elem;
783         }
784         break;
785       }
786       else {
787         theTria1 = elem;
788       }
789     }
790   }
791   return ( theTria1 && theTria2 );
792 }
793
794 //=======================================================================
795 //function : InverseDiag
796 //purpose  : Replace two neighbour triangles sharing theNode1-theNode2 link
797 //           with ones built on the same 4 nodes but having other common link.
798 //           Return false if proper faces not found
799 //=======================================================================
800
801 bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshNode * theNode1,
802                                     const SMDS_MeshNode * theNode2)
803 {
804   myLastCreatedElems.Clear();
805   myLastCreatedNodes.Clear();
806
807   MESSAGE( "::InverseDiag()" );
808
809   const SMDS_MeshElement *tr1, *tr2;
810   if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
811     return false;
812
813   const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( tr1 );
814   if (!F1) return false;
815   const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( tr2 );
816   if (!F2) return false;
817   if ((tr1->GetEntityType() == SMDSEntity_Triangle) &&
818       (tr2->GetEntityType() == SMDSEntity_Triangle)) {
819
820     //  1 +--+ A  tr1: ( 1 A B ) A->2 ( 1 2 B ) 1 +--+ A
821     //    | /|    tr2: ( B A 2 ) B->1 ( 1 A 2 )   |\ |
822     //    |/ |                                    | \|
823     //  B +--+ 2                                B +--+ 2
824
825     // put nodes in array
826     // and find indices of 1,2 and of A in tr1 and of B in tr2
827     int i, iA1 = 0, i1 = 0;
828     const SMDS_MeshNode* aNodes1 [3];
829     SMDS_ElemIteratorPtr it;
830     for (i = 0, it = tr1->nodesIterator(); it->more(); i++ ) {
831       aNodes1[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
832       if ( aNodes1[ i ] == theNode1 )
833         iA1 = i; // node A in tr1
834       else if ( aNodes1[ i ] != theNode2 )
835         i1 = i;  // node 1
836     }
837     int iB2 = 0, i2 = 0;
838     const SMDS_MeshNode* aNodes2 [3];
839     for (i = 0, it = tr2->nodesIterator(); it->more(); i++ ) {
840       aNodes2[ i ] = static_cast<const SMDS_MeshNode*>( it->next() );
841       if ( aNodes2[ i ] == theNode2 )
842         iB2 = i; // node B in tr2
843       else if ( aNodes2[ i ] != theNode1 )
844         i2 = i;  // node 2
845     }
846
847     // nodes 1 and 2 should not be the same
848     if ( aNodes1[ i1 ] == aNodes2[ i2 ] )
849       return false;
850
851     // tr1: A->2
852     aNodes1[ iA1 ] = aNodes2[ i2 ];
853     // tr2: B->1
854     aNodes2[ iB2 ] = aNodes1[ i1 ];
855
856     GetMeshDS()->ChangeElementNodes( tr1, aNodes1, 3 );
857     GetMeshDS()->ChangeElementNodes( tr2, aNodes2, 3 );
858
859     return true;
860   }
861
862   // check case of quadratic faces
863   return InverseDiag(tr1,tr2);
864 }
865
866 //=======================================================================
867 //function : getQuadrangleNodes
868 //purpose  : fill theQuadNodes - nodes of a quadrangle resulting from
869 //           fusion of triangles tr1 and tr2 having shared link on
870 //           theNode1 and theNode2
871 //=======================================================================
872
873 bool getQuadrangleNodes(const SMDS_MeshNode *    theQuadNodes [],
874                         const SMDS_MeshNode *    theNode1,
875                         const SMDS_MeshNode *    theNode2,
876                         const SMDS_MeshElement * tr1,
877                         const SMDS_MeshElement * tr2 )
878 {
879   if( tr1->NbNodes() != tr2->NbNodes() )
880     return false;
881   // find the 4-th node to insert into tr1
882   const SMDS_MeshNode* n4 = 0;
883   SMDS_ElemIteratorPtr it = tr2->nodesIterator();
884   int i=0;
885   while ( !n4 && i<3 ) {
886     const SMDS_MeshNode * n = cast2Node( it->next() );
887     i++;
888     bool isDiag = ( n == theNode1 || n == theNode2 );
889     if ( !isDiag )
890       n4 = n;
891   }
892   // Make an array of nodes to be in a quadrangle
893   int iNode = 0, iFirstDiag = -1;
894   it = tr1->nodesIterator();
895   i=0;
896   while ( i<3 ) {
897     const SMDS_MeshNode * n = cast2Node( it->next() );
898     i++;
899     bool isDiag = ( n == theNode1 || n == theNode2 );
900     if ( isDiag ) {
901       if ( iFirstDiag < 0 )
902         iFirstDiag = iNode;
903       else if ( iNode - iFirstDiag == 1 )
904         theQuadNodes[ iNode++ ] = n4; // insert the 4-th node between diagonal nodes
905     }
906     else if ( n == n4 ) {
907       return false; // tr1 and tr2 should not have all the same nodes
908     }
909     theQuadNodes[ iNode++ ] = n;
910   }
911   if ( iNode == 3 ) // diagonal nodes have 0 and 2 indices
912     theQuadNodes[ iNode ] = n4;
913
914   return true;
915 }
916
917 //=======================================================================
918 //function : DeleteDiag
919 //purpose  : Replace two neighbour triangles sharing theNode1-theNode2 link
920 //           with a quadrangle built on the same 4 nodes.
921 //           Return false if proper faces not found
922 //=======================================================================
923
924 bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1,
925                                    const SMDS_MeshNode * theNode2)
926 {
927   myLastCreatedElems.Clear();
928   myLastCreatedNodes.Clear();
929
930   MESSAGE( "::DeleteDiag()" );
931
932   const SMDS_MeshElement *tr1, *tr2;
933   if ( !findTriangles( theNode1, theNode2, tr1, tr2 ))
934     return false;
935
936   const SMDS_VtkFace* F1 = dynamic_cast<const SMDS_VtkFace*>( tr1 );
937   if (!F1) return false;
938   const SMDS_VtkFace* F2 = dynamic_cast<const SMDS_VtkFace*>( tr2 );
939   if (!F2) return false;
940   SMESHDS_Mesh * aMesh = GetMeshDS();
941
942   if ((tr1->GetEntityType() == SMDSEntity_Triangle) &&
943       (tr2->GetEntityType() == SMDSEntity_Triangle)) {
944
945     const SMDS_MeshNode* aNodes [ 4 ];
946     if ( ! getQuadrangleNodes( aNodes, theNode1, theNode2, tr1, tr2 ))
947       return false;
948
949     const SMDS_MeshElement* newElem = 0;
950     newElem = aMesh->AddFace( aNodes[0], aNodes[1], aNodes[2], aNodes[3] );
951     myLastCreatedElems.Append(newElem);
952     AddToSameGroups( newElem, tr1, aMesh );
953     int aShapeId = tr1->getshapeId();
954     if ( aShapeId )
955       {
956         aMesh->SetMeshElementOnShape( newElem, aShapeId );
957       }
958     aMesh->RemoveElement( tr1 );
959     aMesh->RemoveElement( tr2 );
960
961     return true;
962   }
963
964   // check case of quadratic faces
965   if (tr1->GetEntityType() != SMDSEntity_Quad_Triangle)
966     return false;
967   if (tr2->GetEntityType() != SMDSEntity_Quad_Triangle)
968     return false;
969
970   //       5
971   //  1 +--+--+ 2  tr1: (1 2 4 5 9 7) or (2 4 1 9 7 5) or (4 1 2 7 5 9)
972   //    |    /|    tr2: (2 3 4 6 8 9) or (3 4 2 8 9 6) or (4 2 3 9 6 8)
973   //    |   / |
974   //  7 +  +  + 6
975   //    | /9  |
976   //    |/    |
977   //  4 +--+--+ 3
978   //       8
979
980   const SMDS_MeshNode* N1 [6];
981   const SMDS_MeshNode* N2 [6];
982   if(!GetNodesFromTwoTria(tr1,tr2,N1,N2))
983     return false;
984   // now we receive following N1 and N2 (using numeration as above image)
985   // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
986   // i.e. first nodes from both arrays determ new diagonal
987
988   const SMDS_MeshNode* aNodes[8];
989   aNodes[0] = N1[0];
990   aNodes[1] = N1[1];
991   aNodes[2] = N2[0];
992   aNodes[3] = N2[1];
993   aNodes[4] = N1[3];
994   aNodes[5] = N2[5];
995   aNodes[6] = N2[3];
996   aNodes[7] = N1[5];
997
998   const SMDS_MeshElement* newElem = 0;
999   newElem = aMesh->AddFace( aNodes[0], aNodes[1], aNodes[2], aNodes[3],
1000                             aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
1001   myLastCreatedElems.Append(newElem);
1002   AddToSameGroups( newElem, tr1, aMesh );
1003   int aShapeId = tr1->getshapeId();
1004   if ( aShapeId )
1005     {
1006       aMesh->SetMeshElementOnShape( newElem, aShapeId );
1007     }
1008   aMesh->RemoveElement( tr1 );
1009   aMesh->RemoveElement( tr2 );
1010
1011   // remove middle node (9)
1012   GetMeshDS()->RemoveNode( N1[4] );
1013
1014   return true;
1015 }
1016
1017 //=======================================================================
1018 //function : Reorient
1019 //purpose  : Reverse theElement orientation
1020 //=======================================================================
1021
1022 bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem)
1023 {
1024   MESSAGE("Reorient");
1025   myLastCreatedElems.Clear();
1026   myLastCreatedNodes.Clear();
1027
1028   if (!theElem)
1029     return false;
1030   SMDS_ElemIteratorPtr it = theElem->nodesIterator();
1031   if ( !it || !it->more() )
1032     return false;
1033
1034   switch ( theElem->GetType() ) {
1035
1036   case SMDSAbs_Edge:
1037   case SMDSAbs_Face: {
1038     if(!theElem->IsQuadratic()) {
1039       int i = theElem->NbNodes();
1040       vector<const SMDS_MeshNode*> aNodes( i );
1041       while ( it->more() )
1042         aNodes[ --i ]= static_cast<const SMDS_MeshNode*>( it->next() );
1043       return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], theElem->NbNodes() );
1044     }
1045     else {
1046       // quadratic elements
1047       if(theElem->GetType()==SMDSAbs_Edge) {
1048         vector<const SMDS_MeshNode*> aNodes(3);
1049         aNodes[1]= static_cast<const SMDS_MeshNode*>( it->next() );
1050         aNodes[0]= static_cast<const SMDS_MeshNode*>( it->next() );
1051         aNodes[2]= static_cast<const SMDS_MeshNode*>( it->next() );
1052         return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], 3 );
1053       }
1054       else {
1055         int nbn = theElem->NbNodes();
1056         vector<const SMDS_MeshNode*> aNodes(nbn);
1057         aNodes[0]= static_cast<const SMDS_MeshNode*>( it->next() );
1058         int i=1;
1059         for(; i<nbn/2; i++) {
1060           aNodes[nbn/2-i]= static_cast<const SMDS_MeshNode*>( it->next() );
1061         }
1062         for(i=0; i<nbn/2; i++) {
1063           aNodes[nbn-i-1]= static_cast<const SMDS_MeshNode*>( it->next() );
1064         }
1065         return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], nbn );
1066       }
1067     }
1068   }
1069   case SMDSAbs_Volume: {
1070     if (theElem->IsPoly()) {
1071       // TODO reorient vtk polyhedron
1072       MESSAGE("reorient vtk polyhedron ?");
1073       const SMDS_VtkVolume* aPolyedre =
1074         dynamic_cast<const SMDS_VtkVolume*>( theElem );
1075       if (!aPolyedre) {
1076         MESSAGE("Warning: bad volumic element");
1077         return false;
1078       }
1079
1080       int nbFaces = aPolyedre->NbFaces();
1081       vector<const SMDS_MeshNode *> poly_nodes;
1082       vector<int> quantities (nbFaces);
1083
1084       // reverse each face of the polyedre
1085       for (int iface = 1; iface <= nbFaces; iface++) {
1086         int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface);
1087         quantities[iface - 1] = nbFaceNodes;
1088
1089         for (inode = nbFaceNodes; inode >= 1; inode--) {
1090           const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode);
1091           poly_nodes.push_back(curNode);
1092         }
1093       }
1094
1095       return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities );
1096
1097     }
1098     else {
1099       SMDS_VolumeTool vTool;
1100       if ( !vTool.Set( theElem ))
1101         return false;
1102       vTool.Inverse();
1103       MESSAGE("ChangeElementNodes reorient: check vTool.Inverse");
1104       return GetMeshDS()->ChangeElementNodes( theElem, vTool.GetNodes(), vTool.NbNodes() );
1105     }
1106   }
1107   default:;
1108   }
1109
1110   return false;
1111 }
1112
1113 //================================================================================
1114 /*!
1115  * \brief Reorient faces.
1116  * \param theFaces - the faces to reorient. If empty the whole mesh is meant
1117  * \param theDirection - desired direction of normal of \a theFace
1118  * \param theFace - one of \a theFaces that sould be oriented according to
1119  *        \a theDirection and whose orientation defines orientation of other faces
1120  * \return number of reoriented faces.
1121  */
1122 //================================================================================
1123
1124 int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet &       theFaces,
1125                                   const gp_Dir&            theDirection,
1126                                   const SMDS_MeshElement * theFace)
1127 {
1128   int nbReori = 0;
1129   if ( !theFace || theFace->GetType() != SMDSAbs_Face ) return nbReori;
1130
1131   if ( theFaces.empty() )
1132   {
1133     SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator(/*idInceasingOrder=*/true);
1134     while ( fIt->more() )
1135       theFaces.insert( theFaces.end(), fIt->next() );
1136   }
1137
1138   // orient theFace according to theDirection
1139   gp_XYZ normal;
1140   SMESH_Algo::FaceNormal( theFace, normal, /*normalized=*/false );
1141   if ( normal * theDirection.XYZ() < 0 )
1142     nbReori += Reorient( theFace );
1143
1144   // Orient other faces
1145
1146   set< const SMDS_MeshElement* > startFaces, visitedFaces;
1147   TIDSortedElemSet avoidSet;
1148   set< SMESH_TLink > checkedLinks;
1149   pair< set< SMESH_TLink >::iterator, bool > linkIt_isNew;
1150
1151   if ( theFaces.size() > 1 )// leave 1 face to prevent finding not selected faces
1152     theFaces.erase( theFace );
1153   startFaces.insert( theFace );
1154
1155   int nodeInd1, nodeInd2;
1156   const SMDS_MeshElement*           otherFace;
1157   vector< const SMDS_MeshElement* > facesNearLink;
1158   vector< std::pair< int, int > >   nodeIndsOfFace;
1159
1160   set< const SMDS_MeshElement* >::iterator startFace = startFaces.begin();
1161   while ( !startFaces.empty() )
1162   {
1163     startFace = startFaces.begin();
1164     theFace = *startFace;
1165     startFaces.erase( startFace );
1166     if ( !visitedFaces.insert( theFace ).second )
1167       continue;
1168
1169     avoidSet.clear();
1170     avoidSet.insert(theFace);
1171
1172     NLink link( theFace->GetNode( 0 ), 0 );
1173
1174     const int nbNodes = theFace->NbCornerNodes();
1175     for ( int i = 0; i < nbNodes; ++i ) // loop on links of theFace
1176     {
1177       link.second = theFace->GetNode(( i+1 ) % nbNodes );
1178       linkIt_isNew = checkedLinks.insert( link );
1179       if ( !linkIt_isNew.second )
1180       {
1181         // link has already been checked and won't be encountered more
1182         // if the group (theFaces) is manifold
1183         //checkedLinks.erase( linkIt_isNew.first );
1184       }
1185       else
1186       {
1187         facesNearLink.clear();
1188         nodeIndsOfFace.clear();
1189         while (( otherFace = FindFaceInSet( link.first, link.second,
1190                                             theFaces, avoidSet, &nodeInd1, &nodeInd2 )))
1191           if ( otherFace != theFace)
1192           {
1193             facesNearLink.push_back( otherFace );
1194             nodeIndsOfFace.push_back( make_pair( nodeInd1, nodeInd2 ));
1195             avoidSet.insert( otherFace );
1196           }
1197         if ( facesNearLink.size() > 1 )
1198         {
1199           // NON-MANIFOLD mesh shell !
1200           // select a face most co-directed with theFace,
1201           // other faces won't be visited this time
1202           gp_XYZ NF, NOF;
1203           SMESH_Algo::FaceNormal( theFace, NF, /*normalized=*/false );
1204           double proj, maxProj = -1;
1205           for ( size_t i = 0; i < facesNearLink.size(); ++i ) {
1206             SMESH_Algo::FaceNormal( facesNearLink[i], NOF, /*normalized=*/false );
1207             if (( proj = Abs( NF * NOF )) > maxProj ) {
1208               maxProj = proj;
1209               otherFace = facesNearLink[i];
1210               nodeInd1  = nodeIndsOfFace[i].first;
1211               nodeInd2  = nodeIndsOfFace[i].second;
1212             }
1213           }
1214           // not to visit rejected faces
1215           for ( size_t i = 0; i < facesNearLink.size(); ++i )
1216             if ( facesNearLink[i] != otherFace && theFaces.size() > 1 )
1217               visitedFaces.insert( facesNearLink[i] );
1218         }
1219         else if ( facesNearLink.size() == 1 )
1220         {
1221           otherFace = facesNearLink[0];
1222           nodeInd1  = nodeIndsOfFace.back().first;
1223           nodeInd2  = nodeIndsOfFace.back().second;
1224         }
1225         if ( otherFace && otherFace != theFace)
1226         {
1227           // link must be reverse in otherFace if orientation ot otherFace
1228           // is same as that of theFace
1229           if ( abs(nodeInd2-nodeInd1) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
1230           {
1231             nbReori += Reorient( otherFace );
1232           }
1233           startFaces.insert( otherFace );
1234         }
1235       }
1236       std::swap( link.first, link.second ); // reverse the link
1237     }
1238   }
1239   return nbReori;
1240 }
1241
1242 //=======================================================================
1243 //function : getBadRate
1244 //purpose  :
1245 //=======================================================================
1246
1247 static double getBadRate (const SMDS_MeshElement*               theElem,
1248                           SMESH::Controls::NumericalFunctorPtr& theCrit)
1249 {
1250   SMESH::Controls::TSequenceOfXYZ P;
1251   if ( !theElem || !theCrit->GetPoints( theElem, P ))
1252     return 1e100;
1253   return theCrit->GetBadRate( theCrit->GetValue( P ), theElem->NbNodes() );
1254   //return theCrit->GetBadRate( theCrit->GetValue( theElem->GetID() ), theElem->NbNodes() );
1255 }
1256
1257 //=======================================================================
1258 //function : QuadToTri
1259 //purpose  : Cut quadrangles into triangles.
1260 //           theCrit is used to select a diagonal to cut
1261 //=======================================================================
1262
1263 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet &                   theElems,
1264                                   SMESH::Controls::NumericalFunctorPtr theCrit)
1265 {
1266   myLastCreatedElems.Clear();
1267   myLastCreatedNodes.Clear();
1268
1269   MESSAGE( "::QuadToTri()" );
1270
1271   if ( !theCrit.get() )
1272     return false;
1273
1274   SMESHDS_Mesh * aMesh = GetMeshDS();
1275
1276   Handle(Geom_Surface) surface;
1277   SMESH_MesherHelper   helper( *GetMesh() );
1278
1279   TIDSortedElemSet::iterator itElem;
1280   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
1281     const SMDS_MeshElement* elem = *itElem;
1282     if ( !elem || elem->GetType() != SMDSAbs_Face )
1283       continue;
1284     if ( elem->NbCornerNodes() != 4 )
1285       continue;
1286
1287     // retrieve element nodes
1288     vector< const SMDS_MeshNode* > aNodes( elem->begin_nodes(), elem->end_nodes() );
1289
1290     // compare two sets of possible triangles
1291     double aBadRate1, aBadRate2; // to what extent a set is bad
1292     SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
1293     SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
1294     aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
1295
1296     SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
1297     SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
1298     aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
1299
1300     int aShapeId = FindShape( elem );
1301     const SMDS_MeshElement* newElem1 = 0;
1302     const SMDS_MeshElement* newElem2 = 0;
1303
1304     if( !elem->IsQuadratic() ) {
1305
1306       // split liner quadrangle
1307       // for MaxElementLength2D functor we return minimum diagonal for splitting,
1308       // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
1309       if ( aBadRate1 <= aBadRate2 ) {
1310         // tr1 + tr2 is better
1311         newElem1 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
1312         newElem2 = aMesh->AddFace( aNodes[2], aNodes[0], aNodes[1] );
1313       }
1314       else {
1315         // tr3 + tr4 is better
1316         newElem1 = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
1317         newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] );
1318       }
1319     }
1320     else {
1321
1322       // split quadratic quadrangle
1323
1324       // get surface elem is on
1325       if ( aShapeId != helper.GetSubShapeID() ) {
1326         surface.Nullify();
1327         TopoDS_Shape shape;
1328         if ( aShapeId > 0 )
1329           shape = aMesh->IndexToShape( aShapeId );
1330         if ( !shape.IsNull() && shape.ShapeType() == TopAbs_FACE ) {
1331           TopoDS_Face face = TopoDS::Face( shape );
1332           surface = BRep_Tool::Surface( face );
1333           if ( !surface.IsNull() )
1334             helper.SetSubShape( shape );
1335         }
1336       }
1337       // find middle point for (0,1,2,3)
1338       // and create a node in this point;
1339       const SMDS_MeshNode* newN = 0;
1340       if ( aNodes.size() == 9 )
1341       {
1342         // SMDSEntity_BiQuad_Quadrangle
1343         newN = aNodes.back();
1344       }
1345       else
1346       {
1347         gp_XYZ p( 0,0,0 );
1348         if ( surface.IsNull() )
1349         {
1350           for ( int i = 0; i < 4; i++ )
1351             p += gp_XYZ(aNodes[i]->X(), aNodes[i]->Y(), aNodes[i]->Z() );
1352           p /= 4;
1353         }
1354         else
1355         {
1356           const SMDS_MeshNode* inFaceNode = 0;
1357           if ( helper.GetNodeUVneedInFaceNode() )
1358             for ( size_t i = 0; i < aNodes.size() && !inFaceNode; ++i )
1359               if ( aNodes[ i ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1360                 inFaceNode = aNodes[ i ];
1361
1362           TopoDS_Face face = TopoDS::Face( helper.GetSubShape() );
1363           gp_XY uv( 0,0 );
1364           for ( int i = 0; i < 4; i++ )
1365             uv += helper.GetNodeUV( face, aNodes[i], inFaceNode );
1366           uv /= 4.;
1367           p = surface->Value( uv.X(), uv.Y() ).XYZ();
1368         }
1369         newN = aMesh->AddNode( p.X(), p.Y(), p.Z() );
1370         myLastCreatedNodes.Append(newN);
1371       }
1372       // create a new element
1373       if ( aBadRate1 <= aBadRate2 ) {
1374         newElem1 = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
1375                                   aNodes[6], aNodes[7], newN );
1376         newElem2 = aMesh->AddFace(aNodes[2], aNodes[0], aNodes[1],
1377                                   newN,      aNodes[4], aNodes[5] );
1378       }
1379       else {
1380         newElem1 = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
1381                                   aNodes[7], aNodes[4], newN );
1382         newElem2 = aMesh->AddFace(aNodes[3], aNodes[1], aNodes[2],
1383                                   newN,      aNodes[5], aNodes[6] );
1384       }
1385     } // quadratic case
1386
1387     // care of a new element
1388
1389     myLastCreatedElems.Append(newElem1);
1390     myLastCreatedElems.Append(newElem2);
1391     AddToSameGroups( newElem1, elem, aMesh );
1392     AddToSameGroups( newElem2, elem, aMesh );
1393
1394     // put a new triangle on the same shape
1395     if ( aShapeId )
1396       {
1397         aMesh->SetMeshElementOnShape( newElem1, aShapeId );
1398         aMesh->SetMeshElementOnShape( newElem2, aShapeId );
1399       }
1400     aMesh->RemoveElement( elem );
1401   }
1402   return true;
1403 }
1404
1405 //=======================================================================
1406 //function : BestSplit
1407 //purpose  : Find better diagonal for cutting.
1408 //=======================================================================
1409
1410 int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement*              theQuad,
1411                                  SMESH::Controls::NumericalFunctorPtr theCrit)
1412 {
1413   myLastCreatedElems.Clear();
1414   myLastCreatedNodes.Clear();
1415
1416   if (!theCrit.get())
1417     return -1;
1418
1419   if (!theQuad || theQuad->GetType() != SMDSAbs_Face )
1420     return -1;
1421
1422   if( theQuad->NbNodes()==4 ||
1423       (theQuad->NbNodes()==8 && theQuad->IsQuadratic()) ) {
1424
1425     // retrieve element nodes
1426     const SMDS_MeshNode* aNodes [4];
1427     SMDS_ElemIteratorPtr itN = theQuad->nodesIterator();
1428     int i = 0;
1429     //while (itN->more())
1430     while (i<4) {
1431       aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
1432     }
1433     // compare two sets of possible triangles
1434     double aBadRate1, aBadRate2; // to what extent a set is bad
1435     SMDS_FaceOfNodes tr1 ( aNodes[0], aNodes[1], aNodes[2] );
1436     SMDS_FaceOfNodes tr2 ( aNodes[2], aNodes[3], aNodes[0] );
1437     aBadRate1 = getBadRate( &tr1, theCrit ) + getBadRate( &tr2, theCrit );
1438
1439     SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] );
1440     SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] );
1441     aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit );
1442     // for MaxElementLength2D functor we return minimum diagonal for splitting,
1443     // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4)
1444     if (aBadRate1 <= aBadRate2) // tr1 + tr2 is better
1445       return 1; // diagonal 1-3
1446
1447     return 2; // diagonal 2-4
1448   }
1449   return -1;
1450 }
1451
1452 namespace
1453 {
1454   // Methods of splitting volumes into tetra
1455
1456   const int theHexTo5_1[5*4+1] =
1457     {
1458       0, 1, 2, 5,    0, 4, 5, 7,     0, 2, 3, 7,    2, 5, 6, 7,     0, 5, 2, 7,   -1
1459     };
1460   const int theHexTo5_2[5*4+1] =
1461     {
1462       1, 2, 3, 6,    1, 4, 5, 6,     0, 1, 3, 4,    3, 4, 6, 7,     1, 3, 4, 6,   -1
1463     };
1464   const int* theHexTo5[2] = { theHexTo5_1, theHexTo5_2 };
1465
1466   const int theHexTo6_1[6*4+1] =
1467     {
1468       1, 5, 6, 0,    0, 1, 2, 6,     0, 4, 5, 6,    0, 4, 6, 7,     0, 2, 3, 6,   0, 3, 7, 6,  -1
1469     };
1470   const int theHexTo6_2[6*4+1] =
1471     {
1472       2, 6, 7, 1,    1, 2, 3, 7,     1, 5, 6, 7,    1, 5, 7, 4,     1, 3, 0, 7,   1, 0, 4, 7,  -1
1473     };
1474   const int theHexTo6_3[6*4+1] =
1475     {
1476       3, 7, 4, 2,    2, 3, 0, 4,     2, 6, 7, 4,    2, 6, 4, 5,     2, 0, 1, 4,   2, 1, 5, 4,  -1
1477     };
1478   const int theHexTo6_4[6*4+1] =
1479     {
1480       0, 4, 5, 3,    3, 0, 1, 5,     3, 7, 4, 5,    3, 7, 5, 6,     3, 1, 2, 5,   3, 2, 6, 5,  -1
1481     };
1482   const int* theHexTo6[4] = { theHexTo6_1, theHexTo6_2, theHexTo6_3, theHexTo6_4 };
1483
1484   const int thePyraTo2_1[2*4+1] =
1485     {
1486       0, 1, 2, 4,    0, 2, 3, 4,   -1
1487     };
1488   const int thePyraTo2_2[2*4+1] =
1489     {
1490       1, 2, 3, 4,    1, 3, 0, 4,   -1
1491     };
1492   const int* thePyraTo2[2] = { thePyraTo2_1, thePyraTo2_2 };
1493
1494   const int thePentaTo3_1[3*4+1] =
1495     {
1496       0, 1, 2, 3,    1, 3, 4, 2,     2, 3, 4, 5,    -1
1497     };
1498   const int thePentaTo3_2[3*4+1] =
1499     {
1500       1, 2, 0, 4,    2, 4, 5, 0,     0, 4, 5, 3,    -1
1501     };
1502   const int thePentaTo3_3[3*4+1] =
1503     {
1504       2, 0, 1, 5,    0, 5, 3, 1,     1, 5, 3, 4,    -1
1505     };
1506   const int thePentaTo3_4[3*4+1] =
1507     {
1508       0, 1, 2, 3,    1, 3, 4, 5,     2, 3, 1, 5,    -1
1509     };
1510   const int thePentaTo3_5[3*4+1] =
1511     {
1512       1, 2, 0, 4,    2, 4, 5, 3,     0, 4, 2, 3,    -1
1513     };
1514   const int thePentaTo3_6[3*4+1] =
1515     {
1516       2, 0, 1, 5,    0, 5, 3, 4,     1, 5, 0, 4,    -1
1517     };
1518   const int* thePentaTo3[6] = { thePentaTo3_1, thePentaTo3_2, thePentaTo3_3,
1519                                 thePentaTo3_4, thePentaTo3_5, thePentaTo3_6 };
1520
1521   struct TTriangleFacet //!< stores indices of three nodes of tetra facet
1522   {
1523     int _n1, _n2, _n3;
1524     TTriangleFacet(int n1, int n2, int n3): _n1(n1), _n2(n2), _n3(n3) {}
1525     bool contains(int n) const { return ( n == _n1 || n == _n2 || n == _n3 ); }
1526     bool hasAdjacentTetra( const SMDS_MeshElement* elem ) const;
1527   };
1528   struct TSplitMethod
1529   {
1530     int        _nbTetra;
1531     const int* _connectivity; //!< foursomes of tetra connectivy finished by -1
1532     bool       _baryNode;     //!< additional node is to be created at cell barycenter
1533     bool       _ownConn;      //!< to delete _connectivity in destructor
1534     map<int, const SMDS_MeshNode*> _faceBaryNode; //!< map face index to node at BC of face
1535
1536     TSplitMethod( int nbTet=0, const int* conn=0, bool addNode=false)
1537       : _nbTetra(nbTet), _connectivity(conn), _baryNode(addNode), _ownConn(false) {}
1538     ~TSplitMethod() { if ( _ownConn ) delete [] _connectivity; _connectivity = 0; }
1539     bool hasFacet( const TTriangleFacet& facet ) const
1540     {
1541       const int* tetConn = _connectivity;
1542       for ( ; tetConn[0] >= 0; tetConn += 4 )
1543         if (( facet.contains( tetConn[0] ) +
1544               facet.contains( tetConn[1] ) +
1545               facet.contains( tetConn[2] ) +
1546               facet.contains( tetConn[3] )) == 3 )
1547           return true;
1548       return false;
1549     }
1550   };
1551
1552   //=======================================================================
1553   /*!
1554    * \brief return TSplitMethod for the given element
1555    */
1556   //=======================================================================
1557
1558   TSplitMethod getSplitMethod( SMDS_VolumeTool& vol, const int theMethodFlags)
1559   {
1560     const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
1561
1562     // at HEXA_TO_24 method, each face of volume is split into triangles each based on
1563     // an edge and a face barycenter; tertaherdons are based on triangles and
1564     // a volume barycenter
1565     const bool is24TetMode = ( theMethodFlags == SMESH_MeshEditor::HEXA_TO_24 );
1566
1567     // Find out how adjacent volumes are split
1568
1569     vector < list< TTriangleFacet > > triaSplitsByFace( vol.NbFaces() ); // splits of each side
1570     int hasAdjacentSplits = 0, maxTetConnSize = 0;
1571     for ( int iF = 0; iF < vol.NbFaces(); ++iF )
1572     {
1573       int nbNodes = vol.NbFaceNodes( iF ) / iQ;
1574       maxTetConnSize += 4 * ( nbNodes - (is24TetMode ? 0 : 2));
1575       if ( nbNodes < 4 ) continue;
1576
1577       list< TTriangleFacet >& triaSplits = triaSplitsByFace[ iF ];
1578       const int* nInd = vol.GetFaceNodesIndices( iF );
1579       if ( nbNodes == 4 )
1580       {
1581         TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] );
1582         TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] );
1583         if      ( t012.hasAdjacentTetra( vol.Element() )) triaSplits.push_back( t012 );
1584         else if ( t123.hasAdjacentTetra( vol.Element() )) triaSplits.push_back( t123 );
1585       }
1586       else
1587       {
1588         int iCom = 0; // common node of triangle faces to split into
1589         for ( int iVar = 0; iVar < nbNodes; ++iVar, ++iCom )
1590         {
1591           TTriangleFacet t012( nInd[ iQ * ( iCom             )],
1592                                nInd[ iQ * ( (iCom+1)%nbNodes )],
1593                                nInd[ iQ * ( (iCom+2)%nbNodes )]);
1594           TTriangleFacet t023( nInd[ iQ * ( iCom             )],
1595                                nInd[ iQ * ( (iCom+2)%nbNodes )],
1596                                nInd[ iQ * ( (iCom+3)%nbNodes )]);
1597           if ( t012.hasAdjacentTetra( vol.Element() ) && t023.hasAdjacentTetra( vol.Element() ))
1598           {
1599             triaSplits.push_back( t012 );
1600             triaSplits.push_back( t023 );
1601             break;
1602           }
1603         }
1604       }
1605       if ( !triaSplits.empty() )
1606         hasAdjacentSplits = true;
1607     }
1608
1609     // Among variants of split method select one compliant with adjacent volumes
1610
1611     TSplitMethod method;
1612     if ( !vol.Element()->IsPoly() && !is24TetMode )
1613     {
1614       int nbVariants = 2, nbTet = 0;
1615       const int** connVariants = 0;
1616       switch ( vol.Element()->GetEntityType() )
1617       {
1618       case SMDSEntity_Hexa:
1619       case SMDSEntity_Quad_Hexa:
1620       case SMDSEntity_TriQuad_Hexa:
1621         if ( theMethodFlags == SMESH_MeshEditor::HEXA_TO_5 )
1622           connVariants = theHexTo5, nbTet = 5;
1623         else
1624           connVariants = theHexTo6, nbTet = 6, nbVariants = 4;
1625         break;
1626       case SMDSEntity_Pyramid:
1627       case SMDSEntity_Quad_Pyramid:
1628         connVariants = thePyraTo2;  nbTet = 2;
1629         break;
1630       case SMDSEntity_Penta:
1631       case SMDSEntity_Quad_Penta:
1632         connVariants = thePentaTo3; nbTet = 3; nbVariants = 6;
1633         break;
1634       default:
1635         nbVariants = 0;
1636       }
1637       for ( int variant = 0; variant < nbVariants && method._nbTetra == 0; ++variant )
1638       {
1639         // check method compliancy with adjacent tetras,
1640         // all found splits must be among facets of tetras described by this method
1641         method = TSplitMethod( nbTet, connVariants[variant] );
1642         if ( hasAdjacentSplits && method._nbTetra > 0 )
1643         {
1644           bool facetCreated = true;
1645           for ( int iF = 0; facetCreated && iF < triaSplitsByFace.size(); ++iF )
1646           {
1647             list< TTriangleFacet >::const_iterator facet = triaSplitsByFace[iF].begin();
1648             for ( ; facetCreated && facet != triaSplitsByFace[iF].end(); ++facet )
1649               facetCreated = method.hasFacet( *facet );
1650           }
1651           if ( !facetCreated )
1652             method = TSplitMethod(0); // incompatible method
1653         }
1654       }
1655     }
1656     if ( method._nbTetra < 1 )
1657     {
1658       // No standard method is applicable, use a generic solution:
1659       // each facet of a volume is split into triangles and
1660       // each of triangles and a volume barycenter form a tetrahedron.
1661
1662       const bool isHex27 = ( vol.Element()->GetEntityType() == SMDSEntity_TriQuad_Hexa );
1663
1664       int* connectivity = new int[ maxTetConnSize + 1 ];
1665       method._connectivity = connectivity;
1666       method._ownConn = true;
1667       method._baryNode = !isHex27; // to create central node or not
1668
1669       int connSize = 0;
1670       int baryCenInd = vol.NbNodes() - int( isHex27 );
1671       for ( int iF = 0; iF < vol.NbFaces(); ++iF )
1672       {
1673         const int nbNodes = vol.NbFaceNodes( iF ) / iQ;
1674         const int*   nInd = vol.GetFaceNodesIndices( iF );
1675         // find common node of triangle facets of tetra to create
1676         int iCommon = 0; // index in linear numeration
1677         const list< TTriangleFacet >& triaSplits = triaSplitsByFace[ iF ];
1678         if ( !triaSplits.empty() )
1679         {
1680           // by found facets
1681           const TTriangleFacet* facet = &triaSplits.front();
1682           for ( ; iCommon < nbNodes-1 ; ++iCommon )
1683             if ( facet->contains( nInd[ iQ * iCommon ]) &&
1684                  facet->contains( nInd[ iQ * ((iCommon+2)%nbNodes) ]))
1685               break;
1686         }
1687         else if ( nbNodes > 3 && !is24TetMode )
1688         {
1689           // find the best method of splitting into triangles by aspect ratio
1690           SMESH::Controls::NumericalFunctorPtr aspectRatio( new SMESH::Controls::AspectRatio);
1691           map< double, int > badness2iCommon;
1692           const SMDS_MeshNode** nodes = vol.GetFaceNodes( iF );
1693           int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
1694           for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCommon )
1695           {
1696             double badness = 0;
1697             for ( int iLast = iCommon+2; iLast < iCommon+nbNodes; ++iLast )
1698             {
1699               SMDS_FaceOfNodes tria ( nodes[ iQ*( iCommon         )],
1700                                       nodes[ iQ*((iLast-1)%nbNodes)],
1701                                       nodes[ iQ*((iLast  )%nbNodes)]);
1702               badness += getBadRate( &tria, aspectRatio );
1703             }
1704             badness2iCommon.insert( make_pair( badness, iCommon ));
1705           }
1706           // use iCommon with lowest badness
1707           iCommon = badness2iCommon.begin()->second;
1708         }
1709         if ( iCommon >= nbNodes )
1710           iCommon = 0; // something wrong
1711
1712         // fill connectivity of tetrahedra based on a current face
1713         int nbTet = nbNodes - 2;
1714         if ( is24TetMode && nbNodes > 3 && triaSplits.empty())
1715         {
1716           int faceBaryCenInd;
1717           if ( isHex27 )
1718           {
1719             faceBaryCenInd = vol.GetCenterNodeIndex( iF );
1720             method._faceBaryNode[ iF ] = vol.GetNodes()[ faceBaryCenInd ];
1721           }
1722           else
1723           {
1724             method._faceBaryNode[ iF ] = 0;
1725             faceBaryCenInd = baryCenInd + method._faceBaryNode.size();
1726           }
1727           nbTet = nbNodes;
1728           for ( int i = 0; i < nbTet; ++i )
1729           {
1730             int i1 = i, i2 = (i+1) % nbNodes;
1731             if ( !vol.IsFaceExternal( iF )) swap( i1, i2 );
1732             connectivity[ connSize++ ] = nInd[ iQ * i1 ];
1733             connectivity[ connSize++ ] = nInd[ iQ * i2 ];
1734             connectivity[ connSize++ ] = faceBaryCenInd;
1735             connectivity[ connSize++ ] = baryCenInd;
1736           }
1737         }
1738         else
1739         {
1740           for ( int i = 0; i < nbTet; ++i )
1741           {
1742             int i1 = (iCommon+1+i) % nbNodes, i2 = (iCommon+2+i) % nbNodes;
1743             if ( !vol.IsFaceExternal( iF )) swap( i1, i2 );
1744             connectivity[ connSize++ ] = nInd[ iQ * iCommon ];
1745             connectivity[ connSize++ ] = nInd[ iQ * i1 ];
1746             connectivity[ connSize++ ] = nInd[ iQ * i2 ];
1747             connectivity[ connSize++ ] = baryCenInd;
1748           }
1749         }
1750         method._nbTetra += nbTet;
1751
1752       } // loop on volume faces
1753
1754       connectivity[ connSize++ ] = -1;
1755
1756     } // end of generic solution
1757
1758     return method;
1759   }
1760   //================================================================================
1761   /*!
1762    * \brief Check if there is a tetraherdon adjacent to the given element via this facet
1763    */
1764   //================================================================================
1765
1766   bool TTriangleFacet::hasAdjacentTetra( const SMDS_MeshElement* elem ) const
1767   {
1768     // find the tetrahedron including the three nodes of facet
1769     const SMDS_MeshNode* n1 = elem->GetNode(_n1);
1770     const SMDS_MeshNode* n2 = elem->GetNode(_n2);
1771     const SMDS_MeshNode* n3 = elem->GetNode(_n3);
1772     SMDS_ElemIteratorPtr volIt1 = n1->GetInverseElementIterator(SMDSAbs_Volume);
1773     while ( volIt1->more() )
1774     {
1775       const SMDS_MeshElement* v = volIt1->next();
1776       SMDSAbs_EntityType type = v->GetEntityType();
1777       if ( type != SMDSEntity_Tetra && type != SMDSEntity_Quad_Tetra )
1778         continue;
1779       if ( type == SMDSEntity_Quad_Tetra && v->GetNodeIndex( n1 ) > 3 )
1780         continue; // medium node not allowed
1781       const int ind2 = v->GetNodeIndex( n2 );
1782       if ( ind2 < 0 || 3 < ind2 )
1783         continue;
1784       const int ind3 = v->GetNodeIndex( n3 );
1785       if ( ind3 < 0 || 3 < ind3 )
1786         continue;
1787       return true;
1788     }
1789     return false;
1790   }
1791
1792   //=======================================================================
1793   /*!
1794    * \brief A key of a face of volume
1795    */
1796   //=======================================================================
1797
1798   struct TVolumeFaceKey: pair< pair< int, int>, pair< int, int> >
1799   {
1800     TVolumeFaceKey( SMDS_VolumeTool& vol, int iF )
1801     {
1802       TIDSortedNodeSet sortedNodes;
1803       const int iQ = vol.Element()->IsQuadratic() ? 2 : 1;
1804       int nbNodes = vol.NbFaceNodes( iF );
1805       const SMDS_MeshNode** fNodes = vol.GetFaceNodes( iF );
1806       for ( int i = 0; i < nbNodes; i += iQ )
1807         sortedNodes.insert( fNodes[i] );
1808       TIDSortedNodeSet::iterator n = sortedNodes.begin();
1809       first.first   = (*(n++))->GetID();
1810       first.second  = (*(n++))->GetID();
1811       second.first  = (*(n++))->GetID();
1812       second.second = ( sortedNodes.size() > 3 ) ? (*(n++))->GetID() : 0;
1813     }
1814   };
1815 } // namespace
1816
1817 //=======================================================================
1818 //function : SplitVolumesIntoTetra
1819 //purpose  : Split volume elements into tetrahedra.
1820 //=======================================================================
1821
1822 void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
1823                                               const int                theMethodFlags)
1824 {
1825   // std-like iterator on coordinates of nodes of mesh element
1826   typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > NXyzIterator;
1827   NXyzIterator xyzEnd;
1828
1829   SMDS_VolumeTool    volTool;
1830   SMESH_MesherHelper helper( *GetMesh());
1831
1832   SMESHDS_SubMesh* subMesh = 0;//GetMeshDS()->MeshElements(1);
1833   SMESHDS_SubMesh* fSubMesh = 0;//subMesh;
1834
1835   SMESH_SequenceOfElemPtr newNodes, newElems;
1836
1837   // map face of volume to it's baricenrtic node
1838   map< TVolumeFaceKey, const SMDS_MeshNode* > volFace2BaryNode;
1839   double bc[3];
1840
1841   TIDSortedElemSet::const_iterator elem = theElems.begin();
1842   for ( ; elem != theElems.end(); ++elem )
1843   {
1844     if ( (*elem)->GetType() != SMDSAbs_Volume )
1845       continue;
1846     SMDSAbs_EntityType geomType = (*elem)->GetEntityType();
1847     if ( geomType == SMDSEntity_Tetra || geomType == SMDSEntity_Quad_Tetra )
1848       continue;
1849
1850     if ( !volTool.Set( *elem, /*ignoreCentralNodes=*/false )) continue; // strange...
1851
1852     TSplitMethod splitMethod = getSplitMethod( volTool, theMethodFlags );
1853     if ( splitMethod._nbTetra < 1 ) continue;
1854
1855     // find submesh to add new tetras to
1856     if ( !subMesh || !subMesh->Contains( *elem ))
1857     {
1858       int shapeID = FindShape( *elem );
1859       helper.SetSubShape( shapeID ); // helper will add tetras to the found submesh
1860       subMesh = GetMeshDS()->MeshElements( shapeID );
1861     }
1862     int iQ;
1863     if ( (*elem)->IsQuadratic() )
1864     {
1865       iQ = 2;
1866       // add quadratic links to the helper
1867       for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
1868       {
1869         const SMDS_MeshNode** fNodes = volTool.GetFaceNodes( iF );
1870         int nbN = volTool.NbFaceNodes( iF ) - bool( volTool.GetCenterNodeIndex(iF) > 0 );
1871         for ( int iN = 0; iN < nbN; iN += iQ )
1872           helper.AddTLinkNode( fNodes[iN], fNodes[iN+2], fNodes[iN+1] );
1873       }
1874       helper.SetIsQuadratic( true );
1875     }
1876     else
1877     {
1878       iQ = 1;
1879       helper.SetIsQuadratic( false );
1880     }
1881     vector<const SMDS_MeshNode*> nodes( (*elem)->begin_nodes(), (*elem)->end_nodes() );
1882     helper.SetElementsOnShape( true );
1883     if ( splitMethod._baryNode )
1884     {
1885       // make a node at barycenter
1886       volTool.GetBaryCenter( bc[0], bc[1], bc[2] );
1887       SMDS_MeshNode* gcNode = helper.AddNode( bc[0], bc[1], bc[2] );
1888       nodes.push_back( gcNode );
1889       newNodes.Append( gcNode );
1890     }
1891     if ( !splitMethod._faceBaryNode.empty() )
1892     {
1893       // make or find baricentric nodes of faces
1894       map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.begin();
1895       for ( ; iF_n != splitMethod._faceBaryNode.end(); ++iF_n )
1896       {
1897         map< TVolumeFaceKey, const SMDS_MeshNode* >::iterator f_n =
1898           volFace2BaryNode.insert
1899           ( make_pair( TVolumeFaceKey( volTool,iF_n->first ), iF_n->second )).first;
1900         if ( !f_n->second )
1901         {
1902           volTool.GetFaceBaryCenter( iF_n->first, bc[0], bc[1], bc[2] );
1903           newNodes.Append( f_n->second = helper.AddNode( bc[0], bc[1], bc[2] ));
1904         }
1905         nodes.push_back( iF_n->second = f_n->second );
1906       }
1907     }
1908
1909     // make tetras
1910     vector<const SMDS_MeshElement* > tetras( splitMethod._nbTetra ); // splits of a volume
1911     const int* tetConn = splitMethod._connectivity;
1912     for ( int i = 0; i < splitMethod._nbTetra; ++i, tetConn += 4 )
1913       newElems.Append( tetras[ i ] = helper.AddVolume( nodes[ tetConn[0] ],
1914                                                        nodes[ tetConn[1] ],
1915                                                        nodes[ tetConn[2] ],
1916                                                        nodes[ tetConn[3] ]));
1917
1918     ReplaceElemInGroups( *elem, tetras, GetMeshDS() );
1919
1920     // Split faces on sides of the split volume
1921
1922     const SMDS_MeshNode** volNodes = volTool.GetNodes();
1923     for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
1924     {
1925       const int nbNodes = volTool.NbFaceNodes( iF ) / iQ;
1926       if ( nbNodes < 4 ) continue;
1927
1928       // find an existing face
1929       vector<const SMDS_MeshNode*> fNodes( volTool.GetFaceNodes( iF ),
1930                                            volTool.GetFaceNodes( iF ) + volTool.NbFaceNodes( iF ));
1931       while ( const SMDS_MeshElement* face = GetMeshDS()->FindElement( fNodes, SMDSAbs_Face,
1932                                                                        /*noMedium=*/false))
1933       {
1934         // make triangles
1935         helper.SetElementsOnShape( false );
1936         vector< const SMDS_MeshElement* > triangles;
1937
1938         // find submesh to add new triangles in
1939         if ( !fSubMesh || !fSubMesh->Contains( face ))
1940         {
1941           int shapeID = FindShape( face );
1942           fSubMesh = GetMeshDS()->MeshElements( shapeID );
1943         }
1944         map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.find(iF);
1945         if ( iF_n != splitMethod._faceBaryNode.end() )
1946         {
1947           for ( int iN = 0; iN < nbNodes*iQ; iN += iQ )
1948           {
1949             const SMDS_MeshNode* n1 = fNodes[iN];
1950             const SMDS_MeshNode *n2 = fNodes[(iN+iQ)%(nbNodes*iQ)];
1951             const SMDS_MeshNode *n3 = iF_n->second;
1952             if ( !volTool.IsFaceExternal( iF ))
1953               swap( n2, n3 );
1954             triangles.push_back( helper.AddFace( n1,n2,n3 ));
1955
1956             if ( fSubMesh && n3->getshapeId() < 1 )
1957               fSubMesh->AddNode( n3 );
1958           }
1959         }
1960         else
1961         {
1962           // among possible triangles create ones discribed by split method
1963           const int* nInd = volTool.GetFaceNodesIndices( iF );
1964           int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
1965           int iCom = 0; // common node of triangle faces to split into
1966           list< TTriangleFacet > facets;
1967           for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCom )
1968           {
1969             TTriangleFacet t012( nInd[ iQ * ( iCom                )],
1970                                  nInd[ iQ * ( (iCom+1)%nbNodes )],
1971                                  nInd[ iQ * ( (iCom+2)%nbNodes )]);
1972             TTriangleFacet t023( nInd[ iQ * ( iCom                )],
1973                                  nInd[ iQ * ( (iCom+2)%nbNodes )],
1974                                  nInd[ iQ * ( (iCom+3)%nbNodes )]);
1975             if ( splitMethod.hasFacet( t012 ) && splitMethod.hasFacet( t023 ))
1976             {
1977               facets.push_back( t012 );
1978               facets.push_back( t023 );
1979               for ( int iLast = iCom+4; iLast < iCom+nbNodes; ++iLast )
1980                 facets.push_back( TTriangleFacet( nInd[ iQ * ( iCom             )],
1981                                                   nInd[ iQ * ((iLast-1)%nbNodes )],
1982                                                   nInd[ iQ * ((iLast  )%nbNodes )]));
1983               break;
1984             }
1985           }
1986           list< TTriangleFacet >::iterator facet = facets.begin();
1987           for ( ; facet != facets.end(); ++facet )
1988           {
1989             if ( !volTool.IsFaceExternal( iF ))
1990               swap( facet->_n2, facet->_n3 );
1991             triangles.push_back( helper.AddFace( volNodes[ facet->_n1 ],
1992                                                  volNodes[ facet->_n2 ],
1993                                                  volNodes[ facet->_n3 ]));
1994           }
1995         }
1996         for ( int i = 0; i < triangles.size(); ++i )
1997         {
1998           if ( !triangles[i] ) continue;
1999           if ( fSubMesh )
2000             fSubMesh->AddElement( triangles[i]);
2001           newElems.Append( triangles[i] );
2002         }
2003         ReplaceElemInGroups( face, triangles, GetMeshDS() );
2004         GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false );
2005       }
2006
2007     } // loop on volume faces to split them into triangles
2008
2009     GetMeshDS()->RemoveFreeElement( *elem, subMesh, /*fromGroups=*/false );
2010
2011     if ( geomType == SMDSEntity_TriQuad_Hexa )
2012     {
2013       // remove medium nodes that could become free
2014       for ( int i = 20; i < volTool.NbNodes(); ++i )
2015         if ( volNodes[i]->NbInverseElements() == 0 )
2016           GetMeshDS()->RemoveNode( volNodes[i] );
2017     }
2018   } // loop on volumes to split
2019
2020   myLastCreatedNodes = newNodes;
2021   myLastCreatedElems = newElems;
2022 }
2023
2024 //=======================================================================
2025 //function : AddToSameGroups
2026 //purpose  : add elemToAdd to the groups the elemInGroups belongs to
2027 //=======================================================================
2028
2029 void SMESH_MeshEditor::AddToSameGroups (const SMDS_MeshElement* elemToAdd,
2030                                         const SMDS_MeshElement* elemInGroups,
2031                                         SMESHDS_Mesh *          aMesh)
2032 {
2033   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2034   if (!groups.empty()) {
2035     set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2036     for ( ; grIt != groups.end(); grIt++ ) {
2037       SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2038       if ( group && group->Contains( elemInGroups ))
2039         group->SMDSGroup().Add( elemToAdd );
2040     }
2041   }
2042 }
2043
2044
2045 //=======================================================================
2046 //function : RemoveElemFromGroups
2047 //purpose  : Remove removeelem to the groups the elemInGroups belongs to
2048 //=======================================================================
2049 void SMESH_MeshEditor::RemoveElemFromGroups (const SMDS_MeshElement* removeelem,
2050                                              SMESHDS_Mesh *          aMesh)
2051 {
2052   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2053   if (!groups.empty())
2054   {
2055     set<SMESHDS_GroupBase*>::const_iterator GrIt = groups.begin();
2056     for (; GrIt != groups.end(); GrIt++)
2057     {
2058       SMESHDS_Group* grp = dynamic_cast<SMESHDS_Group*>(*GrIt);
2059       if (!grp || grp->IsEmpty()) continue;
2060       grp->SMDSGroup().Remove(removeelem);
2061     }
2062   }
2063 }
2064
2065 //================================================================================
2066 /*!
2067  * \brief Replace elemToRm by elemToAdd in the all groups
2068  */
2069 //================================================================================
2070
2071 void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
2072                                             const SMDS_MeshElement* elemToAdd,
2073                                             SMESHDS_Mesh *          aMesh)
2074 {
2075   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2076   if (!groups.empty()) {
2077     set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2078     for ( ; grIt != groups.end(); grIt++ ) {
2079       SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2080       if ( group && group->SMDSGroup().Remove( elemToRm ) && elemToAdd )
2081         group->SMDSGroup().Add( elemToAdd );
2082     }
2083   }
2084 }
2085
2086 //================================================================================
2087 /*!
2088  * \brief Replace elemToRm by elemToAdd in the all groups
2089  */
2090 //================================================================================
2091
2092 void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement*                elemToRm,
2093                                             const vector<const SMDS_MeshElement*>& elemToAdd,
2094                                             SMESHDS_Mesh *                         aMesh)
2095 {
2096   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2097   if (!groups.empty())
2098   {
2099     set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2100     for ( ; grIt != groups.end(); grIt++ ) {
2101       SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2102       if ( group && group->SMDSGroup().Remove( elemToRm ) )
2103         for ( int i = 0; i < elemToAdd.size(); ++i )
2104           group->SMDSGroup().Add( elemToAdd[ i ] );
2105     }
2106   }
2107 }
2108
2109 //=======================================================================
2110 //function : QuadToTri
2111 //purpose  : Cut quadrangles into triangles.
2112 //           theCrit is used to select a diagonal to cut
2113 //=======================================================================
2114
2115 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
2116                                   const bool         the13Diag)
2117 {
2118   myLastCreatedElems.Clear();
2119   myLastCreatedNodes.Clear();
2120
2121   MESSAGE( "::QuadToTri()" );
2122
2123   SMESHDS_Mesh * aMesh = GetMeshDS();
2124
2125   Handle(Geom_Surface) surface;
2126   SMESH_MesherHelper   helper( *GetMesh() );
2127
2128   TIDSortedElemSet::iterator itElem;
2129   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
2130     const SMDS_MeshElement* elem = *itElem;
2131     if ( !elem || elem->GetType() != SMDSAbs_Face )
2132       continue;
2133     bool isquad = elem->NbNodes()==4 || elem->NbNodes()==8;
2134     if(!isquad) continue;
2135
2136     if(elem->NbNodes()==4) {
2137       // retrieve element nodes
2138       const SMDS_MeshNode* aNodes [4];
2139       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2140       int i = 0;
2141       while ( itN->more() )
2142         aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
2143
2144       int aShapeId = FindShape( elem );
2145       const SMDS_MeshElement* newElem1 = 0;
2146       const SMDS_MeshElement* newElem2 = 0;
2147       if ( the13Diag ) {
2148         newElem1 = aMesh->AddFace( aNodes[2], aNodes[0], aNodes[1] );
2149         newElem2 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
2150       }
2151       else {
2152         newElem1 = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
2153         newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] );
2154       }
2155       myLastCreatedElems.Append(newElem1);
2156       myLastCreatedElems.Append(newElem2);
2157       // put a new triangle on the same shape and add to the same groups
2158       if ( aShapeId )
2159         {
2160           aMesh->SetMeshElementOnShape( newElem1, aShapeId );
2161           aMesh->SetMeshElementOnShape( newElem2, aShapeId );
2162         }
2163       AddToSameGroups( newElem1, elem, aMesh );
2164       AddToSameGroups( newElem2, elem, aMesh );
2165       //aMesh->RemoveFreeElement(elem, aMesh->MeshElements(aShapeId), true);
2166       aMesh->RemoveElement( elem );
2167     }
2168
2169     // Quadratic quadrangle
2170
2171     if( elem->NbNodes()==8 && elem->IsQuadratic() ) {
2172
2173       // get surface elem is on
2174       int aShapeId = FindShape( elem );
2175       if ( aShapeId != helper.GetSubShapeID() ) {
2176         surface.Nullify();
2177         TopoDS_Shape shape;
2178         if ( aShapeId > 0 )
2179           shape = aMesh->IndexToShape( aShapeId );
2180         if ( !shape.IsNull() && shape.ShapeType() == TopAbs_FACE ) {
2181           TopoDS_Face face = TopoDS::Face( shape );
2182           surface = BRep_Tool::Surface( face );
2183           if ( !surface.IsNull() )
2184             helper.SetSubShape( shape );
2185         }
2186       }
2187
2188       const SMDS_MeshNode* aNodes [8];
2189       const SMDS_MeshNode* inFaceNode = 0;
2190       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2191       int i = 0;
2192       while ( itN->more() ) {
2193         aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
2194         if ( !inFaceNode && helper.GetNodeUVneedInFaceNode() &&
2195              aNodes[ i-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
2196         {
2197           inFaceNode = aNodes[ i-1 ];
2198         }
2199       }
2200
2201       // find middle point for (0,1,2,3)
2202       // and create a node in this point;
2203       gp_XYZ p( 0,0,0 );
2204       if ( surface.IsNull() ) {
2205         for(i=0; i<4; i++)
2206           p += gp_XYZ(aNodes[i]->X(), aNodes[i]->Y(), aNodes[i]->Z() );
2207         p /= 4;
2208       }
2209       else {
2210         TopoDS_Face geomFace = TopoDS::Face( helper.GetSubShape() );
2211         gp_XY uv( 0,0 );
2212         for(i=0; i<4; i++)
2213           uv += helper.GetNodeUV( geomFace, aNodes[i], inFaceNode );
2214         uv /= 4.;
2215         p = surface->Value( uv.X(), uv.Y() ).XYZ();
2216       }
2217       const SMDS_MeshNode* newN = aMesh->AddNode( p.X(), p.Y(), p.Z() );
2218       myLastCreatedNodes.Append(newN);
2219
2220       // create a new element
2221       const SMDS_MeshElement* newElem1 = 0;
2222       const SMDS_MeshElement* newElem2 = 0;
2223       if ( the13Diag ) {
2224         newElem1 = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
2225                                   aNodes[6], aNodes[7], newN );
2226         newElem2 = aMesh->AddFace(aNodes[2], aNodes[0], aNodes[1],
2227                                   newN,      aNodes[4], aNodes[5] );
2228       }
2229       else {
2230         newElem1 = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
2231                                   aNodes[7], aNodes[4], newN );
2232         newElem2 = aMesh->AddFace(aNodes[3], aNodes[1], aNodes[2],
2233                                   newN,      aNodes[5], aNodes[6] );
2234       }
2235       myLastCreatedElems.Append(newElem1);
2236       myLastCreatedElems.Append(newElem2);
2237       // put a new triangle on the same shape and add to the same groups
2238       if ( aShapeId )
2239         {
2240           aMesh->SetMeshElementOnShape( newElem1, aShapeId );
2241           aMesh->SetMeshElementOnShape( newElem2, aShapeId );
2242         }
2243       AddToSameGroups( newElem1, elem, aMesh );
2244       AddToSameGroups( newElem2, elem, aMesh );
2245       aMesh->RemoveElement( elem );
2246     }
2247   }
2248
2249   return true;
2250 }
2251
2252 //=======================================================================
2253 //function : getAngle
2254 //purpose  :
2255 //=======================================================================
2256
2257 double getAngle(const SMDS_MeshElement * tr1,
2258                 const SMDS_MeshElement * tr2,
2259                 const SMDS_MeshNode *    n1,
2260                 const SMDS_MeshNode *    n2)
2261 {
2262   double angle = 2. * M_PI; // bad angle
2263
2264   // get normals
2265   SMESH::Controls::TSequenceOfXYZ P1, P2;
2266   if ( !SMESH::Controls::NumericalFunctor::GetPoints( tr1, P1 ) ||
2267        !SMESH::Controls::NumericalFunctor::GetPoints( tr2, P2 ))
2268     return angle;
2269   gp_Vec N1,N2;
2270   if(!tr1->IsQuadratic())
2271     N1 = gp_Vec( P1(2) - P1(1) ) ^ gp_Vec( P1(3) - P1(1) );
2272   else
2273     N1 = gp_Vec( P1(3) - P1(1) ) ^ gp_Vec( P1(5) - P1(1) );
2274   if ( N1.SquareMagnitude() <= gp::Resolution() )
2275     return angle;
2276   if(!tr2->IsQuadratic())
2277     N2 = gp_Vec( P2(2) - P2(1) ) ^ gp_Vec( P2(3) - P2(1) );
2278   else
2279     N2 = gp_Vec( P2(3) - P2(1) ) ^ gp_Vec( P2(5) - P2(1) );
2280   if ( N2.SquareMagnitude() <= gp::Resolution() )
2281     return angle;
2282
2283   // find the first diagonal node n1 in the triangles:
2284   // take in account a diagonal link orientation
2285   const SMDS_MeshElement *nFirst[2], *tr[] = { tr1, tr2 };
2286   for ( int t = 0; t < 2; t++ ) {
2287     SMDS_ElemIteratorPtr it = tr[ t ]->nodesIterator();
2288     int i = 0, iDiag = -1;
2289     while ( it->more()) {
2290       const SMDS_MeshElement *n = it->next();
2291       if ( n == n1 || n == n2 ) {
2292         if ( iDiag < 0)
2293           iDiag = i;
2294         else {
2295           if ( i - iDiag == 1 )
2296             nFirst[ t ] = ( n == n1 ? n2 : n1 );
2297           else
2298             nFirst[ t ] = n;
2299           break;
2300         }
2301       }
2302       i++;
2303     }
2304   }
2305   if ( nFirst[ 0 ] == nFirst[ 1 ] )
2306     N2.Reverse();
2307
2308   angle = N1.Angle( N2 );
2309   //SCRUTE( angle );
2310   return angle;
2311 }
2312
2313 // =================================================
2314 // class generating a unique ID for a pair of nodes
2315 // and able to return nodes by that ID
2316 // =================================================
2317 class LinkID_Gen {
2318 public:
2319
2320   LinkID_Gen( const SMESHDS_Mesh* theMesh )
2321     :myMesh( theMesh ), myMaxID( theMesh->MaxNodeID() + 1)
2322   {}
2323
2324   long GetLinkID (const SMDS_MeshNode * n1,
2325                   const SMDS_MeshNode * n2) const
2326   {
2327     return ( Min(n1->GetID(),n2->GetID()) * myMaxID + Max(n1->GetID(),n2->GetID()));
2328   }
2329
2330   bool GetNodes (const long             theLinkID,
2331                  const SMDS_MeshNode* & theNode1,
2332                  const SMDS_MeshNode* & theNode2) const
2333   {
2334     theNode1 = myMesh->FindNode( theLinkID / myMaxID );
2335     if ( !theNode1 ) return false;
2336     theNode2 = myMesh->FindNode( theLinkID % myMaxID );
2337     if ( !theNode2 ) return false;
2338     return true;
2339   }
2340
2341 private:
2342   LinkID_Gen();
2343   const SMESHDS_Mesh* myMesh;
2344   long                myMaxID;
2345 };
2346
2347
2348 //=======================================================================
2349 //function : TriToQuad
2350 //purpose  : Fuse neighbour triangles into quadrangles.
2351 //           theCrit is used to select a neighbour to fuse with.
2352 //           theMaxAngle is a max angle between element normals at which
2353 //           fusion is still performed.
2354 //=======================================================================
2355
2356 bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
2357                                   SMESH::Controls::NumericalFunctorPtr theCrit,
2358                                   const double                         theMaxAngle)
2359 {
2360   myLastCreatedElems.Clear();
2361   myLastCreatedNodes.Clear();
2362
2363   MESSAGE( "::TriToQuad()" );
2364
2365   if ( !theCrit.get() )
2366     return false;
2367
2368   SMESHDS_Mesh * aMesh = GetMeshDS();
2369
2370   // Prepare data for algo: build
2371   // 1. map of elements with their linkIDs
2372   // 2. map of linkIDs with their elements
2373
2374   map< SMESH_TLink, list< const SMDS_MeshElement* > > mapLi_listEl;
2375   map< SMESH_TLink, list< const SMDS_MeshElement* > >::iterator itLE;
2376   map< const SMDS_MeshElement*, set< SMESH_TLink > >  mapEl_setLi;
2377   map< const SMDS_MeshElement*, set< SMESH_TLink > >::iterator itEL;
2378
2379   TIDSortedElemSet::iterator itElem;
2380   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
2381     const SMDS_MeshElement* elem = *itElem;
2382     if(!elem || elem->GetType() != SMDSAbs_Face ) continue;
2383     bool IsTria = elem->NbNodes()==3 || (elem->NbNodes()==6 && elem->IsQuadratic());
2384     if(!IsTria) continue;
2385
2386     // retrieve element nodes
2387     const SMDS_MeshNode* aNodes [4];
2388     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2389     int i = 0;
2390     while ( i<3 )
2391       aNodes[ i++ ] = cast2Node( itN->next() );
2392     aNodes[ 3 ] = aNodes[ 0 ];
2393
2394     // fill maps
2395     for ( i = 0; i < 3; i++ ) {
2396       SMESH_TLink link( aNodes[i], aNodes[i+1] );
2397       // check if elements sharing a link can be fused
2398       itLE = mapLi_listEl.find( link );
2399       if ( itLE != mapLi_listEl.end() ) {
2400         if ((*itLE).second.size() > 1 ) // consider only 2 elems adjacent by a link
2401           continue;
2402         const SMDS_MeshElement* elem2 = (*itLE).second.front();
2403         //if ( FindShape( elem ) != FindShape( elem2 ))
2404         //  continue; // do not fuse triangles laying on different shapes
2405         if ( getAngle( elem, elem2, aNodes[i], aNodes[i+1] ) > theMaxAngle )
2406           continue; // avoid making badly shaped quads
2407         (*itLE).second.push_back( elem );
2408       }
2409       else {
2410         mapLi_listEl[ link ].push_back( elem );
2411       }
2412       mapEl_setLi [ elem ].insert( link );
2413     }
2414   }
2415   // Clean the maps from the links shared by a sole element, ie
2416   // links to which only one element is bound in mapLi_listEl
2417
2418   for ( itLE = mapLi_listEl.begin(); itLE != mapLi_listEl.end(); itLE++ ) {
2419     int nbElems = (*itLE).second.size();
2420     if ( nbElems < 2  ) {
2421       const SMDS_MeshElement* elem = (*itLE).second.front();
2422       SMESH_TLink link = (*itLE).first;
2423       mapEl_setLi[ elem ].erase( link );
2424       if ( mapEl_setLi[ elem ].empty() )
2425         mapEl_setLi.erase( elem );
2426     }
2427   }
2428
2429   // Algo: fuse triangles into quadrangles
2430
2431   while ( ! mapEl_setLi.empty() ) {
2432     // Look for the start element:
2433     // the element having the least nb of shared links
2434     const SMDS_MeshElement* startElem = 0;
2435     int minNbLinks = 4;
2436     for ( itEL = mapEl_setLi.begin(); itEL != mapEl_setLi.end(); itEL++ ) {
2437       int nbLinks = (*itEL).second.size();
2438       if ( nbLinks < minNbLinks ) {
2439         startElem = (*itEL).first;
2440         minNbLinks = nbLinks;
2441         if ( minNbLinks == 1 )
2442           break;
2443       }
2444     }
2445
2446     // search elements to fuse starting from startElem or links of elements
2447     // fused earlyer - startLinks
2448     list< SMESH_TLink > startLinks;
2449     while ( startElem || !startLinks.empty() ) {
2450       while ( !startElem && !startLinks.empty() ) {
2451         // Get an element to start, by a link
2452         SMESH_TLink linkId = startLinks.front();
2453         startLinks.pop_front();
2454         itLE = mapLi_listEl.find( linkId );
2455         if ( itLE != mapLi_listEl.end() ) {
2456           list< const SMDS_MeshElement* > & listElem = (*itLE).second;
2457           list< const SMDS_MeshElement* >::iterator itE = listElem.begin();
2458           for ( ; itE != listElem.end() ; itE++ )
2459             if ( mapEl_setLi.find( (*itE) ) != mapEl_setLi.end() )
2460               startElem = (*itE);
2461           mapLi_listEl.erase( itLE );
2462         }
2463       }
2464
2465       if ( startElem ) {
2466         // Get candidates to be fused
2467         const SMDS_MeshElement *tr1 = startElem, *tr2 = 0, *tr3 = 0;
2468         const SMESH_TLink *link12, *link13;
2469         startElem = 0;
2470         ASSERT( mapEl_setLi.find( tr1 ) != mapEl_setLi.end() );
2471         set< SMESH_TLink >& setLi = mapEl_setLi[ tr1 ];
2472         ASSERT( !setLi.empty() );
2473         set< SMESH_TLink >::iterator itLi;
2474         for ( itLi = setLi.begin(); itLi != setLi.end(); itLi++ )
2475         {
2476           const SMESH_TLink & link = (*itLi);
2477           itLE = mapLi_listEl.find( link );
2478           if ( itLE == mapLi_listEl.end() )
2479             continue;
2480
2481           const SMDS_MeshElement* elem = (*itLE).second.front();
2482           if ( elem == tr1 )
2483             elem = (*itLE).second.back();
2484           mapLi_listEl.erase( itLE );
2485           if ( mapEl_setLi.find( elem ) == mapEl_setLi.end())
2486             continue;
2487           if ( tr2 ) {
2488             tr3 = elem;
2489             link13 = &link;
2490           }
2491           else {
2492             tr2 = elem;
2493             link12 = &link;
2494           }
2495
2496           // add other links of elem to list of links to re-start from
2497           set< SMESH_TLink >& links = mapEl_setLi[ elem ];
2498           set< SMESH_TLink >::iterator it;
2499           for ( it = links.begin(); it != links.end(); it++ ) {
2500             const SMESH_TLink& link2 = (*it);
2501             if ( link2 != link )
2502               startLinks.push_back( link2 );
2503           }
2504         }
2505
2506         // Get nodes of possible quadrangles
2507         const SMDS_MeshNode *n12 [4], *n13 [4];
2508         bool Ok12 = false, Ok13 = false;
2509         const SMDS_MeshNode *linkNode1, *linkNode2;
2510         if(tr2) {
2511           linkNode1 = link12->first;
2512           linkNode2 = link12->second;
2513           if ( tr2 && getQuadrangleNodes( n12, linkNode1, linkNode2, tr1, tr2 ))
2514             Ok12 = true;
2515         }
2516         if(tr3) {
2517           linkNode1 = link13->first;
2518           linkNode2 = link13->second;
2519           if ( tr3 && getQuadrangleNodes( n13, linkNode1, linkNode2, tr1, tr3 ))
2520             Ok13 = true;
2521         }
2522
2523         // Choose a pair to fuse
2524         if ( Ok12 && Ok13 ) {
2525           SMDS_FaceOfNodes quad12 ( n12[ 0 ], n12[ 1 ], n12[ 2 ], n12[ 3 ] );
2526           SMDS_FaceOfNodes quad13 ( n13[ 0 ], n13[ 1 ], n13[ 2 ], n13[ 3 ] );
2527           double aBadRate12 = getBadRate( &quad12, theCrit );
2528           double aBadRate13 = getBadRate( &quad13, theCrit );
2529           if (  aBadRate13 < aBadRate12 )
2530             Ok12 = false;
2531           else
2532             Ok13 = false;
2533         }
2534
2535         // Make quadrangles
2536         // and remove fused elems and removed links from the maps
2537         mapEl_setLi.erase( tr1 );
2538         if ( Ok12 ) {
2539           mapEl_setLi.erase( tr2 );
2540           mapLi_listEl.erase( *link12 );
2541           if(tr1->NbNodes()==3) {
2542             const SMDS_MeshElement* newElem = 0;
2543             newElem = aMesh->AddFace(n12[0], n12[1], n12[2], n12[3] );
2544             myLastCreatedElems.Append(newElem);
2545             AddToSameGroups( newElem, tr1, aMesh );
2546             int aShapeId = tr1->getshapeId();
2547             if ( aShapeId )
2548               {
2549                 aMesh->SetMeshElementOnShape( newElem, aShapeId );
2550               }
2551             aMesh->RemoveElement( tr1 );
2552             aMesh->RemoveElement( tr2 );
2553           }
2554           else {
2555             const SMDS_MeshNode* N1 [6];
2556             const SMDS_MeshNode* N2 [6];
2557             GetNodesFromTwoTria(tr1,tr2,N1,N2);
2558             // now we receive following N1 and N2 (using numeration as above image)
2559             // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
2560             // i.e. first nodes from both arrays determ new diagonal
2561             const SMDS_MeshNode* aNodes[8];
2562             aNodes[0] = N1[0];
2563             aNodes[1] = N1[1];
2564             aNodes[2] = N2[0];
2565             aNodes[3] = N2[1];
2566             aNodes[4] = N1[3];
2567             aNodes[5] = N2[5];
2568             aNodes[6] = N2[3];
2569             aNodes[7] = N1[5];
2570             const SMDS_MeshElement* newElem = 0;
2571             newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3],
2572                                      aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
2573             myLastCreatedElems.Append(newElem);
2574             AddToSameGroups( newElem, tr1, aMesh );
2575             int aShapeId = tr1->getshapeId();
2576             if ( aShapeId )
2577               {
2578                 aMesh->SetMeshElementOnShape( newElem, aShapeId );
2579               }
2580             aMesh->RemoveElement( tr1 );
2581             aMesh->RemoveElement( tr2 );
2582             // remove middle node (9)
2583             GetMeshDS()->RemoveNode( N1[4] );
2584           }
2585         }
2586         else if ( Ok13 ) {
2587           mapEl_setLi.erase( tr3 );
2588           mapLi_listEl.erase( *link13 );
2589           if(tr1->NbNodes()==3) {
2590             const SMDS_MeshElement* newElem = 0;
2591             newElem = aMesh->AddFace(n13[0], n13[1], n13[2], n13[3] );
2592             myLastCreatedElems.Append(newElem);
2593             AddToSameGroups( newElem, tr1, aMesh );
2594             int aShapeId = tr1->getshapeId();
2595             if ( aShapeId )
2596               {
2597                 aMesh->SetMeshElementOnShape( newElem, aShapeId );
2598               }
2599             aMesh->RemoveElement( tr1 );
2600             aMesh->RemoveElement( tr3 );
2601           }
2602           else {
2603             const SMDS_MeshNode* N1 [6];
2604             const SMDS_MeshNode* N2 [6];
2605             GetNodesFromTwoTria(tr1,tr3,N1,N2);
2606             // now we receive following N1 and N2 (using numeration as above image)
2607             // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
2608             // i.e. first nodes from both arrays determ new diagonal
2609             const SMDS_MeshNode* aNodes[8];
2610             aNodes[0] = N1[0];
2611             aNodes[1] = N1[1];
2612             aNodes[2] = N2[0];
2613             aNodes[3] = N2[1];
2614             aNodes[4] = N1[3];
2615             aNodes[5] = N2[5];
2616             aNodes[6] = N2[3];
2617             aNodes[7] = N1[5];
2618             const SMDS_MeshElement* newElem = 0;
2619             newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3],
2620                                      aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
2621             myLastCreatedElems.Append(newElem);
2622             AddToSameGroups( newElem, tr1, aMesh );
2623             int aShapeId = tr1->getshapeId();
2624             if ( aShapeId )
2625               {
2626                 aMesh->SetMeshElementOnShape( newElem, aShapeId );
2627               }
2628             aMesh->RemoveElement( tr1 );
2629             aMesh->RemoveElement( tr3 );
2630             // remove middle node (9)
2631             GetMeshDS()->RemoveNode( N1[4] );
2632           }
2633         }
2634
2635         // Next element to fuse: the rejected one
2636         if ( tr3 )
2637           startElem = Ok12 ? tr3 : tr2;
2638
2639       } // if ( startElem )
2640     } // while ( startElem || !startLinks.empty() )
2641   } // while ( ! mapEl_setLi.empty() )
2642
2643   return true;
2644 }
2645
2646
2647 /*#define DUMPSO(txt) \
2648 //  cout << txt << endl;
2649 //=============================================================================
2650 //
2651 //
2652 //
2653 //=============================================================================
2654 static void swap( int i1, int i2, int idNodes[], gp_Pnt P[] )
2655 {
2656 if ( i1 == i2 )
2657 return;
2658 int tmp = idNodes[ i1 ];
2659 idNodes[ i1 ] = idNodes[ i2 ];
2660 idNodes[ i2 ] = tmp;
2661 gp_Pnt Ptmp = P[ i1 ];
2662 P[ i1 ] = P[ i2 ];
2663 P[ i2 ] = Ptmp;
2664 DUMPSO( i1 << "(" << idNodes[ i2 ] << ") <-> " << i2 << "(" << idNodes[ i1 ] << ")");
2665 }
2666
2667 //=======================================================================
2668 //function : SortQuadNodes
2669 //purpose  : Set 4 nodes of a quadrangle face in a good order.
2670 //           Swap 1<->2 or 2<->3 nodes and correspondingly return
2671 //           1 or 2 else 0.
2672 //=======================================================================
2673
2674 int SMESH_MeshEditor::SortQuadNodes (const SMDS_Mesh * theMesh,
2675 int               idNodes[] )
2676 {
2677   gp_Pnt P[4];
2678   int i;
2679   for ( i = 0; i < 4; i++ ) {
2680     const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
2681     if ( !n ) return 0;
2682     P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
2683   }
2684
2685   gp_Vec V1(P[0], P[1]);
2686   gp_Vec V2(P[0], P[2]);
2687   gp_Vec V3(P[0], P[3]);
2688
2689   gp_Vec Cross1 = V1 ^ V2;
2690   gp_Vec Cross2 = V2 ^ V3;
2691
2692   i = 0;
2693   if (Cross1.Dot(Cross2) < 0)
2694   {
2695     Cross1 = V2 ^ V1;
2696     Cross2 = V1 ^ V3;
2697
2698     if (Cross1.Dot(Cross2) < 0)
2699       i = 2;
2700     else
2701       i = 1;
2702     swap ( i, i + 1, idNodes, P );
2703
2704     //     for ( int ii = 0; ii < 4; ii++ ) {
2705     //       const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
2706     //       DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
2707     //     }
2708   }
2709   return i;
2710 }
2711
2712 //=======================================================================
2713 //function : SortHexaNodes
2714 //purpose  : Set 8 nodes of a hexahedron in a good order.
2715 //           Return success status
2716 //=======================================================================
2717
2718 bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh,
2719                                       int               idNodes[] )
2720 {
2721   gp_Pnt P[8];
2722   int i;
2723   DUMPSO( "INPUT: ========================================");
2724   for ( i = 0; i < 8; i++ ) {
2725     const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
2726     if ( !n ) return false;
2727     P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
2728     DUMPSO( i << "(" << idNodes[i] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
2729   }
2730   DUMPSO( "========================================");
2731
2732
2733   set<int> faceNodes;  // ids of bottom face nodes, to be found
2734   set<int> checkedId1; // ids of tried 2-nd nodes
2735   Standard_Real leastDist = DBL_MAX; // dist of the 4-th node from 123 plane
2736   const Standard_Real tol = 1.e-6;   // tolerance to find nodes in plane
2737   int iMin, iLoop1 = 0;
2738
2739   // Loop to try the 2-nd nodes
2740
2741   while ( leastDist > DBL_MIN && ++iLoop1 < 8 )
2742   {
2743     // Find not checked 2-nd node
2744     for ( i = 1; i < 8; i++ )
2745       if ( checkedId1.find( idNodes[i] ) == checkedId1.end() ) {
2746         int id1 = idNodes[i];
2747         swap ( 1, i, idNodes, P );
2748         checkedId1.insert ( id1 );
2749         break;
2750       }
2751
2752     // Find the 3-d node so that 1-2-3 triangle to be on a hexa face,
2753     // ie that all but meybe one (id3 which is on the same face) nodes
2754     // lay on the same side from the triangle plane.
2755
2756     bool manyInPlane = false; // more than 4 nodes lay in plane
2757     int iLoop2 = 0;
2758     while ( ++iLoop2 < 6 ) {
2759
2760       // get 1-2-3 plane coeffs
2761       Standard_Real A, B, C, D;
2762       gp_Vec N = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
2763       if ( N.SquareMagnitude() > gp::Resolution() )
2764       {
2765         gp_Pln pln ( P[0], N );
2766         pln.Coefficients( A, B, C, D );
2767
2768         // find the node (iMin) closest to pln
2769         Standard_Real dist[ 8 ], minDist = DBL_MAX;
2770         set<int> idInPln;
2771         for ( i = 3; i < 8; i++ ) {
2772           dist[i] = A * P[i].X() + B * P[i].Y() + C * P[i].Z() + D;
2773           if ( fabs( dist[i] ) < minDist ) {
2774             minDist = fabs( dist[i] );
2775             iMin = i;
2776           }
2777           if ( fabs( dist[i] ) <= tol )
2778             idInPln.insert( idNodes[i] );
2779         }
2780
2781         // there should not be more than 4 nodes in bottom plane
2782         if ( idInPln.size() > 1 )
2783         {
2784           DUMPSO( "### idInPln.size() = " << idInPln.size());
2785           // idInPlane does not contain the first 3 nodes
2786           if ( manyInPlane || idInPln.size() == 5)
2787             return false; // all nodes in one plane
2788           manyInPlane = true;
2789
2790           // set the 1-st node to be not in plane
2791           for ( i = 3; i < 8; i++ ) {
2792             if ( idInPln.find( idNodes[ i ] ) == idInPln.end() ) {
2793               DUMPSO( "### Reset 0-th node");
2794               swap( 0, i, idNodes, P );
2795               break;
2796             }
2797           }
2798
2799           // reset to re-check second nodes
2800           leastDist = DBL_MAX;
2801           faceNodes.clear();
2802           checkedId1.clear();
2803           iLoop1 = 0;
2804           break; // from iLoop2;
2805         }
2806
2807         // check that the other 4 nodes are on the same side
2808         bool sameSide = true;
2809         bool isNeg = dist[ iMin == 3 ? 4 : 3 ] <= 0.;
2810         for ( i = 3; sameSide && i < 8; i++ ) {
2811           if ( i != iMin )
2812             sameSide = ( isNeg == dist[i] <= 0.);
2813         }
2814
2815         // keep best solution
2816         if ( sameSide && minDist < leastDist ) {
2817           leastDist = minDist;
2818           faceNodes.clear();
2819           faceNodes.insert( idNodes[ 1 ] );
2820           faceNodes.insert( idNodes[ 2 ] );
2821           faceNodes.insert( idNodes[ iMin ] );
2822           DUMPSO( "loop " << iLoop2 << " id2 " << idNodes[ 1 ] << " id3 " << idNodes[ 2 ]
2823                   << " leastDist = " << leastDist);
2824           if ( leastDist <= DBL_MIN )
2825             break;
2826         }
2827       }
2828
2829       // set next 3-d node to check
2830       int iNext = 2 + iLoop2;
2831       if ( iNext < 8 ) {
2832         DUMPSO( "Try 2-nd");
2833         swap ( 2, iNext, idNodes, P );
2834       }
2835     } // while ( iLoop2 < 6 )
2836   } // iLoop1
2837
2838   if ( faceNodes.empty() ) return false;
2839
2840   // Put the faceNodes in proper places
2841   for ( i = 4; i < 8; i++ ) {
2842     if ( faceNodes.find( idNodes[ i ] ) != faceNodes.end() ) {
2843       // find a place to put
2844       int iTo = 1;
2845       while ( faceNodes.find( idNodes[ iTo ] ) != faceNodes.end() )
2846         iTo++;
2847       DUMPSO( "Set faceNodes");
2848       swap ( iTo, i, idNodes, P );
2849     }
2850   }
2851
2852
2853   // Set nodes of the found bottom face in good order
2854   DUMPSO( " Found bottom face: ");
2855   i = SortQuadNodes( theMesh, idNodes );
2856   if ( i ) {
2857     gp_Pnt Ptmp = P[ i ];
2858     P[ i ] = P[ i+1 ];
2859     P[ i+1 ] = Ptmp;
2860   }
2861   //   else
2862   //     for ( int ii = 0; ii < 4; ii++ ) {
2863   //       const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
2864   //       DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
2865   //    }
2866
2867   // Gravity center of the top and bottom faces
2868   gp_Pnt aGCb = ( P[0].XYZ() + P[1].XYZ() + P[2].XYZ() + P[3].XYZ() ) / 4.;
2869   gp_Pnt aGCt = ( P[4].XYZ() + P[5].XYZ() + P[6].XYZ() + P[7].XYZ() ) / 4.;
2870
2871   // Get direction from the bottom to the top face
2872   gp_Vec upDir ( aGCb, aGCt );
2873   Standard_Real upDirSize = upDir.Magnitude();
2874   if ( upDirSize <= gp::Resolution() ) return false;
2875   upDir / upDirSize;
2876
2877   // Assure that the bottom face normal points up
2878   gp_Vec Nb = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
2879   Nb += gp_Vec (P[0], P[2]).Crossed( gp_Vec (P[0], P[3]) );
2880   if ( Nb.Dot( upDir ) < 0 ) {
2881     DUMPSO( "Reverse bottom face");
2882     swap( 1, 3, idNodes, P );
2883   }
2884
2885   // Find 5-th node - the one closest to the 1-st among the last 4 nodes.
2886   Standard_Real minDist = DBL_MAX;
2887   for ( i = 4; i < 8; i++ ) {
2888     // projection of P[i] to the plane defined by P[0] and upDir
2889     gp_Pnt Pp = P[i].Translated( upDir * ( upDir.Dot( gp_Vec( P[i], P[0] ))));
2890     Standard_Real sqDist = P[0].SquareDistance( Pp );
2891     if ( sqDist < minDist ) {
2892       minDist = sqDist;
2893       iMin = i;
2894     }
2895   }
2896   DUMPSO( "Set 4-th");
2897   swap ( 4, iMin, idNodes, P );
2898
2899   // Set nodes of the top face in good order
2900   DUMPSO( "Sort top face");
2901   i = SortQuadNodes( theMesh, &idNodes[4] );
2902   if ( i ) {
2903     i += 4;
2904     gp_Pnt Ptmp = P[ i ];
2905     P[ i ] = P[ i+1 ];
2906     P[ i+1 ] = Ptmp;
2907   }
2908
2909   // Assure that direction of the top face normal is from the bottom face
2910   gp_Vec Nt = gp_Vec (P[4], P[5]).Crossed( gp_Vec (P[4], P[6]) );
2911   Nt += gp_Vec (P[4], P[6]).Crossed( gp_Vec (P[4], P[7]) );
2912   if ( Nt.Dot( upDir ) < 0 ) {
2913     DUMPSO( "Reverse top face");
2914     swap( 5, 7, idNodes, P );
2915   }
2916
2917   //   DUMPSO( "OUTPUT: ========================================");
2918   //   for ( i = 0; i < 8; i++ ) {
2919   //     float *p = ugrid->GetPoint(idNodes[i]);
2920   //     DUMPSO( i << "(" << idNodes[i] << ") : " << p[0] << " " << p[1] << " " << p[2]);
2921   //   }
2922
2923   return true;
2924 }*/
2925
2926 //================================================================================
2927 /*!
2928  * \brief Return nodes linked to the given one
2929  * \param theNode - the node
2930  * \param linkedNodes - the found nodes
2931  * \param type - the type of elements to check
2932  *
2933  * Medium nodes are ignored
2934  */
2935 //================================================================================
2936
2937 void SMESH_MeshEditor::GetLinkedNodes( const SMDS_MeshNode* theNode,
2938                                        TIDSortedElemSet &   linkedNodes,
2939                                        SMDSAbs_ElementType  type )
2940 {
2941   SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(type);
2942   while ( elemIt->more() )
2943   {
2944     const SMDS_MeshElement* elem = elemIt->next();
2945     if(elem->GetType() == SMDSAbs_0DElement)
2946       continue;
2947
2948     SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
2949     if ( elem->GetType() == SMDSAbs_Volume )
2950     {
2951       SMDS_VolumeTool vol( elem );
2952       while ( nodeIt->more() ) {
2953         const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
2954         if ( theNode != n && vol.IsLinked( theNode, n ))
2955           linkedNodes.insert( n );
2956       }
2957     }
2958     else
2959     {
2960       for ( int i = 0; nodeIt->more(); ++i ) {
2961         const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
2962         if ( n == theNode ) {
2963           int iBefore = i - 1;
2964           int iAfter  = i + 1;
2965           if ( elem->IsQuadratic() ) {
2966             int nb = elem->NbNodes() / 2;
2967             iAfter  = SMESH_MesherHelper::WrapIndex( iAfter, nb );
2968             iBefore = SMESH_MesherHelper::WrapIndex( iBefore, nb );
2969           }
2970           linkedNodes.insert( elem->GetNodeWrap( iAfter ));
2971           linkedNodes.insert( elem->GetNodeWrap( iBefore ));
2972         }
2973       }
2974     }
2975   }
2976 }
2977
2978 //=======================================================================
2979 //function : laplacianSmooth
2980 //purpose  : pulls theNode toward the center of surrounding nodes directly
2981 //           connected to that node along an element edge
2982 //=======================================================================
2983
2984 void laplacianSmooth(const SMDS_MeshNode*                 theNode,
2985                      const Handle(Geom_Surface)&          theSurface,
2986                      map< const SMDS_MeshNode*, gp_XY* >& theUVMap)
2987 {
2988   // find surrounding nodes
2989
2990   TIDSortedElemSet nodeSet;
2991   SMESH_MeshEditor::GetLinkedNodes( theNode, nodeSet, SMDSAbs_Face );
2992
2993   // compute new coodrs
2994
2995   double coord[] = { 0., 0., 0. };
2996   TIDSortedElemSet::iterator nodeSetIt = nodeSet.begin();
2997   for ( ; nodeSetIt != nodeSet.end(); nodeSetIt++ ) {
2998     const SMDS_MeshNode* node = cast2Node(*nodeSetIt);
2999     if ( theSurface.IsNull() ) { // smooth in 3D
3000       coord[0] += node->X();
3001       coord[1] += node->Y();
3002       coord[2] += node->Z();
3003     }