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