]> SALOME platform Git repositories - modules/smesh.git/blob - src/SMDS/SMDS_UnstructuredGrid.cxx
Salome HOME
untabify
[modules/smesh.git] / src / SMDS / SMDS_UnstructuredGrid.cxx
1 // Copyright (C) 2010-2011  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19
20 #define CHRONODEF
21 #include "SMDS_UnstructuredGrid.hxx"
22 #include "SMDS_Mesh.hxx"
23 #include "SMDS_MeshInfo.hxx"
24 #include "SMDS_Downward.hxx"
25 #include "SMDS_MeshVolume.hxx"
26
27 #include "utilities.h"
28
29 #include <vtkCellArray.h>
30 #include <vtkCellLinks.h>
31 #include <vtkIdTypeArray.h>
32 #include <vtkUnsignedCharArray.h>
33
34 #include <list>
35 #include <climits>
36
37 using namespace std;
38
39 SMDS_CellLinks* SMDS_CellLinks::New()
40 {
41   MESSAGE("SMDS_CellLinks::New");
42   return new SMDS_CellLinks();
43 }
44
45 vtkCellLinks::Link* SMDS_CellLinks::ResizeL(vtkIdType sz)
46 {
47   return vtkCellLinks::Resize(sz);
48 }
49
50 vtkIdType SMDS_CellLinks::GetLinksSize()
51 {
52   return this->Size;
53 }
54
55 SMDS_CellLinks::SMDS_CellLinks() :
56   vtkCellLinks()
57 {
58 }
59
60 SMDS_CellLinks::~SMDS_CellLinks()
61 {
62 }
63
64 SMDS_UnstructuredGrid* SMDS_UnstructuredGrid::New()
65 {
66   MESSAGE("SMDS_UnstructuredGrid::New");
67   return new SMDS_UnstructuredGrid();
68 }
69
70 SMDS_UnstructuredGrid::SMDS_UnstructuredGrid() :
71   vtkUnstructuredGrid()
72 {
73   _cellIdToDownId.clear();
74   _downTypes.clear();
75   _downArray.clear();
76   _mesh = 0;
77 }
78
79 SMDS_UnstructuredGrid::~SMDS_UnstructuredGrid()
80 {
81 }
82
83 unsigned long SMDS_UnstructuredGrid::GetMTime()
84 {
85   unsigned long mtime = vtkUnstructuredGrid::GetMTime();
86   MESSAGE("vtkUnstructuredGrid::GetMTime: " << mtime);
87   return mtime;
88 }
89
90 void SMDS_UnstructuredGrid::Update()
91 {
92   MESSAGE("SMDS_UnstructuredGrid::Update");
93   return vtkUnstructuredGrid::Update();
94 }
95
96 void SMDS_UnstructuredGrid::UpdateInformation()
97 {
98   MESSAGE("SMDS_UnstructuredGrid::UpdateInformation");
99   return vtkUnstructuredGrid::UpdateInformation();
100 }
101
102 vtkPoints* SMDS_UnstructuredGrid::GetPoints()
103 {
104   // TODO erreur incomprehensible de la macro vtk GetPoints apparue avec la version paraview de fin aout 2010
105   //MESSAGE("*********************** SMDS_UnstructuredGrid::GetPoints " << this->Points << " " << vtkUnstructuredGrid::GetPoints());
106   return this->Points;
107 }
108
109 //#ifdef VTK_HAVE_POLYHEDRON
110 int SMDS_UnstructuredGrid::InsertNextLinkedCell(int type, int npts, vtkIdType *pts)
111 {
112   if (type != VTK_POLYHEDRON)
113     return vtkUnstructuredGrid::InsertNextLinkedCell(type, npts, pts);
114
115   // --- type = VTK_POLYHEDRON
116   //MESSAGE("InsertNextLinkedCell VTK_POLYHEDRON");
117   int cellid = this->InsertNextCell(type, npts, pts);
118
119   set<vtkIdType> setOfNodes;
120   setOfNodes.clear();
121   int nbfaces = npts;
122   int i = 0;
123   for (int nf = 0; nf < nbfaces; nf++)
124     {
125       int nbnodes = pts[i];
126       i++;
127       for (int k = 0; k < nbnodes; k++)
128         {
129           //MESSAGE(" cell " << cellid << " face " << nf << " node " << pts[i]);
130           setOfNodes.insert(pts[i]);
131           i++;
132         }
133     }
134
135   set<vtkIdType>::iterator it = setOfNodes.begin();
136   for (; it != setOfNodes.end(); ++it)
137     {
138       //MESSAGE("reverse link for node " << *it << " cell " << cellid);
139       this->Links->ResizeCellList(*it, 1);
140       this->Links->AddCellReference(cellid, *it);
141     }
142
143   return cellid;
144 }
145 //#endif
146
147 void SMDS_UnstructuredGrid::setSMDS_mesh(SMDS_Mesh *mesh)
148 {
149   _mesh = mesh;
150 }
151
152 void SMDS_UnstructuredGrid::compactGrid(std::vector<int>& idNodesOldToNew, int newNodeSize,
153                                         std::vector<int>& idCellsOldToNew, int newCellSize)
154 {
155   MESSAGE("------------------------- SMDS_UnstructuredGrid::compactGrid " << newNodeSize << " " << newCellSize);CHRONO(1);
156   int alreadyCopied = 0;
157
158   // --- if newNodeSize, create a new compacted vtkPoints
159
160   vtkPoints *newPoints = vtkPoints::New();
161   newPoints->SetDataType(VTK_DOUBLE);
162   newPoints->SetNumberOfPoints(newNodeSize);
163   if (newNodeSize)
164     {
165       MESSAGE("-------------- compactGrid, newNodeSize " << newNodeSize);
166       // rnv: to fix bug "21125: EDF 1233 SMESH: Degradation of precision in a test case for quadratic conversion"
167       // using double type for storing coordinates of nodes instead float.
168       int oldNodeSize = idNodesOldToNew.size();
169
170       int i = 0;
171       while ( i < oldNodeSize )
172       {
173         // skip a hole if any
174         while ( i < oldNodeSize && idNodesOldToNew[i] < 0 )
175           ++i;
176         int startBloc = i;
177         // look for a block end
178         while ( i < oldNodeSize && idNodesOldToNew[i] >= 0 )
179           ++i;
180         int endBloc = i;
181         copyNodes(newPoints, idNodesOldToNew, alreadyCopied, startBloc, endBloc);
182       }
183       newPoints->Squeeze();
184     }
185
186   // --- create new compacted Connectivity, Locations and Types
187
188   int oldCellSize = this->Types->GetNumberOfTuples();
189
190   vtkCellArray *newConnectivity = vtkCellArray::New();
191   newConnectivity->Initialize();
192   int oldCellDataSize = this->Connectivity->GetData()->GetSize();
193   newConnectivity->Allocate(oldCellDataSize);
194   MESSAGE("oldCellSize="<< oldCellSize << " oldCellDataSize=" << oldCellDataSize);
195
196   vtkUnsignedCharArray *newTypes = vtkUnsignedCharArray::New();
197   newTypes->Initialize();
198   newTypes->SetNumberOfValues(newCellSize);
199
200   vtkIdTypeArray *newLocations = vtkIdTypeArray::New();
201   newLocations->Initialize();
202   newLocations->SetNumberOfValues(newCellSize);
203
204   // TODO some polyhedron may be huge (only in some tests)
205   vtkIdType tmpid[NBMAXNODESINCELL];
206   vtkIdType *pointsCell = &tmpid[0]; // --- points id to fill a new cell
207
208   alreadyCopied = 0;
209   int i = 0;
210   while ( i < oldCellSize )
211   {
212     // skip a hole if any
213     while ( i < oldCellSize && this->Types->GetValue(i) == VTK_EMPTY_CELL )
214       ++i;
215     int startBloc = i;
216     // look for a block end
217     while ( i < oldCellSize && this->Types->GetValue(i) != VTK_EMPTY_CELL )
218       ++i;
219     int endBloc = i;
220     if ( endBloc > startBloc )
221       copyBloc(newTypes, idCellsOldToNew, idNodesOldToNew, newConnectivity, newLocations, pointsCell,
222                alreadyCopied, startBloc, endBloc);
223   }
224
225   newConnectivity->Squeeze();
226
227   if (1/*newNodeSize*/)
228     {
229       MESSAGE("------- newNodeSize, setPoints");
230       this->SetPoints(newPoints);
231       MESSAGE("NumberOfPoints: " << this->GetNumberOfPoints());
232     }
233
234   if (this->FaceLocations)
235     {
236       vtkIdTypeArray *newFaceLocations = vtkIdTypeArray::New();
237       newFaceLocations->Initialize();
238       newFaceLocations->Allocate(newTypes->GetSize());
239       vtkIdTypeArray *newFaces = vtkIdTypeArray::New();
240       newFaces->Initialize();
241       newFaces->Allocate(this->Faces->GetSize());
242       for (int i = 0; i < oldCellSize; i++)
243         {
244           if (this->Types->GetValue(i) == VTK_EMPTY_CELL)
245             continue;
246           int newCellId = idCellsOldToNew[i];
247           if (newTypes->GetValue(newCellId) == VTK_POLYHEDRON)
248             {
249                newFaceLocations->InsertNextValue(newFaces->GetMaxId()+1);
250                int oldFaceLoc = this->FaceLocations->GetValue(i);
251                int nCellFaces = this->Faces->GetValue(oldFaceLoc++);
252                newFaces->InsertNextValue(nCellFaces);
253                for (int n=0; n<nCellFaces; n++)
254                  {
255                    int nptsInFace = this->Faces->GetValue(oldFaceLoc++);
256                    newFaces->InsertNextValue(nptsInFace);
257                    for (int k=0; k<nptsInFace; k++)
258                      {
259                        int oldpt = this->Faces->GetValue(oldFaceLoc++);
260                        newFaces->InsertNextValue(idNodesOldToNew[oldpt]);
261                      }
262                  }
263             }
264           else
265             {
266                newFaceLocations->InsertNextValue(-1);
267             }
268         }
269       newFaceLocations->Squeeze();
270       newFaces->Squeeze();
271       newFaceLocations->Register(this);
272       newFaces->Register(this);
273       this->SetCells(newTypes, newLocations, newConnectivity, newFaceLocations, newFaces);
274       newFaceLocations->Delete();
275       newFaces->Delete();
276     }
277   else
278     this->SetCells(newTypes, newLocations, newConnectivity, FaceLocations, Faces);
279
280   newPoints->Delete();
281   newTypes->Delete();
282   newLocations->Delete();
283   newConnectivity->Delete();
284   this->BuildLinks();
285 }
286
287 void SMDS_UnstructuredGrid::copyNodes(vtkPoints *newPoints, std::vector<int>& idNodesOldToNew, int& alreadyCopied,
288                                       int start, int end)
289 {
290   MESSAGE("copyNodes " << alreadyCopied << " " << start << " " << end << " size: " << end - start << " total: " << alreadyCopied + end - start);
291   void *target = newPoints->GetVoidPointer(3 * alreadyCopied);
292   void *source = this->Points->GetVoidPointer(3 * start);
293   int nbPoints = end - start;
294   if (nbPoints > 0)
295     {
296       memcpy(target, source, 3 * sizeof(double) * nbPoints);
297       for (int j = start; j < end; j++)
298         idNodesOldToNew[j] = alreadyCopied++; // old vtkId --> new vtkId
299     }
300 }
301
302 void SMDS_UnstructuredGrid::copyBloc(vtkUnsignedCharArray *newTypes, std::vector<int>& idCellsOldToNew,
303                                      std::vector<int>& idNodesOldToNew, vtkCellArray* newConnectivity,
304                                      vtkIdTypeArray* newLocations, vtkIdType* pointsCell, int& alreadyCopied,
305                                      int start, int end)
306 {
307   MESSAGE("copyBloc " << alreadyCopied << " " << start << " " << end << " size: " << end - start << " total: " << alreadyCopied + end - start);
308   for (int j = start; j < end; j++)
309     {
310       newTypes->SetValue(alreadyCopied, this->Types->GetValue(j));
311       idCellsOldToNew[j] = alreadyCopied; // old vtkId --> new vtkId
312       vtkIdType oldLoc = this->Locations->GetValue(j);
313       vtkIdType nbpts;
314       vtkIdType *oldPtsCell = 0;
315       this->Connectivity->GetCell(oldLoc, nbpts, oldPtsCell);
316       assert(nbpts < NBMAXNODESINCELL);
317       //MESSAGE(j << " " << alreadyCopied << " " << (int)this->Types->GetValue(j) << " " << oldLoc << " " << nbpts );
318       for (int l = 0; l < nbpts; l++)
319         {
320           int oldval = oldPtsCell[l];
321           pointsCell[l] = idNodesOldToNew[oldval];
322           //MESSAGE("   " << oldval << " " << pointsCell[l]);
323         }
324       /*int newcnt = */newConnectivity->InsertNextCell(nbpts, pointsCell);
325       int newLoc = newConnectivity->GetInsertLocation(nbpts);
326       //MESSAGE(newcnt << " " << newLoc);
327       newLocations->SetValue(alreadyCopied, newLoc);
328       alreadyCopied++;
329     }
330 }
331
332 int SMDS_UnstructuredGrid::CellIdToDownId(int vtkCellId)
333 {
334   if((vtkCellId < 0) || (vtkCellId >= _cellIdToDownId.size()))
335     {
336       //MESSAGE("SMDS_UnstructuredGrid::CellIdToDownId structure not up to date: vtkCellId="
337       //    << vtkCellId << " max="<< _cellIdToDownId.size());
338       return -1;
339     }
340   return _cellIdToDownId[vtkCellId];
341 }
342
343 void SMDS_UnstructuredGrid::setCellIdToDownId(int vtkCellId, int downId)
344 {
345   // ASSERT((vtkCellId >= 0) && (vtkCellId < _cellIdToDownId.size()));
346   _cellIdToDownId[vtkCellId] = downId;
347 }
348
349 void SMDS_UnstructuredGrid::CleanDownwardConnectivity()
350 {
351   for (int i = 0; i < _downArray.size(); i++)
352     {
353       if (_downArray[i])
354         delete _downArray[i];
355       _downArray[i] = 0;
356     }
357   _cellIdToDownId.clear();
358 }
359
360 /*! Build downward connectivity: to do only when needed because heavy memory load.
361  *  Downward connectivity is no more valid if vtkUnstructuredGrid is modified.
362  *
363  */
364 void SMDS_UnstructuredGrid::BuildDownwardConnectivity(bool withEdges)
365 {
366   MESSAGE("SMDS_UnstructuredGrid::BuildDownwardConnectivity");CHRONO(2);
367   // TODO calcul partiel sans edges
368
369   // --- erase previous data if any
370
371   this->CleanDownwardConnectivity();
372
373   // --- create SMDS_Downward structures (in _downArray vector[vtkCellType])
374
375   _downArray.resize(VTK_MAXTYPE + 1, 0); // --- max. type value = VTK_QUADRATIC_PYRAMID
376
377   _downArray[VTK_LINE] = new SMDS_DownEdge(this);
378   _downArray[VTK_QUADRATIC_EDGE] = new SMDS_DownQuadEdge(this);
379   _downArray[VTK_TRIANGLE] = new SMDS_DownTriangle(this);
380   _downArray[VTK_QUADRATIC_TRIANGLE] = new SMDS_DownQuadTriangle(this);
381   _downArray[VTK_QUAD] = new SMDS_DownQuadrangle(this);
382   _downArray[VTK_QUADRATIC_QUAD] = new SMDS_DownQuadQuadrangle(this);
383   _downArray[VTK_TETRA] = new SMDS_DownTetra(this);
384   _downArray[VTK_QUADRATIC_TETRA] = new SMDS_DownQuadTetra(this);
385   _downArray[VTK_PYRAMID] = new SMDS_DownPyramid(this);
386   _downArray[VTK_QUADRATIC_PYRAMID] = new SMDS_DownQuadPyramid(this);
387   _downArray[VTK_WEDGE] = new SMDS_DownPenta(this);
388   _downArray[VTK_QUADRATIC_WEDGE] = new SMDS_DownQuadPenta(this);
389   _downArray[VTK_HEXAHEDRON] = new SMDS_DownHexa(this);
390   _downArray[VTK_QUADRATIC_HEXAHEDRON] = new SMDS_DownQuadHexa(this);
391
392   // --- get detailed info of number of cells of each type, allocate SMDS_downward structures
393
394   const SMDS_MeshInfo &meshInfo = _mesh->GetMeshInfo();
395
396   int nbLinTetra = meshInfo.NbTetras(ORDER_LINEAR);
397   int nbQuadTetra = meshInfo.NbTetras(ORDER_QUADRATIC);
398   int nbLinPyra = meshInfo.NbPyramids(ORDER_LINEAR);
399   int nbQuadPyra = meshInfo.NbPyramids(ORDER_QUADRATIC);
400   int nbLinPrism = meshInfo.NbPrisms(ORDER_LINEAR);
401   int nbQuadPrism = meshInfo.NbPrisms(ORDER_QUADRATIC);
402   int nbLinHexa = meshInfo.NbHexas(ORDER_LINEAR);
403   int nbQuadHexa = meshInfo.NbHexas(ORDER_QUADRATIC);
404
405   int nbLineGuess = int((4.0 / 3.0) * nbLinTetra + 2 * nbLinPrism + 2.5 * nbLinPyra + 3 * nbLinHexa);
406   int nbQuadEdgeGuess = int((4.0 / 3.0) * nbQuadTetra + 2 * nbQuadPrism + 2.5 * nbQuadPyra + 3 * nbQuadHexa);
407   int nbLinTriaGuess = 2 * nbLinTetra + nbLinPrism + 2 * nbLinPyra;
408   int nbQuadTriaGuess = 2 * nbQuadTetra + nbQuadPrism + 2 * nbQuadPyra;
409   int nbLinQuadGuess = int((2.0 / 3.0) * nbLinPrism + (1.0 / 2.0) * nbLinPyra + 3 * nbLinHexa);
410   int nbQuadQuadGuess = int((2.0 / 3.0) * nbQuadPrism + (1.0 / 2.0) * nbQuadPyra + 3 * nbQuadHexa);
411
412   int GuessSize[VTK_QUADRATIC_TETRA];
413   GuessSize[VTK_LINE] = nbLineGuess;
414   GuessSize[VTK_QUADRATIC_EDGE] = nbQuadEdgeGuess;
415   GuessSize[VTK_TRIANGLE] = nbLinTriaGuess;
416   GuessSize[VTK_QUADRATIC_TRIANGLE] = nbQuadTriaGuess;
417   GuessSize[VTK_QUAD] = nbLinQuadGuess;
418   GuessSize[VTK_QUADRATIC_QUAD] = nbQuadQuadGuess;
419   GuessSize[VTK_TETRA] = nbLinTetra;
420   GuessSize[VTK_QUADRATIC_TETRA] = nbQuadTetra;
421   GuessSize[VTK_PYRAMID] = nbLinPyra;
422   GuessSize[VTK_QUADRATIC_PYRAMID] = nbQuadPyra;
423   GuessSize[VTK_WEDGE] = nbLinPrism;
424   GuessSize[VTK_QUADRATIC_WEDGE] = nbQuadPrism;
425   GuessSize[VTK_HEXAHEDRON] = nbLinHexa;
426   GuessSize[VTK_QUADRATIC_HEXAHEDRON] = nbQuadHexa;
427
428   _downArray[VTK_LINE]->allocate(nbLineGuess);
429   _downArray[VTK_QUADRATIC_EDGE]->allocate(nbQuadEdgeGuess);
430   _downArray[VTK_TRIANGLE]->allocate(nbLinTriaGuess);
431   _downArray[VTK_QUADRATIC_TRIANGLE]->allocate(nbQuadTriaGuess);
432   _downArray[VTK_QUAD]->allocate(nbLinQuadGuess);
433   _downArray[VTK_QUADRATIC_QUAD]->allocate(nbQuadQuadGuess);
434   _downArray[VTK_TETRA]->allocate(nbLinTetra);
435   _downArray[VTK_QUADRATIC_TETRA]->allocate(nbQuadTetra);
436   _downArray[VTK_PYRAMID]->allocate(nbLinPyra);
437   _downArray[VTK_QUADRATIC_PYRAMID]->allocate(nbQuadPyra);
438   _downArray[VTK_WEDGE]->allocate(nbLinPrism);
439   _downArray[VTK_QUADRATIC_WEDGE]->allocate(nbQuadPrism);
440   _downArray[VTK_HEXAHEDRON]->allocate(nbLinHexa);
441   _downArray[VTK_QUADRATIC_HEXAHEDRON]->allocate(nbQuadHexa);
442
443   // --- iteration on vtkUnstructuredGrid cells, only faces
444   //     for each vtk face:
445   //       create a downward face entry with its downward id.
446   //       compute vtk volumes, create downward volumes entry.
447   //       mark face in downward volumes
448   //       mark volumes in downward face
449
450   MESSAGE("--- iteration on vtkUnstructuredGrid cells, only faces");CHRONO(20);
451   int cellSize = this->Types->GetNumberOfTuples();
452   _cellIdToDownId.resize(cellSize, -1);
453
454   for (int i = 0; i < cellSize; i++)
455     {
456       int vtkFaceType = this->GetCellType(i);
457       if (SMDS_Downward::getCellDimension(vtkFaceType) == 2)
458         {
459           int vtkFaceId = i;
460           //ASSERT(_downArray[vtkFaceType]);
461           int connFaceId = _downArray[vtkFaceType]->addCell(vtkFaceId);
462           SMDS_Down2D* downFace = static_cast<SMDS_Down2D*> (_downArray[vtkFaceType]);
463           downFace->setTempNodes(connFaceId, vtkFaceId);
464           int vols[2] = { -1, -1 };
465           int nbVolumes = downFace->computeVolumeIds(vtkFaceId, vols);
466           //MESSAGE("nbVolumes="<< nbVolumes);
467           for (int ivol = 0; ivol < nbVolumes; ivol++)
468             {
469               int vtkVolId = vols[ivol];
470               int vtkVolType = this->GetCellType(vtkVolId);
471               //ASSERT(_downArray[vtkVolType]);
472               int connVolId = _downArray[vtkVolType]->addCell(vtkVolId);
473               _downArray[vtkVolType]->addDownCell(connVolId, connFaceId, vtkFaceType);
474               _downArray[vtkFaceType]->addUpCell(connFaceId, connVolId, vtkVolType);
475               // MESSAGE("Face " << vtkFaceId << " belongs to volume " << vtkVolId);
476             }
477         }
478     }
479
480   // --- iteration on vtkUnstructuredGrid cells, only volumes
481   //     for each vtk volume:
482   //       create downward volumes entry if not already done
483   //       build a temporary list of faces described with their nodes
484   //       for each face
485   //         compute the vtk volumes containing this face
486   //         check if the face is already listed in the volumes (comparison of ordered list of nodes)
487   //         if not, create a downward face entry (resizing of structure required)
488   //         (the downward faces store a temporary list of nodes to ease the comparison)
489   //         create downward volumes entry if not already done
490   //         mark volumes in downward face
491   //         mark face in downward volumes
492
493   CHRONOSTOP(20);
494   MESSAGE("--- iteration on vtkUnstructuredGrid cells, only volumes");CHRONO(21);
495
496   for (int i = 0; i < cellSize; i++)
497     {
498       int vtkType = this->GetCellType(i);
499       if (SMDS_Downward::getCellDimension(vtkType) == 3)
500         {
501           //CHRONO(31);
502           int vtkVolId = i;
503           // MESSAGE("vtk volume " << vtkVolId);
504           //ASSERT(_downArray[vtkType]);
505           /*int connVolId = */_downArray[vtkType]->addCell(vtkVolId);
506
507           // --- find all the faces of the volume, describe the faces by their nodes
508
509           SMDS_Down3D* downVol = static_cast<SMDS_Down3D*> (_downArray[vtkType]);
510           ListElemByNodesType facesWithNodes;
511           downVol->computeFacesWithNodes(vtkVolId, facesWithNodes);
512           // MESSAGE("vtk volume " << vtkVolId << " contains " << facesWithNodes.nbElems << " faces");
513           //CHRONOSTOP(31);
514           for (int iface = 0; iface < facesWithNodes.nbElems; iface++)
515             {
516               // --- find the volumes containing the face
517
518               //CHRONO(32);
519               int vtkFaceType = facesWithNodes.elems[iface].vtkType;
520               SMDS_Down2D* downFace = static_cast<SMDS_Down2D*> (_downArray[vtkFaceType]);
521               int vols[2] = { -1, -1 };
522               int *nodes = &facesWithNodes.elems[iface].nodeIds[0];
523               int lg = facesWithNodes.elems[iface].nbNodes;
524               int nbVolumes = downFace->computeVolumeIdsFromNodesFace(nodes, lg, vols);
525               // MESSAGE("vtk volume " << vtkVolId << " face " << iface << " belongs to " << nbVolumes << " volumes");
526
527               // --- check if face is registered in the volumes
528               //CHRONOSTOP(32);
529
530               //CHRONO(33);
531               int connFaceId = -1;
532               for (int ivol = 0; ivol < nbVolumes; ivol++)
533                 {
534                   int vtkVolId2 = vols[ivol];
535                   int vtkVolType = this->GetCellType(vtkVolId2);
536                   //ASSERT(_downArray[vtkVolType]);
537                   int connVolId2 = _downArray[vtkVolType]->addCell(vtkVolId2);
538                   SMDS_Down3D* downVol2 = static_cast<SMDS_Down3D*> (_downArray[vtkVolType]);
539                   connFaceId = downVol2->FindFaceByNodes(connVolId2, facesWithNodes.elems[iface]);
540                   if (connFaceId >= 0)
541                     break; // --- face already created
542                 }//CHRONOSTOP(33);
543
544               // --- if face is not registered in the volumes, create face
545
546               //CHRONO(34);
547               if (connFaceId < 0)
548                 {
549                   connFaceId = _downArray[vtkFaceType]->addCell();
550                   downFace->setTempNodes(connFaceId, facesWithNodes.elems[iface]);
551                 }//CHRONOSTOP(34);
552
553               // --- mark volumes in downward face and mark face in downward volumes
554
555               //CHRONO(35);
556               for (int ivol = 0; ivol < nbVolumes; ivol++)
557                 {
558                   int vtkVolId2 = vols[ivol];
559                   int vtkVolType = this->GetCellType(vtkVolId2);
560                   //ASSERT(_downArray[vtkVolType]);
561                   int connVolId2 = _downArray[vtkVolType]->addCell(vtkVolId2);
562                   _downArray[vtkVolType]->addDownCell(connVolId2, connFaceId, vtkFaceType);
563                   _downArray[vtkFaceType]->addUpCell(connFaceId, connVolId2, vtkVolType);
564                   // MESSAGE(" From volume " << vtkVolId << " face " << connFaceId << " belongs to volume " << vtkVolId2);
565                 }//CHRONOSTOP(35);
566             }
567         }
568     }
569
570   // --- iteration on vtkUnstructuredGrid cells, only edges
571   //     for each vtk edge:
572   //       create downward edge entry
573   //       store the nodes id's in downward edge (redundant with vtkUnstructuredGrid)
574   //       find downward faces
575   //       (from vtk faces or volumes, get downward faces, they have a temporary list of nodes)
576   //       mark edge in downward faces
577   //       mark faces in downward edge
578
579   CHRONOSTOP(21);
580   MESSAGE("--- iteration on vtkUnstructuredGrid cells, only edges");CHRONO(22);
581
582   for (int i = 0; i < cellSize; i++)
583     {
584       int vtkEdgeType = this->GetCellType(i);
585       if (SMDS_Downward::getCellDimension(vtkEdgeType) == 1)
586         {
587           int vtkEdgeId = i;
588           //ASSERT(_downArray[vtkEdgeType]);
589           int connEdgeId = _downArray[vtkEdgeType]->addCell(vtkEdgeId);
590           SMDS_Down1D* downEdge = static_cast<SMDS_Down1D*> (_downArray[vtkEdgeType]);
591           downEdge->setNodes(connEdgeId, vtkEdgeId);
592           vector<int> vtkIds;
593           int nbVtkCells = downEdge->computeVtkCells(connEdgeId, vtkIds);
594           int downFaces[1000];
595           unsigned char downTypes[1000];
596           int nbDownFaces = downEdge->computeFaces(connEdgeId, &vtkIds[0], nbVtkCells, downFaces, downTypes);
597           for (int n = 0; n < nbDownFaces; n++)
598             {
599               _downArray[downTypes[n]]->addDownCell(downFaces[n], connEdgeId, vtkEdgeType);
600               _downArray[vtkEdgeType]->addUpCell(connEdgeId, downFaces[n], downTypes[n]);
601             }
602         }
603     }
604
605   // --- iteration on downward faces (they are all listed now)
606   //     for each downward face:
607   //       build a temporary list of edges with their ordered list of nodes
608   //       for each edge:
609   //         find all the vtk cells containing this edge
610   //         then identify all the downward faces containing the edge, from the vtk cells
611   //         check if the edge is already listed in the faces (comparison of ordered list of nodes)
612   //         if not, create a downward edge entry with the node id's
613   //         mark edge in downward faces
614   //         mark downward faces in edge (size of list unknown, to be allocated)
615
616   CHRONOSTOP(22);CHRONO(23);
617
618   for (int vtkFaceType = 0; vtkFaceType < VTK_QUADRATIC_PYRAMID; vtkFaceType++)
619     {
620       if (SMDS_Downward::getCellDimension(vtkFaceType) != 2)
621         continue;
622
623       // --- find all the edges of the face, describe the edges by their nodes
624
625       SMDS_Down2D* downFace = static_cast<SMDS_Down2D*> (_downArray[vtkFaceType]);
626       int maxId = downFace->getMaxId();
627       for (int faceId = 0; faceId < maxId; faceId++)
628         {
629           //CHRONO(40);
630           ListElemByNodesType edgesWithNodes;
631           downFace->computeEdgesWithNodes(faceId, edgesWithNodes);
632           // MESSAGE("downward face type " << vtkFaceType << " num " << faceId << " contains " << edgesWithNodes.nbElems << " edges");
633
634           //CHRONOSTOP(40);
635           for (int iedge = 0; iedge < edgesWithNodes.nbElems; iedge++)
636             {
637
638               // --- check if the edge is already registered by exploration of the faces
639
640               //CHRONO(41);
641               vector<int> vtkIds;
642               unsigned char vtkEdgeType = edgesWithNodes.elems[iedge].vtkType;
643               int *pts = &edgesWithNodes.elems[iedge].nodeIds[0];
644               SMDS_Down1D* downEdge = static_cast<SMDS_Down1D*> (_downArray[vtkEdgeType]);
645               int nbVtkCells = downEdge->computeVtkCells(pts, vtkIds);
646               //CHRONOSTOP(41);CHRONO(42);
647               int downFaces[1000];
648               unsigned char downTypes[1000];
649               int nbDownFaces = downEdge->computeFaces(pts, &vtkIds[0], nbVtkCells, downFaces, downTypes);
650               //CHRONOSTOP(42);
651
652               //CHRONO(43);
653               int connEdgeId = -1;
654               for (int idf = 0; idf < nbDownFaces; idf++)
655                 {
656                   int faceId2 = downFaces[idf];
657                   int faceType = downTypes[idf];
658                   //ASSERT(_downArray[faceType]);
659                   SMDS_Down2D* downFace2 = static_cast<SMDS_Down2D*> (_downArray[faceType]);
660                   connEdgeId = downFace2->FindEdgeByNodes(faceId2, edgesWithNodes.elems[iedge]);
661                   if (connEdgeId >= 0)
662                     break; // --- edge already created
663                 }//CHRONOSTOP(43);
664
665               // --- if edge is not registered in the faces, create edge
666
667               if (connEdgeId < 0)
668                 {
669                   //CHRONO(44);
670                   connEdgeId = _downArray[vtkEdgeType]->addCell();
671                   downEdge->setNodes(connEdgeId, edgesWithNodes.elems[iedge].nodeIds);
672                   //CHRONOSTOP(44);
673                 }
674
675               // --- mark faces in downward edge and mark edge in downward faces
676
677               //CHRONO(45);
678               for (int idf = 0; idf < nbDownFaces; idf++)
679                 {
680                   int faceId2 = downFaces[idf];
681                   int faceType = downTypes[idf];
682                   //ASSERT(_downArray[faceType]);
683                   //SMDS_Down2D* downFace2 = static_cast<SMDS_Down2D*> (_downArray[faceType]);
684                   _downArray[vtkEdgeType]->addUpCell(connEdgeId, faceId2, faceType);
685                   _downArray[faceType]->addDownCell(faceId2, connEdgeId, vtkEdgeType);
686                   // MESSAGE(" From face t:" << vtkFaceType << " " << faceId <<
687                   //  " edge " << connEdgeId << " belongs to face t:" << faceType << " " << faceId2);
688                 }//CHRONOSTOP(45);
689             }
690         }
691     }
692
693   CHRONOSTOP(23);CHRONO(24);
694
695   // compact downward connectivity structure: adjust downward arrays size, replace vector<vector int>> by a single vector<int>
696   // 3D first then 2D and last 1D to release memory before edge upCells reorganization, (temporary memory use)
697
698   for (int vtkType = VTK_QUADRATIC_PYRAMID; vtkType >= 0; vtkType--)
699     {
700       if (SMDS_Downward *down = _downArray[vtkType])
701         {
702           down->compactStorage();
703         }
704     }
705
706   // --- Statistics
707
708   for (int vtkType = 0; vtkType <= VTK_QUADRATIC_PYRAMID; vtkType++)
709     {
710       if (SMDS_Downward *down = _downArray[vtkType])
711         {
712           if (down->getMaxId())
713             {
714               MESSAGE("Cells of Type " << vtkType << " : number of entities, est: "
715                   << GuessSize[vtkType] << " real: " << down->getMaxId());
716             }
717         }
718     }CHRONOSTOP(24);CHRONOSTOP(2);
719   counters::stats();
720 }
721
722 /*! Get the neighbors of a cell.
723  * Only the neighbors having the dimension of the cell are taken into account
724  * (neighbors of a volume are the volumes sharing a face with this volume,
725  *  neighbors of a face are the faces sharing an edge with this face...).
726  * @param neighborsVtkIds vector of neighbors vtk id's to fill (reserve enough space).
727  * @param downIds downward id's of cells of dimension n-1, to fill (reserve enough space).
728  * @param downTypes vtk types of cells of dimension n-1, to fill (reserve enough space).
729  * @param vtkId the vtk id of the cell
730  * @return number of neighbors
731  */
732 int SMDS_UnstructuredGrid::GetNeighbors(int* neighborsVtkIds, int* downIds, unsigned char* downTypes, int vtkId)
733 {
734   int vtkType = this->GetCellType(vtkId);
735   int cellDim = SMDS_Downward::getCellDimension(vtkType);
736   if (cellDim <2)
737     return 0; // TODO voisins des edges = edges connectees
738   int cellId = this->_cellIdToDownId[vtkId];
739
740   int nbCells = _downArray[vtkType]->getNumberOfDownCells(cellId);
741   const int *downCells = _downArray[vtkType]->getDownCells(cellId);
742   const unsigned char* downTyp = _downArray[vtkType]->getDownTypes(cellId);
743
744   // --- iteration on faces of the 3D cell (or edges on the 2D cell).
745
746   int nb = 0;
747   for (int i = 0; i < nbCells; i++)
748     {
749       int downId = downCells[i];
750       int cellType = downTyp[i];
751       int nbUp = _downArray[cellType]->getNumberOfUpCells(downId);
752       const int *upCells = _downArray[cellType]->getUpCells(downId);
753       const unsigned char* upTypes = _downArray[cellType]->getUpTypes(downId);
754
755       // ---for a volume, max 2 upCells, one is this cell, the other is a neighbor
756       //    for a face, number of neighbors (connected faces) not known
757
758       for (int j = 0; j < nbUp; j++)
759         {
760           if ((upCells[j] == cellId) && (upTypes[j] == vtkType))
761             continue;
762           int vtkNeighbor = _downArray[upTypes[j]]->getVtkCellId(upCells[j]);
763           neighborsVtkIds[nb] = vtkNeighbor;
764           downIds[nb] = downId;
765           downTypes[nb] = cellType;
766           nb++;
767         }
768       if (nb >= NBMAXNEIGHBORS)
769         assert(0);
770     }
771   return nb;
772 }
773
774 /*! get the volumes containing a face or an edge of the grid
775  * The edge or face belongs to the vtkUnstructuredGrid
776  * @param volVtkIds vector of parent volume ids to fill (reserve enough space!)
777  * @param vtkId vtk id of the face or edge
778  */
779 int SMDS_UnstructuredGrid::GetParentVolumes(int* volVtkIds, int vtkId)
780 {
781   int vtkType = this->GetCellType(vtkId);
782   int dim = SMDS_Downward::getCellDimension(vtkType);
783   int nbFaces = 0;
784   unsigned char cellTypes[1000];
785   int downCellId[1000];
786   if (dim == 1)
787     {
788       int downId = this->CellIdToDownId(vtkId);
789       if (downId < 0)
790         {
791           MESSAGE("Downward structure not up to date: new edge not taken into account");
792           return 0;
793         }
794       nbFaces = _downArray[vtkType]->getNumberOfUpCells(downId);
795       const int *upCells = _downArray[vtkType]->getUpCells(downId);
796       const unsigned char* upTypes = _downArray[vtkType]->getUpTypes(downId);
797       for (int i=0; i< nbFaces; i++)
798         {
799           cellTypes[i] = upTypes[i];
800           downCellId[i] = upCells[i];
801         }
802     }
803   else if (dim == 2)
804     {
805       nbFaces = 1;
806       cellTypes[0] = this->GetCellType(vtkId);
807       int downId = this->CellIdToDownId(vtkId);
808       if (downId < 0)
809         {
810           MESSAGE("Downward structure not up to date: new face not taken into account");
811           return 0;
812         }
813       downCellId[0] = downId;
814     }
815
816   int nbvol =0;
817   for (int i=0; i<nbFaces; i++)
818     {
819       int vtkTypeFace = cellTypes[i];
820       int downId = downCellId[i];
821       int nv = _downArray[vtkTypeFace]->getNumberOfUpCells(downId);
822       const int *upCells = _downArray[vtkTypeFace]->getUpCells(downId);
823       const unsigned char* upTypes = _downArray[vtkTypeFace]->getUpTypes(downId);
824        for (int j=0; j<nv; j++)
825         {
826           int vtkVolId = _downArray[upTypes[j]]->getVtkCellId(upCells[j]);
827           if (vtkVolId >= 0)
828             volVtkIds[nbvol++] = vtkVolId;
829         }
830     }
831   return nbvol;
832 }
833
834 /*! get the volumes containing a face or an edge of the downward structure
835  * The edge or face does not necessary belong to the vtkUnstructuredGrid
836  * @param volVtkIds vector of parent volume ids to fill (reserve enough space!)
837  * @param downId id in the downward structure
838  * @param downType type of cell
839  */
840 int SMDS_UnstructuredGrid::GetParentVolumes(int* volVtkIds, int downId, unsigned char downType)
841 {
842   int vtkType = downType;
843   int dim = SMDS_Downward::getCellDimension(vtkType);
844   int nbFaces = 0;
845   unsigned char cellTypes[1000];
846   int downCellId[1000];
847   if (dim == 1)
848     {
849       nbFaces = _downArray[vtkType]->getNumberOfUpCells(downId);
850       const int *upCells = _downArray[vtkType]->getUpCells(downId);
851       const unsigned char* upTypes = _downArray[vtkType]->getUpTypes(downId);
852       for (int i=0; i< nbFaces; i++)
853         {
854           cellTypes[i] = upTypes[i];
855           downCellId[i] = upCells[i];
856         }
857     }
858   else if (dim == 2)
859     {
860       nbFaces = 1;
861       cellTypes[0] = vtkType;
862       downCellId[0] = downId;
863     }
864
865   int nbvol =0;
866   for (int i=0; i<nbFaces; i++)
867     {
868       int vtkTypeFace = cellTypes[i];
869       int downId = downCellId[i];
870       int nv = _downArray[vtkTypeFace]->getNumberOfUpCells(downId);
871       const int *upCells = _downArray[vtkTypeFace]->getUpCells(downId);
872       const unsigned char* upTypes = _downArray[vtkTypeFace]->getUpTypes(downId);
873        for (int j=0; j<nv; j++)
874         {
875           int vtkVolId = _downArray[upTypes[j]]->getVtkCellId(upCells[j]);
876           if (vtkVolId >= 0)
877             volVtkIds[nbvol++] = vtkVolId;
878         }
879     }
880   return nbvol;
881 }
882
883 /*! get the node id's of a cell.
884  * The cell is defined by it's downward connectivity id and type.
885  * @param nodeSet set of of vtk node id's to fill.
886  * @param downId downward connectivity id of the cell.
887  * @param downType type of cell.
888  */
889 void SMDS_UnstructuredGrid::GetNodeIds(std::set<int>& nodeSet, int downId, unsigned char downType)
890 {
891   _downArray[downType]->getNodeIds(downId, nodeSet);
892 }
893
894 /*! change some nodes in cell without modifying type or internal connectivity.
895  * Nodes inverse connectivity is maintained up to date.
896  * @param vtkVolId vtk id of the cell
897  * @param localClonedNodeIds map old node id to new node id.
898  */
899 void SMDS_UnstructuredGrid::ModifyCellNodes(int vtkVolId, std::map<int, int> localClonedNodeIds)
900 {
901   vtkIdType npts = 0;
902   vtkIdType *pts; // will refer to the point id's of the face
903   this->GetCellPoints(vtkVolId, npts, pts);
904   for (int i = 0; i < npts; i++)
905     {
906       if (localClonedNodeIds.count(pts[i]))
907         {
908           vtkIdType oldpt = pts[i];
909           pts[i] = localClonedNodeIds[oldpt];
910           //MESSAGE(oldpt << " --> " << pts[i]);
911           //this->RemoveReferenceToCell(oldpt, vtkVolId);
912           //this->AddReferenceToCell(pts[i], vtkVolId);
913         }
914     }
915 }
916
917 /*! reorder the nodes of a face
918  * @param vtkVolId vtk id of a volume containing the face, to get an orientation for the face.
919  * @param orderedNodes list of nodes to reorder (in out)
920  * @return size of the list
921  */
922 int SMDS_UnstructuredGrid::getOrderedNodesOfFace(int vtkVolId, std::vector<vtkIdType>& orderedNodes)
923 {
924   int vtkType = this->GetCellType(vtkVolId);
925   int cellDim = SMDS_Downward::getCellDimension(vtkType);
926   if (cellDim != 3)
927     return 0;
928   SMDS_Down3D *downvol = static_cast<SMDS_Down3D*> (_downArray[vtkType]);
929   int downVolId = this->_cellIdToDownId[vtkVolId];
930   downvol->getOrderedNodesOfFace(downVolId, orderedNodes);
931   return orderedNodes.size();
932 }
933
934 void SMDS_UnstructuredGrid::BuildLinks()
935 {
936   // Remove the old links if they are already built
937   if (this->Links)
938     {
939     this->Links->UnRegister(this);
940     }
941
942   this->Links = SMDS_CellLinks::New();
943   this->Links->Allocate(this->GetNumberOfPoints());
944   this->Links->Register(this);
945   this->Links->BuildLinks(this, this->Connectivity);
946   this->Links->Delete();
947 }
948
949 /*! Create a volume (prism or hexahedron) by duplication of a face.
950  * Designed for use in creation of flat elements separating volume domains.
951  * A face separating two domains is shared by two volume cells.
952  * All the nodes are already created (for the two faces).
953  * Each original Node is associated to corresponding nodes in the domains.
954  * Some nodes may be duplicated for more than two domains, when domain separations intersect.
955  * In that case, even some of the nodes to use for the original face may be changed.
956  * @param vtkVolId: vtk id of a volume containing the face, to get an orientation for the face.
957  * @param domain1: domain of the original face
958  * @param domain2: domain of the duplicated face
959  * @param originalNodes: the vtk node ids of the original face
960  * @param nodeDomains: map(original id --> map(domain --> duplicated node id))
961  * @return ok if success.
962  */
963 SMDS_MeshVolume* SMDS_UnstructuredGrid::extrudeVolumeFromFace(int vtkVolId,
964                                                   int domain1,
965                                                   int domain2,
966                                                   std::set<int>& originalNodes,
967                                                   std::map<int, std::map<int, int> >& nodeDomains,
968                                                   std::map<int, std::map<long, int> >& nodeQuadDomains)
969 {
970   //MESSAGE("extrudeVolumeFromFace " << vtkVolId);
971   vector<vtkIdType> orderedOriginals;
972   orderedOriginals.clear();
973   set<int>::const_iterator it = originalNodes.begin();
974   for (; it != originalNodes.end(); ++it)
975     orderedOriginals.push_back(*it);
976
977   int nbNodes = this->getOrderedNodesOfFace(vtkVolId, orderedOriginals);
978   vector<vtkIdType> orderedNodes;
979
980   switch (orderedOriginals.size())
981   {
982     case 3:
983     case 4:
984       for (int i = 0; i < nbNodes; i++)
985         orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
986       for (int i = 0; i < nbNodes; i++)
987         orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
988       break;
989     case 6:
990     case 8:
991       {
992         long dom1 = domain1;
993         long dom2 = domain2;
994         long dom1_2; // for nodeQuadDomains
995         if (domain1 < domain2)
996           dom1_2 = dom1 + INT_MAX * dom2;
997         else
998           dom1_2 = dom2 + INT_MAX * dom1;
999         //cerr << "dom1=" << dom1 << " dom2=" << dom2 << " dom1_2=" << dom1_2 << endl;
1000         int ima = orderedOriginals.size();
1001         int mid = orderedOriginals.size() / 2;
1002         //cerr << "ima=" << ima << " mid=" << mid << endl;
1003         for (int i = 0; i < mid; i++)
1004           orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
1005         for (int i = 0; i < mid; i++)
1006           orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
1007         for (int i = mid; i < ima; i++)
1008           orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain1]);
1009         for (int i = mid; i < ima; i++)
1010           orderedNodes.push_back(nodeDomains[orderedOriginals[i]][domain2]);
1011         for (int i = 0; i < mid; i++)
1012           {
1013             int oldId = orderedOriginals[i];
1014             int newId;
1015             if (nodeQuadDomains.count(oldId) && nodeQuadDomains[oldId].count(dom1_2))
1016               newId = nodeQuadDomains[oldId][dom1_2];
1017             else
1018               {
1019                 double *coords = this->GetPoint(oldId);
1020                 SMDS_MeshNode *newNode = _mesh->AddNode(coords[0], coords[1], coords[2]);
1021                 newId = newNode->getVtkId();
1022                 std::map<long, int> emptyMap;
1023                 nodeQuadDomains[oldId] = emptyMap;
1024                 nodeQuadDomains[oldId][dom1_2] = newId;
1025               }
1026             orderedNodes.push_back(newId);
1027           }
1028       }
1029       break;
1030     default:
1031       ASSERT(0);
1032   }
1033
1034   SMDS_MeshVolume *vol = _mesh->AddVolumeFromVtkIds(orderedNodes);
1035
1036   // TODO update subshape list of elements and nodes
1037   return vol;
1038 }