]> SALOME platform Git repositories - modules/smesh.git/blob - src/OBJECT/SMESH_ExtractGeometry.cxx
Salome HOME
0021825: Error in the example of "Projection Algorithms" in the user's guide
[modules/smesh.git] / src / OBJECT / SMESH_ExtractGeometry.cxx
1 // Copyright (C) 2007-2012  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.
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 #include "SMESH_ExtractGeometry.h"
24
25 #include <vtkCell.h>
26 #include <vtkCellData.h>
27 #include <vtkFloatArray.h>
28 #include <vtkIdList.h>
29 #include <vtkImplicitFunction.h>
30 #include <vtkObjectFactory.h>
31 #include <vtkPointData.h>
32 #include <vtkUnstructuredGrid.h>
33 #include <vtkInformation.h>
34 #include <vtkInformationVector.h>
35
36 using namespace std;
37
38 //#ifdef _DEBUG_
39 //static int MYDEBUG = 0;
40 //#else
41 //static int MYDEBUG = 0;
42 //#endif
43
44 #if defined __GNUC__
45   #if __GNUC__ == 2
46     #define __GNUC_2__
47   #endif
48 #endif
49
50
51 vtkStandardNewMacro(SMESH_ExtractGeometry);
52
53
54 SMESH_ExtractGeometry::SMESH_ExtractGeometry()
55 {}
56
57
58 SMESH_ExtractGeometry::~SMESH_ExtractGeometry(){}
59
60
61 vtkIdType SMESH_ExtractGeometry::GetElemObjId(int theVtkID){
62   if( theVtkID < 0 || theVtkID >= myElemVTK2ObjIds.size()) return -1;
63   return myElemVTK2ObjIds[theVtkID];
64 }
65
66
67 vtkIdType SMESH_ExtractGeometry::GetNodeObjId(int theVtkID){
68   if ( theVtkID < 0 || theVtkID >= myNodeVTK2ObjIds.size()) return -1;
69   return myNodeVTK2ObjIds[theVtkID];
70 }
71
72
73 int SMESH_ExtractGeometry::RequestData(
74   vtkInformation *vtkNotUsed(request),
75   vtkInformationVector **inputVector,
76   vtkInformationVector *outputVector)
77 {
78   // get the info objects
79   vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
80   vtkInformation *outInfo = outputVector->GetInformationObject(0);
81
82   // get the input and ouptut
83   vtkDataSet *input = vtkDataSet::SafeDownCast(
84     inInfo->Get(vtkDataObject::DATA_OBJECT()));
85   vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
86     outInfo->Get(vtkDataObject::DATA_OBJECT()));
87
88   vtkIdType ptId, numPts, numCells, i, cellId, newCellId, newId, *pointMap;
89   vtkIdList *cellPts;
90   vtkCell *cell;
91   int numCellPts;
92   vtkFloatingPointType *x;
93   vtkFloatingPointType multiplier;
94   vtkPoints *newPts;
95   vtkIdList *newCellPts;
96   vtkPointData *pd = input->GetPointData();
97   vtkCellData *cd = input->GetCellData();
98   vtkPointData *outputPD = output->GetPointData();
99   vtkCellData *outputCD = output->GetCellData();
100   int npts;
101   numCells = input->GetNumberOfCells();
102   numPts = input->GetNumberOfPoints();
103   
104   vtkDebugMacro(<< "Extracting geometry");
105
106   if ( ! this->ImplicitFunction )
107     {
108     vtkErrorMacro(<<"No implicit function specified");
109     return 0;
110     }
111
112   newCellPts = vtkIdList::New();
113   newCellPts->Allocate(VTK_CELL_SIZE);
114
115   if ( this->ExtractInside )
116     {
117     multiplier = 1.0;
118     }
119   else 
120     {
121     multiplier = -1.0;
122     }
123
124   // Loop over all points determining whether they are inside the
125   // implicit function. Copy the points and point data if they are.
126   //
127   pointMap = new vtkIdType[numPts]; // maps old point ids into new
128   for (i=0; i < numPts; i++)
129     {
130     pointMap[i] = -1;
131     }
132
133   output->Allocate(numCells/4); //allocate storage for geometry/topology
134   newPts = vtkPoints::New();
135   newPts->Allocate(numPts/4,numPts);
136   outputPD->CopyAllocate(pd);
137   outputCD->CopyAllocate(cd);
138   vtkFloatArray *newScalars = NULL;
139   
140   if(myStoreMapping){
141     myElemVTK2ObjIds.clear();
142     myElemVTK2ObjIds.reserve(numCells);
143     myNodeVTK2ObjIds.clear();
144     myNodeVTK2ObjIds.reserve(numPts);
145   }
146
147   if ( ! this->ExtractBoundaryCells )
148     {
149     for ( ptId=0; ptId < numPts; ptId++ )
150       {
151       x = input->GetPoint(ptId);
152       if ( (this->ImplicitFunction->FunctionValue(x)*multiplier) < 0.0 )
153         {
154         newId = newPts->InsertNextPoint(x);
155         pointMap[ptId] = newId;
156         myNodeVTK2ObjIds.push_back(ptId);
157         outputPD->CopyData(pd,ptId,newId);
158         }
159       }
160     }
161   else
162     {
163     // To extract boundary cells, we have to create supplemental information
164     if ( this->ExtractBoundaryCells )
165       {
166       vtkFloatingPointType val;
167       newScalars = vtkFloatArray::New();
168       newScalars->SetNumberOfValues(numPts);
169
170       for (ptId=0; ptId < numPts; ptId++ )
171         {
172         x = input->GetPoint(ptId);
173         val = this->ImplicitFunction->FunctionValue(x) * multiplier;
174         newScalars->SetValue(ptId, val);
175         if ( val < 0.0 )
176           {
177           newId = newPts->InsertNextPoint(x);
178           pointMap[ptId] = newId;
179           myNodeVTK2ObjIds.push_back(ptId);
180           outputPD->CopyData(pd,ptId,newId);
181           }
182         }
183       }
184     }
185
186   // Now loop over all cells to see whether they are inside implicit
187   // function (or on boundary if ExtractBoundaryCells is on).
188   //
189   for (cellId=0; cellId < numCells; cellId++)
190     {
191     cell = input->GetCell(cellId);
192     cellPts = cell->GetPointIds();
193     numCellPts = cell->GetNumberOfPoints();
194
195     newCellPts->Reset();
196     if ( ! this->ExtractBoundaryCells ) //requires less work
197       {
198       for ( npts=0, i=0; i < numCellPts; i++, npts++)
199         {
200         ptId = cellPts->GetId(i);
201         if ( pointMap[ptId] < 0 )
202           {
203           break; //this cell won't be inserted
204           }
205         else
206           {
207           newCellPts->InsertId(i,pointMap[ptId]);
208           }
209         }
210       } //if don't want to extract boundary cells
211     
212     else //want boundary cells
213       {
214       for ( npts=0, i=0; i < numCellPts; i++ )
215         {
216         ptId = cellPts->GetId(i);
217         if ( newScalars->GetValue(ptId) <= 0.0 )
218           {
219           npts++;
220           }
221         }
222       if ( npts > 0 )
223         {
224         for ( i=0; i < numCellPts; i++ )
225           {
226           ptId = cellPts->GetId(i);
227           if ( pointMap[ptId] < 0 )
228             {
229             x = input->GetPoint(ptId);
230             newId = newPts->InsertNextPoint(x);
231             pointMap[ptId] = newId;
232             myNodeVTK2ObjIds.push_back(ptId);
233             outputPD->CopyData(pd,ptId,newId);
234             }
235           newCellPts->InsertId(i,pointMap[ptId]);
236           }
237         }//a boundary or interior cell
238       }//if mapping boundary cells
239       
240     if ( npts >= numCellPts || (this->ExtractBoundaryCells && npts > 0) )
241       {
242         if(cell->GetCellType() == VTK_POLYHEDRON) {
243           newCellPts->Reset();
244           vtkUnstructuredGrid::SafeDownCast(input)->GetFaceStream( cellId ,newCellPts );        
245           vtkUnstructuredGrid::ConvertFaceStreamPointIds(newCellPts, pointMap);
246         }
247           newCellId = output->InsertNextCell(cell->GetCellType(),newCellPts);
248           myElemVTK2ObjIds.push_back(cellId);
249           outputCD->CopyData(cd,cellId,newCellId);
250       }
251     }//for all cells
252
253   // Update ourselves and release memory
254   //
255   delete [] pointMap;
256   newCellPts->Delete();
257   output->SetPoints(newPts);
258   newPts->Delete();
259   
260   if ( this->ExtractBoundaryCells )
261     {
262     newScalars->Delete();
263     }
264
265   output->Squeeze();
266   return 1;
267 }