1 // Copyright (C) 2007-2008 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
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.
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.
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
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 #include "SMESH_ExtractGeometry.h"
25 #include <vtkCellData.h>
26 #include <vtkFloatArray.h>
27 #include <vtkIdList.h>
28 #include <vtkImplicitFunction.h>
29 #include <vtkObjectFactory.h>
30 #include <vtkPointData.h>
31 #include <vtkUnstructuredGrid.h>
32 #include <vtkInformation.h>
33 #include <vtkInformationVector.h>
38 static int MYDEBUG = 0;
40 static int MYDEBUG = 0;
50 vtkStandardNewMacro(SMESH_ExtractGeometry);
53 SMESH_ExtractGeometry::SMESH_ExtractGeometry()
57 SMESH_ExtractGeometry::~SMESH_ExtractGeometry(){}
60 vtkIdType SMESH_ExtractGeometry::GetElemObjId(int theVtkID){
61 if(myElemVTK2ObjIds.empty() || theVtkID > myElemVTK2ObjIds.size()) return -1;
62 #if defined __GNUC_2__
63 return myElemVTK2ObjIds[theVtkID];
65 return myElemVTK2ObjIds.at(theVtkID);
70 vtkIdType SMESH_ExtractGeometry::GetNodeObjId(int theVtkID){
71 if(myNodeVTK2ObjIds.empty() || theVtkID > myNodeVTK2ObjIds.size()) return -1;
72 #if defined __GNUC_2__
73 return myNodeVTK2ObjIds[theVtkID];
75 return myNodeVTK2ObjIds.at(theVtkID);
80 int SMESH_ExtractGeometry::RequestData(
81 vtkInformation *vtkNotUsed(request),
82 vtkInformationVector **inputVector,
83 vtkInformationVector *outputVector)
85 // get the info objects
86 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
87 vtkInformation *outInfo = outputVector->GetInformationObject(0);
89 // get the input and ouptut
90 vtkDataSet *input = vtkDataSet::SafeDownCast(
91 inInfo->Get(vtkDataObject::DATA_OBJECT()));
92 vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
93 outInfo->Get(vtkDataObject::DATA_OBJECT()));
95 vtkIdType ptId, numPts, numCells, i, cellId, newCellId, newId, *pointMap;
99 vtkFloatingPointType *x;
100 vtkFloatingPointType multiplier;
102 vtkIdList *newCellPts;
103 vtkPointData *pd = input->GetPointData();
104 vtkCellData *cd = input->GetCellData();
105 vtkPointData *outputPD = output->GetPointData();
106 vtkCellData *outputCD = output->GetCellData();
108 numCells = input->GetNumberOfCells();
109 numPts = input->GetNumberOfPoints();
111 vtkDebugMacro(<< "Extracting geometry");
113 if ( ! this->ImplicitFunction )
115 vtkErrorMacro(<<"No implicit function specified");
119 newCellPts = vtkIdList::New();
120 newCellPts->Allocate(VTK_CELL_SIZE);
122 if ( this->ExtractInside )
131 // Loop over all points determining whether they are inside the
132 // implicit function. Copy the points and point data if they are.
134 pointMap = new vtkIdType[numPts]; // maps old point ids into new
135 for (i=0; i < numPts; i++)
140 output->Allocate(numCells/4); //allocate storage for geometry/topology
141 newPts = vtkPoints::New();
142 newPts->Allocate(numPts/4,numPts);
143 outputPD->CopyAllocate(pd);
144 outputCD->CopyAllocate(cd);
145 vtkFloatArray *newScalars = NULL;
148 myElemVTK2ObjIds.clear();
149 myElemVTK2ObjIds.reserve(numCells);
150 myNodeVTK2ObjIds.clear();
151 myNodeVTK2ObjIds.reserve(numPts);
154 if ( ! this->ExtractBoundaryCells )
156 for ( ptId=0; ptId < numPts; ptId++ )
158 x = input->GetPoint(ptId);
159 if ( (this->ImplicitFunction->FunctionValue(x)*multiplier) < 0.0 )
161 newId = newPts->InsertNextPoint(x);
162 pointMap[ptId] = newId;
163 myNodeVTK2ObjIds.push_back(ptId);
164 outputPD->CopyData(pd,ptId,newId);
170 // To extract boundary cells, we have to create supplemental information
171 if ( this->ExtractBoundaryCells )
173 vtkFloatingPointType val;
174 newScalars = vtkFloatArray::New();
175 newScalars->SetNumberOfValues(numPts);
177 for (ptId=0; ptId < numPts; ptId++ )
179 x = input->GetPoint(ptId);
180 val = this->ImplicitFunction->FunctionValue(x) * multiplier;
181 newScalars->SetValue(ptId, val);
184 newId = newPts->InsertNextPoint(x);
185 pointMap[ptId] = newId;
186 myNodeVTK2ObjIds.push_back(ptId);
187 outputPD->CopyData(pd,ptId,newId);
193 // Now loop over all cells to see whether they are inside implicit
194 // function (or on boundary if ExtractBoundaryCells is on).
196 for (cellId=0; cellId < numCells; cellId++)
198 cell = input->GetCell(cellId);
199 cellPts = cell->GetPointIds();
200 numCellPts = cell->GetNumberOfPoints();
203 if ( ! this->ExtractBoundaryCells ) //requires less work
205 for ( npts=0, i=0; i < numCellPts; i++, npts++)
207 ptId = cellPts->GetId(i);
208 if ( pointMap[ptId] < 0 )
210 break; //this cell won't be inserted
214 newCellPts->InsertId(i,pointMap[ptId]);
217 } //if don't want to extract boundary cells
219 else //want boundary cells
221 for ( npts=0, i=0; i < numCellPts; i++ )
223 ptId = cellPts->GetId(i);
224 if ( newScalars->GetValue(ptId) <= 0.0 )
231 for ( i=0; i < numCellPts; i++ )
233 ptId = cellPts->GetId(i);
234 if ( pointMap[ptId] < 0 )
236 x = input->GetPoint(ptId);
237 newId = newPts->InsertNextPoint(x);
238 pointMap[ptId] = newId;
239 myNodeVTK2ObjIds.push_back(ptId);
240 outputPD->CopyData(pd,ptId,newId);
242 newCellPts->InsertId(i,pointMap[ptId]);
244 }//a boundary or interior cell
245 }//if mapping boundary cells
247 if ( npts >= numCellPts || (this->ExtractBoundaryCells && npts > 0) )
249 newCellId = output->InsertNextCell(cell->GetCellType(),newCellPts);
250 myElemVTK2ObjIds.push_back(cellId);
251 outputCD->CopyData(cd,cellId,newCellId);
255 // Update ourselves and release memory
258 newCellPts->Delete();
259 output->SetPoints(newPts);
262 if ( this->ExtractBoundaryCells )
264 newScalars->Delete();