Salome HOME
1d5698d249bc67e59ed207f38ab86ebc735811ca
[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 class TElemToDelete
1903 {
1904 public:
1905   TElemToDelete(const SMDS_MeshElement* theElem, SMESHDS_SubMesh* theSubMesh)
1906   {
1907     elem = theElem;
1908     subMesh = theSubMesh;
1909   }
1910   const SMDS_MeshElement* Elem() const {return elem;}
1911   SMESHDS_SubMesh* Submesh() {return subMesh;}
1912   const SMDS_MeshElement* elem;
1913   SMESHDS_SubMesh* subMesh;
1914 };
1915
1916 //=======================================================================
1917 //function : SplitVolumesIntoTetra
1918 //purpose  : Split volume elements into tetrahedra.
1919 //=======================================================================
1920
1921 void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems,
1922                                               const int                theMethodFlags)
1923 {
1924   // std-like iterator on coordinates of nodes of mesh element
1925   typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > NXyzIterator;
1926   NXyzIterator xyzEnd;
1927
1928   SMDS_VolumeTool    volTool;
1929   SMESH_MesherHelper helper( *GetMesh());
1930
1931   SMESHDS_SubMesh* subMesh = 0;//GetMeshDS()->MeshElements(1);
1932   SMESHDS_SubMesh* fSubMesh = 0;//subMesh;
1933
1934   SMESH_SequenceOfElemPtr newNodes, newElems;
1935
1936   // map face of volume to it's baricenrtic node
1937   map< TVolumeFaceKey, const SMDS_MeshNode* > volFace2BaryNode;
1938   double bc[3];
1939
1940   TIDSortedElemSet::const_iterator elem = theElems.begin();
1941   std::vector<TElemToDelete> elem_to_delete;
1942   for ( ; elem != theElems.end(); ++elem )
1943   {
1944     if ( (*elem)->GetType() != SMDSAbs_Volume )
1945       continue;
1946     SMDSAbs_EntityType geomType = (*elem)->GetEntityType();
1947     if ( geomType == SMDSEntity_Tetra || geomType == SMDSEntity_Quad_Tetra )
1948       continue;
1949
1950     if ( !volTool.Set( *elem, /*ignoreCentralNodes=*/false )) continue; // strange...
1951
1952     TSplitMethod splitMethod = getSplitMethod( volTool, theMethodFlags );
1953     if ( splitMethod._nbTetra < 1 ) continue;
1954
1955     // find submesh to add new tetras to
1956     if ( !subMesh || !subMesh->Contains( *elem ))
1957     {
1958       int shapeID = FindShape( *elem );
1959       helper.SetSubShape( shapeID ); // helper will add tetras to the found submesh
1960       subMesh = GetMeshDS()->MeshElements( shapeID );
1961     }
1962     int iQ;
1963     if ( (*elem)->IsQuadratic() )
1964     {
1965       iQ = 2;
1966       // add quadratic links to the helper
1967       for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
1968       {
1969         const SMDS_MeshNode** fNodes = volTool.GetFaceNodes( iF );
1970         int nbN = volTool.NbFaceNodes( iF ) - bool( volTool.GetCenterNodeIndex(iF) > 0 );
1971         for ( int iN = 0; iN < nbN; iN += iQ )
1972           helper.AddTLinkNode( fNodes[iN], fNodes[iN+2], fNodes[iN+1] );
1973       }
1974       helper.SetIsQuadratic( true );
1975     }
1976     else
1977     {
1978       iQ = 1;
1979       helper.SetIsQuadratic( false );
1980     }
1981     vector<const SMDS_MeshNode*> nodes( (*elem)->begin_nodes(), (*elem)->end_nodes() );
1982     helper.SetElementsOnShape( true );
1983     if ( splitMethod._baryNode )
1984     {
1985       // make a node at barycenter
1986       volTool.GetBaryCenter( bc[0], bc[1], bc[2] );
1987       SMDS_MeshNode* gcNode = helper.AddNode( bc[0], bc[1], bc[2] );
1988       nodes.push_back( gcNode );
1989       newNodes.Append( gcNode );
1990     }
1991     if ( !splitMethod._faceBaryNode.empty() )
1992     {
1993       // make or find baricentric nodes of faces
1994       map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.begin();
1995       for ( ; iF_n != splitMethod._faceBaryNode.end(); ++iF_n )
1996       {
1997         map< TVolumeFaceKey, const SMDS_MeshNode* >::iterator f_n =
1998           volFace2BaryNode.insert
1999           ( make_pair( TVolumeFaceKey( volTool,iF_n->first ), iF_n->second )).first;
2000         if ( !f_n->second )
2001         {
2002           volTool.GetFaceBaryCenter( iF_n->first, bc[0], bc[1], bc[2] );
2003           newNodes.Append( f_n->second = helper.AddNode( bc[0], bc[1], bc[2] ));
2004         }
2005         nodes.push_back( iF_n->second = f_n->second );
2006       }
2007     }
2008
2009     // make tetras
2010     vector<const SMDS_MeshElement* > tetras( splitMethod._nbTetra ); // splits of a volume
2011     const int* tetConn = splitMethod._connectivity;
2012     for ( int i = 0; i < splitMethod._nbTetra; ++i, tetConn += 4 )
2013       newElems.Append( tetras[ i ] = helper.AddVolume( nodes[ tetConn[0] ],
2014                                                        nodes[ tetConn[1] ],
2015                                                        nodes[ tetConn[2] ],
2016                                                        nodes[ tetConn[3] ]));
2017
2018     ReplaceElemInGroups( *elem, tetras, GetMeshDS() );
2019
2020     // Split faces on sides of the split volume
2021
2022     const SMDS_MeshNode** volNodes = volTool.GetNodes();
2023     for ( int iF = 0; iF < volTool.NbFaces(); ++iF )
2024     {
2025       const int nbNodes = volTool.NbFaceNodes( iF ) / iQ;
2026       if ( nbNodes < 4 ) continue;
2027
2028       // find an existing face
2029       vector<const SMDS_MeshNode*> fNodes( volTool.GetFaceNodes( iF ),
2030                                            volTool.GetFaceNodes( iF ) + volTool.NbFaceNodes( iF ));
2031       while ( const SMDS_MeshElement* face = GetMeshDS()->FindElement( fNodes, SMDSAbs_Face,
2032                                                                        /*noMedium=*/false))
2033       {
2034         // make triangles
2035         helper.SetElementsOnShape( false );
2036         vector< const SMDS_MeshElement* > triangles;
2037
2038         // find submesh to add new triangles in
2039         if ( !fSubMesh || !fSubMesh->Contains( face ))
2040         {
2041           int shapeID = FindShape( face );
2042           fSubMesh = GetMeshDS()->MeshElements( shapeID );
2043         }
2044         map<int, const SMDS_MeshNode*>::iterator iF_n = splitMethod._faceBaryNode.find(iF);
2045         if ( iF_n != splitMethod._faceBaryNode.end() )
2046         {
2047           for ( int iN = 0; iN < nbNodes*iQ; iN += iQ )
2048           {
2049             const SMDS_MeshNode* n1 = fNodes[iN];
2050             const SMDS_MeshNode *n2 = fNodes[(iN+iQ)%(nbNodes*iQ)];
2051             const SMDS_MeshNode *n3 = iF_n->second;
2052             if ( !volTool.IsFaceExternal( iF ))
2053               swap( n2, n3 );
2054             triangles.push_back( helper.AddFace( n1,n2,n3 ));
2055
2056             if ( fSubMesh && n3->getshapeId() < 1 )
2057               fSubMesh->AddNode( n3 );
2058           }
2059         }
2060         else
2061         {
2062           // among possible triangles create ones discribed by split method
2063           const int* nInd = volTool.GetFaceNodesIndices( iF );
2064           int nbVariants = ( nbNodes == 4 ? 2 : nbNodes );
2065           int iCom = 0; // common node of triangle faces to split into
2066           list< TTriangleFacet > facets;
2067           for ( int iVar = 0; iVar < nbVariants; ++iVar, ++iCom )
2068           {
2069             TTriangleFacet t012( nInd[ iQ * ( iCom                )],
2070                                  nInd[ iQ * ( (iCom+1)%nbNodes )],
2071                                  nInd[ iQ * ( (iCom+2)%nbNodes )]);
2072             TTriangleFacet t023( nInd[ iQ * ( iCom                )],
2073                                  nInd[ iQ * ( (iCom+2)%nbNodes )],
2074                                  nInd[ iQ * ( (iCom+3)%nbNodes )]);
2075             if ( splitMethod.hasFacet( t012 ) && splitMethod.hasFacet( t023 ))
2076             {
2077               facets.push_back( t012 );
2078               facets.push_back( t023 );
2079               for ( int iLast = iCom+4; iLast < iCom+nbNodes; ++iLast )
2080                 facets.push_back( TTriangleFacet( nInd[ iQ * ( iCom             )],
2081                                                   nInd[ iQ * ((iLast-1)%nbNodes )],
2082                                                   nInd[ iQ * ((iLast  )%nbNodes )]));
2083               break;
2084             }
2085           }
2086           list< TTriangleFacet >::iterator facet = facets.begin();
2087           for ( ; facet != facets.end(); ++facet )
2088           {
2089             if ( !volTool.IsFaceExternal( iF ))
2090               swap( facet->_n2, facet->_n3 );
2091             triangles.push_back( helper.AddFace( volNodes[ facet->_n1 ],
2092                                                  volNodes[ facet->_n2 ],
2093                                                  volNodes[ facet->_n3 ]));
2094           }
2095         }
2096         for ( int i = 0; i < triangles.size(); ++i )
2097         {
2098           if ( !triangles[i] ) continue;
2099           if ( fSubMesh )
2100             fSubMesh->AddElement( triangles[i]);
2101           newElems.Append( triangles[i] );
2102         }
2103         ReplaceElemInGroups( face, triangles, GetMeshDS() );
2104         GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false );
2105 //         TElemToDelete faceToDelete(face, fSubMesh);
2106 //         elem_to_delete.push_back(faceToDelete); 
2107       }
2108
2109     } // loop on volume faces to split them into triangles
2110
2111 //     GetMeshDS()->RemoveFreeElement( *elem, subMesh, /*fromGroups=*/false );
2112
2113     // rnc : don't delete the elem here because it results in a mesh with a free
2114     // ID at the beginning of the ID list. The first tetra is then inserted in O(1)
2115     // but the second one is inserted in O(n), then the whole procedure has almost a O(n^2)
2116     // complexity. If all elements to remove are stored and removed after tetra creation 
2117     // we get a O(n) complexity for the whole procedure. 
2118     // The memory cost is at worst a 6*n*constant memory occupation (where n is the number of elements) 
2119     // before deletion of the hexas and then 5*n*constant instead of a maximum of 5*n*constant.
2120     // So there is a transient 1/5*(memory occupation) additional cost.
2121
2122     // Store the elements to delete
2123     TElemToDelete elemToDelete(*elem, subMesh);
2124     elem_to_delete.push_back(elemToDelete);
2125
2126     if ( geomType == SMDSEntity_TriQuad_Hexa )
2127     {
2128       // remove medium nodes that could become free
2129       for ( int i = 20; i < volTool.NbNodes(); ++i )
2130         if ( volNodes[i]->NbInverseElements() == 0 )
2131           GetMeshDS()->RemoveNode( volNodes[i] );
2132     }
2133   } // loop on volumes to split
2134
2135   // Delete stored elements
2136   std::vector<TElemToDelete>::iterator it;
2137   for( it = elem_to_delete.begin(); it!= elem_to_delete.end(); it++)
2138   {
2139     GetMeshDS()->RemoveFreeElement( it->Elem(), it->Submesh(), /*fromGroups=*/false );
2140   }
2141   
2142   myLastCreatedNodes = newNodes;
2143   myLastCreatedElems = newElems;
2144 }
2145
2146 //=======================================================================
2147 //function : AddToSameGroups
2148 //purpose  : add elemToAdd to the groups the elemInGroups belongs to
2149 //=======================================================================
2150
2151 void SMESH_MeshEditor::AddToSameGroups (const SMDS_MeshElement* elemToAdd,
2152                                         const SMDS_MeshElement* elemInGroups,
2153                                         SMESHDS_Mesh *          aMesh)
2154 {
2155   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2156   if (!groups.empty()) {
2157     set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2158     for ( ; grIt != groups.end(); grIt++ ) {
2159       SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2160       if ( group && group->Contains( elemInGroups ))
2161         group->SMDSGroup().Add( elemToAdd );
2162     }
2163   }
2164 }
2165
2166
2167 //=======================================================================
2168 //function : RemoveElemFromGroups
2169 //purpose  : Remove removeelem to the groups the elemInGroups belongs to
2170 //=======================================================================
2171 void SMESH_MeshEditor::RemoveElemFromGroups (const SMDS_MeshElement* removeelem,
2172                                              SMESHDS_Mesh *          aMesh)
2173 {
2174   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2175   if (!groups.empty())
2176   {
2177     set<SMESHDS_GroupBase*>::const_iterator GrIt = groups.begin();
2178     for (; GrIt != groups.end(); GrIt++)
2179     {
2180       SMESHDS_Group* grp = dynamic_cast<SMESHDS_Group*>(*GrIt);
2181       if (!grp || grp->IsEmpty()) continue;
2182       grp->SMDSGroup().Remove(removeelem);
2183     }
2184   }
2185 }
2186
2187 //================================================================================
2188 /*!
2189  * \brief Replace elemToRm by elemToAdd in the all groups
2190  */
2191 //================================================================================
2192
2193 void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement* elemToRm,
2194                                             const SMDS_MeshElement* elemToAdd,
2195                                             SMESHDS_Mesh *          aMesh)
2196 {
2197   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2198   if (!groups.empty()) {
2199     set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2200     for ( ; grIt != groups.end(); grIt++ ) {
2201       SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2202       if ( group && group->SMDSGroup().Remove( elemToRm ) && elemToAdd )
2203         group->SMDSGroup().Add( elemToAdd );
2204     }
2205   }
2206 }
2207
2208 //================================================================================
2209 /*!
2210  * \brief Replace elemToRm by elemToAdd in the all groups
2211  */
2212 //================================================================================
2213
2214 void SMESH_MeshEditor::ReplaceElemInGroups (const SMDS_MeshElement*                elemToRm,
2215                                             const vector<const SMDS_MeshElement*>& elemToAdd,
2216                                             SMESHDS_Mesh *                         aMesh)
2217 {
2218   const set<SMESHDS_GroupBase*>& groups = aMesh->GetGroups();
2219   if (!groups.empty())
2220   {
2221     set<SMESHDS_GroupBase*>::const_iterator grIt = groups.begin();
2222     for ( ; grIt != groups.end(); grIt++ ) {
2223       SMESHDS_Group* group = dynamic_cast<SMESHDS_Group*>( *grIt );
2224       if ( group && group->SMDSGroup().Remove( elemToRm ) )
2225         for ( int i = 0; i < elemToAdd.size(); ++i )
2226           group->SMDSGroup().Add( elemToAdd[ i ] );
2227     }
2228   }
2229 }
2230
2231 //=======================================================================
2232 //function : QuadToTri
2233 //purpose  : Cut quadrangles into triangles.
2234 //           theCrit is used to select a diagonal to cut
2235 //=======================================================================
2236
2237 bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems,
2238                                   const bool         the13Diag)
2239 {
2240   myLastCreatedElems.Clear();
2241   myLastCreatedNodes.Clear();
2242
2243   MESSAGE( "::QuadToTri()" );
2244
2245   SMESHDS_Mesh * aMesh = GetMeshDS();
2246
2247   Handle(Geom_Surface) surface;
2248   SMESH_MesherHelper   helper( *GetMesh() );
2249
2250   TIDSortedElemSet::iterator itElem;
2251   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
2252     const SMDS_MeshElement* elem = *itElem;
2253     if ( !elem || elem->GetType() != SMDSAbs_Face )
2254       continue;
2255     bool isquad = elem->NbNodes()==4 || elem->NbNodes()==8;
2256     if(!isquad) continue;
2257
2258     if(elem->NbNodes()==4) {
2259       // retrieve element nodes
2260       const SMDS_MeshNode* aNodes [4];
2261       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2262       int i = 0;
2263       while ( itN->more() )
2264         aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
2265
2266       int aShapeId = FindShape( elem );
2267       const SMDS_MeshElement* newElem1 = 0;
2268       const SMDS_MeshElement* newElem2 = 0;
2269       if ( the13Diag ) {
2270         newElem1 = aMesh->AddFace( aNodes[2], aNodes[0], aNodes[1] );
2271         newElem2 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] );
2272       }
2273       else {
2274         newElem1 = aMesh->AddFace( aNodes[3], aNodes[0], aNodes[1] );
2275         newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] );
2276       }
2277       myLastCreatedElems.Append(newElem1);
2278       myLastCreatedElems.Append(newElem2);
2279       // put a new triangle on the same shape and add to the same groups
2280       if ( aShapeId )
2281         {
2282           aMesh->SetMeshElementOnShape( newElem1, aShapeId );
2283           aMesh->SetMeshElementOnShape( newElem2, aShapeId );
2284         }
2285       AddToSameGroups( newElem1, elem, aMesh );
2286       AddToSameGroups( newElem2, elem, aMesh );
2287       //aMesh->RemoveFreeElement(elem, aMesh->MeshElements(aShapeId), true);
2288       aMesh->RemoveElement( elem );
2289     }
2290
2291     // Quadratic quadrangle
2292
2293     if( elem->NbNodes()==8 && elem->IsQuadratic() ) {
2294
2295       // get surface elem is on
2296       int aShapeId = FindShape( elem );
2297       if ( aShapeId != helper.GetSubShapeID() ) {
2298         surface.Nullify();
2299         TopoDS_Shape shape;
2300         if ( aShapeId > 0 )
2301           shape = aMesh->IndexToShape( aShapeId );
2302         if ( !shape.IsNull() && shape.ShapeType() == TopAbs_FACE ) {
2303           TopoDS_Face face = TopoDS::Face( shape );
2304           surface = BRep_Tool::Surface( face );
2305           if ( !surface.IsNull() )
2306             helper.SetSubShape( shape );
2307         }
2308       }
2309
2310       const SMDS_MeshNode* aNodes [8];
2311       const SMDS_MeshNode* inFaceNode = 0;
2312       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
2313       int i = 0;
2314       while ( itN->more() ) {
2315         aNodes[ i++ ] = static_cast<const SMDS_MeshNode*>( itN->next() );
2316         if ( !inFaceNode && helper.GetNodeUVneedInFaceNode() &&
2317              aNodes[ i-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
2318         {
2319           inFaceNode = aNodes[ i-1 ];
2320         }
2321       }
2322
2323       // find middle point for (0,1,2,3)
2324       // and create a node in this point;
2325       gp_XYZ p( 0,0,0 );
2326       if ( surface.IsNull() ) {
2327         for(i=0; i<4; i++)
2328           p += gp_XYZ(aNodes[i]->X(), aNodes[i]->Y(), aNodes[i]->Z() );
2329         p /= 4;
2330       }
2331       else {
2332         TopoDS_Face geomFace = TopoDS::Face( helper.GetSubShape() );
2333         gp_XY uv( 0,0 );
2334         for(i=0; i<4; i++)
2335           uv += helper.GetNodeUV( geomFace, aNodes[i], inFaceNode );
2336         uv /= 4.;
2337         p = surface->Value( uv.X(), uv.Y() ).XYZ();
2338       }
2339       const SMDS_MeshNode* newN = aMesh->AddNode( p.X(), p.Y(), p.Z() );
2340       myLastCreatedNodes.Append(newN);
2341
2342       // create a new element
2343       const SMDS_MeshElement* newElem1 = 0;
2344       const SMDS_MeshElement* newElem2 = 0;
2345       if ( the13Diag ) {
2346         newElem1 = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0],
2347                                   aNodes[6], aNodes[7], newN );
2348         newElem2 = aMesh->AddFace(aNodes[2], aNodes[0], aNodes[1],
2349                                   newN,      aNodes[4], aNodes[5] );
2350       }
2351       else {
2352         newElem1 = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1],
2353                                   aNodes[7], aNodes[4], newN );
2354         newElem2 = aMesh->AddFace(aNodes[3], aNodes[1], aNodes[2],
2355                                   newN,      aNodes[5], aNodes[6] );
2356       }
2357       myLastCreatedElems.Append(newElem1);
2358       myLastCreatedElems.Append(newElem2);
2359       // put a new triangle on the same shape and add to the same groups
2360       if ( aShapeId )
2361         {
2362           aMesh->SetMeshElementOnShape( newElem1, aShapeId );
2363           aMesh->SetMeshElementOnShape( newElem2, aShapeId );
2364         }
2365       AddToSameGroups( newElem1, elem, aMesh );
2366       AddToSameGroups( newElem2, elem, aMesh );
2367       aMesh->RemoveElement( elem );
2368     }
2369   }
2370
2371   return true;
2372 }
2373
2374 //=======================================================================
2375 //function : getAngle
2376 //purpose  :
2377 //=======================================================================
2378
2379 double getAngle(const SMDS_MeshElement * tr1,
2380                 const SMDS_MeshElement * tr2,
2381                 const SMDS_MeshNode *    n1,
2382                 const SMDS_MeshNode *    n2)
2383 {
2384   double angle = 2. * M_PI; // bad angle
2385
2386   // get normals
2387   SMESH::Controls::TSequenceOfXYZ P1, P2;
2388   if ( !SMESH::Controls::NumericalFunctor::GetPoints( tr1, P1 ) ||
2389        !SMESH::Controls::NumericalFunctor::GetPoints( tr2, P2 ))
2390     return angle;
2391   gp_Vec N1,N2;
2392   if(!tr1->IsQuadratic())
2393     N1 = gp_Vec( P1(2) - P1(1) ) ^ gp_Vec( P1(3) - P1(1) );
2394   else
2395     N1 = gp_Vec( P1(3) - P1(1) ) ^ gp_Vec( P1(5) - P1(1) );
2396   if ( N1.SquareMagnitude() <= gp::Resolution() )
2397     return angle;
2398   if(!tr2->IsQuadratic())
2399     N2 = gp_Vec( P2(2) - P2(1) ) ^ gp_Vec( P2(3) - P2(1) );
2400   else
2401     N2 = gp_Vec( P2(3) - P2(1) ) ^ gp_Vec( P2(5) - P2(1) );
2402   if ( N2.SquareMagnitude() <= gp::Resolution() )
2403     return angle;
2404
2405   // find the first diagonal node n1 in the triangles:
2406   // take in account a diagonal link orientation
2407   const SMDS_MeshElement *nFirst[2], *tr[] = { tr1, tr2 };
2408   for ( int t = 0; t < 2; t++ ) {
2409     SMDS_ElemIteratorPtr it = tr[ t ]->nodesIterator();
2410     int i = 0, iDiag = -1;
2411     while ( it->more()) {
2412       const SMDS_MeshElement *n = it->next();
2413       if ( n == n1 || n == n2 ) {
2414         if ( iDiag < 0)
2415           iDiag = i;
2416         else {
2417           if ( i - iDiag == 1 )
2418             nFirst[ t ] = ( n == n1 ? n2 : n1 );
2419           else
2420             nFirst[ t ] = n;
2421           break;
2422         }
2423       }
2424       i++;
2425     }
2426   }
2427   if ( nFirst[ 0 ] == nFirst[ 1 ] )
2428     N2.Reverse();
2429
2430   angle = N1.Angle( N2 );
2431   //SCRUTE( angle );
2432   return angle;
2433 }
2434
2435 // =================================================
2436 // class generating a unique ID for a pair of nodes
2437 // and able to return nodes by that ID
2438 // =================================================
2439 class LinkID_Gen {
2440 public:
2441
2442   LinkID_Gen( const SMESHDS_Mesh* theMesh )
2443     :myMesh( theMesh ), myMaxID( theMesh->MaxNodeID() + 1)
2444   {}
2445
2446   long GetLinkID (const SMDS_MeshNode * n1,
2447                   const SMDS_MeshNode * n2) const
2448   {
2449     return ( Min(n1->GetID(),n2->GetID()) * myMaxID + Max(n1->GetID(),n2->GetID()));
2450   }
2451
2452   bool GetNodes (const long             theLinkID,
2453                  const SMDS_MeshNode* & theNode1,
2454                  const SMDS_MeshNode* & theNode2) const
2455   {
2456     theNode1 = myMesh->FindNode( theLinkID / myMaxID );
2457     if ( !theNode1 ) return false;
2458     theNode2 = myMesh->FindNode( theLinkID % myMaxID );
2459     if ( !theNode2 ) return false;
2460     return true;
2461   }
2462
2463 private:
2464   LinkID_Gen();
2465   const SMESHDS_Mesh* myMesh;
2466   long                myMaxID;
2467 };
2468
2469
2470 //=======================================================================
2471 //function : TriToQuad
2472 //purpose  : Fuse neighbour triangles into quadrangles.
2473 //           theCrit is used to select a neighbour to fuse with.
2474 //           theMaxAngle is a max angle between element normals at which
2475 //           fusion is still performed.
2476 //=======================================================================
2477
2478 bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet &                   theElems,
2479                                   SMESH::Controls::NumericalFunctorPtr theCrit,
2480                                   const double                         theMaxAngle)
2481 {
2482   myLastCreatedElems.Clear();
2483   myLastCreatedNodes.Clear();
2484
2485   MESSAGE( "::TriToQuad()" );
2486
2487   if ( !theCrit.get() )
2488     return false;
2489
2490   SMESHDS_Mesh * aMesh = GetMeshDS();
2491
2492   // Prepare data for algo: build
2493   // 1. map of elements with their linkIDs
2494   // 2. map of linkIDs with their elements
2495
2496   map< SMESH_TLink, list< const SMDS_MeshElement* > > mapLi_listEl;
2497   map< SMESH_TLink, list< const SMDS_MeshElement* > >::iterator itLE;
2498   map< const SMDS_MeshElement*, set< SMESH_TLink > >  mapEl_setLi;
2499   map< const SMDS_MeshElement*, set< SMESH_TLink > >::iterator itEL;
2500
2501   TIDSortedElemSet::iterator itElem;
2502   for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ )
2503   {
2504     const SMDS_MeshElement* elem = *itElem;
2505     if(!elem || elem->GetType() != SMDSAbs_Face ) continue;
2506     bool IsTria = ( elem->NbCornerNodes()==3 );
2507     if (!IsTria) continue;
2508
2509     // retrieve element nodes
2510     const SMDS_MeshNode* aNodes [4];
2511     SMDS_NodeIteratorPtr itN = elem->nodeIterator();
2512     int i = 0;
2513     while ( i < 3 )
2514       aNodes[ i++ ] = itN->next();
2515     aNodes[ 3 ] = aNodes[ 0 ];
2516
2517     // fill maps
2518     for ( i = 0; i < 3; i++ ) {
2519       SMESH_TLink link( aNodes[i], aNodes[i+1] );
2520       // check if elements sharing a link can be fused
2521       itLE = mapLi_listEl.find( link );
2522       if ( itLE != mapLi_listEl.end() ) {
2523         if ((*itLE).second.size() > 1 ) // consider only 2 elems adjacent by a link
2524           continue;
2525         const SMDS_MeshElement* elem2 = (*itLE).second.front();
2526         //if ( FindShape( elem ) != FindShape( elem2 ))
2527         //  continue; // do not fuse triangles laying on different shapes
2528         if ( getAngle( elem, elem2, aNodes[i], aNodes[i+1] ) > theMaxAngle )
2529           continue; // avoid making badly shaped quads
2530         (*itLE).second.push_back( elem );
2531       }
2532       else {
2533         mapLi_listEl[ link ].push_back( elem );
2534       }
2535       mapEl_setLi [ elem ].insert( link );
2536     }
2537   }
2538   // Clean the maps from the links shared by a sole element, ie
2539   // links to which only one element is bound in mapLi_listEl
2540
2541   for ( itLE = mapLi_listEl.begin(); itLE != mapLi_listEl.end(); itLE++ ) {
2542     int nbElems = (*itLE).second.size();
2543     if ( nbElems < 2  ) {
2544       const SMDS_MeshElement* elem = (*itLE).second.front();
2545       SMESH_TLink link = (*itLE).first;
2546       mapEl_setLi[ elem ].erase( link );
2547       if ( mapEl_setLi[ elem ].empty() )
2548         mapEl_setLi.erase( elem );
2549     }
2550   }
2551
2552   // Algo: fuse triangles into quadrangles
2553
2554   while ( ! mapEl_setLi.empty() ) {
2555     // Look for the start element:
2556     // the element having the least nb of shared links
2557     const SMDS_MeshElement* startElem = 0;
2558     int minNbLinks = 4;
2559     for ( itEL = mapEl_setLi.begin(); itEL != mapEl_setLi.end(); itEL++ ) {
2560       int nbLinks = (*itEL).second.size();
2561       if ( nbLinks < minNbLinks ) {
2562         startElem = (*itEL).first;
2563         minNbLinks = nbLinks;
2564         if ( minNbLinks == 1 )
2565           break;
2566       }
2567     }
2568
2569     // search elements to fuse starting from startElem or links of elements
2570     // fused earlyer - startLinks
2571     list< SMESH_TLink > startLinks;
2572     while ( startElem || !startLinks.empty() ) {
2573       while ( !startElem && !startLinks.empty() ) {
2574         // Get an element to start, by a link
2575         SMESH_TLink linkId = startLinks.front();
2576         startLinks.pop_front();
2577         itLE = mapLi_listEl.find( linkId );
2578         if ( itLE != mapLi_listEl.end() ) {
2579           list< const SMDS_MeshElement* > & listElem = (*itLE).second;
2580           list< const SMDS_MeshElement* >::iterator itE = listElem.begin();
2581           for ( ; itE != listElem.end() ; itE++ )
2582             if ( mapEl_setLi.find( (*itE) ) != mapEl_setLi.end() )
2583               startElem = (*itE);
2584           mapLi_listEl.erase( itLE );
2585         }
2586       }
2587
2588       if ( startElem ) {
2589         // Get candidates to be fused
2590         const SMDS_MeshElement *tr1 = startElem, *tr2 = 0, *tr3 = 0;
2591         const SMESH_TLink *link12, *link13;
2592         startElem = 0;
2593         ASSERT( mapEl_setLi.find( tr1 ) != mapEl_setLi.end() );
2594         set< SMESH_TLink >& setLi = mapEl_setLi[ tr1 ];
2595         ASSERT( !setLi.empty() );
2596         set< SMESH_TLink >::iterator itLi;
2597         for ( itLi = setLi.begin(); itLi != setLi.end(); itLi++ )
2598         {
2599           const SMESH_TLink & link = (*itLi);
2600           itLE = mapLi_listEl.find( link );
2601           if ( itLE == mapLi_listEl.end() )
2602             continue;
2603
2604           const SMDS_MeshElement* elem = (*itLE).second.front();
2605           if ( elem == tr1 )
2606             elem = (*itLE).second.back();
2607           mapLi_listEl.erase( itLE );
2608           if ( mapEl_setLi.find( elem ) == mapEl_setLi.end())
2609             continue;
2610           if ( tr2 ) {
2611             tr3 = elem;
2612             link13 = &link;
2613           }
2614           else {
2615             tr2 = elem;
2616             link12 = &link;
2617           }
2618
2619           // add other links of elem to list of links to re-start from
2620           set< SMESH_TLink >& links = mapEl_setLi[ elem ];
2621           set< SMESH_TLink >::iterator it;
2622           for ( it = links.begin(); it != links.end(); it++ ) {
2623             const SMESH_TLink& link2 = (*it);
2624             if ( link2 != link )
2625               startLinks.push_back( link2 );
2626           }
2627         }
2628
2629         // Get nodes of possible quadrangles
2630         const SMDS_MeshNode *n12 [4], *n13 [4];
2631         bool Ok12 = false, Ok13 = false;
2632         const SMDS_MeshNode *linkNode1, *linkNode2;
2633         if(tr2) {
2634           linkNode1 = link12->first;
2635           linkNode2 = link12->second;
2636           if ( tr2 && getQuadrangleNodes( n12, linkNode1, linkNode2, tr1, tr2 ))
2637             Ok12 = true;
2638         }
2639         if(tr3) {
2640           linkNode1 = link13->first;
2641           linkNode2 = link13->second;
2642           if ( tr3 && getQuadrangleNodes( n13, linkNode1, linkNode2, tr1, tr3 ))
2643             Ok13 = true;
2644         }
2645
2646         // Choose a pair to fuse
2647         if ( Ok12 && Ok13 ) {
2648           SMDS_FaceOfNodes quad12 ( n12[ 0 ], n12[ 1 ], n12[ 2 ], n12[ 3 ] );
2649           SMDS_FaceOfNodes quad13 ( n13[ 0 ], n13[ 1 ], n13[ 2 ], n13[ 3 ] );
2650           double aBadRate12 = getBadRate( &quad12, theCrit );
2651           double aBadRate13 = getBadRate( &quad13, theCrit );
2652           if (  aBadRate13 < aBadRate12 )
2653             Ok12 = false;
2654           else
2655             Ok13 = false;
2656         }
2657
2658         // Make quadrangles
2659         // and remove fused elems and remove links from the maps
2660         mapEl_setLi.erase( tr1 );
2661         if ( Ok12 )
2662         {
2663           mapEl_setLi.erase( tr2 );
2664           mapLi_listEl.erase( *link12 );
2665           if ( tr1->NbNodes() == 3 )
2666           {
2667             const SMDS_MeshElement* newElem = 0;
2668             newElem = aMesh->AddFace(n12[0], n12[1], n12[2], n12[3] );
2669             myLastCreatedElems.Append(newElem);
2670             AddToSameGroups( newElem, tr1, aMesh );
2671             int aShapeId = tr1->getshapeId();
2672             if ( aShapeId )
2673               aMesh->SetMeshElementOnShape( newElem, aShapeId );
2674             aMesh->RemoveElement( tr1 );
2675             aMesh->RemoveElement( tr2 );
2676           }
2677           else {
2678             vector< const SMDS_MeshNode* > N1;
2679             vector< const SMDS_MeshNode* > N2;
2680             getNodesFromTwoTria(tr1,tr2,N1,N2);
2681             // now we receive following N1 and N2 (using numeration as in image in InverseDiag())
2682             // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
2683             // i.e. first nodes from both arrays form a new diagonal
2684             const SMDS_MeshNode* aNodes[8];
2685             aNodes[0] = N1[0];
2686             aNodes[1] = N1[1];
2687             aNodes[2] = N2[0];
2688             aNodes[3] = N2[1];
2689             aNodes[4] = N1[3];
2690             aNodes[5] = N2[5];
2691             aNodes[6] = N2[3];
2692             aNodes[7] = N1[5];
2693             const SMDS_MeshElement* newElem = 0;
2694             if ( N1.size() == 7 || N2.size() == 7 ) // biquadratic
2695               newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3],
2696                                        aNodes[4], aNodes[5], aNodes[6], aNodes[7], N1[4]);
2697             else
2698               newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3],
2699                                        aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
2700             myLastCreatedElems.Append(newElem);
2701             AddToSameGroups( newElem, tr1, aMesh );
2702             int aShapeId = tr1->getshapeId();
2703             if ( aShapeId )
2704               aMesh->SetMeshElementOnShape( newElem, aShapeId );
2705             aMesh->RemoveElement( tr1 );
2706             aMesh->RemoveElement( tr2 );
2707             // remove middle node (9)
2708             if ( N1[4]->NbInverseElements() == 0 )
2709               aMesh->RemoveNode( N1[4] );
2710             if ( N1.size() == 7 && N1[6]->NbInverseElements() == 0 )
2711               aMesh->RemoveNode( N1[6] );
2712             if ( N2.size() == 7 && N2[6]->NbInverseElements() == 0 )
2713               aMesh->RemoveNode( N2[6] );
2714           }
2715         }
2716         else if ( Ok13 )
2717         {
2718           mapEl_setLi.erase( tr3 );
2719           mapLi_listEl.erase( *link13 );
2720           if ( tr1->NbNodes() == 3 ) {
2721             const SMDS_MeshElement* newElem = 0;
2722             newElem = aMesh->AddFace(n13[0], n13[1], n13[2], n13[3] );
2723             myLastCreatedElems.Append(newElem);
2724             AddToSameGroups( newElem, tr1, aMesh );
2725             int aShapeId = tr1->getshapeId();
2726             if ( aShapeId )
2727               aMesh->SetMeshElementOnShape( newElem, aShapeId );
2728             aMesh->RemoveElement( tr1 );
2729             aMesh->RemoveElement( tr3 );
2730           }
2731           else {
2732             vector< const SMDS_MeshNode* > N1;
2733             vector< const SMDS_MeshNode* > N2;
2734             getNodesFromTwoTria(tr1,tr3,N1,N2);
2735             // now we receive following N1 and N2 (using numeration as above image)
2736             // tria1 : (1 2 4 5 9 7)  and  tria2 : (3 4 2 8 9 6)
2737             // i.e. first nodes from both arrays form a new diagonal
2738             const SMDS_MeshNode* aNodes[8];
2739             aNodes[0] = N1[0];
2740             aNodes[1] = N1[1];
2741             aNodes[2] = N2[0];
2742             aNodes[3] = N2[1];
2743             aNodes[4] = N1[3];
2744             aNodes[5] = N2[5];
2745             aNodes[6] = N2[3];
2746             aNodes[7] = N1[5];
2747             const SMDS_MeshElement* newElem = 0;
2748             if ( N1.size() == 7 || N2.size() == 7 ) // biquadratic
2749               newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3],
2750                                        aNodes[4], aNodes[5], aNodes[6], aNodes[7], N1[4]);
2751             else
2752               newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3],
2753                                        aNodes[4], aNodes[5], aNodes[6], aNodes[7]);
2754             myLastCreatedElems.Append(newElem);
2755             AddToSameGroups( newElem, tr1, aMesh );
2756             int aShapeId = tr1->getshapeId();
2757             if ( aShapeId )
2758               aMesh->SetMeshElementOnShape( newElem, aShapeId );
2759             aMesh->RemoveElement( tr1 );
2760             aMesh->RemoveElement( tr3 );
2761             // remove middle node (9)
2762             if ( N1[4]->NbInverseElements() == 0 )
2763               aMesh->RemoveNode( N1[4] );
2764             if ( N1.size() == 7 && N1[6]->NbInverseElements() == 0 )
2765               aMesh->RemoveNode( N1[6] );
2766             if ( N2.size() == 7 && N2[6]->NbInverseElements() == 0 )
2767               aMesh->RemoveNode( N2[6] );
2768           }
2769         }
2770
2771         // Next element to fuse: the rejected one
2772         if ( tr3 )
2773           startElem = Ok12 ? tr3 : tr2;
2774
2775       } // if ( startElem )
2776     } // while ( startElem || !startLinks.empty() )
2777   } // while ( ! mapEl_setLi.empty() )
2778
2779   return true;
2780 }
2781
2782
2783 /*#define DUMPSO(txt) \
2784 //  cout << txt << endl;
2785 //=============================================================================
2786 //
2787 //
2788 //
2789 //=============================================================================
2790 static void swap( int i1, int i2, int idNodes[], gp_Pnt P[] )
2791 {
2792 if ( i1 == i2 )
2793 return;
2794 int tmp = idNodes[ i1 ];
2795 idNodes[ i1 ] = idNodes[ i2 ];
2796 idNodes[ i2 ] = tmp;
2797 gp_Pnt Ptmp = P[ i1 ];
2798 P[ i1 ] = P[ i2 ];
2799 P[ i2 ] = Ptmp;
2800 DUMPSO( i1 << "(" << idNodes[ i2 ] << ") <-> " << i2 << "(" << idNodes[ i1 ] << ")");
2801 }
2802
2803 //=======================================================================
2804 //function : SortQuadNodes
2805 //purpose  : Set 4 nodes of a quadrangle face in a good order.
2806 //           Swap 1<->2 or 2<->3 nodes and correspondingly return
2807 //           1 or 2 else 0.
2808 //=======================================================================
2809
2810 int SMESH_MeshEditor::SortQuadNodes (const SMDS_Mesh * theMesh,
2811 int               idNodes[] )
2812 {
2813   gp_Pnt P[4];
2814   int i;
2815   for ( i = 0; i < 4; i++ ) {
2816     const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
2817     if ( !n ) return 0;
2818     P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
2819   }
2820
2821   gp_Vec V1(P[0], P[1]);
2822   gp_Vec V2(P[0], P[2]);
2823   gp_Vec V3(P[0], P[3]);
2824
2825   gp_Vec Cross1 = V1 ^ V2;
2826   gp_Vec Cross2 = V2 ^ V3;
2827
2828   i = 0;
2829   if (Cross1.Dot(Cross2) < 0)
2830   {
2831     Cross1 = V2 ^ V1;
2832     Cross2 = V1 ^ V3;
2833
2834     if (Cross1.Dot(Cross2) < 0)
2835       i = 2;
2836     else
2837       i = 1;
2838     swap ( i, i + 1, idNodes, P );
2839
2840     //     for ( int ii = 0; ii < 4; ii++ ) {
2841     //       const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
2842     //       DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
2843     //     }
2844   }
2845   return i;
2846 }
2847
2848 //=======================================================================
2849 //function : SortHexaNodes
2850 //purpose  : Set 8 nodes of a hexahedron in a good order.
2851 //           Return success status
2852 //=======================================================================
2853
2854 bool SMESH_MeshEditor::SortHexaNodes (const SMDS_Mesh * theMesh,
2855                                       int               idNodes[] )
2856 {
2857   gp_Pnt P[8];
2858   int i;
2859   DUMPSO( "INPUT: ========================================");
2860   for ( i = 0; i < 8; i++ ) {
2861     const SMDS_MeshNode *n = theMesh->FindNode( idNodes[i] );
2862     if ( !n ) return false;
2863     P[ i ].SetCoord( n->X(), n->Y(), n->Z() );
2864     DUMPSO( i << "(" << idNodes[i] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
2865   }
2866   DUMPSO( "========================================");
2867
2868
2869   set<int> faceNodes;  // ids of bottom face nodes, to be found
2870   set<int> checkedId1; // ids of tried 2-nd nodes
2871   Standard_Real leastDist = DBL_MAX; // dist of the 4-th node from 123 plane
2872   const Standard_Real tol = 1.e-6;   // tolerance to find nodes in plane
2873   int iMin, iLoop1 = 0;
2874
2875   // Loop to try the 2-nd nodes
2876
2877   while ( leastDist > DBL_MIN && ++iLoop1 < 8 )
2878   {
2879     // Find not checked 2-nd node
2880     for ( i = 1; i < 8; i++ )
2881       if ( checkedId1.find( idNodes[i] ) == checkedId1.end() ) {
2882         int id1 = idNodes[i];
2883         swap ( 1, i, idNodes, P );
2884         checkedId1.insert ( id1 );
2885         break;
2886       }
2887
2888     // Find the 3-d node so that 1-2-3 triangle to be on a hexa face,
2889     // ie that all but meybe one (id3 which is on the same face) nodes
2890     // lay on the same side from the triangle plane.
2891
2892     bool manyInPlane = false; // more than 4 nodes lay in plane
2893     int iLoop2 = 0;
2894     while ( ++iLoop2 < 6 ) {
2895
2896       // get 1-2-3 plane coeffs
2897       Standard_Real A, B, C, D;
2898       gp_Vec N = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
2899       if ( N.SquareMagnitude() > gp::Resolution() )
2900       {
2901         gp_Pln pln ( P[0], N );
2902         pln.Coefficients( A, B, C, D );
2903
2904         // find the node (iMin) closest to pln
2905         Standard_Real dist[ 8 ], minDist = DBL_MAX;
2906         set<int> idInPln;
2907         for ( i = 3; i < 8; i++ ) {
2908           dist[i] = A * P[i].X() + B * P[i].Y() + C * P[i].Z() + D;
2909           if ( fabs( dist[i] ) < minDist ) {
2910             minDist = fabs( dist[i] );
2911             iMin = i;
2912           }
2913           if ( fabs( dist[i] ) <= tol )
2914             idInPln.insert( idNodes[i] );
2915         }
2916
2917         // there should not be more than 4 nodes in bottom plane
2918         if ( idInPln.size() > 1 )
2919         {
2920           DUMPSO( "### idInPln.size() = " << idInPln.size());
2921           // idInPlane does not contain the first 3 nodes
2922           if ( manyInPlane || idInPln.size() == 5)
2923             return false; // all nodes in one plane
2924           manyInPlane = true;
2925
2926           // set the 1-st node to be not in plane
2927           for ( i = 3; i < 8; i++ ) {
2928             if ( idInPln.find( idNodes[ i ] ) == idInPln.end() ) {
2929               DUMPSO( "### Reset 0-th node");
2930               swap( 0, i, idNodes, P );
2931               break;
2932             }
2933           }
2934
2935           // reset to re-check second nodes
2936           leastDist = DBL_MAX;
2937           faceNodes.clear();
2938           checkedId1.clear();
2939           iLoop1 = 0;
2940           break; // from iLoop2;
2941         }
2942
2943         // check that the other 4 nodes are on the same side
2944         bool sameSide = true;
2945         bool isNeg = dist[ iMin == 3 ? 4 : 3 ] <= 0.;
2946         for ( i = 3; sameSide && i < 8; i++ ) {
2947           if ( i != iMin )
2948             sameSide = ( isNeg == dist[i] <= 0.);
2949         }
2950
2951         // keep best solution
2952         if ( sameSide && minDist < leastDist ) {
2953           leastDist = minDist;
2954           faceNodes.clear();
2955           faceNodes.insert( idNodes[ 1 ] );
2956           faceNodes.insert( idNodes[ 2 ] );
2957           faceNodes.insert( idNodes[ iMin ] );
2958           DUMPSO( "loop " << iLoop2 << " id2 " << idNodes[ 1 ] << " id3 " << idNodes[ 2 ]
2959                   << " leastDist = " << leastDist);
2960           if ( leastDist <= DBL_MIN )
2961             break;
2962         }
2963       }
2964
2965       // set next 3-d node to check
2966       int iNext = 2 + iLoop2;
2967       if ( iNext < 8 ) {
2968         DUMPSO( "Try 2-nd");
2969         swap ( 2, iNext, idNodes, P );
2970       }
2971     } // while ( iLoop2 < 6 )
2972   } // iLoop1
2973
2974   if ( faceNodes.empty() ) return false;
2975
2976   // Put the faceNodes in proper places
2977   for ( i = 4; i < 8; i++ ) {
2978     if ( faceNodes.find( idNodes[ i ] ) != faceNodes.end() ) {
2979       // find a place to put
2980       int iTo = 1;
2981       while ( faceNodes.find( idNodes[ iTo ] ) != faceNodes.end() )
2982         iTo++;
2983       DUMPSO( "Set faceNodes");
2984       swap ( iTo, i, idNodes, P );
2985     }
2986   }
2987
2988
2989   // Set nodes of the found bottom face in good order
2990   DUMPSO( " Found bottom face: ");
2991   i = SortQuadNodes( theMesh, idNodes );
2992   if ( i ) {
2993     gp_Pnt Ptmp = P[ i ];
2994     P[ i ] = P[ i+1 ];
2995     P[ i+1 ] = Ptmp;
2996   }
2997   //   else
2998   //     for ( int ii = 0; ii < 4; ii++ ) {
2999   //       const SMDS_MeshNode *n = theMesh->FindNode( idNodes[ii] );
3000   //       DUMPSO( ii << "(" << idNodes[ii] <<") : "<<n->X()<<" "<<n->Y()<<" "<<n->Z());
3001   //    }
3002
3003   // Gravity center of the top and bottom faces
3004   gp_Pnt aGCb = ( P[0].XYZ() + P[1].XYZ() + P[2].XYZ() + P[3].XYZ() ) / 4.;
3005   gp_Pnt aGCt = ( P[4].XYZ() + P[5].XYZ() + P[6].XYZ() + P[7].XYZ() ) / 4.;
3006
3007   // Get direction from the bottom to the top face
3008   gp_Vec upDir ( aGCb, aGCt );
3009   Standard_Real upDirSize = upDir.Magnitude();
3010   if ( upDirSize <= gp::Resolution() ) return false;
3011   upDir / upDirSize;
3012
3013   // Assure that the bottom face normal points up
3014   gp_Vec Nb = gp_Vec (P[0], P[1]).Crossed( gp_Vec (P[0], P[2]) );
3015   Nb += gp_Vec (P[0], P[2]).Crossed( gp_Vec (P[0], P[3]) );
3016   if ( Nb.Dot( upDir ) < 0 ) {
3017     DUMPSO( "Reverse bottom face");
3018     swap( 1, 3, idNodes, P );
3019   }
3020
3021   // Find 5-th node - the one closest to the 1-st among the last 4 nodes.
3022   Standard_Real minDist = DBL_MAX;
3023   for ( i = 4; i < 8; i++ ) {
3024     // projection of P[i] to the plane defined by P[0] and upDir
3025     gp_Pnt Pp = P[i].Translated( upDir * ( upDir.Dot( gp_Vec( P[i], P[0] ))));
3026     Standard_Real sqDist = P[0].SquareDistance( Pp );
3027     if ( sqDist < minDist ) {
3028       minDist = sqDist;
3029       iMin = i;
3030     }
3031   }
3032   DUMPSO( "Set 4-th");
3033   swap ( 4, iMin, idNodes, P );
3034
3035   // Set nodes of the top face in good order
3036   DUMPSO( "Sort top face");
3037   i = SortQuadNodes( theMesh, &idNodes[4] );
3038   if ( i ) {
3039     i += 4;
3040     gp_Pnt Ptmp = P[ i ];
3041     P[ i ] = P[ i+1 ];
3042     P[ i+1 ] = Ptmp;
3043   }
3044
3045   // Assure that direction of the top face normal is from the bottom face
3046   gp_Vec Nt = gp_Vec (P[4], P[5]).Crossed( gp_Vec (P[4], P[6]) );
3047   Nt += gp_Vec (P[4], P[6]).Crossed( gp_Vec (P[4], P[7]) );
3048   if ( Nt.Dot( upDir ) < 0 ) {
3049     DUMPSO( "Reverse top face");
3050     swap( 5, 7, idNodes, P );
3051   }
3052
3053   //   DUMPSO( "OUTPUT: ========================================");
3054   //   for ( i = 0; i < 8; i++ ) {
3055   //     float *p = ugrid->GetPoint(idNodes[i]);
3056   //     DUMPSO( i << "(" << idNodes[i] << ") : " << p[0] << " " << p[1] << " " << p[2]);
3057   //   }
3058
3059   return true;
3060 }*/
3061
3062 //================================================================================
3063 /*!
3064  * \brief Return nodes linked to the given one
3065  * \param theNode - the node
3066  * \param linkedNodes - the found nodes
3067  * \param type - the type of elements to check
3068  *
3069  * Medium nodes are ignored
3070  */
3071 //================================================================================
3072
3073 void SMESH_MeshEditor::GetLinkedNodes( const SMDS_MeshNode* theNode,
3074                                        TIDSortedElemSet &   linkedNodes,
3075                                        SMDSAbs_ElementType  type )
3076 {
3077   SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(type);
3078   while ( elemIt->more() )
3079   {
3080     const SMDS_MeshElement* elem = elemIt->next();
3081     if(elem->GetType() == SMDSAbs_0DElement)
3082       continue;
3083
3084     SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
3085     if ( elem->GetType() == SMDSAbs_Volume )
3086     {
3087       SMDS_VolumeTool vol( elem );
3088       while ( nodeIt->more() ) {
3089         const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
3090         if ( theNode != n && vol.IsLinked( theNode, n ))
3091           linkedNodes.insert( n );
3092       }
3093     }
3094     else
3095     {
3096       for ( int i = 0; nodeIt->more(); ++i ) {
3097         const SMDS_MeshNode* n = cast2Node( nodeIt->next() );
3098         if ( n == theNode ) {
3099           int iBefore = i - 1;
3100           int iAfter  = i + 1;
3101           if ( elem->IsQuadratic() ) {
3102             int nb = elem->NbNodes() / 2;
3103             iAfter  = SMESH_MesherHelper::WrapIndex( iAfter, nb );
3104             iBefore = SMESH_MesherHelper::WrapIndex( iBefore, nb );
3105           }
3106           linkedNodes.insert( elem->GetNodeWrap( iAfter ));
3107           linkedNodes.insert( elem->GetNodeWrap( iBefore ));
3108         }
3109       }
3110     }
3111   }
3112 }
3113
3114 //=======================================================================
3115 //function : laplacianSmooth
3116 //purpose  : pulls theNode toward the center of surrounding nodes directly
3117 //           connected to that node along an element edge
3118 //=======================================================================
3119
3120 void laplacianSmooth(const SMDS_MeshNode*                 theNode,
3121                      const Handle(Geom_Surface)&          theSurface,
3122                      map< const SMDS_MeshNode*, gp_XY* >& theUVMap)
3123 {
3124   // find surrounding nodes
3125
3126   TIDSortedElemSet nodeSet;
3127   SMESH_MeshEditor::GetLinkedNodes( theNode, nodeSet, SMDSAbs_Face );
3128
3129   // compute new coodrs
3130
3131   double coord[] = { 0., 0., 0. };
3132   TIDSortedElemSet::iterator nodeSetIt = nodeSet.begin();
3133   for ( ; nodeSetIt != nodeSet.end(); nodeSetIt++ ) {
3134     const SMDS_MeshNode* node = cast2Node(*nodeSetIt);
3135     if ( theSurface.IsNull() ) { // smooth in 3D
3136       coord[0] += node->X();
3137       coord[1] += node->Y();
3138       coord[2] += node->Z();
3139     }
3140     else { // smooth in 2D
3141       ASSERT( theUVMap.find( node ) != theUVMap.end() );
3142       gp_XY* uv = theUVMap[ node ];
3143       coord[0] += uv->X();
3144       coord[1] += uv->Y();
3145     }
3146   }
3147   int nbNodes = nodeSet.size();
3148   if ( !nbNodes )
3149     return;
3150   coord[0] /= nbNodes;
3151   coord[1] /= nbNodes;
3152
3153   if ( !theSurface.IsNull() ) {
3154     ASSERT( theUVMap.find( theNode ) != theUVMap.end() );
3155     theUVMap[ theNode ]->SetCoord( coord[0], coord[1] );
3156     gp_Pnt p3d = theSurface->Value( coord[0], coord[1] );
3157     coord[0] = p3d.X();
3158     coord[1] = p3d.Y();
3159     coord[2] = p3d.Z();
3160   }
3161   else
3162     coord[2] /= nbNodes;
3163
3164   // move node
3165
3166   const_cast< SMDS_MeshNode* >( theNode )->setXYZ(coord[0],coord[1],coord[2]);
3167 }
3168
3169 //=======================================================================
3170 //function : centroidalSmooth
3171 //purpose  : pulls theNode toward the element-area-weighted centroid of the
3172 //           surrounding elements
3173 //=======================================================================
3174
3175 void centroidalSmooth(const SMDS_MeshNode*                 theNode,
3176                       const Handle(Geom_Surface)&          theSurface,
3177                       map< const SMDS_MeshNode*, gp_XY* >& theUVMap)
3178 {
3179   gp_XYZ aNewXYZ(0.,0.,0.);
3180   SMESH::Controls::Area anAreaFunc;
3181   double totalArea = 0.;
3182   int nbElems = 0;
3183
3184   // compute new XYZ
3185
3186   SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(SMDSAbs_Face);
3187   while ( elemIt->more() )
3188   {
3189     const SMDS_MeshElement* elem = elemIt->next();
3190     nbElems++;
3191
3192     gp_XYZ elemCenter(0.,0.,0.);
3193     SMESH::Controls::TSequenceOfXYZ aNodePoints;
3194     SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3195     int nn = elem->NbNodes();
3196     if(elem->IsQuadratic()) nn = nn/2;
3197     int i=0;
3198     //while ( itN->more() ) {
3199     while ( i<nn ) {
3200       const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( itN->next() );
3201       i++;
3202       gp_XYZ aP( aNode->X(), aNode->Y(), aNode->Z() );
3203       aNodePoints.push_back( aP );
3204       if ( !theSurface.IsNull() ) { // smooth in 2D
3205         ASSERT( theUVMap.find( aNode ) != theUVMap.end() );
3206         gp_XY* uv = theUVMap[ aNode ];
3207         aP.SetCoord( uv->X(), uv->Y(), 0. );
3208       }
3209       elemCenter += aP;
3210     }
3211     double elemArea = anAreaFunc.GetValue( aNodePoints );
3212     totalArea += elemArea;
3213     elemCenter /= nn;
3214     aNewXYZ += elemCenter * elemArea;
3215   }
3216   aNewXYZ /= totalArea;
3217   if ( !theSurface.IsNull() ) {
3218     theUVMap[ theNode ]->SetCoord( aNewXYZ.X(), aNewXYZ.Y() );
3219     aNewXYZ = theSurface->Value( aNewXYZ.X(), aNewXYZ.Y() ).XYZ();
3220   }
3221
3222   // move node
3223
3224   const_cast< SMDS_MeshNode* >( theNode )->setXYZ(aNewXYZ.X(),aNewXYZ.Y(),aNewXYZ.Z());
3225 }
3226
3227 //=======================================================================
3228 //function : getClosestUV
3229 //purpose  : return UV of closest projection
3230 //=======================================================================
3231
3232 static bool getClosestUV (Extrema_GenExtPS& projector,
3233                           const gp_Pnt&     point,
3234                           gp_XY &           result)
3235 {
3236   projector.Perform( point );
3237   if ( projector.IsDone() ) {
3238     double u, v, minVal = DBL_MAX;
3239     for ( int i = projector.NbExt(); i > 0; i-- )
3240 #if OCC_VERSION_LARGE > 0x06040000 // Porting to OCCT6.5.1
3241       if ( projector.SquareDistance( i ) < minVal ) {
3242         minVal = projector.SquareDistance( i );
3243 #else
3244       if ( projector.Value( i ) < minVal ) {
3245         minVal = projector.Value( i );
3246 #endif
3247         projector.Point( i ).Parameter( u, v );
3248       }
3249     result.SetCoord( u, v );
3250     return true;
3251   }
3252   return false;
3253 }
3254
3255 //=======================================================================
3256 //function : Smooth
3257 //purpose  : Smooth theElements during theNbIterations or until a worst
3258 //           element has aspect ratio <= theTgtAspectRatio.
3259 //           Aspect Ratio varies in range [1.0, inf].
3260 //           If theElements is empty, the whole mesh is smoothed.
3261 //           theFixedNodes contains additionally fixed nodes. Nodes built
3262 //           on edges and boundary nodes are always fixed.
3263 //=======================================================================
3264
3265 void SMESH_MeshEditor::Smooth (TIDSortedElemSet &          theElems,
3266                                set<const SMDS_MeshNode*> & theFixedNodes,
3267                                const SmoothMethod          theSmoothMethod,
3268                                const int                   theNbIterations,
3269                                double                      theTgtAspectRatio,
3270                                const bool                  the2D)
3271 {
3272   myLastCreatedElems.Clear();
3273   myLastCreatedNodes.Clear();
3274
3275   MESSAGE((theSmoothMethod==LAPLACIAN ? "LAPLACIAN" : "CENTROIDAL") << "--::Smooth()");
3276
3277   if ( theTgtAspectRatio < 1.0 )
3278     theTgtAspectRatio = 1.0;
3279
3280   const double disttol = 1.e-16;
3281
3282   SMESH::Controls::AspectRatio aQualityFunc;
3283
3284   SMESHDS_Mesh* aMesh = GetMeshDS();
3285
3286   if ( theElems.empty() ) {
3287     // add all faces to theElems
3288     SMDS_FaceIteratorPtr fIt = aMesh->facesIterator();
3289     while ( fIt->more() ) {
3290       const SMDS_MeshElement* face = fIt->next();
3291       theElems.insert( theElems.end(), face );
3292     }
3293   }
3294   // get all face ids theElems are on
3295   set< int > faceIdSet;
3296   TIDSortedElemSet::iterator itElem;
3297   if ( the2D )
3298     for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) {
3299       int fId = FindShape( *itElem );
3300       // check that corresponding submesh exists and a shape is face
3301       if (fId &&
3302           faceIdSet.find( fId ) == faceIdSet.end() &&
3303           aMesh->MeshElements( fId )) {
3304         TopoDS_Shape F = aMesh->IndexToShape( fId );
3305         if ( !F.IsNull() && F.ShapeType() == TopAbs_FACE )
3306           faceIdSet.insert( fId );
3307       }
3308     }
3309   faceIdSet.insert( 0 ); // to smooth elements that are not on any TopoDS_Face
3310
3311   // ===============================================
3312   // smooth elements on each TopoDS_Face separately
3313   // ===============================================
3314
3315   set< int >::reverse_iterator fId = faceIdSet.rbegin(); // treate 0 fId at the end
3316   for ( ; fId != faceIdSet.rend(); ++fId ) {
3317     // get face surface and submesh
3318     Handle(Geom_Surface) surface;
3319     SMESHDS_SubMesh* faceSubMesh = 0;
3320     TopoDS_Face face;
3321     double fToler2 = 0, f,l;
3322     double u1 = 0, u2 = 0, v1 = 0, v2 = 0;
3323     bool isUPeriodic = false, isVPeriodic = false;
3324     if ( *fId ) {
3325       face = TopoDS::Face( aMesh->IndexToShape( *fId ));
3326       surface = BRep_Tool::Surface( face );
3327       faceSubMesh = aMesh->MeshElements( *fId );
3328       fToler2 = BRep_Tool::Tolerance( face );
3329       fToler2 *= fToler2 * 10.;
3330       isUPeriodic = surface->IsUPeriodic();
3331       if ( isUPeriodic )
3332         surface->UPeriod();
3333       isVPeriodic = surface->IsVPeriodic();
3334       if ( isVPeriodic )
3335         surface->VPeriod();
3336       surface->Bounds( u1, u2, v1, v2 );
3337     }
3338     // ---------------------------------------------------------
3339     // for elements on a face, find movable and fixed nodes and
3340     // compute UV for them
3341     // ---------------------------------------------------------
3342     bool checkBoundaryNodes = false;
3343     bool isQuadratic = false;
3344     set<const SMDS_MeshNode*> setMovableNodes;
3345     map< const SMDS_MeshNode*, gp_XY* > uvMap, uvMap2;
3346     list< gp_XY > listUV; // uvs the 2 uvMaps refer to
3347     list< const SMDS_MeshElement* > elemsOnFace;
3348
3349     Extrema_GenExtPS projector;
3350     GeomAdaptor_Surface surfAdaptor;
3351     if ( !surface.IsNull() ) {
3352       surfAdaptor.Load( surface );
3353       projector.Initialize( surfAdaptor, 20,20, 1e-5,1e-5 );
3354     }
3355     int nbElemOnFace = 0;
3356     itElem = theElems.begin();
3357     // loop on not yet smoothed elements: look for elems on a face
3358     while ( itElem != theElems.end() ) {
3359       if ( faceSubMesh && nbElemOnFace == faceSubMesh->NbElements() )
3360         break; // all elements found
3361
3362       const SMDS_MeshElement* elem = *itElem;
3363       if ( !elem || elem->GetType() != SMDSAbs_Face || elem->NbNodes() < 3 ||
3364            ( faceSubMesh && !faceSubMesh->Contains( elem ))) {
3365         ++itElem;
3366         continue;
3367       }
3368       elemsOnFace.push_back( elem );
3369       theElems.erase( itElem++ );
3370       nbElemOnFace++;
3371
3372       if ( !isQuadratic )
3373         isQuadratic = elem->IsQuadratic();
3374
3375       // get movable nodes of elem
3376       const SMDS_MeshNode* node;
3377       SMDS_TypeOfPosition posType;
3378       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
3379       int nn = 0, nbn =  elem->NbNodes();
3380       if(elem->IsQuadratic())
3381         nbn = nbn/2;
3382       while ( nn++ < nbn ) {
3383         node = static_cast<const SMDS_MeshNode*>( itN->next() );
3384         const SMDS_PositionPtr& pos = node->GetPosition();
3385         posType = pos ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
3386         if (posType != SMDS_TOP_EDGE &&
3387             posType != SMDS_TOP_VERTEX &&
3388             theFixedNodes.find( node ) == theFixedNodes.end())
3389         {
3390           // check if all faces around the node are on faceSubMesh
3391           // because a node on edge may be bound to face
3392           SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Face);
3393           bool all = true;
3394           if ( faceSubMesh ) {
3395             while ( eIt->more() && all ) {
3396               const SMDS_MeshElement* e = eIt->next();
3397               all = faceSubMesh->Contains( e );
3398             }
3399           }
3400           if ( all )
3401             setMovableNodes.insert( node );
3402           else
3403             checkBoundaryNodes = true;
3404         }
3405         if ( posType == SMDS_TOP_3DSPACE )
3406           checkBoundaryNodes = true;
3407       }
3408
3409       if ( surface.IsNull() )
3410         continue;
3411
3412       // get nodes to check UV
3413       list< const SMDS_MeshNode* > uvCheckNodes;
3414       itN = elem->nodesIterator();
3415       nn = 0; nbn =  elem->NbNodes();
3416       if(elem->IsQuadratic())
3417         nbn = nbn/2;
3418       while ( nn++ < nbn ) {
3419         node = static_cast<const SMDS_MeshNode*>( itN->next() );
3420         if ( uvMap.find( node ) == uvMap.end() )
3421           uvCheckNodes.push_back( node );
3422         // add nodes of elems sharing node
3423         //         SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Face);
3424         //         while ( eIt->more() ) {
3425         //           const SMDS_MeshElement* e = eIt->next();
3426         //           if ( e != elem ) {
3427         //             SMDS_ElemIteratorPtr nIt = e->nodesIterator();
3428         //             while ( nIt->more() ) {
3429         //               const SMDS_MeshNode* n =
3430         //                 static_cast<const SMDS_MeshNode*>( nIt->next() );
3431         //               if ( uvMap.find( n ) == uvMap.end() )
3432         //                 uvCheckNodes.push_back( n );
3433         //             }
3434         //           }
3435         //         }
3436       }
3437       // check UV on face
3438       list< const SMDS_MeshNode* >::iterator n = uvCheckNodes.begin();
3439       for ( ; n != uvCheckNodes.end(); ++n ) {
3440         node = *n;
3441         gp_XY uv( 0, 0 );
3442         const SMDS_PositionPtr& pos = node->GetPosition();
3443         posType = pos ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
3444         // get existing&nb