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