Salome HOME
53225: Viscous Layers hyp: problem with selection of input/output faces
[modules/smesh.git] / src / SMDS / SMDS_Mesh.cxx
1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
5 //
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
10 //
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14 // Lesser General Public License for more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
19 //
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 //
22
23 //  SMESH SMDS : implementation of Salome mesh data structure
24 //
25 #ifdef _MSC_VER
26 #pragma warning(disable:4786)
27 #endif
28
29 #include "SMDS_FaceOfEdges.hxx"
30 #include "SMDS_FaceOfNodes.hxx"
31 #include "SMDS_Mesh.hxx"
32 #include "SMDS_PolygonalFaceOfNodes.hxx"
33 #include "SMDS_PolyhedralVolumeOfNodes.hxx"
34 #include "SMDS_QuadraticEdge.hxx"
35 #include "SMDS_QuadraticFaceOfNodes.hxx"
36 #include "SMDS_QuadraticVolumeOfNodes.hxx"
37 #include "SMDS_SetIterator.hxx"
38 #include "SMDS_SpacePosition.hxx"
39 #include "SMDS_UnstructuredGrid.hxx"
40 #include "SMDS_VolumeOfFaces.hxx"
41 #include "SMDS_VolumeOfNodes.hxx"
42
43 #include "utilities.h"
44
45 #include <vtkUnstructuredGrid.h>
46 #include <vtkUnstructuredGridWriter.h>
47 #include <vtkUnsignedCharArray.h>
48 #include <vtkCell.h>
49 #include <vtkCellLinks.h>
50 #include <vtkIdList.h>
51
52 #include <algorithm>
53 #include <map>
54 #include <iostream>
55 #include <fstream>
56 #include <iterator>
57 using namespace std;
58
59 #ifndef WIN32
60 #include <sys/sysinfo.h>
61 #endif
62
63 // number of added entities to check memory after
64 #define CHECKMEMORY_INTERVAL 100000
65
66 vector<SMDS_Mesh*> SMDS_Mesh::_meshList = vector<SMDS_Mesh*>();
67 int SMDS_Mesh::chunkSize = 1024;
68
69
70 //================================================================================
71 /*!
72  * \brief Raise an exception if free memory (ram+swap) too low
73  * \param doNotRaise - if true, suppres exception, just return free memory size
74  * \retval int - amount of available memory in MB or negative number in failure case
75  */
76 //================================================================================
77
78 int SMDS_Mesh::CheckMemory(const bool doNotRaise) throw (std::bad_alloc)
79 {
80 #ifndef WIN32
81   struct sysinfo si;
82   int err = sysinfo( &si );
83   if ( err )
84     return -1;
85
86   const unsigned long Mbyte = 1024 * 1024;
87
88   static int limit = -1;
89   if ( limit < 0 ) {
90     int status = system("SMDS_MemoryLimit"); // it returns lower limit of free RAM
91     if (status >= 0 ) {
92       limit = WEXITSTATUS(status);
93     }
94     else {
95       double factor = ( si.totalswap == 0 ) ? 0.1 : 0.2;
96       limit = int(( factor * si.totalram * si.mem_unit ) / Mbyte );
97     }
98     if ( limit < 20 )
99       limit = 20;
100     else
101       limit = int ( limit * 1.5 );
102     MESSAGE ( "SMDS_Mesh::CheckMemory() memory limit = " << limit << " MB" );
103   }
104
105   // compute separately to avoid overflow
106   int freeMb =
107     ( si.freeram  * si.mem_unit ) / Mbyte +
108     ( si.freeswap * si.mem_unit ) / Mbyte;
109   //cout << "freeMb = " << freeMb << " limit = " << limit << endl;
110
111   if ( freeMb > limit )
112     return freeMb - limit;
113
114   if ( doNotRaise )
115     return 0;
116
117   MESSAGE ("SMDS_Mesh::CheckMemory() throws as free memory too low: " << freeMb <<" MB" );
118   throw std::bad_alloc();
119 #else
120   return -1;
121 #endif
122 }
123
124 ///////////////////////////////////////////////////////////////////////////////
125 /// Create a new mesh object
126 ///////////////////////////////////////////////////////////////////////////////
127 SMDS_Mesh::SMDS_Mesh():
128   myNodePool(0), myVolumePool(0), myFacePool(0), myEdgePool(0), myBallPool(0),
129   myParent(NULL),
130   myNodeIDFactory(new SMDS_MeshNodeIDFactory()),
131   myElementIDFactory(new SMDS_MeshElementIDFactory()),
132   myModified(false), myModifTime(0), myCompactTime(0),
133   myNodeMin(0), myNodeMax(0),
134   myHasConstructionEdges(false), myHasConstructionFaces(false),
135   myHasInverseElements(true),
136   xmin(0), xmax(0), ymin(0), ymax(0), zmin(0), zmax(0)
137 {
138   myMeshId = _meshList.size();         // --- index of the mesh to push back in the vector
139   myNodeIDFactory->SetMesh(this);
140   myElementIDFactory->SetMesh(this);
141   _meshList.push_back(this);
142   myNodePool = new ObjectPool<SMDS_MeshNode>(SMDS_Mesh::chunkSize);
143   myEdgePool = new ObjectPool<SMDS_VtkEdge>(SMDS_Mesh::chunkSize);
144   myFacePool = new ObjectPool<SMDS_VtkFace>(SMDS_Mesh::chunkSize);
145   myVolumePool = new ObjectPool<SMDS_VtkVolume>(SMDS_Mesh::chunkSize);
146   myBallPool = new ObjectPool<SMDS_BallElement>(SMDS_Mesh::chunkSize);
147
148   myNodes.clear();
149   myCells.clear();
150   //myCellIdSmdsToVtk.clear();
151   myCellIdVtkToSmds.clear();
152   myGrid = SMDS_UnstructuredGrid::New();
153   myGrid->setSMDS_mesh(this);
154   myGrid->Initialize();
155   myGrid->Allocate();
156   vtkPoints* points = vtkPoints::New();
157   // bug "21125: EDF 1233 SMESH: Degrardation of precision in a test case for quadratic conversion"
158   // Use double type for storing coordinates of nodes instead of float.
159   points->SetDataType(VTK_DOUBLE);
160   points->SetNumberOfPoints(0 /*SMDS_Mesh::chunkSize*/);
161   myGrid->SetPoints( points );
162   points->Delete();
163   myGrid->BuildLinks();
164   this->Modified();
165 }
166
167 ///////////////////////////////////////////////////////////////////////////////
168 /// Create a new child mesh
169 /// Note that the tree structure of SMDS_Mesh seems to be unused in this version
170 /// (2003-09-08) of SMESH
171 ///////////////////////////////////////////////////////////////////////////////
172 SMDS_Mesh::SMDS_Mesh(SMDS_Mesh * parent):
173   myNodePool(parent->myNodePool),
174   myVolumePool(parent->myVolumePool),
175   myFacePool(parent->myFacePool),
176   myEdgePool(parent->myEdgePool),
177   myBallPool(parent->myBallPool),
178   myParent(parent), myNodeIDFactory(parent->myNodeIDFactory),
179   myElementIDFactory(parent->myElementIDFactory),
180   myHasConstructionEdges(false), myHasConstructionFaces(false),
181   myHasInverseElements(true)
182 {
183 }
184
185 ///////////////////////////////////////////////////////////////////////////////
186 ///Create a submesh and add it to the current mesh
187 ///////////////////////////////////////////////////////////////////////////////
188
189 SMDS_Mesh *SMDS_Mesh::AddSubMesh()
190 {
191   SMDS_Mesh *submesh = new SMDS_Mesh(this);
192   myChildren.insert(myChildren.end(), submesh);
193   return submesh;
194 }
195
196 ///////////////////////////////////////////////////////////////////////////////
197 ///create a MeshNode and add it to the current Mesh
198 ///An ID is automatically assigned to the node.
199 ///@return : The created node
200 ///////////////////////////////////////////////////////////////////////////////
201
202 SMDS_MeshNode * SMDS_Mesh::AddNode(double x, double y, double z)
203 {
204   return SMDS_Mesh::AddNodeWithID(x,y,z,myNodeIDFactory->GetFreeID());
205 }
206
207 ///////////////////////////////////////////////////////////////////////////////
208 ///create a MeshNode and add it to the current Mesh
209 ///@param ID : The ID of the MeshNode to create
210 ///@return : The created node or NULL if a node with this ID already exists
211 ///////////////////////////////////////////////////////////////////////////////
212 SMDS_MeshNode * SMDS_Mesh::AddNodeWithID(double x, double y, double z, int ID)
213 {
214   // find the MeshNode corresponding to ID
215   const SMDS_MeshElement *node = myNodeIDFactory->MeshElement(ID);
216   if(!node){
217     if (ID < 1)
218       {
219         MESSAGE("=============>  Bad Node Id: " << ID);
220         ID = myNodeIDFactory->GetFreeID();
221       }
222     myNodeIDFactory->adjustMaxId(ID);
223     SMDS_MeshNode * node = myNodePool->getNew();
224     node->init(ID, myMeshId, 0, x, y, z);
225
226     if (ID >= (int)myNodes.size())
227     {
228         myNodes.resize(ID+SMDS_Mesh::chunkSize, 0);
229 //         MESSAGE(" ------------------ myNodes resize " << ID << " --> " << ID+SMDS_Mesh::chunkSize);
230     }
231     myNodes[ID] = node;
232     myNodeIDFactory->BindID(ID,node);
233     myInfo.myNbNodes++;
234     myModified = true;
235     this->adjustBoundingBox(x, y, z);
236     return node;
237   }else
238     return NULL;
239 }
240
241 ///////////////////////////////////////////////////////////////////////////////
242 /// create a Mesh0DElement and add it to the current Mesh
243 /// @return : The created Mesh0DElement
244 ///////////////////////////////////////////////////////////////////////////////
245 SMDS_Mesh0DElement* SMDS_Mesh::Add0DElementWithID(int idnode, int ID)
246 {
247   SMDS_MeshNode * node = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode);
248   if (!node) return NULL;
249   return SMDS_Mesh::Add0DElementWithID(node, ID);
250 }
251
252 ///////////////////////////////////////////////////////////////////////////////
253 /// create a Mesh0DElement and add it to the current Mesh
254 /// @return : The created Mesh0DElement
255 ///////////////////////////////////////////////////////////////////////////////
256 SMDS_Mesh0DElement* SMDS_Mesh::Add0DElement(const SMDS_MeshNode * node)
257 {
258   return SMDS_Mesh::Add0DElementWithID(node, myElementIDFactory->GetFreeID());
259 }
260
261 ///////////////////////////////////////////////////////////////////////////////
262 /// Create a new Mesh0DElement and at it to the mesh
263 /// @param idnode ID of the node
264 /// @param ID ID of the 0D element to create
265 /// @return The created 0D element or NULL if an element with this
266 ///         ID already exists or if input node is not found.
267 ///////////////////////////////////////////////////////////////////////////////
268 SMDS_Mesh0DElement* SMDS_Mesh::Add0DElementWithID(const SMDS_MeshNode * n, int ID)
269 {
270   if (!n) return 0;
271
272   if (Nb0DElements() % CHECKMEMORY_INTERVAL == 0) CheckMemory();
273   //MESSAGE("Add0DElementWithID" << ID)
274   SMDS_Mesh0DElement * el0d = new SMDS_Mesh0DElement(n);
275   if (myElementIDFactory->BindID(ID, el0d)) {
276     //SMDS_MeshNode *node = const_cast<SMDS_MeshNode*>(n);
277     //node->AddInverseElement(el0d);// --- fait avec BindID
278     adjustmyCellsCapacity(ID);
279     myCells[ID] = el0d;
280     myInfo.myNb0DElements++;
281     return el0d;
282   }
283
284   delete el0d;
285   return NULL;
286 }
287
288 ///////////////////////////////////////////////////////////////////////////////
289 /// create a Ball and add it to the current Mesh
290 /// @return : The created Ball
291 ///////////////////////////////////////////////////////////////////////////////
292 SMDS_BallElement* SMDS_Mesh::AddBallWithID(int idnode, double diameter, int ID)
293 {
294   SMDS_MeshNode * node = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode);
295   if (!node) return NULL;
296   return SMDS_Mesh::AddBallWithID(node, diameter, ID);
297 }
298
299 ///////////////////////////////////////////////////////////////////////////////
300 /// create a Ball and add it to the current Mesh
301 /// @return : The created Ball
302 ///////////////////////////////////////////////////////////////////////////////
303 SMDS_BallElement* SMDS_Mesh::AddBall(const SMDS_MeshNode * node, double diameter)
304 {
305   return SMDS_Mesh::AddBallWithID(node, diameter, myElementIDFactory->GetFreeID());
306 }
307
308 ///////////////////////////////////////////////////////////////////////////////
309 /// Create a new Ball and at it to the mesh
310 /// @param idnode ID of the node
311 //  @param diameter ball diameter
312 /// @param ID ID of the 0D element to create
313 /// @return The created 0D element or NULL if an element with this
314 ///         ID already exists or if input node is not found.
315 ///////////////////////////////////////////////////////////////////////////////
316 SMDS_BallElement* SMDS_Mesh::AddBallWithID(const SMDS_MeshNode * n, double diameter, int ID)
317 {
318   if (!n) return 0;
319
320   if (NbBalls() % CHECKMEMORY_INTERVAL == 0) CheckMemory();
321
322   SMDS_BallElement *ball = myBallPool->getNew();
323   ball->init(n->getVtkId(), diameter, this);
324   if (!this->registerElement(ID,ball))
325   {
326     this->myGrid->GetCellTypesArray()->SetValue(ball->getVtkId(), VTK_EMPTY_CELL);
327     myBallPool->destroy(ball);
328     return 0;
329   }
330   adjustmyCellsCapacity(ID);
331   myCells[ID] = ball;
332   myInfo.myNbBalls++;
333   return ball;
334 }
335
336 ///////////////////////////////////////////////////////////////////////////////
337 /// create a MeshEdge and add it to the current Mesh
338 /// @return : The created MeshEdge
339 ///////////////////////////////////////////////////////////////////////////////
340
341 SMDS_MeshEdge* SMDS_Mesh::AddEdgeWithID(int idnode1, int idnode2, int ID)
342 {
343   SMDS_MeshNode * node1 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode1);
344   SMDS_MeshNode * node2 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode2);
345   if(!node1 || !node2) return NULL;
346   return SMDS_Mesh::AddEdgeWithID(node1, node2, ID);
347 }
348
349 ///////////////////////////////////////////////////////////////////////////////
350 /// create a MeshEdge and add it to the current Mesh
351 /// @return : The created MeshEdge
352 ///////////////////////////////////////////////////////////////////////////////
353
354 SMDS_MeshEdge* SMDS_Mesh::AddEdge(const SMDS_MeshNode * node1,
355                                   const SMDS_MeshNode * node2)
356 {
357   return SMDS_Mesh::AddEdgeWithID(node1, node2, myElementIDFactory->GetFreeID());
358 }
359
360 ///////////////////////////////////////////////////////////////////////////////
361 /// Create a new edge and at it to the mesh
362 /// @param idnode1 ID of the first node
363 /// @param idnode2 ID of the second node
364 /// @param ID ID of the edge to create
365 /// @return The created edge or NULL if an element with this ID already exists or
366 /// if input nodes are not found.
367 ///////////////////////////////////////////////////////////////////////////////
368
369 SMDS_MeshEdge* SMDS_Mesh::AddEdgeWithID(const SMDS_MeshNode * n1,
370                                         const SMDS_MeshNode * n2,
371                                         int ID)
372 {
373   if ( !n1 || !n2 ) return 0;
374   SMDS_MeshEdge * edge = 0;
375
376   // --- retreive nodes ID
377   vector<vtkIdType> nodeIds;
378   nodeIds.clear();
379   nodeIds.push_back(n1->getVtkId());
380   nodeIds.push_back(n2->getVtkId());
381
382   SMDS_VtkEdge *edgevtk = myEdgePool->getNew();
383   edgevtk->init(nodeIds, this);
384   if (!this->registerElement(ID,edgevtk))
385     {
386       this->myGrid->GetCellTypesArray()->SetValue(edgevtk->getVtkId(), VTK_EMPTY_CELL);
387       myEdgePool->destroy(edgevtk);
388       return 0;
389     }
390   edge = edgevtk;
391   adjustmyCellsCapacity(ID);
392   myCells[ID] = edge;
393   myInfo.myNbEdges++;
394
395 //  if (edge && !registerElement(ID, edge))
396 //  {
397 //    RemoveElement(edge, false);
398 //    edge = NULL;
399 //  }
400   return edge;
401 }
402
403 ///////////////////////////////////////////////////////////////////////////////
404 /// Add a triangle defined by its nodes. An ID is automatically affected to the
405 /// Created face
406 ///////////////////////////////////////////////////////////////////////////////
407
408 SMDS_MeshFace* SMDS_Mesh::AddFace(const SMDS_MeshNode * n1,
409                                   const SMDS_MeshNode * n2,
410                                   const SMDS_MeshNode * n3)
411 {
412   return SMDS_Mesh::AddFaceWithID(n1,n2,n3, myElementIDFactory->GetFreeID());
413 }
414
415 ///////////////////////////////////////////////////////////////////////////////
416 /// Add a triangle defined by its nodes IDs
417 ///////////////////////////////////////////////////////////////////////////////
418
419 SMDS_MeshFace* SMDS_Mesh::AddFaceWithID(int idnode1, int idnode2, int idnode3, int ID)
420 {
421   SMDS_MeshNode * node1 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode1);
422   SMDS_MeshNode * node2 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode2);
423   SMDS_MeshNode * node3 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode3);
424   if(!node1 || !node2 || !node3) return NULL;
425   return SMDS_Mesh::AddFaceWithID(node1, node2, node3, ID);
426 }
427
428 ///////////////////////////////////////////////////////////////////////////////
429 /// Add a triangle defined by its nodes
430 ///////////////////////////////////////////////////////////////////////////////
431
432 SMDS_MeshFace* SMDS_Mesh::AddFaceWithID(const SMDS_MeshNode * n1,
433                                         const SMDS_MeshNode * n2,
434                                         const SMDS_MeshNode * n3,
435                                         int ID)
436 {
437   //MESSAGE("AddFaceWithID " << ID)
438   SMDS_MeshFace * face=createTriangle(n1, n2, n3, ID);
439
440 //  if (face && !registerElement(ID, face)) {
441 //    RemoveElement(face, false);
442 //    face = NULL;
443 //  }
444   return face;
445 }
446
447 ///////////////////////////////////////////////////////////////////////////////
448 /// Add a quadrangle defined by its nodes. An ID is automatically affected to the
449 /// created face
450 ///////////////////////////////////////////////////////////////////////////////
451
452 SMDS_MeshFace* SMDS_Mesh::AddFace(const SMDS_MeshNode * n1,
453                                   const SMDS_MeshNode * n2,
454                                   const SMDS_MeshNode * n3,
455                                   const SMDS_MeshNode * n4)
456 {
457   return SMDS_Mesh::AddFaceWithID(n1,n2,n3, n4, myElementIDFactory->GetFreeID());
458 }
459
460 ///////////////////////////////////////////////////////////////////////////////
461 /// Add a quadrangle defined by its nodes IDs
462 ///////////////////////////////////////////////////////////////////////////////
463
464 SMDS_MeshFace* SMDS_Mesh::AddFaceWithID(int idnode1,
465                                         int idnode2,
466                                         int idnode3,
467                                         int idnode4,
468                                         int ID)
469 {
470   SMDS_MeshNode *node1, *node2, *node3, *node4;
471   node1 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode1);
472   node2 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode2);
473   node3 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode3);
474   node4 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode4);
475   if(!node1 || !node2 || !node3 || !node4) return NULL;
476   return SMDS_Mesh::AddFaceWithID(node1, node2, node3, node4, ID);
477 }
478
479 ///////////////////////////////////////////////////////////////////////////////
480 /// Add a quadrangle defined by its nodes
481 ///////////////////////////////////////////////////////////////////////////////
482
483 SMDS_MeshFace* SMDS_Mesh::AddFaceWithID(const SMDS_MeshNode * n1,
484                                         const SMDS_MeshNode * n2,
485                                         const SMDS_MeshNode * n3,
486                                         const SMDS_MeshNode * n4,
487                                         int ID)
488 {
489   //MESSAGE("AddFaceWithID " << ID);
490   SMDS_MeshFace * face=createQuadrangle(n1, n2, n3, n4, ID);
491
492 //  if (face && !registerElement(ID, face)) {
493 //    RemoveElement(face, false);
494 //    face = NULL;
495 //  }
496   return face;
497 }
498
499 ///////////////////////////////////////////////////////////////////////////////
500 /// Add a triangle defined by its edges. An ID is automatically assigned to the
501 /// Created face
502 ///////////////////////////////////////////////////////////////////////////////
503
504 SMDS_MeshFace* SMDS_Mesh::AddFace(const SMDS_MeshEdge * e1,
505                                   const SMDS_MeshEdge * e2,
506                                   const SMDS_MeshEdge * e3)
507 {
508   if (!hasConstructionEdges())
509     return NULL;
510      //MESSAGE("AddFaceWithID");
511  return AddFaceWithID(e1,e2,e3, myElementIDFactory->GetFreeID());
512 }
513
514 ///////////////////////////////////////////////////////////////////////////////
515 /// Add a triangle defined by its edges
516 ///////////////////////////////////////////////////////////////////////////////
517
518 SMDS_MeshFace* SMDS_Mesh::AddFaceWithID(const SMDS_MeshEdge * e1,
519                                         const SMDS_MeshEdge * e2,
520                                         const SMDS_MeshEdge * e3,
521                                         int ID)
522 {
523   if (!hasConstructionEdges())
524     return NULL;
525   if ( !e1 || !e2 || !e3 ) return 0;
526
527   if ( NbFaces() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
528
529   SMDS_MeshFace * face = new SMDS_FaceOfEdges(e1,e2,e3);
530   adjustmyCellsCapacity(ID);
531   myCells[ID] = face;
532   myInfo.myNbTriangles++;
533
534   if (!registerElement(ID, face)) {
535     registerElement(myElementIDFactory->GetFreeID(), face);
536     //RemoveElement(face, false);
537     //face = NULL;
538   }
539   return face;
540 }
541
542 ///////////////////////////////////////////////////////////////////////////////
543 /// Add a quadrangle defined by its edges. An ID is automatically assigned to the
544 /// Created face
545 ///////////////////////////////////////////////////////////////////////////////
546
547 SMDS_MeshFace* SMDS_Mesh::AddFace(const SMDS_MeshEdge * e1,
548                                   const SMDS_MeshEdge * e2,
549                                   const SMDS_MeshEdge * e3,
550                                   const SMDS_MeshEdge * e4)
551 {
552   if (!hasConstructionEdges())
553     return NULL;
554  return AddFaceWithID(e1,e2,e3,e4, myElementIDFactory->GetFreeID());
555 }
556
557 ///////////////////////////////////////////////////////////////////////////////
558 /// Add a quadrangle defined by its edges
559 ///////////////////////////////////////////////////////////////////////////////
560
561 SMDS_MeshFace* SMDS_Mesh::AddFaceWithID(const SMDS_MeshEdge * e1,
562                                         const SMDS_MeshEdge * e2,
563                                         const SMDS_MeshEdge * e3,
564                                         const SMDS_MeshEdge * e4,
565                                         int ID)
566 {
567   if (!hasConstructionEdges())
568     return NULL;
569   if ( !e1 || !e2 || !e3 || !e4 ) return 0;
570   if ( NbFaces() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
571   SMDS_MeshFace * face = new SMDS_FaceOfEdges(e1,e2,e3,e4);
572   adjustmyCellsCapacity(ID);
573   myCells[ID] = face;
574   myInfo.myNbQuadrangles++;
575
576   if (!registerElement(ID, face))
577   {
578     registerElement(myElementIDFactory->GetFreeID(), face);
579     //RemoveElement(face, false);
580     //face = NULL;
581   }
582   return face;
583 }
584
585 ///////////////////////////////////////////////////////////////////////////////
586 ///Create a new tetrahedron and add it to the mesh.
587 ///@return The created tetrahedron
588 ///////////////////////////////////////////////////////////////////////////////
589
590 SMDS_MeshVolume* SMDS_Mesh::AddVolume(const SMDS_MeshNode * n1,
591                                       const SMDS_MeshNode * n2,
592                                       const SMDS_MeshNode * n3,
593                                       const SMDS_MeshNode * n4)
594 {
595   int ID = myElementIDFactory->GetFreeID();
596   SMDS_MeshVolume * v = SMDS_Mesh::AddVolumeWithID(n1, n2, n3, n4, ID);
597   if(v==NULL) myElementIDFactory->ReleaseID(ID);
598   return v;
599 }
600
601 ///////////////////////////////////////////////////////////////////////////////
602 ///Create a new tetrahedron and add it to the mesh.
603 ///@param ID The ID of the new volume
604 ///@return The created tetrahedron or NULL if an element with this ID already exists
605 ///or if input nodes are not found.
606 ///////////////////////////////////////////////////////////////////////////////
607
608 SMDS_MeshVolume * SMDS_Mesh::AddVolumeWithID(int idnode1,
609                                              int idnode2,
610                                              int idnode3,
611                                              int idnode4,
612                                              int ID)
613 {
614   SMDS_MeshNode *node1, *node2, *node3, *node4;
615   node1 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode1);
616   node2 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode2);
617   node3 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode3);
618   node4 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode4);
619   if(!node1 || !node2 || !node3 || !node4) return NULL;
620   return SMDS_Mesh::AddVolumeWithID(node1, node2, node3, node4, ID);
621 }
622
623 ///////////////////////////////////////////////////////////////////////////////
624 ///Create a new tetrahedron and add it to the mesh.
625 ///@param ID The ID of the new volume
626 ///@return The created tetrahedron
627 ///////////////////////////////////////////////////////////////////////////////
628
629 SMDS_MeshVolume* SMDS_Mesh::AddVolumeWithID(const SMDS_MeshNode * n1,
630                                             const SMDS_MeshNode * n2,
631                                             const SMDS_MeshNode * n3,
632                                             const SMDS_MeshNode * n4,
633                                             int ID)
634 {
635   SMDS_MeshVolume* volume = 0;
636   if ( !n1 || !n2 || !n3 || !n4) return volume;
637   if ( NbVolumes() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
638   if(hasConstructionFaces()) {
639     SMDS_MeshFace * f1=FindFaceOrCreate(n1,n2,n3);
640     SMDS_MeshFace * f2=FindFaceOrCreate(n1,n2,n4);
641     SMDS_MeshFace * f3=FindFaceOrCreate(n1,n3,n4);
642     SMDS_MeshFace * f4=FindFaceOrCreate(n2,n3,n4);
643     volume=new SMDS_VolumeOfFaces(f1,f2,f3,f4);
644     adjustmyCellsCapacity(ID);
645     myCells[ID] = volume;
646     myInfo.myNbTetras++;
647   }
648   else if(hasConstructionEdges()) {
649     return NULL;
650   }
651   else {
652     // --- retrieve nodes ID
653     myNodeIds.resize(4);
654     myNodeIds[0] = n1->getVtkId();
655     myNodeIds[1] = n3->getVtkId(); // order SMDS-->VTK
656     myNodeIds[2] = n2->getVtkId();
657     myNodeIds[3] = n4->getVtkId();
658
659     SMDS_VtkVolume *volvtk = myVolumePool->getNew();
660     volvtk->init(myNodeIds, this);
661     if (!this->registerElement(ID,volvtk))
662     {
663       this->myGrid->GetCellTypesArray()->SetValue(volvtk->getVtkId(), VTK_EMPTY_CELL);
664       myVolumePool->destroy(volvtk);
665       return 0;
666     }
667     volume = volvtk;
668     adjustmyCellsCapacity(ID);
669     myCells[ID] = volume;
670     myInfo.myNbTetras++;
671   }
672
673   //  if (!registerElement(ID, volume)) {
674   //    RemoveElement(volume, false);
675   //    volume = NULL;
676   //  }
677   return volume;
678 }
679
680 ///////////////////////////////////////////////////////////////////////////////
681 ///Create a new pyramid and add it to the mesh.
682 ///Nodes 1,2,3 and 4 define the base of the pyramid
683 ///@return The created pyramid
684 ///////////////////////////////////////////////////////////////////////////////
685
686 SMDS_MeshVolume* SMDS_Mesh::AddVolume(const SMDS_MeshNode * n1,
687                                       const SMDS_MeshNode * n2,
688                                       const SMDS_MeshNode * n3,
689                                       const SMDS_MeshNode * n4,
690                                       const SMDS_MeshNode * n5)
691 {
692   int ID = myElementIDFactory->GetFreeID();
693   SMDS_MeshVolume * v = SMDS_Mesh::AddVolumeWithID(n1, n2, n3, n4, n5, ID);
694   if(v==NULL) myElementIDFactory->ReleaseID(ID);
695   return v;
696 }
697
698 ///////////////////////////////////////////////////////////////////////////////
699 ///Create a new pyramid and add it to the mesh.
700 ///Nodes 1,2,3 and 4 define the base of the pyramid
701 ///@param ID The ID of the new volume
702 ///@return The created pyramid or NULL if an element with this ID already exists
703 ///or if input nodes are not found.
704 ///////////////////////////////////////////////////////////////////////////////
705
706 SMDS_MeshVolume * SMDS_Mesh::AddVolumeWithID(int idnode1,
707                                              int idnode2,
708                                              int idnode3,
709                                              int idnode4,
710                                              int idnode5,
711                                              int ID)
712 {
713   SMDS_MeshNode *node1, *node2, *node3, *node4, *node5;
714   node1 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode1);
715   node2 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode2);
716   node3 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode3);
717   node4 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode4);
718   node5 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode5);
719   if(!node1 || !node2 || !node3 || !node4 || !node5) return NULL;
720   return SMDS_Mesh::AddVolumeWithID(node1, node2, node3, node4, node5, ID);
721 }
722
723 ///////////////////////////////////////////////////////////////////////////////
724 ///Create a new pyramid and add it to the mesh.
725 ///Nodes 1,2,3 and 4 define the base of the pyramid
726 ///@param ID The ID of the new volume
727 ///@return The created pyramid
728 ///////////////////////////////////////////////////////////////////////////////
729
730 SMDS_MeshVolume* SMDS_Mesh::AddVolumeWithID(const SMDS_MeshNode * n1,
731                                             const SMDS_MeshNode * n2,
732                                             const SMDS_MeshNode * n3,
733                                             const SMDS_MeshNode * n4,
734                                             const SMDS_MeshNode * n5,
735                                             int ID)
736 {
737   SMDS_MeshVolume* volume = 0;
738   if ( !n1 || !n2 || !n3 || !n4 || !n5) return volume;
739   if ( NbVolumes() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
740   if(hasConstructionFaces()) {
741     SMDS_MeshFace * f1=FindFaceOrCreate(n1,n2,n3,n4);
742     SMDS_MeshFace * f2=FindFaceOrCreate(n1,n2,n5);
743     SMDS_MeshFace * f3=FindFaceOrCreate(n2,n3,n5);
744     SMDS_MeshFace * f4=FindFaceOrCreate(n3,n4,n5);
745     volume=new SMDS_VolumeOfFaces(f1,f2,f3,f4);
746     adjustmyCellsCapacity(ID);
747     myCells[ID] = volume;
748     myInfo.myNbPyramids++;
749   }
750   else if(hasConstructionEdges()) {
751     return NULL;
752   }
753   else {
754     // --- retrieve nodes ID
755     myNodeIds.resize(5);
756     myNodeIds[0] = n1->getVtkId();
757     myNodeIds[1] = n4->getVtkId();
758     myNodeIds[2] = n3->getVtkId();
759     myNodeIds[3] = n2->getVtkId();
760     myNodeIds[4] = n5->getVtkId();
761
762     SMDS_VtkVolume *volvtk = myVolumePool->getNew();
763     volvtk->init(myNodeIds, this);
764     if (!this->registerElement(ID,volvtk))
765     {
766       this->myGrid->GetCellTypesArray()->SetValue(volvtk->getVtkId(), VTK_EMPTY_CELL);
767       myVolumePool->destroy(volvtk);
768       return 0;
769     }
770     volume = volvtk;
771     adjustmyCellsCapacity(ID);
772     myCells[ID] = volume;
773     myInfo.myNbPyramids++;
774   }
775
776   //  if (!registerElement(ID, volume)) {
777   //    RemoveElement(volume, false);
778   //    volume = NULL;
779   //  }
780   return volume;
781 }
782
783 ///////////////////////////////////////////////////////////////////////////////
784 ///Create a new prism and add it to the mesh.
785 ///Nodes 1,2,3 is a triangle and 1,2,5,4 a quadrangle.
786 ///@return The created prism
787 ///////////////////////////////////////////////////////////////////////////////
788
789 SMDS_MeshVolume* SMDS_Mesh::AddVolume(const SMDS_MeshNode * n1,
790                                       const SMDS_MeshNode * n2,
791                                       const SMDS_MeshNode * n3,
792                                       const SMDS_MeshNode * n4,
793                                       const SMDS_MeshNode * n5,
794                                       const SMDS_MeshNode * n6)
795 {
796   int ID = myElementIDFactory->GetFreeID();
797   SMDS_MeshVolume * v = SMDS_Mesh::AddVolumeWithID(n1, n2, n3, n4, n5, n6, ID);
798   if(v==NULL) myElementIDFactory->ReleaseID(ID);
799   return v;
800 }
801
802 ///////////////////////////////////////////////////////////////////////////////
803 ///Create a new prism and add it to the mesh.
804 ///Nodes 1,2,3 is a triangle and 1,2,5,4 a quadrangle.
805 ///@param ID The ID of the new volume
806 ///@return The created prism or NULL if an element with this ID already exists
807 ///or if input nodes are not found.
808 ///////////////////////////////////////////////////////////////////////////////
809
810 SMDS_MeshVolume * SMDS_Mesh::AddVolumeWithID(int idnode1,
811                                              int idnode2,
812                                              int idnode3,
813                                              int idnode4,
814                                              int idnode5,
815                                              int idnode6,
816                                              int ID)
817 {
818   SMDS_MeshNode *node1, *node2, *node3, *node4, *node5, *node6;
819   node1 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode1);
820   node2 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode2);
821   node3 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode3);
822   node4 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode4);
823   node5 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode5);
824   node6 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode6);
825   if(!node1 || !node2 || !node3 || !node4 || !node5 || !node6) return NULL;
826   return SMDS_Mesh::AddVolumeWithID(node1, node2, node3, node4, node5, node6, ID);
827 }
828
829 ///////////////////////////////////////////////////////////////////////////////
830 ///Create a new prism and add it to the mesh.
831 ///Nodes 1,2,3 is a triangle and 1,2,5,4 a quadrangle.
832 ///@param ID The ID of the new volume
833 ///@return The created prism
834 ///////////////////////////////////////////////////////////////////////////////
835
836 SMDS_MeshVolume* SMDS_Mesh::AddVolumeWithID(const SMDS_MeshNode * n1,
837                                             const SMDS_MeshNode * n2,
838                                             const SMDS_MeshNode * n3,
839                                             const SMDS_MeshNode * n4,
840                                             const SMDS_MeshNode * n5,
841                                             const SMDS_MeshNode * n6,
842                                             int ID)
843 {
844   SMDS_MeshVolume* volume = 0;
845   if ( !n1 || !n2 || !n3 || !n4 || !n5 || !n6) return volume;
846   if ( NbVolumes() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
847   if(hasConstructionFaces()) {
848     SMDS_MeshFace * f1=FindFaceOrCreate(n1,n2,n3);
849     SMDS_MeshFace * f2=FindFaceOrCreate(n4,n5,n6);
850     SMDS_MeshFace * f3=FindFaceOrCreate(n1,n4,n5,n2);
851     SMDS_MeshFace * f4=FindFaceOrCreate(n2,n5,n6,n3);
852     SMDS_MeshFace * f5=FindFaceOrCreate(n3,n6,n4,n1);
853     volume=new SMDS_VolumeOfFaces(f1,f2,f3,f4,f5);
854     adjustmyCellsCapacity(ID);
855     myCells[ID] = volume;
856     myInfo.myNbPrisms++;
857   }
858   else if(hasConstructionEdges()) {
859     return NULL;
860   }
861   else {
862     // --- retrieve nodes ID
863     myNodeIds.resize(6);
864     myNodeIds[0] = n1->getVtkId();
865     myNodeIds[1] = n2->getVtkId();
866     myNodeIds[2] = n3->getVtkId();
867     myNodeIds[3] = n4->getVtkId();
868     myNodeIds[4] = n5->getVtkId();
869     myNodeIds[5] = n6->getVtkId();
870
871     SMDS_VtkVolume *volvtk = myVolumePool->getNew();
872     volvtk->init(myNodeIds, this);
873     if (!this->registerElement(ID,volvtk))
874     {
875       this->myGrid->GetCellTypesArray()->SetValue(volvtk->getVtkId(), VTK_EMPTY_CELL);
876       myVolumePool->destroy(volvtk);
877       return 0;
878     }
879     volume = volvtk;
880     adjustmyCellsCapacity(ID);
881     myCells[ID] = volume;
882     myInfo.myNbPrisms++;
883   }
884
885   //  if (!registerElement(ID, volume)) {
886   //    RemoveElement(volume, false);
887   //    volume = NULL;
888   //  }
889   return volume;
890 }
891
892 ///////////////////////////////////////////////////////////////////////////////
893 ///Create a new hexagonal prism and add it to the mesh.
894 ///@return The created prism
895 ///////////////////////////////////////////////////////////////////////////////
896
897 SMDS_MeshVolume* SMDS_Mesh::AddVolume(const SMDS_MeshNode * n1,
898                                       const SMDS_MeshNode * n2,
899                                       const SMDS_MeshNode * n3,
900                                       const SMDS_MeshNode * n4,
901                                       const SMDS_MeshNode * n5,
902                                       const SMDS_MeshNode * n6,
903                                       const SMDS_MeshNode * n7,
904                                       const SMDS_MeshNode * n8,
905                                       const SMDS_MeshNode * n9,
906                                       const SMDS_MeshNode * n10,
907                                       const SMDS_MeshNode * n11,
908                                       const SMDS_MeshNode * n12)
909 {
910   int ID = myElementIDFactory->GetFreeID();
911   SMDS_MeshVolume * v = SMDS_Mesh::AddVolumeWithID(n1, n2, n3, n4, n5, n6,
912                                                    n7, n8, n9, n10, n11, n12,
913                                                    ID);
914   if(v==NULL) myElementIDFactory->ReleaseID(ID);
915   return v;
916 }
917
918 ///////////////////////////////////////////////////////////////////////////////
919 ///Create a new hexagonal prism and add it to the mesh.
920 ///@param ID The ID of the new volume
921 ///@return The created prism or NULL if an element with this ID already exists
922 ///or if input nodes are not found.
923 ///////////////////////////////////////////////////////////////////////////////
924
925 SMDS_MeshVolume * SMDS_Mesh::AddVolumeWithID(int idnode1,
926                                              int idnode2,
927                                              int idnode3,
928                                              int idnode4,
929                                              int idnode5,
930                                              int idnode6,
931                                              int idnode7,
932                                              int idnode8,
933                                              int idnode9,
934                                              int idnode10,
935                                              int idnode11,
936                                              int idnode12,
937                                              int ID)
938 {
939   SMDS_MeshNode *node1 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode1);
940   SMDS_MeshNode *node2 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode2);
941   SMDS_MeshNode *node3 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode3);
942   SMDS_MeshNode *node4 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode4);
943   SMDS_MeshNode *node5 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode5);
944   SMDS_MeshNode *node6 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode6);
945   SMDS_MeshNode *node7 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode7);
946   SMDS_MeshNode *node8 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode8);
947   SMDS_MeshNode *node9 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode9);
948   SMDS_MeshNode *node10 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode10);
949   SMDS_MeshNode *node11 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode11);
950   SMDS_MeshNode *node12 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode12);
951   return SMDS_Mesh::AddVolumeWithID(node1, node2, node3, node4, node5, node6,
952                                     node7, node8, node9, node10, node11, node12,
953                                     ID);
954 }
955
956 ///////////////////////////////////////////////////////////////////////////////
957 ///Create a new hexagonal prism and add it to the mesh.
958 ///@param ID The ID of the new volume
959 ///@return The created prism
960 ///////////////////////////////////////////////////////////////////////////////
961
962 SMDS_MeshVolume* SMDS_Mesh::AddVolumeWithID(const SMDS_MeshNode * n1,
963                                             const SMDS_MeshNode * n2,
964                                             const SMDS_MeshNode * n3,
965                                             const SMDS_MeshNode * n4,
966                                             const SMDS_MeshNode * n5,
967                                             const SMDS_MeshNode * n6,
968                                             const SMDS_MeshNode * n7,
969                                             const SMDS_MeshNode * n8,
970                                             const SMDS_MeshNode * n9,
971                                             const SMDS_MeshNode * n10,
972                                             const SMDS_MeshNode * n11,
973                                             const SMDS_MeshNode * n12,
974                                             int ID)
975 {
976   SMDS_MeshVolume* volume = 0;
977   if(!n1 || !n2 || !n3 || !n4 || !n5 || !n6 ||
978      !n7 || !n8 || !n9 || !n10 || !n11 || !n12 )
979     return volume;
980   if ( NbVolumes() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
981   if(hasConstructionFaces()) {
982     return NULL;
983   }
984   else if(hasConstructionEdges()) {
985     return NULL;
986   }
987   else {
988     // --- retrieve nodes ID
989     myNodeIds.resize(12);
990     myNodeIds[0] = n1->getVtkId();
991     myNodeIds[1] = n6->getVtkId();
992     myNodeIds[2] = n5->getVtkId();
993     myNodeIds[3] = n4->getVtkId();
994     myNodeIds[4] = n3->getVtkId();
995     myNodeIds[5] = n2->getVtkId();
996
997     myNodeIds[6] = n7->getVtkId();
998     myNodeIds[7] = n12->getVtkId();
999     myNodeIds[8] = n11->getVtkId();
1000     myNodeIds[9] = n10->getVtkId();
1001     myNodeIds[10] = n9->getVtkId();
1002     myNodeIds[11] = n8->getVtkId();
1003
1004     SMDS_VtkVolume *volvtk = myVolumePool->getNew();
1005     volvtk->init(myNodeIds, this);
1006     if (!this->registerElement(ID,volvtk))
1007     {
1008       this->myGrid->GetCellTypesArray()->SetValue(volvtk->getVtkId(), VTK_EMPTY_CELL);
1009       myVolumePool->destroy(volvtk);
1010       return 0;
1011     }
1012     volume = volvtk;
1013     adjustmyCellsCapacity(ID);
1014     myCells[ID] = volume;
1015     myInfo.myNbHexPrism++;
1016   }
1017
1018   return volume;
1019 }
1020
1021 ///////////////////////////////////////////////////////////////////////////////
1022 ///Create a new hexahedron and add it to the mesh.
1023 ///Nodes 1,2,3,4 and 5,6,7,8 are quadrangle and 5,1 and 7,3 are an edges.
1024 ///@return The created hexahedron
1025 ///////////////////////////////////////////////////////////////////////////////
1026
1027 SMDS_MeshVolume* SMDS_Mesh::AddVolume(const SMDS_MeshNode * n1,
1028                                       const SMDS_MeshNode * n2,
1029                                       const SMDS_MeshNode * n3,
1030                                       const SMDS_MeshNode * n4,
1031                                       const SMDS_MeshNode * n5,
1032                                       const SMDS_MeshNode * n6,
1033                                       const SMDS_MeshNode * n7,
1034                                       const SMDS_MeshNode * n8)
1035 {
1036   int ID = myElementIDFactory->GetFreeID();
1037  SMDS_MeshVolume * v = SMDS_Mesh::AddVolumeWithID(n1, n2, n3, n4, n5, n6, n7, n8, ID);
1038   if(v==NULL) myElementIDFactory->ReleaseID(ID);
1039   return v;
1040 }
1041
1042 ///////////////////////////////////////////////////////////////////////////////
1043 ///Create a new hexahedron and add it to the mesh.
1044 ///Nodes 1,2,3,4 and 5,6,7,8 are quadrangle and 5,1 and 7,3 are an edges.
1045 ///@param ID The ID of the new volume
1046 ///@return The created hexahedron or NULL if an element with this ID already
1047 ///exists or if input nodes are not found.
1048 ///////////////////////////////////////////////////////////////////////////////
1049
1050 SMDS_MeshVolume * SMDS_Mesh::AddVolumeWithID(int idnode1,
1051                                              int idnode2,
1052                                              int idnode3,
1053                                              int idnode4,
1054                                              int idnode5,
1055                                              int idnode6,
1056                                              int idnode7,
1057                                              int idnode8,
1058                                              int ID)
1059 {
1060   SMDS_MeshNode *node1, *node2, *node3, *node4, *node5, *node6, *node7, *node8;
1061   node1 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode1);
1062   node2 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode2);
1063   node3 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode3);
1064   node4 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode4);
1065   node5 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode5);
1066   node6 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode6);
1067   node7 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode7);
1068   node8 = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(idnode8);
1069   if(!node1 || !node2 || !node3 || !node4 || !node5 || !node6 || !node7 || !node8)
1070     return NULL;
1071   return SMDS_Mesh::AddVolumeWithID(node1, node2, node3, node4, node5, node6,
1072                                     node7, node8, ID);
1073 }
1074
1075 ///////////////////////////////////////////////////////////////////////////////
1076 ///Create a new hexahedron and add it to the mesh.
1077 ///Nodes 1,2,3,4 and 5,6,7,8 are quadrangle and 5,1 and 7,3 are an edges.
1078 ///@param ID The ID of the new volume
1079 ///@return The created prism or NULL if an element with this ID already exists
1080 ///or if input nodes are not found.
1081 ///////////////////////////////////////////////////////////////////////////////
1082
1083 SMDS_MeshVolume* SMDS_Mesh::AddVolumeWithID(const SMDS_MeshNode * n1,
1084                                             const SMDS_MeshNode * n2,
1085                                             const SMDS_MeshNode * n3,
1086                                             const SMDS_MeshNode * n4,
1087                                             const SMDS_MeshNode * n5,
1088                                             const SMDS_MeshNode * n6,
1089                                             const SMDS_MeshNode * n7,
1090                                             const SMDS_MeshNode * n8,
1091                                             int ID)
1092 {
1093   SMDS_MeshVolume* volume = 0;
1094   if ( !n1 || !n2 || !n3 || !n4 || !n5 || !n6 || !n7 || !n8) return volume;
1095   if ( NbVolumes() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
1096   if(hasConstructionFaces()) {
1097     SMDS_MeshFace * f1=FindFaceOrCreate(n1,n2,n3,n4);
1098     SMDS_MeshFace * f2=FindFaceOrCreate(n5,n6,n7,n8);
1099     SMDS_MeshFace * f3=FindFaceOrCreate(n1,n4,n8,n5);
1100     SMDS_MeshFace * f4=FindFaceOrCreate(n1,n2,n6,n5);
1101     SMDS_MeshFace * f5=FindFaceOrCreate(n2,n3,n7,n6);
1102     SMDS_MeshFace * f6=FindFaceOrCreate(n3,n4,n8,n7);
1103     volume=new SMDS_VolumeOfFaces(f1,f2,f3,f4,f5,f6);
1104     adjustmyCellsCapacity(ID);
1105     myCells[ID] = volume;
1106     myInfo.myNbHexas++;
1107   }
1108   else if(hasConstructionEdges()) {
1109     return NULL;
1110   }
1111   else {
1112     // --- retrieve nodes ID
1113     myNodeIds.resize(8);
1114     myNodeIds[0] = n1->getVtkId();
1115     myNodeIds[1] = n4->getVtkId();
1116     myNodeIds[2] = n3->getVtkId();
1117     myNodeIds[3] = n2->getVtkId();
1118     myNodeIds[4] = n5->getVtkId();
1119     myNodeIds[5] = n8->getVtkId();
1120     myNodeIds[6] = n7->getVtkId();
1121     myNodeIds[7] = n6->getVtkId();
1122
1123     SMDS_VtkVolume *volvtk = myVolumePool->getNew();
1124     volvtk->init(myNodeIds, this);
1125     if (!this->registerElement(ID,volvtk))
1126     {
1127       this->myGrid->GetCellTypesArray()->SetValue(volvtk->getVtkId(), VTK_EMPTY_CELL);
1128       myVolumePool->destroy(volvtk);
1129       return 0;
1130     }
1131     volume = volvtk;
1132     adjustmyCellsCapacity(ID);
1133     myCells[ID] = volume;
1134     myInfo.myNbHexas++;
1135   }
1136
1137   //  if (!registerElement(ID, volume)) {
1138   //    RemoveElement(volume, false);
1139   //    volume = NULL;
1140   //  }
1141   return volume;
1142 }
1143
1144 ///////////////////////////////////////////////////////////////////////////////
1145 ///Create a new tetrahedron defined by its faces and add it to the mesh.
1146 ///@return The created tetrahedron
1147 ///////////////////////////////////////////////////////////////////////////////
1148
1149 SMDS_MeshVolume* SMDS_Mesh::AddVolume(const SMDS_MeshFace * f1,
1150                                       const SMDS_MeshFace * f2,
1151                                       const SMDS_MeshFace * f3,
1152                                       const SMDS_MeshFace * f4)
1153 {
1154   if (!hasConstructionFaces())
1155     return NULL;
1156   return AddVolumeWithID(f1,f2,f3,f4, myElementIDFactory->GetFreeID());
1157 }
1158
1159 ///////////////////////////////////////////////////////////////////////////////
1160 ///Create a new tetrahedron defined by its faces and add it to the mesh.
1161 ///@param ID The ID of the new volume
1162 ///@return The created tetrahedron
1163 ///////////////////////////////////////////////////////////////////////////////
1164
1165 SMDS_MeshVolume* SMDS_Mesh::AddVolumeWithID(const SMDS_MeshFace * f1,
1166                                             const SMDS_MeshFace * f2,
1167                                             const SMDS_MeshFace * f3,
1168                                             const SMDS_MeshFace * f4,
1169                                             int ID)
1170 {
1171   if (!hasConstructionFaces())
1172     return NULL;
1173   if ( !f1 || !f2 || !f3 || !f4) return 0;
1174   if ( NbVolumes() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
1175   SMDS_MeshVolume * volume = new SMDS_VolumeOfFaces(f1,f2,f3,f4);
1176   adjustmyCellsCapacity(ID);
1177   myCells[ID] = volume;
1178   myInfo.myNbTetras++;
1179
1180   if (!registerElement(ID, volume)) {
1181     registerElement(myElementIDFactory->GetFreeID(), volume);
1182     //RemoveElement(volume, false);
1183     //volume = NULL;
1184   }
1185   return volume;
1186 }
1187
1188 ///////////////////////////////////////////////////////////////////////////////
1189 ///Create a new pyramid defined by its faces and add it to the mesh.
1190 ///@return The created pyramid
1191 ///////////////////////////////////////////////////////////////////////////////
1192
1193 SMDS_MeshVolume* SMDS_Mesh::AddVolume(const SMDS_MeshFace * f1,
1194                                       const SMDS_MeshFace * f2,
1195                                       const SMDS_MeshFace * f3,
1196                                       const SMDS_MeshFace * f4,
1197                                       const SMDS_MeshFace * f5)
1198 {
1199  if (!hasConstructionFaces())
1200     return NULL;
1201   return AddVolumeWithID(f1,f2,f3,f4,f5, myElementIDFactory->GetFreeID());
1202 }
1203
1204 ///////////////////////////////////////////////////////////////////////////////
1205 ///Create a new pyramid defined by its faces and add it to the mesh.
1206 ///@param ID The ID of the new volume
1207 ///@return The created pyramid
1208 ///////////////////////////////////////////////////////////////////////////////
1209
1210 SMDS_MeshVolume* SMDS_Mesh::AddVolumeWithID(const SMDS_MeshFace * f1,
1211                                             const SMDS_MeshFace * f2,
1212                                             const SMDS_MeshFace * f3,
1213                                             const SMDS_MeshFace * f4,
1214                                             const SMDS_MeshFace * f5,
1215                                             int ID)
1216 {
1217   if (!hasConstructionFaces())
1218     return NULL;
1219   if ( !f1 || !f2 || !f3 || !f4 || !f5) return 0;
1220   if ( NbVolumes() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
1221   SMDS_MeshVolume * volume = new SMDS_VolumeOfFaces(f1,f2,f3,f4,f5);
1222   adjustmyCellsCapacity(ID);
1223   myCells[ID] = volume;
1224   myInfo.myNbPyramids++;
1225
1226   if (!registerElement(ID, volume)) {
1227     registerElement(myElementIDFactory->GetFreeID(), volume);
1228     //RemoveElement(volume, false);
1229     //volume = NULL;
1230   }
1231   return volume;
1232 }
1233
1234 ///////////////////////////////////////////////////////////////////////////////
1235 ///Create a new prism defined by its faces and add it to the mesh.
1236 ///@return The created prism
1237 ///////////////////////////////////////////////////////////////////////////////
1238
1239 SMDS_MeshVolume* SMDS_Mesh::AddVolume(const SMDS_MeshFace * f1,
1240                                       const SMDS_MeshFace * f2,
1241                                       const SMDS_MeshFace * f3,
1242                                       const SMDS_MeshFace * f4,
1243                                       const SMDS_MeshFace * f5,
1244                                       const SMDS_MeshFace * f6)
1245 {
1246  if (!hasConstructionFaces())
1247     return NULL;
1248   return AddVolumeWithID(f1,f2,f3,f4,f5,f6, myElementIDFactory->GetFreeID());
1249 }
1250
1251 ///////////////////////////////////////////////////////////////////////////////
1252 ///Create a new prism defined by its faces and add it to the mesh.
1253 ///@param ID The ID of the new volume
1254 ///@return The created prism
1255 ///////////////////////////////////////////////////////////////////////////////
1256
1257 SMDS_MeshVolume* SMDS_Mesh::AddVolumeWithID(const SMDS_MeshFace * f1,
1258                                             const SMDS_MeshFace * f2,
1259                                             const SMDS_MeshFace * f3,
1260                                             const SMDS_MeshFace * f4,
1261                                             const SMDS_MeshFace * f5,
1262                                             const SMDS_MeshFace * f6,
1263                                             int ID)
1264 {
1265   if (!hasConstructionFaces())
1266     return NULL;
1267   if ( !f1 || !f2 || !f3 || !f4 || !f5 || !f6) return 0;
1268   if ( NbVolumes() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
1269   SMDS_MeshVolume * volume = new SMDS_VolumeOfFaces(f1,f2,f3,f4,f5,f6);
1270   adjustmyCellsCapacity(ID);
1271   myCells[ID] = volume;
1272   myInfo.myNbPrisms++;
1273
1274   if (!registerElement(ID, volume)) {
1275     registerElement(myElementIDFactory->GetFreeID(), volume);
1276     //RemoveElement(volume, false);
1277     //volume = NULL;
1278   }
1279   return volume;
1280 }
1281
1282 ///////////////////////////////////////////////////////////////////////////////
1283 /// Add a polygon defined by its nodes IDs
1284 ///////////////////////////////////////////////////////////////////////////////
1285
1286 SMDS_MeshFace* SMDS_Mesh::AddPolygonalFaceWithID (const vector<int> & nodes_ids,
1287                                                   const int           ID)
1288 {
1289   int nbNodes = nodes_ids.size();
1290   vector<const SMDS_MeshNode*> nodes (nbNodes);
1291   for (int i = 0; i < nbNodes; i++) {
1292     nodes[i] = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(nodes_ids[i]);
1293     if (!nodes[i]) return NULL;
1294   }
1295   return SMDS_Mesh::AddPolygonalFaceWithID(nodes, ID);
1296 }
1297
1298 ///////////////////////////////////////////////////////////////////////////////
1299 /// Add a polygon defined by its nodes
1300 ///////////////////////////////////////////////////////////////////////////////
1301
1302 SMDS_MeshFace*
1303 SMDS_Mesh::AddPolygonalFaceWithID (const vector<const SMDS_MeshNode*> & nodes,
1304                                    const int                            ID)
1305 {
1306   SMDS_MeshFace * face;
1307
1308   if ( NbFaces() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
1309   if (hasConstructionEdges())
1310   {
1311     return NULL;
1312   }
1313   else
1314   {
1315     myNodeIds.resize( nodes.size() );
1316     for ( size_t i = 0; i < nodes.size(); ++i )
1317       myNodeIds[i] = nodes[i]->getVtkId();
1318
1319     SMDS_VtkFace *facevtk = myFacePool->getNew();
1320     facevtk->initPoly(myNodeIds, this);
1321     if (!this->registerElement(ID,facevtk))
1322     {
1323       this->myGrid->GetCellTypesArray()->SetValue(facevtk->getVtkId(), VTK_EMPTY_CELL);
1324       myFacePool->destroy(facevtk);
1325       return 0;
1326     }
1327     face = facevtk;
1328
1329     adjustmyCellsCapacity(ID);
1330     myCells[ID] = face;
1331     myInfo.myNbPolygons++;
1332   }
1333
1334   return face;
1335 }
1336
1337 ///////////////////////////////////////////////////////////////////////////////
1338 /// Add a polygon defined by its nodes.
1339 /// An ID is automatically affected to the created face.
1340 ///////////////////////////////////////////////////////////////////////////////
1341
1342 SMDS_MeshFace* SMDS_Mesh::AddPolygonalFace (const vector<const SMDS_MeshNode*> & nodes)
1343 {
1344   return SMDS_Mesh::AddPolygonalFaceWithID(nodes, myElementIDFactory->GetFreeID());
1345 }
1346
1347 ///////////////////////////////////////////////////////////////////////////////
1348 /// Add a quadratic polygon defined by its nodes IDs
1349 ///////////////////////////////////////////////////////////////////////////////
1350
1351 SMDS_MeshFace* SMDS_Mesh::AddQuadPolygonalFaceWithID (const vector<int> & nodes_ids,
1352                                                       const int           ID)
1353 {
1354   vector<const SMDS_MeshNode*> nodes( nodes_ids.size() );
1355   for ( size_t i = 0; i < nodes.size(); i++) {
1356     nodes[i] = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(nodes_ids[i]);
1357     if (!nodes[i]) return NULL;
1358   }
1359   return SMDS_Mesh::AddQuadPolygonalFaceWithID(nodes, ID);
1360 }
1361
1362 ///////////////////////////////////////////////////////////////////////////////
1363 /// Add a quadratic polygon defined by its nodes
1364 ///////////////////////////////////////////////////////////////////////////////
1365
1366 SMDS_MeshFace*
1367 SMDS_Mesh::AddQuadPolygonalFaceWithID (const vector<const SMDS_MeshNode*> & nodes,
1368                                        const int                            ID)
1369 {
1370   SMDS_MeshFace * face;
1371
1372   if ( NbFaces() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
1373   if (hasConstructionEdges())
1374   {
1375     return NULL;
1376   }
1377   else
1378   {
1379     myNodeIds.resize( nodes.size() );
1380     for ( size_t i = 0; i < nodes.size(); ++i )
1381       myNodeIds[i] = nodes[i]->getVtkId();
1382
1383     SMDS_VtkFace *facevtk = myFacePool->getNew();
1384     facevtk->initQuadPoly(myNodeIds, this);
1385     if (!this->registerElement(ID,facevtk))
1386     {
1387       this->myGrid->GetCellTypesArray()->SetValue(facevtk->getVtkId(), VTK_EMPTY_CELL);
1388       myFacePool->destroy(facevtk);
1389       return 0;
1390     }
1391     face = facevtk;
1392     adjustmyCellsCapacity(ID);
1393     myCells[ID] = face;
1394     myInfo.myNbQuadPolygons++;
1395   }
1396   return face;
1397 }
1398
1399 ///////////////////////////////////////////////////////////////////////////////
1400 /// Add a quadratic polygon defined by its nodes.
1401 /// An ID is automatically affected to the created face.
1402 ///////////////////////////////////////////////////////////////////////////////
1403
1404 SMDS_MeshFace* SMDS_Mesh::AddQuadPolygonalFace (const vector<const SMDS_MeshNode*> & nodes)
1405 {
1406   return SMDS_Mesh::AddQuadPolygonalFaceWithID(nodes, myElementIDFactory->GetFreeID());
1407 }
1408
1409 ///////////////////////////////////////////////////////////////////////////////
1410 /// Create a new polyhedral volume and add it to the mesh.
1411 /// @param ID The ID of the new volume
1412 /// @return The created volume or NULL if an element with this ID already exists
1413 /// or if input nodes are not found.
1414 ///////////////////////////////////////////////////////////////////////////////
1415
1416 SMDS_MeshVolume * SMDS_Mesh::AddPolyhedralVolumeWithID
1417                              (const vector<int> & nodes_ids,
1418                               const vector<int> & quantities,
1419                               const int           ID)
1420 {
1421   int nbNodes = nodes_ids.size();
1422   vector<const SMDS_MeshNode*> nodes (nbNodes);
1423   for (int i = 0; i < nbNodes; i++) {
1424     nodes[i] = (SMDS_MeshNode *)myNodeIDFactory->MeshElement(nodes_ids[i]);
1425     if (!nodes[i]) return NULL;
1426   }
1427   return SMDS_Mesh::AddPolyhedralVolumeWithID(nodes, quantities, ID);
1428 }
1429
1430 ///////////////////////////////////////////////////////////////////////////////
1431 /// Create a new polyhedral volume and add it to the mesh.
1432 /// @param ID The ID of the new volume
1433 /// @return The created  volume
1434 ///////////////////////////////////////////////////////////////////////////////
1435
1436 SMDS_MeshVolume*
1437 SMDS_Mesh::AddPolyhedralVolumeWithID (const vector<const SMDS_MeshNode*>& nodes,
1438                                       const vector<int>                 & quantities,
1439                                       const int                           ID)
1440 {
1441   SMDS_MeshVolume* volume = 0;
1442   if ( nodes.empty() || quantities.empty() )
1443     return NULL;
1444   if ( NbVolumes() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
1445   if (hasConstructionFaces())
1446   {
1447     return NULL;
1448   }
1449   else if (hasConstructionEdges())
1450   {
1451     return NULL;
1452   }
1453   else
1454   {
1455     //#ifdef VTK_HAVE_POLYHEDRON
1456     myNodeIds.resize( nodes.size() );
1457     for ( size_t i = 0; i < nodes.size(); ++i )
1458       myNodeIds[i] = nodes[i]->getVtkId();
1459
1460     SMDS_VtkVolume *volvtk = myVolumePool->getNew();
1461     volvtk->initPoly(myNodeIds, quantities, this);
1462     if (!this->registerElement(ID, volvtk))
1463     {
1464       this->myGrid->GetCellTypesArray()->SetValue(volvtk->getVtkId(), VTK_EMPTY_CELL);
1465       myVolumePool->destroy(volvtk);
1466       return 0;
1467     }
1468     volume = volvtk;
1469     //#else
1470     //      for ( int i = 0; i < nodes.size(); ++i )
1471     //      if ( !nodes[ i ] ) return 0;
1472     //      volume = new SMDS_PolyhedralVolumeOfNodes(nodes, quantities);
1473     //#endif
1474     adjustmyCellsCapacity(ID);
1475     myCells[ID] = volume;
1476     myInfo.myNbPolyhedrons++;
1477   }
1478
1479   //#ifndef VTK_HAVE_POLYHEDRON
1480   //  if (!registerElement(ID, volume))
1481   //    {
1482   //      registerElement(myElementIDFactory->GetFreeID(), volume);
1483   //      //RemoveElement(volume, false);
1484   //      //volume = NULL;
1485   //    }
1486   //#endif
1487   return volume;
1488 }
1489
1490 ///////////////////////////////////////////////////////////////////////////////
1491 /// Create a new polyhedral volume and add it to the mesh.
1492 /// @return The created  volume
1493 ///////////////////////////////////////////////////////////////////////////////
1494
1495 SMDS_MeshVolume* SMDS_Mesh::AddPolyhedralVolume
1496                             (const vector<const SMDS_MeshNode*> & nodes,
1497                              const vector<int>                  & quantities)
1498 {
1499   int ID = myElementIDFactory->GetFreeID();
1500   SMDS_MeshVolume * v = SMDS_Mesh::AddPolyhedralVolumeWithID(nodes, quantities, ID);
1501   if (v == NULL) myElementIDFactory->ReleaseID(ID);
1502   return v;
1503 }
1504
1505 SMDS_MeshVolume* SMDS_Mesh::AddVolumeFromVtkIds(const std::vector<vtkIdType>& vtkNodeIds)
1506 {
1507   int ID = myElementIDFactory->GetFreeID();
1508   SMDS_MeshVolume * v = SMDS_Mesh::AddVolumeFromVtkIdsWithID(vtkNodeIds, ID);
1509   if (v == NULL) myElementIDFactory->ReleaseID(ID);
1510   return v;
1511 }
1512
1513 SMDS_MeshVolume* SMDS_Mesh::AddVolumeFromVtkIdsWithID(const std::vector<vtkIdType>& vtkNodeIds, const int ID)
1514 {
1515   SMDS_VtkVolume *volvtk = myVolumePool->getNew();
1516   volvtk->init(vtkNodeIds, this);
1517   if (!this->registerElement(ID,volvtk))
1518     {
1519       this->myGrid->GetCellTypesArray()->SetValue(volvtk->getVtkId(), VTK_EMPTY_CELL);
1520       myVolumePool->destroy(volvtk);
1521       return 0;
1522     }
1523   adjustmyCellsCapacity(ID);
1524   myCells[ID] = volvtk;
1525   vtkIdType aVtkType = volvtk->GetVtkType();
1526   switch (aVtkType)
1527   {
1528     case VTK_TETRA:
1529       myInfo.myNbTetras++;
1530       break;
1531     case VTK_PYRAMID:
1532       myInfo.myNbPyramids++;
1533       break;
1534     case VTK_WEDGE:
1535       myInfo.myNbPrisms++;
1536       break;
1537     case VTK_HEXAHEDRON:
1538       myInfo.myNbHexas++;
1539       break;
1540     case VTK_QUADRATIC_TETRA:
1541       myInfo.myNbQuadTetras++;
1542       break;
1543     case VTK_QUADRATIC_PYRAMID:
1544       myInfo.myNbQuadPyramids++;
1545       break;
1546     case VTK_QUADRATIC_WEDGE:
1547       myInfo.myNbQuadPrisms++;
1548       break;
1549     case VTK_QUADRATIC_HEXAHEDRON:
1550       myInfo.myNbQuadHexas++;
1551       break;
1552 //#ifdef VTK_HAVE_POLYHEDRON
1553     case VTK_POLYHEDRON:
1554       myInfo.myNbPolyhedrons++;
1555       break;
1556 //#endif
1557     default:
1558       myInfo.myNbPolyhedrons++;
1559       break;
1560   }
1561   return volvtk;
1562 }
1563
1564 SMDS_MeshFace* SMDS_Mesh::AddFaceFromVtkIds(const std::vector<vtkIdType>& vtkNodeIds)
1565 {
1566   int ID = myElementIDFactory->GetFreeID();
1567   SMDS_MeshFace * f = SMDS_Mesh::AddFaceFromVtkIdsWithID(vtkNodeIds, ID);
1568   if (f == NULL) myElementIDFactory->ReleaseID(ID);
1569   return f;
1570 }
1571
1572 SMDS_MeshFace* SMDS_Mesh::AddFaceFromVtkIdsWithID(const std::vector<vtkIdType>& vtkNodeIds, const int ID)
1573 {
1574   SMDS_VtkFace *facevtk = myFacePool->getNew();
1575   facevtk->init(vtkNodeIds, this);
1576   if (!this->registerElement(ID,facevtk))
1577     {
1578       this->myGrid->GetCellTypesArray()->SetValue(facevtk->getVtkId(), VTK_EMPTY_CELL);
1579       myFacePool->destroy(facevtk);
1580       return 0;
1581     }
1582   adjustmyCellsCapacity(ID);
1583   myCells[ID] = facevtk;
1584   vtkIdType aVtkType = facevtk->GetVtkType();
1585   switch (aVtkType)
1586   {
1587     case VTK_TRIANGLE:
1588       myInfo.myNbTriangles++;
1589       break;
1590     case VTK_QUAD:
1591       myInfo.myNbQuadrangles++;
1592       break;
1593     case VTK_QUADRATIC_TRIANGLE:
1594       myInfo.myNbQuadTriangles++;
1595       break;
1596     case VTK_QUADRATIC_QUAD:
1597       myInfo.myNbQuadQuadrangles++;
1598       break;
1599     case VTK_BIQUADRATIC_QUAD:
1600       myInfo.myNbBiQuadQuadrangles++;
1601       break;
1602     case VTK_BIQUADRATIC_TRIANGLE:
1603       myInfo.myNbBiQuadTriangles++;
1604       break;
1605     case VTK_POLYGON:
1606       myInfo.myNbPolygons++;
1607       break;
1608      default:
1609       myInfo.myNbPolygons++;
1610   }
1611   return facevtk;
1612 }
1613
1614 ///////////////////////////////////////////////////////////////////////////////
1615 /// Registers element with the given ID, maintains inverse connections
1616 ///////////////////////////////////////////////////////////////////////////////
1617 bool SMDS_Mesh::registerElement(int ID, SMDS_MeshElement* element)
1618 {
1619   if ((ID >=0) && (ID < (int)myCells.size()) && myCells[ID]) // --- already bound
1620   {
1621     MESSAGE(" ------------------ already bound "<< ID << " " << myCells[ID]->getVtkId());
1622     return false;
1623   }
1624
1625   element->myID = ID;
1626   element->myMeshId = myMeshId;
1627
1628   SMDS_MeshCell *cell = dynamic_cast<SMDS_MeshCell*>(element);
1629   MYASSERT(cell);
1630   int vtkId = cell->getVtkId();  
1631   if (vtkId == -1)
1632     vtkId = myElementIDFactory->SetInVtkGrid(element);
1633
1634   if (vtkId >= (int)myCellIdVtkToSmds.size()) // --- resize local vector
1635   {
1636     myCellIdVtkToSmds.resize(vtkId + SMDS_Mesh::chunkSize, -1);
1637   }
1638   myCellIdVtkToSmds[vtkId] = ID;
1639
1640   myElementIDFactory->updateMinMax(ID);
1641   return true;
1642 }
1643
1644 //=======================================================================
1645 //function : MoveNode
1646 //purpose  : 
1647 //=======================================================================
1648
1649 void SMDS_Mesh::MoveNode(const SMDS_MeshNode *n, double x, double y, double z)
1650 {
1651   SMDS_MeshNode * node=const_cast<SMDS_MeshNode*>(n);
1652   node->setXYZ(x,y,z);
1653 }
1654
1655 ///////////////////////////////////////////////////////////////////////////////
1656 /// Return the node whose SMDS ID is 'ID'.
1657 ///////////////////////////////////////////////////////////////////////////////
1658 const SMDS_MeshNode * SMDS_Mesh::FindNode(int ID) const
1659 {
1660   if (ID < 1 || ID >= (int)myNodes.size())
1661   {
1662     return 0;
1663   }
1664   return (const SMDS_MeshNode *)myNodes[ID];
1665 }
1666
1667 ///////////////////////////////////////////////////////////////////////////////
1668 /// Return the node whose VTK ID is 'vtkId'.
1669 ///////////////////////////////////////////////////////////////////////////////
1670 const SMDS_MeshNode * SMDS_Mesh::FindNodeVtk(int vtkId) const
1671 {
1672   // TODO if needed use mesh->nodeIdFromVtkToSmds
1673   if ( vtkId < 0 || vtkId+1 >= (int) myNodes.size() )
1674   {
1675     MESSAGE("------------------------------------------------------------------------- ");
1676     MESSAGE("---------------------------- bad VTK ID " << vtkId << " " << myNodes.size());
1677     MESSAGE("------------------------------------------------------------------------- ");
1678     return 0;
1679   }
1680   return (const SMDS_MeshNode *)myNodes[vtkId+1];
1681 }
1682
1683 ///////////////////////////////////////////////////////////////////////////////
1684 ///Create a triangle and add it to the current mesh. This method do not bind an
1685 ///ID to the create triangle.
1686 ///////////////////////////////////////////////////////////////////////////////
1687 SMDS_MeshFace * SMDS_Mesh::createTriangle(const SMDS_MeshNode * node1,
1688                                           const SMDS_MeshNode * node2,
1689                                           const SMDS_MeshNode * node3,
1690                                           int ID)
1691 {
1692   if ( !node1 || !node2 || !node3) return 0;
1693   if ( NbFaces() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
1694   if(hasConstructionEdges())
1695   {
1696     SMDS_MeshEdge *edge1, *edge2, *edge3;
1697     edge1=FindEdgeOrCreate(node1,node2);
1698     edge2=FindEdgeOrCreate(node2,node3);
1699     edge3=FindEdgeOrCreate(node3,node1);
1700
1701     //int ID = myElementIDFactory->GetFreeID(); // -PR- voir si on range cet element
1702     SMDS_MeshFace * face = new SMDS_FaceOfEdges(edge1,edge2,edge3);
1703     adjustmyCellsCapacity(ID);
1704     myCells[ID] = face;
1705     myInfo.myNbTriangles++;
1706     return face;
1707   }
1708   else
1709   {
1710     // --- retrieve nodes ID
1711     myNodeIds.resize(3);
1712     myNodeIds[0] = node1->getVtkId();
1713     myNodeIds[1] = node2->getVtkId();
1714     myNodeIds[2] = node3->getVtkId();
1715
1716     SMDS_MeshFace * face = 0;
1717     SMDS_VtkFace *facevtk = myFacePool->getNew();
1718     facevtk->init(myNodeIds, this); // put in vtkUnstructuredGrid
1719     if (!this->registerElement(ID,facevtk))
1720     {
1721       this->myGrid->GetCellTypesArray()->SetValue(facevtk->getVtkId(), VTK_EMPTY_CELL);
1722       myFacePool->destroy(facevtk);
1723       return 0;
1724     }
1725     face = facevtk;
1726     adjustmyCellsCapacity(ID);
1727     myCells[ID] = face;
1728     myInfo.myNbTriangles++;
1729     return face;
1730   }
1731 }
1732
1733 ///////////////////////////////////////////////////////////////////////////////
1734 ///Create a quadrangle and add it to the current mesh. This methode do not bind
1735 ///a ID to the create triangle.
1736 ///////////////////////////////////////////////////////////////////////////////
1737 SMDS_MeshFace * SMDS_Mesh::createQuadrangle(const SMDS_MeshNode * node1,
1738                                             const SMDS_MeshNode * node2,
1739                                             const SMDS_MeshNode * node3,
1740                                             const SMDS_MeshNode * node4,
1741                                             int ID)
1742 {
1743   if ( !node1 || !node2 || !node3 || !node4 ) return 0;
1744   if ( NbFaces() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
1745   if(hasConstructionEdges())
1746   {
1747     SMDS_MeshEdge *edge1, *edge2, *edge3, *edge4;
1748     edge1=FindEdgeOrCreate(node1,node2);
1749     edge2=FindEdgeOrCreate(node2,node3);
1750     edge3=FindEdgeOrCreate(node3,node4);
1751     edge4=FindEdgeOrCreate(node4,node1);
1752
1753     SMDS_MeshFace * face = new SMDS_FaceOfEdges(edge1,edge2,edge3,edge4);
1754     adjustmyCellsCapacity(ID);
1755     myCells[ID] = face;
1756     myInfo.myNbQuadrangles++;
1757     return face;
1758   }
1759   else
1760   {
1761     // --- retrieve nodes ID
1762     myNodeIds.resize(4);
1763     myNodeIds[0] = node1->getVtkId();
1764     myNodeIds[1] = node2->getVtkId();
1765     myNodeIds[2] = node3->getVtkId();
1766     myNodeIds[3] = node4->getVtkId();
1767
1768     SMDS_MeshFace * face = 0;
1769     SMDS_VtkFace *facevtk = myFacePool->getNew();
1770     facevtk->init(myNodeIds, this);
1771     if (!this->registerElement(ID,facevtk))
1772     {
1773       this->myGrid->GetCellTypesArray()->SetValue(facevtk->getVtkId(), VTK_EMPTY_CELL);
1774       myFacePool->destroy(facevtk);
1775       return 0;
1776     }
1777     face = facevtk;
1778     adjustmyCellsCapacity(ID);
1779     myCells[ID] = face;
1780     myInfo.myNbQuadrangles++;
1781     return face;
1782   }
1783 }
1784
1785 ///////////////////////////////////////////////////////////////////////////////
1786 /// Remove a node and all the elements which own this node
1787 ///////////////////////////////////////////////////////////////////////////////
1788
1789 void SMDS_Mesh::RemoveNode(const SMDS_MeshNode * node)
1790 {
1791   RemoveElement(node, true);
1792 }
1793
1794 ///////////////////////////////////////////////////////////////////////////////
1795 /// Remove an edge and all the elements which own this edge
1796 ///////////////////////////////////////////////////////////////////////////////
1797
1798 void SMDS_Mesh::Remove0DElement(const SMDS_Mesh0DElement * elem0d)
1799 {
1800   RemoveElement(elem0d,true);
1801 }
1802
1803 ///////////////////////////////////////////////////////////////////////////////
1804 /// Remove an edge and all the elements which own this edge
1805 ///////////////////////////////////////////////////////////////////////////////
1806
1807 void SMDS_Mesh::RemoveEdge(const SMDS_MeshEdge * edge)
1808 {
1809   RemoveElement(edge,true);
1810 }
1811
1812 ///////////////////////////////////////////////////////////////////////////////
1813 /// Remove an face and all the elements which own this face
1814 ///////////////////////////////////////////////////////////////////////////////
1815
1816 void SMDS_Mesh::RemoveFace(const SMDS_MeshFace * face)
1817 {
1818   RemoveElement(face, true);
1819 }
1820
1821 ///////////////////////////////////////////////////////////////////////////////
1822 /// Remove a volume
1823 ///////////////////////////////////////////////////////////////////////////////
1824
1825 void SMDS_Mesh::RemoveVolume(const SMDS_MeshVolume * volume)
1826 {
1827   RemoveElement(volume, true);
1828 }
1829
1830 //=======================================================================
1831 //function : RemoveFromParent
1832 //purpose  :
1833 //=======================================================================
1834
1835 bool SMDS_Mesh::RemoveFromParent()
1836 {
1837   if (myParent==NULL) return false;
1838   else return (myParent->RemoveSubMesh(this));
1839 }
1840
1841 //=======================================================================
1842 //function : RemoveSubMesh
1843 //purpose  :
1844 //=======================================================================
1845
1846 bool SMDS_Mesh::RemoveSubMesh(const SMDS_Mesh * aMesh)
1847 {
1848   bool found = false;
1849
1850   list<SMDS_Mesh *>::iterator itmsh=myChildren.begin();
1851   for (; itmsh!=myChildren.end() && !found; itmsh++)
1852   {
1853     SMDS_Mesh * submesh = *itmsh;
1854     if (submesh == aMesh)
1855     {
1856       found = true;
1857       myChildren.erase(itmsh);
1858     }
1859   }
1860
1861   return found;
1862 }
1863
1864 //=======================================================================
1865 //function : ChangeElementNodes
1866 //purpose  :
1867 //=======================================================================
1868
1869 bool SMDS_Mesh::ChangeElementNodes(const SMDS_MeshElement * element,
1870                                    const SMDS_MeshNode    * nodes[],
1871                                    const int                nbnodes)
1872 {
1873   // keep current nodes of elem
1874   set<const SMDS_MeshNode*> oldNodes( element->begin_nodes(), element->end_nodes() );
1875
1876   // change nodes
1877   bool Ok = false;
1878   SMDS_MeshCell* cell = dynamic_cast<SMDS_MeshCell*>((SMDS_MeshElement*) element);
1879   if (cell)
1880   {
1881     Ok = cell->vtkOrder(nodes, nbnodes);
1882     Ok = cell->ChangeNodes(nodes, nbnodes);
1883   }
1884
1885   if ( Ok ) { // update InverseElements
1886
1887     set<const SMDS_MeshNode*>::iterator it;
1888
1889     // AddInverseElement to new nodes
1890     for ( int i = 0; i < nbnodes; i++ ) {
1891       it = oldNodes.find( nodes[i] );
1892       if ( it == oldNodes.end() )
1893         // new node
1894         const_cast<SMDS_MeshNode*>( nodes[i] )->AddInverseElement( cell );
1895       else
1896         // remove from oldNodes a node that remains in elem
1897         oldNodes.erase( it );
1898     }
1899     // RemoveInverseElement from the nodes removed from elem
1900     for ( it = oldNodes.begin(); it != oldNodes.end(); it++ )
1901     {
1902       SMDS_MeshNode * n = const_cast<SMDS_MeshNode *>( *it );
1903       n->RemoveInverseElement( cell );
1904     }
1905   }
1906
1907   return Ok;
1908 }
1909
1910 //=======================================================================
1911 //function : ChangePolyhedronNodes
1912 //purpose  : to change nodes of polyhedral volume
1913 //=======================================================================
1914 bool SMDS_Mesh::ChangePolyhedronNodes (const SMDS_MeshElement *            elem,
1915                                        const vector<const SMDS_MeshNode*>& nodes,
1916                                        const vector<int>                 & quantities)
1917 {
1918   if (elem->GetType() != SMDSAbs_Volume) {
1919     MESSAGE("WRONG ELEM TYPE");
1920     return false;
1921   }
1922
1923   const SMDS_VtkVolume* vol = dynamic_cast<const SMDS_VtkVolume*>(elem);
1924   if (!vol) {
1925     return false;
1926   }
1927
1928   // keep current nodes of elem
1929   set<const SMDS_MeshElement*> oldNodes;
1930   SMDS_ElemIteratorPtr itn = elem->nodesIterator();
1931   while (itn->more()) {
1932     oldNodes.insert(itn->next());
1933   }
1934
1935   // change nodes
1936   // TODO remove this function
1937   //bool Ok = const_cast<SMDS_VtkVolume*>(vol)->ChangeNodes(nodes, quantities);
1938   bool Ok = false;
1939   if (!Ok) {
1940     return false;
1941   }
1942
1943   // update InverseElements
1944
1945   // AddInverseElement to new nodes
1946   int nbnodes = nodes.size();
1947   set<const SMDS_MeshElement*>::iterator it;
1948   for (int i = 0; i < nbnodes; i++) {
1949     it = oldNodes.find(nodes[i]);
1950     if (it == oldNodes.end()) {
1951       // new node
1952       const_cast<SMDS_MeshNode*>(nodes[i])->AddInverseElement(elem);
1953     } else {
1954       // remove from oldNodes a node that remains in elem
1955       oldNodes.erase(it);
1956     }
1957   }
1958
1959   // RemoveInverseElement from the nodes removed from elem
1960   for (it = oldNodes.begin(); it != oldNodes.end(); it++) {
1961     SMDS_MeshNode * n = static_cast<SMDS_MeshNode *>
1962       (const_cast<SMDS_MeshElement *>( *it ));
1963     n->RemoveInverseElement(elem);
1964   }
1965
1966   return Ok;
1967 }
1968
1969
1970 //=======================================================================
1971 //function : Find0DElement
1972 //purpose  :
1973 //=======================================================================
1974 const SMDS_Mesh0DElement* SMDS_Mesh::Find0DElement(int idnode) const
1975 {
1976   const SMDS_MeshNode * node = FindNode(idnode);
1977   if(node == NULL) return NULL;
1978   return Find0DElement(node);
1979 }
1980
1981 const SMDS_Mesh0DElement* SMDS_Mesh::Find0DElement(const SMDS_MeshNode * node)
1982 {
1983   if (!node) return 0;
1984   const SMDS_Mesh0DElement* toReturn = NULL;
1985   SMDS_ElemIteratorPtr it1 = node->GetInverseElementIterator(SMDSAbs_0DElement);
1986   while (it1->more() && (toReturn == NULL)) {
1987     const SMDS_MeshElement* e = it1->next();
1988     if (e->NbNodes() == 1) {
1989       toReturn = static_cast<const SMDS_Mesh0DElement*>(e);
1990     }
1991   }
1992   return toReturn;
1993 }
1994
1995 //=======================================================================
1996 //function : FindBall
1997 //purpose  :
1998 //=======================================================================
1999
2000 const SMDS_BallElement* SMDS_Mesh::FindBall(int idnode) const
2001 {
2002   const SMDS_MeshNode * node = FindNode(idnode);
2003   if(node == NULL) return NULL;
2004   return FindBall(node);
2005 }
2006
2007 const SMDS_BallElement* SMDS_Mesh::FindBall(const SMDS_MeshNode * node)
2008 {
2009   if (!node) return 0;
2010   const SMDS_BallElement* toReturn = NULL;
2011   SMDS_ElemIteratorPtr it1 = node->GetInverseElementIterator(SMDSAbs_Ball);
2012   while (it1->more() && (toReturn == NULL)) {
2013     const SMDS_MeshElement* e = it1->next();
2014     if (e->GetGeomType() == SMDSGeom_BALL)
2015       toReturn = static_cast<const SMDS_BallElement*>(e);
2016   }
2017   return toReturn;
2018 }
2019
2020 //=======================================================================
2021 //function : Find0DElementOrCreate
2022 //purpose  :
2023 //=======================================================================
2024 //SMDS_Mesh0DElement* SMDS_Mesh::Find0DElementOrCreate(const SMDS_MeshNode * node)
2025 //{
2026 //  if (!node) return 0;
2027 //  SMDS_Mesh0DElement * toReturn = NULL;
2028 //  toReturn = const_cast<SMDS_Mesh0DElement*>(Find0DElement(node));
2029 //  if (toReturn == NULL) {
2030 //    //if (my0DElements.Extent() % CHECKMEMORY_INTERVAL == 0) CheckMemory();
2031 //    toReturn = new SMDS_Mesh0DElement(node);
2032 //    my0DElements.Add(toReturn);
2033 //    myInfo.myNb0DElements++;
2034 //  }
2035 //  return toReturn;
2036 //}
2037
2038
2039 //=======================================================================
2040 //function : FindEdge
2041 //purpose  :
2042 //=======================================================================
2043
2044 const SMDS_MeshEdge* SMDS_Mesh::FindEdge(int idnode1, int idnode2) const
2045 {
2046   const SMDS_MeshNode * node1=FindNode(idnode1);
2047   const SMDS_MeshNode * node2=FindNode(idnode2);
2048   if((node1==NULL)||(node2==NULL)) return NULL;
2049   return FindEdge(node1,node2);
2050 }
2051
2052 //#include "Profiler.h"
2053 const SMDS_MeshEdge* SMDS_Mesh::FindEdge(const SMDS_MeshNode * node1,
2054                                          const SMDS_MeshNode * node2)
2055 {
2056   if ( !node1 ) return 0;
2057   const SMDS_MeshEdge * toReturn=NULL;
2058   //PROFILER_Init();
2059   //PROFILER_Set();
2060   SMDS_ElemIteratorPtr it1=node1->GetInverseElementIterator(SMDSAbs_Edge);
2061   //PROFILER_Get(0);
2062   //PROFILER_Set();
2063   while(it1->more()) {
2064     const SMDS_MeshElement * e = it1->next();
2065     if ( e->NbNodes() == 2 && e->GetNodeIndex( node2 ) >= 0 ) {
2066       toReturn = static_cast<const SMDS_MeshEdge*>( e );
2067       break;
2068     }
2069   }
2070   //PROFILER_Get(1);
2071   return toReturn;
2072 }
2073
2074
2075 //=======================================================================
2076 //function : FindEdgeOrCreate
2077 //purpose  :
2078 //=======================================================================
2079
2080 SMDS_MeshEdge* SMDS_Mesh::FindEdgeOrCreate(const SMDS_MeshNode * node1,
2081                                            const SMDS_MeshNode * node2)
2082 {
2083   if ( !node1 || !node2) return 0;
2084   SMDS_MeshEdge * toReturn=NULL;
2085   toReturn=const_cast<SMDS_MeshEdge*>(FindEdge(node1,node2));
2086   if(toReturn==NULL) {
2087     if ( NbEdges() % CHECKMEMORY_INTERVAL == 0 ) CheckMemory();
2088     int ID = myElementIDFactory->GetFreeID(); // -PR- voir si on range cet element
2089     adjustmyCellsCapacity(ID);
2090     myNodeIds.resize(2);
2091     myNodeIds[0] = node1->getVtkId();
2092     myNodeIds[1] = node2->getVtkId();
2093
2094     SMDS_VtkEdge *edgevtk = myEdgePool->getNew();
2095     edgevtk->init(myNodeIds, this);
2096     if (!this->registerElement(ID,edgevtk))
2097     {
2098       this->myGrid->GetCellTypesArray()->SetValue(edgevtk->getVtkId(), VTK_EMPTY_CELL);
2099       myEdgePool->destroy(edgevtk);
2100       return 0;
2101     }
2102     toReturn = edgevtk;
2103     myCells[ID] = toReturn;
2104     myInfo.myNbEdges++;
2105   }
2106   return toReturn;
2107 }
2108
2109
2110 //=======================================================================
2111 //function : FindEdge
2112 //purpose  :
2113 //=======================================================================
2114
2115 const SMDS_MeshEdge* SMDS_Mesh::FindEdge(int idnode1, int idnode2,
2116                                          int idnode3) const
2117 {
2118   const SMDS_MeshNode * node1=FindNode(idnode1);
2119   const SMDS_MeshNode * node2=FindNode(idnode2);
2120   const SMDS_MeshNode * node3=FindNode(idnode3);
2121   return FindEdge(node1,node2,node3);
2122 }
2123
2124 const SMDS_MeshEdge* SMDS_Mesh::FindEdge(const SMDS_MeshNode * node1,
2125                                          const SMDS_MeshNode * node2,
2126                                          const SMDS_MeshNode * node3)
2127 {
2128   if ( !node1 ) return 0;
2129   SMDS_ElemIteratorPtr it1 = node1->GetInverseElementIterator(SMDSAbs_Edge);
2130   while(it1->more()) {
2131     const SMDS_MeshElement * e = it1->next();
2132     if ( e->NbNodes() == 3 ) {
2133       SMDS_ElemIteratorPtr it2 = e->nodesIterator();
2134       while(it2->more()) {
2135         const SMDS_MeshElement* n = it2->next();
2136         if( n!=node1 &&
2137             n!=node2 &&
2138             n!=node3 )
2139         {
2140           e = 0;
2141           break;
2142         }
2143       }
2144       if ( e )
2145         return static_cast<const SMDS_MeshEdge *> (e);
2146     }
2147   }
2148   return 0;
2149 }
2150
2151
2152 //=======================================================================
2153 //function : FindFace
2154 //purpose  :
2155 //=======================================================================
2156
2157 const SMDS_MeshFace* SMDS_Mesh::FindFace(int idnode1, int idnode2,
2158         int idnode3) const
2159 {
2160   const SMDS_MeshNode * node1=FindNode(idnode1);
2161   const SMDS_MeshNode * node2=FindNode(idnode2);
2162   const SMDS_MeshNode * node3=FindNode(idnode3);
2163   return FindFace(node1, node2, node3);
2164 }
2165
2166 const SMDS_MeshFace* SMDS_Mesh::FindFace(const SMDS_MeshNode *node1,
2167                                          const SMDS_MeshNode *node2,
2168                                          const SMDS_MeshNode *node3)
2169 {
2170   if ( !node1 ) return 0;
2171   SMDS_ElemIteratorPtr it1 = node1->GetInverseElementIterator(SMDSAbs_Face);
2172   while(it1->more()) {
2173     const SMDS_MeshElement * e = it1->next();
2174     if ( e->NbNodes() == 3 ) {
2175       SMDS_ElemIteratorPtr it2 = e->nodesIterator();
2176       while(it2->more()) {
2177         const SMDS_MeshElement* n = it2->next();
2178         if( n!=node1 &&
2179             n!=node2 &&
2180             n!=node3 )
2181         {
2182           e = 0;
2183           break;
2184         }
2185       }
2186       if ( e )
2187         return static_cast<const SMDS_MeshFace *> (e);
2188     }
2189   }
2190   return 0;
2191 }
2192
2193 SMDS_MeshFace* SMDS_Mesh::FindFaceOrCreate(const SMDS_MeshNode *node1,
2194                                            const SMDS_MeshNode *node2,
2195                                            const SMDS_MeshNode *node3)
2196 {
2197   SMDS_MeshFace * toReturn=NULL;
2198   toReturn = const_cast<SMDS_MeshFace*>(FindFace(node1,node2,node3));
2199   if(toReturn==NULL) {
2200     int ID = myElementIDFactory->GetFreeID();
2201     toReturn = createTriangle(node1,node2,node3, ID);
2202   }
2203   return toReturn;
2204 }
2205
2206
2207 //=======================================================================
2208 //function : FindFace
2209 //purpose  :
2210 //=======================================================================
2211
2212 const SMDS_MeshFace* SMDS_Mesh::FindFace(int idnode1, int idnode2,
2213                                          int idnode3, int idnode4) const
2214 {
2215   const SMDS_MeshNode * node1=FindNode(idnode1);
2216   const SMDS_MeshNode * node2=FindNode(idnode2);
2217   const SMDS_MeshNode * node3=FindNode(idnode3);
2218   const SMDS_MeshNode * node4=FindNode(idnode4);
2219   return FindFace(node1, node2, node3, node4);
2220 }
2221
2222 const SMDS_MeshFace* SMDS_Mesh::FindFace(const SMDS_MeshNode *node1,
2223                                          const SMDS_MeshNode *node2,
2224                                          const SMDS_MeshNode *node3,
2225                                          const SMDS_MeshNode *node4)
2226 {
2227   if ( !node1 ) return 0;
2228   SMDS_ElemIteratorPtr it1 = node1->GetInverseElementIterator(SMDSAbs_Face);
2229   while(it1->more()) {
2230     const SMDS_MeshElement * e = it1->next();
2231     if ( e->NbNodes() == 4 ) {
2232       SMDS_ElemIteratorPtr it2 = e->nodesIterator();
2233       while(it2->more()) {
2234         const SMDS_MeshElement* n = it2->next();
2235         if( n!=node1 &&
2236             n!=node2 &&
2237             n!=node3 &&
2238             n!=node4 )
2239         {
2240           e = 0;
2241           break;
2242         }
2243       }
2244       if ( e )
2245         return static_cast<const SMDS_MeshFace *> (e);
2246     }
2247   }
2248   return 0;
2249 }
2250
2251 SMDS_MeshFace* SMDS_Mesh::FindFaceOrCreate(const SMDS_MeshNode *node1,
2252                                            const SMDS_MeshNode *node2,
2253                                            const SMDS_MeshNode *node3,
2254                                            const SMDS_MeshNode *node4)
2255 {
2256   SMDS_MeshFace * toReturn=NULL;
2257   toReturn=const_cast<SMDS_MeshFace*>(FindFace(node1,node2,node3,node4));
2258   if(toReturn==NULL) {
2259     int ID = myElementIDFactory->GetFreeID();
2260     toReturn=createQuadrangle(node1,node2,node3,node4,ID);
2261   }
2262   return toReturn;
2263 }
2264
2265
2266 //=======================================================================
2267 //function : FindFace
2268 //purpose  :quadratic triangle
2269 //=======================================================================
2270
2271 const SMDS_MeshFace* SMDS_Mesh::FindFace(int idnode1, int idnode2,
2272                                          int idnode3, int idnode4,
2273                                          int idnode5, int idnode6) const
2274 {
2275   const SMDS_MeshNode * node1 = FindNode(idnode1);
2276   const SMDS_MeshNode * node2 = FindNode(idnode2);
2277   const SMDS_MeshNode * node3 = FindNode(idnode3);
2278   const SMDS_MeshNode * node4 = FindNode(idnode4);
2279   const SMDS_MeshNode * node5 = FindNode(idnode5);
2280   const SMDS_MeshNode * node6 = FindNode(idnode6);
2281   return FindFace(node1, node2, node3, node4, node5, node6);
2282 }
2283
2284 const SMDS_MeshFace* SMDS_Mesh::FindFace(const SMDS_MeshNode *node1,
2285                                          const SMDS_MeshNode *node2,
2286                                          const SMDS_MeshNode *node3,
2287                                          const SMDS_MeshNode *node4,
2288                                          const SMDS_MeshNode *node5,
2289                                          const SMDS_MeshNode *node6)
2290 {
2291   if ( !node1 ) return 0;
2292   SMDS_ElemIteratorPtr it1 = node1->GetInverseElementIterator(SMDSAbs_Face);
2293   while(it1->more()) {
2294     const SMDS_MeshElement * e = it1->next();
2295     if ( e->NbNodes() == 6 ) {
2296       SMDS_ElemIteratorPtr it2 = e->nodesIterator();
2297       while(it2->more()) {
2298         const SMDS_MeshElement* n = it2->next();
2299         if( n!=node1 &&
2300             n!=node2 &&
2301             n!=node3 &&
2302             n!=node4 &&
2303             n!=node5 &&
2304             n!=node6 )
2305         {
2306           e = 0;
2307           break;
2308         }
2309       }
2310       if ( e )
2311         return static_cast<const SMDS_MeshFace *> (e);
2312     }
2313   }
2314   return 0;
2315 }
2316
2317
2318 //=======================================================================
2319 //function : FindFace
2320 //purpose  : quadratic quadrangle
2321 //=======================================================================
2322
2323 const SMDS_MeshFace* SMDS_Mesh::FindFace(int idnode1, int idnode2,
2324                                          int idnode3, int idnode4,
2325                                          int idnode5, int idnode6,
2326                                          int idnode7, int idnode8) const
2327 {
2328   const SMDS_MeshNode * node1 = FindNode(idnode1);
2329   const SMDS_MeshNode * node2 = FindNode(idnode2);
2330   const SMDS_MeshNode * node3 = FindNode(idnode3);
2331   const SMDS_MeshNode * node4 = FindNode(idnode4);
2332   const SMDS_MeshNode * node5 = FindNode(idnode5);
2333   const SMDS_MeshNode * node6 = FindNode(idnode6);
2334   const SMDS_MeshNode * node7 = FindNode(idnode7);
2335   const SMDS_MeshNode * node8 = FindNode(idnode8);
2336   return FindFace(node1, node2, node3, node4, node5, node6, node7, node8);
2337 }
2338
2339 const SMDS_MeshFace* SMDS_Mesh::FindFace(const SMDS_MeshNode *node1,
2340                                          const SMDS_MeshNode *node2,
2341                                          const SMDS_MeshNode *node3,
2342                                          const SMDS_MeshNode *node4,
2343                                          const SMDS_MeshNode *node5,
2344                                          const SMDS_MeshNode *node6,
2345                                          const SMDS_MeshNode *node7,
2346                                          const SMDS_MeshNode *node8)
2347 {
2348   if ( !node1 ) return 0;
2349   SMDS_ElemIteratorPtr it1 = node1->GetInverseElementIterator(SMDSAbs_Face);
2350   while(it1->more()) {
2351     const SMDS_MeshElement * e = it1->next();
2352     if ( e->NbNodes() == 8 ) {
2353       SMDS_ElemIteratorPtr it2 = e->nodesIterator();
2354       while(it2->more()) {
2355         const SMDS_MeshElement* n = it2->next();
2356         if( n!=node1 &&
2357             n!=node2 &&
2358             n!=node3 &&
2359             n!=node4 &&
2360             n!=node5 &&
2361             n!=node6 &&
2362             n!=node7 &&
2363             n!=node8 )
2364         {
2365           e = 0;
2366           break;
2367         }
2368       }
2369       if ( e )
2370         return static_cast<const SMDS_MeshFace *> (e);
2371     }
2372   }
2373   return 0;
2374 }
2375
2376
2377 //=======================================================================
2378 //function : FindElement
2379 //purpose  :
2380 //=======================================================================
2381
2382 const SMDS_MeshElement* SMDS_Mesh::FindElement(int IDelem) const
2383 {
2384   if ( IDelem <= 0 || IDelem >= (int)myCells.size() )
2385   {
2386     return 0;
2387   }
2388   return myCells[IDelem];
2389 }
2390
2391 //=======================================================================
2392 //function : FindFace
2393 //purpose  : find polygon
2394 //=======================================================================
2395
2396 const SMDS_MeshFace* SMDS_Mesh::FindFace (const vector<int>& nodes_ids) const
2397 {
2398   int nbnodes = nodes_ids.size();
2399   vector<const SMDS_MeshNode *> poly_nodes (nbnodes);
2400   for (int inode = 0; inode < nbnodes; inode++) {
2401     const SMDS_MeshNode * node = FindNode(nodes_ids[inode]);
2402     if (node == NULL) return NULL;
2403     poly_nodes[inode] = node;
2404   }
2405   return FindFace(poly_nodes);
2406 }
2407
2408 const SMDS_MeshFace* SMDS_Mesh::FindFace (const vector<const SMDS_MeshNode *>& nodes)
2409 {
2410   return (const SMDS_MeshFace*) FindElement( nodes, SMDSAbs_Face );
2411 }
2412
2413
2414 //================================================================================
2415 /*!
2416  * \brief Return element based on all given nodes
2417  *  \param nodes - node of element
2418  *  \param type - type of element
2419  *  \param noMedium - true if medium nodes of quadratic element are not included in <nodes>
2420  *  \retval const SMDS_MeshElement* - found element or NULL
2421  */
2422 //================================================================================
2423
2424 const SMDS_MeshElement* SMDS_Mesh::FindElement (const vector<const SMDS_MeshNode *>& nodes,
2425                                                 const SMDSAbs_ElementType            type,
2426                                                 const bool                           noMedium)
2427 {
2428   if ( nodes.size() > 0 && nodes[0] )
2429   {
2430     SMDS_ElemIteratorPtr itF = nodes[0]->GetInverseElementIterator(type);
2431     while (itF->more())
2432     {
2433       const SMDS_MeshElement* e = itF->next();
2434       int nbNodesToCheck = noMedium ? e->NbCornerNodes() : e->NbNodes();
2435       if ( nbNodesToCheck == (int)nodes.size() )
2436       {
2437         for ( size_t i = 1; e && i < nodes.size(); ++i )
2438         {
2439           int nodeIndex = e->GetNodeIndex( nodes[ i ]);
2440           if ( nodeIndex < 0 || nodeIndex >= nbNodesToCheck )
2441             e = 0;
2442         }
2443         if ( e )
2444           return e;
2445       }
2446     }
2447   }
2448   return NULL;
2449 }
2450
2451 //=======================================================================
2452 //function : DumpNodes
2453 //purpose  :
2454 //=======================================================================
2455
2456 void SMDS_Mesh::DumpNodes() const
2457 {
2458   SMDS_NodeIteratorPtr itnode=nodesIterator();
2459   while(itnode->more()) ; //MESSAGE(itnode->next());
2460 }
2461
2462 //=======================================================================
2463 //function : Dump0DElements
2464 //purpose  :
2465 //=======================================================================
2466 void SMDS_Mesh::Dump0DElements() const
2467 {
2468   SMDS_ElemIteratorPtr it0d = elementsIterator(SMDSAbs_0DElement);
2469   while(it0d->more()) ; //MESSAGE(it0d->next());
2470 }
2471
2472 //=======================================================================
2473 //function : DumpEdges
2474 //purpose  :
2475 //=======================================================================
2476
2477 void SMDS_Mesh::DumpEdges() const
2478 {
2479   SMDS_EdgeIteratorPtr itedge=edgesIterator();
2480   while(itedge->more()) ; //MESSAGE(itedge->next());
2481 }
2482
2483 //=======================================================================
2484 //function : DumpFaces
2485 //purpose  :
2486 //=======================================================================
2487
2488 void SMDS_Mesh::DumpFaces() const
2489 {
2490   SMDS_FaceIteratorPtr itface=facesIterator();
2491   while(itface->more()) ; //MESSAGE(itface->next());
2492 }
2493
2494 //=======================================================================
2495 //function : DumpVolumes
2496 //purpose  :
2497 //=======================================================================
2498
2499 void SMDS_Mesh::DumpVolumes() const
2500 {
2501   SMDS_VolumeIteratorPtr itvol=volumesIterator();
2502   while(itvol->more()) ; //MESSAGE(itvol->next());
2503 }
2504
2505 //=======================================================================
2506 //function : DebugStats
2507 //purpose  :
2508 //=======================================================================
2509
2510 void SMDS_Mesh::DebugStats() const
2511 {
2512   MESSAGE("Debug stats of mesh : ");
2513
2514   MESSAGE("===== NODES ====="<<NbNodes());
2515   MESSAGE("===== 0DELEMS ====="<<Nb0DElements());
2516   MESSAGE("===== EDGES ====="<<NbEdges());
2517   MESSAGE("===== FACES ====="<<NbFaces());
2518   MESSAGE("===== VOLUMES ====="<<NbVolumes());
2519
2520   MESSAGE("End Debug stats of mesh ");
2521
2522   //#ifdef DEB
2523
2524   SMDS_NodeIteratorPtr itnode=nodesIterator();
2525   int sizeofnodes = 0;
2526   int sizeoffaces = 0;
2527
2528   while(itnode->more())
2529   {
2530     const SMDS_MeshNode *node = itnode->next();
2531
2532     sizeofnodes += sizeof(*node);
2533
2534     SMDS_ElemIteratorPtr it = node->GetInverseElementIterator();
2535     while(it->more())
2536     {
2537       const SMDS_MeshElement *me = it->next();
2538       sizeofnodes += sizeof(me);
2539     }
2540   }
2541
2542   SMDS_FaceIteratorPtr itface=facesIterator();
2543   while(itface->more())
2544   {
2545     const SMDS_MeshElement *face = itface->next();
2546     sizeoffaces += sizeof(*face);
2547   }
2548
2549   MESSAGE("total size of node elements = " << sizeofnodes);;
2550   MESSAGE("total size of face elements = " << sizeoffaces);;
2551
2552   //#endif
2553 }
2554
2555 ///////////////////////////////////////////////////////////////////////////////
2556 /// Return the number of nodes
2557 ///////////////////////////////////////////////////////////////////////////////
2558 int SMDS_Mesh::NbNodes() const
2559 {
2560   return myInfo.NbNodes();
2561 }
2562
2563 ///////////////////////////////////////////////////////////////////////////////
2564 /// Return the number of 0D elements
2565 ///////////////////////////////////////////////////////////////////////////////
2566 int SMDS_Mesh::Nb0DElements() const
2567 {
2568   return myInfo.Nb0DElements();
2569 }
2570
2571 ///////////////////////////////////////////////////////////////////////////////
2572 /// Return the number of 0D elements
2573 ///////////////////////////////////////////////////////////////////////////////
2574 int SMDS_Mesh::NbBalls() const
2575 {
2576   return myInfo.NbBalls();
2577 }
2578
2579 ///////////////////////////////////////////////////////////////////////////////
2580 /// Return the number of edges (including construction edges)
2581 ///////////////////////////////////////////////////////////////////////////////
2582 int SMDS_Mesh::NbEdges() const
2583 {
2584   return myInfo.NbEdges();
2585 }
2586
2587 ///////////////////////////////////////////////////////////////////////////////
2588 /// Return the number of faces (including construction faces)
2589 ///////////////////////////////////////////////////////////////////////////////
2590 int SMDS_Mesh::NbFaces() const
2591 {
2592   return myInfo.NbFaces();
2593 }
2594
2595 ///////////////////////////////////////////////////////////////////////////////
2596 /// Return the number of volumes
2597 ///////////////////////////////////////////////////////////////////////////////
2598 int SMDS_Mesh::NbVolumes() const
2599 {
2600   return myInfo.NbVolumes();
2601 }
2602
2603 ///////////////////////////////////////////////////////////////////////////////
2604 /// Return the number of child mesh of this mesh.
2605 /// Note that the tree structure of SMDS_Mesh is unused in SMESH
2606 ///////////////////////////////////////////////////////////////////////////////
2607 int SMDS_Mesh::NbSubMesh() const
2608 {
2609   return myChildren.size();
2610 }
2611
2612 ///////////////////////////////////////////////////////////////////////////////
2613 /// Destroy the mesh and all its elements
2614 /// All pointer on elements owned by this mesh become illegals.
2615 ///////////////////////////////////////////////////////////////////////////////
2616 SMDS_Mesh::~SMDS_Mesh()
2617 {
2618   list<SMDS_Mesh*>::iterator itc=myChildren.begin();
2619   while(itc!=myChildren.end())
2620   {
2621     delete *itc;
2622     itc++;
2623   }
2624
2625   if(myParent==NULL)
2626   {
2627     delete myNodeIDFactory;
2628     delete myElementIDFactory;
2629   }
2630   else
2631   {
2632     SMDS_ElemIteratorPtr eIt = elementsIterator();
2633     while ( eIt->more() )
2634       {
2635         const SMDS_MeshElement *elem = eIt->next();
2636         myElementIDFactory->ReleaseID(elem->GetID(), elem->getVtkId());
2637       }
2638     SMDS_NodeIteratorPtr itn = nodesIterator();
2639     while (itn->more())
2640       {
2641         const SMDS_MeshNode *node = itn->next();
2642         ((SMDS_MeshNode*)node)->SetPosition(SMDS_SpacePosition::originSpacePosition());
2643         myNodeIDFactory->ReleaseID(node->GetID(), node->getVtkId());
2644       }
2645   }
2646   myGrid->Delete();
2647
2648   delete myNodePool;
2649   delete myVolumePool;
2650   delete myFacePool;
2651   delete myEdgePool;
2652   delete myBallPool;
2653 }
2654
2655 //================================================================================
2656 /*!
2657  * \brief Clear all data
2658  */
2659 //================================================================================
2660
2661 void SMDS_Mesh::Clear()
2662 {
2663   if (myParent!=NULL)
2664     {
2665     SMDS_ElemIteratorPtr eIt = elementsIterator();
2666     while ( eIt->more() )
2667       {
2668         const SMDS_MeshElement *elem = eIt->next();
2669         myElementIDFactory->ReleaseID(elem->GetID(), elem->getVtkId());
2670       }
2671     SMDS_NodeIteratorPtr itn = nodesIterator();
2672     while (itn->more())
2673       {
2674         const SMDS_MeshNode *node = itn->next();
2675         myNodeIDFactory->ReleaseID(node->GetID(), node->getVtkId());
2676       }
2677     }
2678   else
2679     {
2680     myNodeIDFactory->Clear();
2681     myElementIDFactory->Clear();
2682     }
2683
2684   // SMDS_ElemIteratorPtr itv = elementsIterator();
2685   // while (itv->more())
2686   //   {
2687   //     SMDS_MeshElement* elem = (SMDS_MeshElement*)(itv->next());
2688   //     SMDSAbs_ElementType aType = elem->GetType();
2689   //     switch (aType)
2690   //     {
2691   //       case SMDSAbs_0DElement:
2692   //         delete elem;
2693   //         break;
2694   //       case SMDSAbs_Edge:
2695   //          myEdgePool->destroy(static_cast<SMDS_VtkEdge*>(elem));
2696   //         break;
2697   //       case SMDSAbs_Face:
2698   //         myFacePool->destroy(static_cast<SMDS_VtkFace*>(elem));
2699   //         break;
2700   //       case SMDSAbs_Volume:
2701   //         myVolumePool->destroy(static_cast<SMDS_VtkVolume*>(elem));
2702   //         break;
2703   //       case SMDSAbs_Ball:
2704   //         myBallPool->destroy(static_cast<SMDS_BallElement*>(elem));
2705   //         break;
2706   //       default:
2707   //         break;
2708   //     }
2709   //   }
2710   myVolumePool->clear();
2711   myFacePool->clear();
2712   myEdgePool->clear();
2713   myBallPool->clear();
2714
2715   clearVector( myCells );
2716   clearVector( myCellIdVtkToSmds );
2717
2718   SMDS_NodeIteratorPtr itn = nodesIterator();
2719   while (itn->more())
2720     {
2721       SMDS_MeshNode *node = (SMDS_MeshNode*)(itn->next());
2722       node->SetPosition(SMDS_SpacePosition::originSpacePosition());
2723       //myNodePool->destroy(node);
2724     }
2725   myNodePool->clear();
2726   clearVector( myNodes );
2727
2728   list<SMDS_Mesh*>::iterator itc=myChildren.begin();
2729   while(itc!=myChildren.end())
2730     (*itc)->Clear();
2731
2732   myModified = false;
2733   myModifTime++;
2734   xmin = 0; xmax = 0;
2735   ymin = 0; ymax = 0;
2736   zmin = 0; zmax = 0;
2737
2738   myInfo.Clear();
2739
2740   myGrid->Initialize();
2741   myGrid->Allocate();
2742   vtkPoints* points = vtkPoints::New();
2743   // rnv: to fix bug "21125: EDF 1233 SMESH: Degrardation of precision in a test case for quadratic conversion"
2744   // using double type for storing coordinates of nodes instead float.
2745   points->SetDataType(VTK_DOUBLE);
2746   points->SetNumberOfPoints(0 /*SMDS_Mesh::chunkSize*/);
2747   myGrid->SetPoints( points );
2748   points->Delete();
2749   myGrid->BuildLinks();
2750 }
2751
2752 ///////////////////////////////////////////////////////////////////////////////
2753 /// Return true if this mesh create faces with edges.
2754 /// A false returned value mean that faces are created with nodes. A concequence
2755 /// is, iteration on edges (SMDS_Element::edgesIterator) will be unavailable.
2756 ///////////////////////////////////////////////////////////////////////////////
2757 bool SMDS_Mesh::hasConstructionEdges()
2758 {
2759         return myHasConstructionEdges;
2760 }
2761
2762 ///////////////////////////////////////////////////////////////////////////////
2763 /// Return true if this mesh create volumes with faces
2764 /// A false returned value mean that volumes are created with nodes or edges.
2765 /// (see hasConstructionEdges)
2766 /// A concequence is, iteration on faces (SMDS_Element::facesIterator) will be
2767 /// unavailable.
2768 ///////////////////////////////////////////////////////////////////////////////
2769 bool SMDS_Mesh::hasConstructionFaces()
2770 {
2771         return myHasConstructionFaces;
2772 }
2773
2774 ///////////////////////////////////////////////////////////////////////////////
2775 /// Return true if nodes are linked to the finit elements, they are belonging to.
2776 /// Currently, It always return true.
2777 ///////////////////////////////////////////////////////////////////////////////
2778 bool SMDS_Mesh::hasInverseElements()
2779 {
2780         return myHasInverseElements;
2781 }
2782
2783 ///////////////////////////////////////////////////////////////////////////////
2784 /// Make this mesh creating construction edges (see hasConstructionEdges)
2785 /// @param b true to have construction edges, else false.
2786 ///////////////////////////////////////////////////////////////////////////////
2787 void SMDS_Mesh::setConstructionEdges(bool b)
2788 {
2789         myHasConstructionEdges=b;
2790 }
2791
2792 ///////////////////////////////////////////////////////////////////////////////
2793 /// Make this mesh creating construction faces (see hasConstructionFaces)
2794 /// @param b true to have construction faces, else false.
2795 ///////////////////////////////////////////////////////////////////////////////
2796 void SMDS_Mesh::setConstructionFaces(bool b)
2797 {
2798          myHasConstructionFaces=b;
2799 }
2800
2801 ///////////////////////////////////////////////////////////////////////////////
2802 /// Make this mesh creating link from nodes to elements (see hasInverseElements)
2803 /// @param b true to link nodes to elements, else false.
2804 ///////////////////////////////////////////////////////////////////////////////
2805 void SMDS_Mesh::setInverseElements(bool b)
2806 {
2807   if(!b) MESSAGE("Error : inverseElement=false not implemented");
2808   myHasInverseElements=b;
2809 }
2810
2811 namespace {
2812
2813   //================================================================================
2814   /*!
2815    * \brief Iterator on elements in id increasing order
2816    */
2817   //================================================================================
2818
2819   template <typename ELEM=const SMDS_MeshElement*>
2820   class IdSortedIterator : public SMDS_Iterator<ELEM>
2821   {
2822     SMDS_MeshElementIDFactory&       myIDFact;
2823     int                              myID, myMaxID, myNbFound, myTotalNb;
2824     SMDSAbs_ElementType              myType;
2825     ELEM                             myElem;
2826
2827   public:
2828     IdSortedIterator(const SMDS_MeshElementIDFactory& fact,
2829                      const SMDSAbs_ElementType        type, // SMDSAbs_All NOT allowed!!! 
2830                      const int                        totalNb)
2831       :myIDFact( fact ),
2832        myID(1), myMaxID( myIDFact.GetMaxID() ),myNbFound(0), myTotalNb( totalNb ),
2833        myType( type ),
2834        myElem(0)
2835     {
2836       next();
2837     }
2838     bool more()
2839     {
2840       return myElem;
2841     }
2842     ELEM next()
2843     {
2844       ELEM current = myElem;
2845
2846       for ( myElem = 0;  !myElem && myNbFound < myTotalNb && myID <= myMaxID; ++myID )
2847         if ((myElem = (ELEM) myIDFact.MeshElement( myID ))
2848             && myElem->GetType() != myType )
2849           myElem = 0;
2850
2851       myNbFound += bool(myElem);
2852
2853       return current;
2854     }
2855   };
2856
2857   //================================================================================
2858   /*!
2859    * \brief Iterator on vector of elements, possibly being resized while iteration
2860    */
2861   //================================================================================
2862
2863   template<typename RETURN_VALUE,
2864            typename VECTOR_VALUE=SMDS_MeshCell*,
2865            typename VALUE_FILTER=SMDS::NonNullFilter<VECTOR_VALUE> >
2866   class ElemVecIterator: public SMDS_Iterator<RETURN_VALUE>
2867   {
2868     const std::vector<VECTOR_VALUE>& _vector;
2869     size_t                           _index;
2870     bool                             _more;
2871     VALUE_FILTER                     _filter;
2872   public:
2873     ElemVecIterator(const std::vector<VECTOR_VALUE>& vec,
2874                     const VALUE_FILTER&              filter=VALUE_FILTER() )
2875       :_vector( vec ), _index(0), _more( !vec.empty() ), _filter( filter )
2876     {
2877       if ( _more && !_filter( _vector[ _index ]))
2878         next();
2879     }
2880     virtual bool more()
2881     {
2882       return _more;
2883     }
2884     virtual RETURN_VALUE next()
2885     {
2886       if ( !_more ) return NULL;
2887       VECTOR_VALUE current = _vector[ _index ];
2888       _more = 0;
2889       while ( !_more && ++_index < _vector.size() )
2890         _more = _filter( _vector[ _index ]);
2891       return (RETURN_VALUE) current;
2892     }
2893   };
2894 }
2895
2896 ///////////////////////////////////////////////////////////////////////////////
2897 /// Return an iterator on nodes of the current mesh factory
2898 ///////////////////////////////////////////////////////////////////////////////
2899
2900 SMDS_NodeIteratorPtr SMDS_Mesh::nodesIterator(bool idInceasingOrder) const
2901 {
2902   // naturally always sorted by ID
2903   typedef ElemVecIterator<const SMDS_MeshNode*, SMDS_MeshNode*> TIterator;
2904   return SMDS_NodeIteratorPtr( new TIterator(myNodes));
2905 }
2906
2907 SMDS_ElemIteratorPtr SMDS_Mesh::elementGeomIterator(SMDSAbs_GeometryType type) const
2908 {
2909   // naturally always sorted by ID
2910   typedef ElemVecIterator
2911     < const SMDS_MeshElement*, SMDS_MeshCell*, SMDS_MeshElement::GeomFilter > TIterator;
2912   return SMDS_ElemIteratorPtr
2913     (new TIterator(myCells, SMDS_MeshElement::GeomFilter( type )));
2914 }
2915
2916 SMDS_ElemIteratorPtr SMDS_Mesh::elementEntityIterator(SMDSAbs_EntityType type) const
2917 {
2918   if ( type == SMDSEntity_Node )
2919   {
2920     typedef ElemVecIterator<const SMDS_MeshElement*, SMDS_MeshNode*> TIterator;
2921     return SMDS_ElemIteratorPtr( new TIterator(myNodes));
2922   }
2923   // naturally always sorted by ID
2924   typedef ElemVecIterator
2925     < const SMDS_MeshElement*, SMDS_MeshCell*, SMDS_MeshElement::EntityFilter > TIterator;
2926   return SMDS_ElemIteratorPtr
2927     (new TIterator(myCells, SMDS_MeshElement::EntityFilter( type )));
2928 }
2929
2930 ///////////////////////////////////////////////////////////////////////////////
2931 /// Return an iterator on elements of the current mesh factory
2932 ///////////////////////////////////////////////////////////////////////////////
2933 SMDS_ElemIteratorPtr SMDS_Mesh::elementsIterator(SMDSAbs_ElementType type) const
2934 {
2935   // naturally always sorted by ID
2936   switch ( type ) {
2937
2938   case SMDSAbs_All:
2939     return SMDS_ElemIteratorPtr (new ElemVecIterator<const SMDS_MeshElement*>(myCells));
2940
2941   case SMDSAbs_Node:
2942     return SMDS_ElemIteratorPtr
2943       ( new ElemVecIterator<const SMDS_MeshElement*, SMDS_MeshNode*>( myNodes ));
2944
2945   default:
2946     typedef ElemVecIterator
2947       < const SMDS_MeshElement*, SMDS_MeshCell*, SMDS_MeshElement::TypeFilter > TIterator;
2948     return SMDS_ElemIteratorPtr (new TIterator(myCells, SMDS_MeshElement::TypeFilter( type )));
2949   }
2950   return SMDS_ElemIteratorPtr();
2951 }
2952
2953 ///////////////////////////////////////////////////////////////////////////////
2954 ///Return an iterator on edges of the current mesh.
2955 ///////////////////////////////////////////////////////////////////////////////
2956
2957 SMDS_EdgeIteratorPtr SMDS_Mesh::edgesIterator(bool idInceasingOrder) const
2958 {
2959   // naturally always sorted by ID
2960   typedef ElemVecIterator
2961     < const SMDS_MeshEdge*, SMDS_MeshCell*, SMDS_MeshElement::TypeFilter > TIterator;
2962   return SMDS_EdgeIteratorPtr
2963     (new TIterator(myCells, SMDS_MeshElement::TypeFilter( SMDSAbs_Edge )));
2964 }
2965
2966 ///////////////////////////////////////////////////////////////////////////////
2967 ///Return an iterator on faces of the current mesh.
2968 ///////////////////////////////////////////////////////////////////////////////
2969
2970 SMDS_FaceIteratorPtr SMDS_Mesh::facesIterator(bool idInceasingOrder) const
2971 {
2972   // naturally always sorted by ID
2973   typedef ElemVecIterator
2974     < const SMDS_MeshFace*, SMDS_MeshCell*, SMDS_MeshElement::TypeFilter > TIterator;
2975   return SMDS_FaceIteratorPtr
2976     (new TIterator(myCells, SMDS_MeshElement::TypeFilter( SMDSAbs_Face )));
2977 }
2978
2979 ///////////////////////////////////////////////////////////////////////////////
2980 ///Return an iterator on volumes of the current mesh.
2981 ///////////////////////////////////////////////////////////////////////////////
2982
2983 SMDS_VolumeIteratorPtr SMDS_Mesh::volumesIterator(bool idInceasingOrder) const
2984 {
2985   // naturally always sorted by ID
2986   typedef ElemVecIterator
2987     < const SMDS_MeshVolume*, SMDS_MeshCell*, SMDS_MeshElement::TypeFilter > TIterator;
2988   return SMDS_VolumeIteratorPtr
2989     (new TIterator(myCells, SMDS_MeshElement::TypeFilter( SMDSAbs_Volume )));
2990 }
2991
2992 ///////////////////////////////////////////////////////////////////////////////
2993 /// Do intersection of sets (more than 2)
2994 ///////////////////////////////////////////////////////////////////////////////
2995 static set<const SMDS_MeshElement*> * intersectionOfSets(
2996         set<const SMDS_MeshElement*> vs[], int numberOfSets)
2997 {
2998         set<const SMDS_MeshElement*>* rsetA=new set<const SMDS_MeshElement*>(vs[0]);
2999         set<const SMDS_MeshElement*>* rsetB;
3000
3001         for(int i=0; i<numberOfSets-1; i++)
3002         {
3003                 rsetB=new set<const SMDS_MeshElement*>();
3004                 set_intersection(
3005                         rsetA->begin(), rsetA->end(),
3006                         vs[i+1].begin(), vs[i+1].end(),
3007                         inserter(*rsetB, rsetB->begin()));
3008                 delete rsetA;
3009                 rsetA=rsetB;
3010         }
3011         return rsetA;
3012 }
3013
3014 ///////////////////////////////////////////////////////////////////////////////
3015 /// Return the list of finite elements owning the given element: elements
3016 /// containing all the nodes of the given element, for instance faces and
3017 /// volumes containing a given edge.
3018 ///////////////////////////////////////////////////////////////////////////////
3019 static set<const SMDS_MeshElement*> * getFinitElements(const SMDS_MeshElement * element)
3020 {
3021   int numberOfSets=element->NbNodes();
3022   set<const SMDS_MeshElement*> *initSet = new set<const SMDS_MeshElement*>[numberOfSets];
3023
3024   SMDS_ElemIteratorPtr itNodes=element->nodesIterator();
3025
3026   int i = 0;
3027   while ( itNodes->more() )
3028   {
3029     const SMDS_MeshElement* node = itNodes->next();
3030     MYASSERT(node);
3031     const SMDS_MeshNode * n=static_cast<const SMDS_MeshNode*>(node);
3032     SMDS_ElemIteratorPtr itFe = n->GetInverseElementIterator();
3033
3034     while ( itFe->more() )
3035     {
3036       const SMDS_MeshElement* elem = itFe->next();
3037       MYASSERT(elem);
3038       initSet[i].insert(elem);
3039     }
3040
3041     i++;
3042   }
3043   set<const SMDS_MeshElement*> *retSet = intersectionOfSets( initSet, numberOfSets );
3044   delete [] initSet;
3045   return retSet;
3046 }
3047
3048 ///////////////////////////////////////////////////////////////////////////////
3049 /// Return the list of nodes used only by the given elements
3050 ///////////////////////////////////////////////////////////////////////////////
3051 static set<const SMDS_MeshElement*> * getExclusiveNodes(set<const SMDS_MeshElement*>& elements)
3052 {
3053   set<const SMDS_MeshElement*> * toReturn=new set<const SMDS_MeshElement*>();
3054   set<const SMDS_MeshElement*>::iterator itElements=elements.begin();
3055
3056   while(itElements!=elements.end())
3057   {
3058     SMDS_ElemIteratorPtr itNodes = (*itElements)->nodesIterator();
3059     itElements++;
3060
3061     while(itNodes->more())
3062     {
3063       const SMDS_MeshNode * n=static_cast<const SMDS_MeshNode*>(itNodes->next());
3064       SMDS_ElemIteratorPtr itFe = n->GetInverseElementIterator();
3065       set<const SMDS_MeshElement*> s;
3066       while(itFe->more())
3067         s.insert(itFe->next());
3068       if(s==elements) toReturn->insert(n);
3069     }
3070   }
3071   return toReturn;
3072 }
3073
3074 ///////////////////////////////////////////////////////////////////////////////
3075 ///Find the children of an element that are made of given nodes
3076 ///@param setOfChildren The set in which matching children will be inserted
3077 ///@param element The element were to search matching children
3078 ///@param nodes The nodes that the children must have to be selected
3079 ///////////////////////////////////////////////////////////////////////////////
3080 void SMDS_Mesh::addChildrenWithNodes(set<const SMDS_MeshElement*>& setOfChildren,
3081                                      const SMDS_MeshElement *      element,
3082                                      set<const SMDS_MeshElement*>& nodes)
3083 {
3084   switch(element->GetType())
3085   {
3086   case SMDSAbs_Node:
3087     MESSAGE("Internal Error: This should not happen");
3088     break;
3089   case SMDSAbs_0DElement:
3090   case SMDSAbs_Ball:
3091   {
3092   }
3093   break;
3094   case SMDSAbs_Edge:
3095   {
3096     SMDS_ElemIteratorPtr itn=element->nodesIterator();
3097     while(itn->more())
3098     {
3099       const SMDS_MeshElement * e=itn->next();
3100       if(nodes.find(e)!=nodes.end())
3101       {
3102         setOfChildren.insert(element);
3103         break;
3104       }
3105     }
3106   } break;
3107   case SMDSAbs_Face:
3108   {
3109     SMDS_ElemIteratorPtr itn=element->nodesIterator();
3110     while(itn->more())
3111     {
3112       const SMDS_MeshElement * e=itn->next();
3113       if(nodes.find(e)!=nodes.end())
3114       {
3115         setOfChildren.insert(element);
3116         break;
3117       }
3118     }
3119     if(hasConstructionEdges())
3120     {
3121       SMDS_ElemIteratorPtr ite=element->edgesIterator();
3122       while(ite->more())
3123         addChildrenWithNodes(setOfChildren, ite->next(), nodes);
3124     }
3125   } break;
3126   case SMDSAbs_Volume:
3127   {
3128     if(hasConstructionFaces())
3129     {
3130       SMDS_ElemIteratorPtr ite=element->facesIterator();
3131       while(ite->more())
3132         addChildrenWithNodes(setOfChildren, ite->next(), nodes);
3133     }
3134     else if(hasConstructionEdges())
3135     {
3136       SMDS_ElemIteratorPtr ite=element->edgesIterator();
3137       while(ite->more())
3138         addChildrenWithNodes(setOfChildren, ite->next(), nodes);
3139     }
3140   }
3141   case SMDSAbs_NbElementTypes:
3142   case SMDSAbs_All: break;
3143   }
3144 }
3145
3146 ///////////////////////////////////////////////////////////////////////////////
3147 ///@param elem The element to delete
3148 ///@param removenodes if true remaining nodes will be removed
3149 ///////////////////////////////////////////////////////////////////////////////
3150 void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem,
3151                               const bool removenodes)
3152 {
3153   list<const SMDS_MeshElement *> removedElems;
3154   list<const SMDS_MeshElement *> removedNodes;
3155   RemoveElement( elem, removedElems, removedNodes, removenodes );
3156 }
3157
3158 ///////////////////////////////////////////////////////////////////////////////
3159 ///@param elem The element to delete
3160 ///@param removedElems to be filled with all removed elements
3161 ///@param removedNodes to be filled with all removed nodes
3162 ///@param removenodes if true remaining nodes will be removed
3163 ///////////////////////////////////////////////////////////////////////////////
3164 void SMDS_Mesh::RemoveElement(const SMDS_MeshElement *        elem,
3165                               list<const SMDS_MeshElement *>& removedElems,
3166                               list<const SMDS_MeshElement *>& removedNodes,
3167                               bool                            removenodes)
3168 {
3169   // get finite elements built on elem
3170   set<const SMDS_MeshElement*> * s1;
3171   if (    (elem->GetType() == SMDSAbs_0DElement)
3172       || ((elem->GetType() == SMDSAbs_Edge) && !hasConstructionEdges())
3173       || ((elem->GetType() == SMDSAbs_Face) && !hasConstructionFaces())
3174       ||  (elem->GetType() == SMDSAbs_Volume) )
3175     {
3176       s1 = new set<const SMDS_MeshElement*> ();
3177       s1->insert(elem);
3178     }
3179   else
3180     s1 = getFinitElements(elem);
3181
3182   // get exclusive nodes (which would become free afterwards)
3183   set<const SMDS_MeshElement*> * s2;
3184   if (elem->GetType() == SMDSAbs_Node) // a node is removed
3185     {
3186       // do not remove nodes except elem
3187       s2 = new set<const SMDS_MeshElement*> ();
3188       s2->insert(elem);
3189       removenodes = true;
3190     }
3191   else
3192     s2 = getExclusiveNodes(*s1);
3193
3194   // form the set of finite and construction elements to remove
3195   set<const SMDS_MeshElement*> s3;
3196   set<const SMDS_MeshElement*>::iterator it = s1->begin();
3197   while (it != s1->end())
3198     {
3199       addChildrenWithNodes(s3, *it, *s2);
3200       s3.insert(*it);
3201       it++;
3202     }
3203   if (elem->GetType() != SMDSAbs_Node)
3204     s3.insert(elem);
3205
3206   // remove finite and construction elements
3207   it = s3.begin();
3208   while (it != s3.end())
3209     {