src/PIPELINE/SALOME_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 "SALOME_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 #include <vtkImplicitBoolean.h>
  33 #include <vtkImplicitFunctionCollection.h>
  34
  35 using namespace std;
  36
  37 #if defined __GNUC__
  38   #if __GNUC__ == 2
  39     #define __GNUC_2__
  40   #endif
  41 #endif
  42
  43
  44 vtkStandardNewMacro(SALOME_ExtractGeometry);
  45
  46
  47 SALOME_ExtractGeometry::SALOME_ExtractGeometry()
  48 {}
  49
  50 SALOME_ExtractGeometry::~SALOME_ExtractGeometry()
  51 {}
  52
  53 vtkIdType SALOME_ExtractGeometry::GetElemObjId(int theVtkID)
  54 {
  55   if (myElemVTK2ObjIds.empty())
  56     return theVtkID;
  57
  58   if (theVtkID < 0 || myElemVTK2ObjIds.size() <= theVtkID)
  59     return -1;
  60
  61 #if defined __GNUC_2__
  62   return myElemVTK2ObjIds[theVtkID];
  63 #else
  64   return myElemVTK2ObjIds.at(theVtkID);
  65 #endif
  66 }
  67
  68
  69 vtkIdType SALOME_ExtractGeometry::GetNodeObjId(int theVtkID)
  70 {
  71   if (myNodeVTK2ObjIds.empty())
  72     return theVtkID;
  73
  74   if (theVtkID < 0 || myNodeVTK2ObjIds.size() <= theVtkID)
  75     return -1;
  76
  77 #if defined __GNUC_2__
  78   return myNodeVTK2ObjIds[theVtkID];
  79 #else
  80   return myNodeVTK2ObjIds.at(theVtkID);
  81 #endif
  82 }
  83
  84
  85 void SALOME_ExtractGeometry::SetImplicitBoolean(vtkImplicitBoolean* theImplicitBoolean)
  86 {
  87   myImplicitBoolean = theImplicitBoolean;
  88   SetImplicitFunction(theImplicitBoolean);
  89 }
  90
  91
  92 void SALOME_ExtractGeometry::SetStoreMapping(bool theStoreMapping)
  93 {
  94   myStoreMapping = theStoreMapping;
  95   Modified();
  96 }
  97
  98
  99 void SALOME_ExtractGeometry::Execute()
 100 {
 101   if(myImplicitBoolean.GetPointer()){
 102     if(vtkImplicitFunctionCollection* aFunction = myImplicitBoolean->GetFunction()){
 103       if(aFunction->GetNumberOfItems() == 0){
 104         myElemVTK2ObjIds.clear();
 105         myNodeVTK2ObjIds.clear();
 106
 107         vtkDebugMacro(<< "Extracting geometry - ShallowCopy");
 108         GetOutput()->ShallowCopy(GetInput());
 109         Modified();
 110         return;
 111       }
 112     }
 113   }
 114   Execute2();
 115 }
 116
 117 void SALOME_ExtractGeometry::Execute2()
 118 {
 119   vtkIdType ptId, numPts, numCells, i, cellId, newCellId, newId, *pointMap;
 120   vtkIdList *cellPts;
 121   vtkCell *cell;
 122   int numCellPts;
 123   float *x;
 124   float multiplier;
 125   vtkPoints *newPts;
 126   vtkIdList *newCellPts;
 127   vtkDataSet *input = this->GetInput();
 128   vtkPointData *pd = input->GetPointData();
 129   vtkCellData *cd = input->GetCellData();
 130   vtkUnstructuredGrid *output = this->GetOutput();
 131   vtkPointData *outputPD = output->GetPointData();
 132   vtkCellData *outputCD = output->GetCellData();
 133   int npts;
 134   numCells = input->GetNumberOfCells();
 135   numPts = input->GetNumberOfPoints();
 136
 137   vtkDebugMacro(<< "Extracting geometry");
 138
 139   if ( ! this->ImplicitFunction )
 140     {
 141     vtkErrorMacro(<<"No implicit function specified");
 142     return;
 143     }
 144
 145   newCellPts = vtkIdList::New();
 146   newCellPts->Allocate(VTK_CELL_SIZE);
 147
 148   if ( this->ExtractInside )
 149     {
 150     multiplier = 1.0;
 151     }
 152   else
 153     {
 154     multiplier = -1.0;
 155     }
 156
 157   // Loop over all points determining whether they are inside the
 158   // implicit function. Copy the points and point data if they are.
 159   //
 160   pointMap = new vtkIdType[numPts]; // maps old point ids into new
 161   for (i=0; i < numPts; i++)
 162     {
 163     pointMap[i] = -1;
 164     }
 165
 166   output->Allocate(numCells/4); //allocate storage for geometry/topology
 167   newPts = vtkPoints::New();
 168   newPts->Allocate(numPts/4,numPts);
 169   outputPD->CopyAllocate(pd);
 170   outputCD->CopyAllocate(cd);
 171   vtkFloatArray *newScalars = NULL;
 172
 173   if(myStoreMapping){
 174     myElemVTK2ObjIds.clear();
 175     myElemVTK2ObjIds.reserve(numCells);
 176     myNodeVTK2ObjIds.clear();
 177     myNodeVTK2ObjIds.reserve(numPts);
 178   }
 179
 180   if ( ! this->ExtractBoundaryCells )
 181     {
 182     for ( ptId=0; ptId < numPts; ptId++ )
 183       {
 184       x = input->GetPoint(ptId);
 185       if ( (this->ImplicitFunction->FunctionValue(x)*multiplier) < 0.0 )
 186         {
 187         newId = newPts->InsertNextPoint(x);
 188         pointMap[ptId] = newId;
 189         myNodeVTK2ObjIds.push_back(ptId);
 190         outputPD->CopyData(pd,ptId,newId);
 191         }
 192       }
 193     }
 194   else
 195     {
 196     // To extract boundary cells, we have to create supplemental information
 197     if ( this->ExtractBoundaryCells )
 198       {
 199       float val;
 200       newScalars = vtkFloatArray::New();
 201       newScalars->SetNumberOfValues(numPts);
 202
 203       for (ptId=0; ptId < numPts; ptId++ )
 204         {
 205         x = input->GetPoint(ptId);
 206         val = this->ImplicitFunction->FunctionValue(x) * multiplier;
 207         newScalars->SetValue(ptId, val);
 208         if ( val < 0.0 )
 209           {
 210           newId = newPts->InsertNextPoint(x);
 211           pointMap[ptId] = newId;
 212           myNodeVTK2ObjIds.push_back(ptId);
 213           outputPD->CopyData(pd,ptId,newId);
 214           }
 215         }
 216       }
 217     }
 218
 219   // Now loop over all cells to see whether they are inside implicit
 220   // function (or on boundary if ExtractBoundaryCells is on).
 221   //
 222   for (cellId=0; cellId < numCells; cellId++)
 223     {
 224     cell = input->GetCell(cellId);
 225     cellPts = cell->GetPointIds();
 226     numCellPts = cell->GetNumberOfPoints();
 227
 228     newCellPts->Reset();
 229     if ( ! this->ExtractBoundaryCells ) //requires less work
 230       {
 231       for ( npts=0, i=0; i < numCellPts; i++, npts++)
 232         {
 233         ptId = cellPts->GetId(i);
 234         if ( pointMap[ptId] < 0 )
 235           {
 236           break; //this cell won't be inserted
 237           }
 238         else
 239           {
 240           newCellPts->InsertId(i,pointMap[ptId]);
 241           }
 242         }
 243       } //if don't want to extract boundary cells
 244
 245     else //want boundary cells
 246       {
 247       for ( npts=0, i=0; i < numCellPts; i++ )
 248         {
 249         ptId = cellPts->GetId(i);
 250         if ( newScalars->GetValue(ptId) <= 0.0 )
 251           {
 252           npts++;
 253           }
 254         }
 255       if ( npts > 0 )
 256         {
 257         for ( i=0; i < numCellPts; i++ )
 258           {
 259           ptId = cellPts->GetId(i);
 260           if ( pointMap[ptId] < 0 )
 261             {
 262             x = input->GetPoint(ptId);
 263             newId = newPts->InsertNextPoint(x);
 264             pointMap[ptId] = newId;
 265             myNodeVTK2ObjIds.push_back(ptId);
 266             outputPD->CopyData(pd,ptId,newId);
 267             }
 268           newCellPts->InsertId(i,pointMap[ptId]);
 269           }
 270         }//a boundary or interior cell
 271       }//if mapping boundary cells
 272
 273     if ( npts >= numCellPts || (this->ExtractBoundaryCells && npts > 0) )
 274       {
 275       newCellId = output->InsertNextCell(cell->GetCellType(),newCellPts);
 276       myElemVTK2ObjIds.push_back(cellId);
 277       outputCD->CopyData(cd,cellId,newCellId);
 278       }
 279     }//for all cells
 280
 281   // Update ourselves and release memory
 282   //
 283   delete [] pointMap;
 284   newCellPts->Delete();
 285   output->SetPoints(newPts);
 286   newPts->Delete();
 287
 288   if ( this->ExtractBoundaryCells )
 289     {
 290     newScalars->Delete();
 291     }
 292
 293   output->Squeeze();
 294 }