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