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