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