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