]> SALOME platform Git repositories - modules/smesh.git/blob - src/OBJECT/SMESH_ExtractGeometry.cxx
Salome HOME
Porting on Mandriva 64
[modules/smesh.git] / src / OBJECT / SMESH_ExtractGeometry.cxx
1 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
2 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS 
3 // 
4 //  This library is free software; you can redistribute it and/or 
5 //  modify it under the terms of the GNU Lesser General Public 
6 //  License as published by the Free Software Foundation; either 
7 //  version 2.1 of the License. 
8 // 
9 //  This library is distributed in the hope that it will be useful, 
10 //  but WITHOUT ANY WARRANTY; without even the implied warranty of 
11 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
12 //  Lesser General Public License for more details. 
13 // 
14 //  You should have received a copy of the GNU Lesser General Public 
15 //  License along with this library; if not, write to the Free Software 
16 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA 
17 // 
18 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org 
19
20
21 #include "SMESH_ExtractGeometry.h"
22
23 #include "vtkCell.h"
24 #include "vtkCellData.h"
25 #include "vtkFloatArray.h"
26 #include "vtkIdList.h"
27 #include "vtkImplicitFunction.h"
28 #include "vtkObjectFactory.h"
29 #include "vtkPointData.h"
30 #include "vtkUnstructuredGrid.h"
31
32 using namespace std;
33
34 #ifdef _DEBUG_
35 static int MYDEBUG = 0;
36 #else
37 static int MYDEBUG = 0;
38 #endif
39
40 #if defined __GNUC__
41   #if __GNUC__ == 2
42     #define __GNUC_2__
43   #endif
44 #endif
45
46
47 vtkStandardNewMacro(SMESH_ExtractGeometry);
48
49
50 SMESH_ExtractGeometry::SMESH_ExtractGeometry()
51 {}
52
53
54 SMESH_ExtractGeometry::~SMESH_ExtractGeometry(){}
55
56
57 vtkIdType SMESH_ExtractGeometry::GetElemObjId(int theVtkID){
58   if(myElemVTK2ObjIds.empty() || theVtkID > myElemVTK2ObjIds.size()) return -1;
59 #if defined __GNUC_2__
60   return myElemVTK2ObjIds[theVtkID];
61 #else
62   return myElemVTK2ObjIds.at(theVtkID);
63 #endif
64 }
65
66
67 vtkIdType SMESH_ExtractGeometry::GetNodeObjId(int theVtkID){
68   if(myNodeVTK2ObjIds.empty() || theVtkID > myNodeVTK2ObjIds.size()) return -1;
69 #if defined __GNUC_2__
70   return myNodeVTK2ObjIds[theVtkID];
71 #else
72   return myNodeVTK2ObjIds.at(theVtkID);
73 #endif
74 }
75
76
77 void SMESH_ExtractGeometry::Execute()
78 {
79   vtkIdType ptId, numPts, numCells, i, cellId, newCellId, newId, *pointMap;
80   vtkIdList *cellPts;
81   vtkCell *cell;
82   int numCellPts;
83   vtkFloatingPointType *x;
84   vtkFloatingPointType multiplier;
85   vtkPoints *newPts;
86   vtkIdList *newCellPts;
87   vtkDataSet *input = this->GetInput();
88   vtkPointData *pd = input->GetPointData();
89   vtkCellData *cd = input->GetCellData();
90   vtkUnstructuredGrid *output = this->GetOutput();
91   vtkPointData *outputPD = output->GetPointData();
92   vtkCellData *outputCD = output->GetCellData();
93   int npts;
94   numCells = input->GetNumberOfCells();
95   numPts = input->GetNumberOfPoints();
96   
97   vtkDebugMacro(<< "Extracting geometry");
98
99   if ( ! this->ImplicitFunction )
100     {
101     vtkErrorMacro(<<"No implicit function specified");
102     return;
103     }
104
105   newCellPts = vtkIdList::New();
106   newCellPts->Allocate(VTK_CELL_SIZE);
107
108   if ( this->ExtractInside )
109     {
110     multiplier = 1.0;
111     }
112   else 
113     {
114     multiplier = -1.0;
115     }
116
117   // Loop over all points determining whether they are inside the
118   // implicit function. Copy the points and point data if they are.
119   //
120   pointMap = new vtkIdType[numPts]; // maps old point ids into new
121   for (i=0; i < numPts; i++)
122     {
123     pointMap[i] = -1;
124     }
125
126   output->Allocate(numCells/4); //allocate storage for geometry/topology
127   newPts = vtkPoints::New();
128   newPts->Allocate(numPts/4,numPts);
129   outputPD->CopyAllocate(pd);
130   outputCD->CopyAllocate(cd);
131   vtkFloatArray *newScalars = NULL;
132   
133   if(myStoreMapping){
134     myElemVTK2ObjIds.clear();
135     myElemVTK2ObjIds.reserve(numCells);
136     myNodeVTK2ObjIds.clear();
137     myNodeVTK2ObjIds.reserve(numPts);
138   }
139
140   if ( ! this->ExtractBoundaryCells )
141     {
142     for ( ptId=0; ptId < numPts; ptId++ )
143       {
144       x = input->GetPoint(ptId);
145       if ( (this->ImplicitFunction->FunctionValue(x)*multiplier) < 0.0 )
146         {
147         newId = newPts->InsertNextPoint(x);
148         pointMap[ptId] = newId;
149         myNodeVTK2ObjIds.push_back(ptId);
150         outputPD->CopyData(pd,ptId,newId);
151         }
152       }
153     }
154   else
155     {
156     // To extract boundary cells, we have to create supplemental information
157     if ( this->ExtractBoundaryCells )
158       {
159       vtkFloatingPointType val;
160       newScalars = vtkFloatArray::New();
161       newScalars->SetNumberOfValues(numPts);
162
163       for (ptId=0; ptId < numPts; ptId++ )
164         {
165         x = input->GetPoint(ptId);
166         val = this->ImplicitFunction->FunctionValue(x) * multiplier;
167         newScalars->SetValue(ptId, val);
168         if ( val < 0.0 )
169           {
170           newId = newPts->InsertNextPoint(x);
171           pointMap[ptId] = newId;
172           myNodeVTK2ObjIds.push_back(ptId);
173           outputPD->CopyData(pd,ptId,newId);
174           }
175         }
176       }
177     }
178
179   // Now loop over all cells to see whether they are inside implicit
180   // function (or on boundary if ExtractBoundaryCells is on).
181   //
182   for (cellId=0; cellId < numCells; cellId++)
183     {
184     cell = input->GetCell(cellId);
185     cellPts = cell->GetPointIds();
186     numCellPts = cell->GetNumberOfPoints();
187
188     newCellPts->Reset();
189     if ( ! this->ExtractBoundaryCells ) //requires less work
190       {
191       for ( npts=0, i=0; i < numCellPts; i++, npts++)
192         {
193         ptId = cellPts->GetId(i);
194         if ( pointMap[ptId] < 0 )
195           {
196           break; //this cell won't be inserted
197           }
198         else
199           {
200           newCellPts->InsertId(i,pointMap[ptId]);
201           }
202         }
203       } //if don't want to extract boundary cells
204     
205     else //want boundary cells
206       {
207       for ( npts=0, i=0; i < numCellPts; i++ )
208         {
209         ptId = cellPts->GetId(i);
210         if ( newScalars->GetValue(ptId) <= 0.0 )
211           {
212           npts++;
213           }
214         }
215       if ( npts > 0 )
216         {
217         for ( i=0; i < numCellPts; i++ )
218           {
219           ptId = cellPts->GetId(i);
220           if ( pointMap[ptId] < 0 )
221             {
222             x = input->GetPoint(ptId);
223             newId = newPts->InsertNextPoint(x);
224             pointMap[ptId] = newId;
225             myNodeVTK2ObjIds.push_back(ptId);
226             outputPD->CopyData(pd,ptId,newId);
227             }
228           newCellPts->InsertId(i,pointMap[ptId]);
229           }
230         }//a boundary or interior cell
231       }//if mapping boundary cells
232       
233     if ( npts >= numCellPts || (this->ExtractBoundaryCells && npts > 0) )
234       {
235       newCellId = output->InsertNextCell(cell->GetCellType(),newCellPts);
236       myElemVTK2ObjIds.push_back(cellId);
237       outputCD->CopyData(cd,cellId,newCellId);
238       }
239     }//for all cells
240
241   // Update ourselves and release memory
242   //
243   delete [] pointMap;
244   newCellPts->Delete();
245   output->SetPoints(newPts);
246   newPts->Delete();
247   
248   if ( this->ExtractBoundaryCells )
249     {
250     newScalars->Delete();
251     }
252
253   output->Squeeze();
254 }