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