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