Back to index

salome-smesh  6.5.0
SMESH_ExtractGeometry.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 #include "SMESH_ExtractGeometry.h"
00024 
00025 #include <vtkCell.h>
00026 #include <vtkCellData.h>
00027 #include <vtkFloatArray.h>
00028 #include <vtkIdList.h>
00029 #include <vtkImplicitFunction.h>
00030 #include <vtkObjectFactory.h>
00031 #include <vtkPointData.h>
00032 #include <vtkUnstructuredGrid.h>
00033 #include <vtkInformation.h>
00034 #include <vtkInformationVector.h>
00035 
00036 using namespace std;
00037 
00038 //#ifdef _DEBUG_
00039 //static int MYDEBUG = 0;
00040 //#else
00041 //static int MYDEBUG = 0;
00042 //#endif
00043 
00044 #if defined __GNUC__
00045   #if __GNUC__ == 2
00046     #define __GNUC_2__
00047   #endif
00048 #endif
00049 
00050 
00051 vtkStandardNewMacro(SMESH_ExtractGeometry);
00052 
00053 
00054 SMESH_ExtractGeometry::SMESH_ExtractGeometry()
00055 {}
00056 
00057 
00058 SMESH_ExtractGeometry::~SMESH_ExtractGeometry(){}
00059 
00060 
00061 vtkIdType SMESH_ExtractGeometry::GetElemObjId(int theVtkID){
00062   if( theVtkID < 0 || theVtkID >= myElemVTK2ObjIds.size()) return -1;
00063   return myElemVTK2ObjIds[theVtkID];
00064 }
00065 
00066 
00067 vtkIdType SMESH_ExtractGeometry::GetNodeObjId(int theVtkID){
00068   if ( theVtkID < 0 || theVtkID >= myNodeVTK2ObjIds.size()) return -1;
00069   return myNodeVTK2ObjIds[theVtkID];
00070 }
00071 
00072 
00073 int SMESH_ExtractGeometry::RequestData(
00074   vtkInformation *vtkNotUsed(request),
00075   vtkInformationVector **inputVector,
00076   vtkInformationVector *outputVector)
00077 {
00078   // get the info objects
00079   vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
00080   vtkInformation *outInfo = outputVector->GetInformationObject(0);
00081 
00082   // get the input and ouptut
00083   vtkDataSet *input = vtkDataSet::SafeDownCast(
00084     inInfo->Get(vtkDataObject::DATA_OBJECT()));
00085   vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
00086     outInfo->Get(vtkDataObject::DATA_OBJECT()));
00087 
00088   vtkIdType ptId, numPts, numCells, i, cellId, newCellId, newId, *pointMap;
00089   vtkIdList *cellPts;
00090   vtkCell *cell;
00091   int numCellPts;
00092   vtkFloatingPointType *x;
00093   vtkFloatingPointType multiplier;
00094   vtkPoints *newPts;
00095   vtkIdList *newCellPts;
00096   vtkPointData *pd = input->GetPointData();
00097   vtkCellData *cd = input->GetCellData();
00098   vtkPointData *outputPD = output->GetPointData();
00099   vtkCellData *outputCD = output->GetCellData();
00100   int npts;
00101   numCells = input->GetNumberOfCells();
00102   numPts = input->GetNumberOfPoints();
00103   
00104   vtkDebugMacro(<< "Extracting geometry");
00105 
00106   if ( ! this->ImplicitFunction )
00107     {
00108     vtkErrorMacro(<<"No implicit function specified");
00109     return 0;
00110     }
00111 
00112   newCellPts = vtkIdList::New();
00113   newCellPts->Allocate(VTK_CELL_SIZE);
00114 
00115   if ( this->ExtractInside )
00116     {
00117     multiplier = 1.0;
00118     }
00119   else 
00120     {
00121     multiplier = -1.0;
00122     }
00123 
00124   // Loop over all points determining whether they are inside the
00125   // implicit function. Copy the points and point data if they are.
00126   //
00127   pointMap = new vtkIdType[numPts]; // maps old point ids into new
00128   for (i=0; i < numPts; i++)
00129     {
00130     pointMap[i] = -1;
00131     }
00132 
00133   output->Allocate(numCells/4); //allocate storage for geometry/topology
00134   newPts = vtkPoints::New();
00135   newPts->Allocate(numPts/4,numPts);
00136   outputPD->CopyAllocate(pd);
00137   outputCD->CopyAllocate(cd);
00138   vtkFloatArray *newScalars = NULL;
00139   
00140   if(myStoreMapping){
00141     myElemVTK2ObjIds.clear();
00142     myElemVTK2ObjIds.reserve(numCells);
00143     myNodeVTK2ObjIds.clear();
00144     myNodeVTK2ObjIds.reserve(numPts);
00145   }
00146 
00147   if ( ! this->ExtractBoundaryCells )
00148     {
00149     for ( ptId=0; ptId < numPts; ptId++ )
00150       {
00151       x = input->GetPoint(ptId);
00152       if ( (this->ImplicitFunction->FunctionValue(x)*multiplier) < 0.0 )
00153         {
00154         newId = newPts->InsertNextPoint(x);
00155         pointMap[ptId] = newId;
00156         myNodeVTK2ObjIds.push_back(ptId);
00157         outputPD->CopyData(pd,ptId,newId);
00158         }
00159       }
00160     }
00161   else
00162     {
00163     // To extract boundary cells, we have to create supplemental information
00164     if ( this->ExtractBoundaryCells )
00165       {
00166       vtkFloatingPointType val;
00167       newScalars = vtkFloatArray::New();
00168       newScalars->SetNumberOfValues(numPts);
00169 
00170       for (ptId=0; ptId < numPts; ptId++ )
00171         {
00172         x = input->GetPoint(ptId);
00173         val = this->ImplicitFunction->FunctionValue(x) * multiplier;
00174         newScalars->SetValue(ptId, val);
00175         if ( val < 0.0 )
00176           {
00177           newId = newPts->InsertNextPoint(x);
00178           pointMap[ptId] = newId;
00179           myNodeVTK2ObjIds.push_back(ptId);
00180           outputPD->CopyData(pd,ptId,newId);
00181           }
00182         }
00183       }
00184     }
00185 
00186   // Now loop over all cells to see whether they are inside implicit
00187   // function (or on boundary if ExtractBoundaryCells is on).
00188   //
00189   for (cellId=0; cellId < numCells; cellId++)
00190     {
00191     cell = input->GetCell(cellId);
00192     cellPts = cell->GetPointIds();
00193     numCellPts = cell->GetNumberOfPoints();
00194 
00195     newCellPts->Reset();
00196     if ( ! this->ExtractBoundaryCells ) //requires less work
00197       {
00198       for ( npts=0, i=0; i < numCellPts; i++, npts++)
00199         {
00200         ptId = cellPts->GetId(i);
00201         if ( pointMap[ptId] < 0 )
00202           {
00203           break; //this cell won't be inserted
00204           }
00205         else
00206           {
00207           newCellPts->InsertId(i,pointMap[ptId]);
00208           }
00209         }
00210       } //if don't want to extract boundary cells
00211     
00212     else //want boundary cells
00213       {
00214       for ( npts=0, i=0; i < numCellPts; i++ )
00215         {
00216         ptId = cellPts->GetId(i);
00217         if ( newScalars->GetValue(ptId) <= 0.0 )
00218           {
00219           npts++;
00220           }
00221         }
00222       if ( npts > 0 )
00223         {
00224         for ( i=0; i < numCellPts; i++ )
00225           {
00226           ptId = cellPts->GetId(i);
00227           if ( pointMap[ptId] < 0 )
00228             {
00229             x = input->GetPoint(ptId);
00230             newId = newPts->InsertNextPoint(x);
00231             pointMap[ptId] = newId;
00232             myNodeVTK2ObjIds.push_back(ptId);
00233             outputPD->CopyData(pd,ptId,newId);
00234             }
00235           newCellPts->InsertId(i,pointMap[ptId]);
00236           }
00237         }//a boundary or interior cell
00238       }//if mapping boundary cells
00239       
00240     if ( npts >= numCellPts || (this->ExtractBoundaryCells && npts > 0) )
00241       {
00242         if(cell->GetCellType() == VTK_POLYHEDRON) {
00243           newCellPts->Reset();
00244           vtkUnstructuredGrid::SafeDownCast(input)->GetFaceStream( cellId ,newCellPts );        
00245           vtkUnstructuredGrid::ConvertFaceStreamPointIds(newCellPts, pointMap);
00246         }
00247           newCellId = output->InsertNextCell(cell->GetCellType(),newCellPts);
00248           myElemVTK2ObjIds.push_back(cellId);
00249           outputCD->CopyData(cd,cellId,newCellId);
00250       }
00251     }//for all cells
00252 
00253   // Update ourselves and release memory
00254   //
00255   delete [] pointMap;
00256   newCellPts->Delete();
00257   output->SetPoints(newPts);
00258   newPts->Delete();
00259   
00260   if ( this->ExtractBoundaryCells )
00261     {
00262     newScalars->Delete();
00263     }
00264 
00265   output->Squeeze();
00266   return 1;
00267 }