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