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