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