src/PIPELINE/SALOME_ExtractGeometry.cxx

   1 // Copyright (C) 2007-2011  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 "SALOME_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
  34 #include <vtkInformation.h>
  35 #include <vtkInformationVector.h>
  36
  37 #include <vtkImplicitBoolean.h>
  38 #include <vtkImplicitFunctionCollection.h>
  39
  40 #if defined __GNUC__
  41   #if __GNUC__ == 2
  42     #define __GNUC_2__
  43   #endif
  44 #endif
  45
  46
  47 //----------------------------------------------------------------------------
  48 vtkStandardNewMacro(SALOME_ExtractGeometry);
  49
  50
  51 //----------------------------------------------------------------------------
  52 SALOME_ExtractGeometry
  53 ::SALOME_ExtractGeometry():
  54   myStoreMapping(false),
  55   myIsDoneShallowCopy(false)
  56 {}
  57
  58 SALOME_ExtractGeometry
  59 ::~SALOME_ExtractGeometry()
  60 {}
  61
  62
  63 //----------------------------------------------------------------------------
  64 vtkImplicitBoolean*
  65 SALOME_ExtractGeometry
  66 ::GetImplicitBoolean()
  67 {
  68   return myImplicitBoolean.GetPointer();
  69 }
  70
  71
  72 void
  73 SALOME_ExtractGeometry
  74 ::SetImplicitFunction(vtkImplicitFunction* theImplicitFunction)
  75 {
  76   myImplicitBoolean = dynamic_cast<vtkImplicitBoolean*>(theImplicitFunction);
  77   Superclass::SetImplicitFunction(theImplicitFunction);
  78 }
  79
  80
  81 //----------------------------------------------------------------------------
  82 void
  83 SALOME_ExtractGeometry
  84 ::SetStoreMapping(bool theStoreMapping)
  85 {
  86   if(myStoreMapping != theStoreMapping){
  87     myStoreMapping = theStoreMapping;
  88     Modified();
  89   }
  90 }
  91
  92 bool
  93 SALOME_ExtractGeometry
  94 ::GetStoreMapping() const
  95 {
  96   return myStoreMapping;
  97 }
  98
  99
 100 //----------------------------------------------------------------------------
 101 vtkIdType
 102 SALOME_ExtractGeometry
 103 ::GetElemVTKId(vtkIdType theID)
 104 {
 105   if(!myStoreMapping || myIsDoneShallowCopy)
 106     return theID;
 107
 108   vtkIdType iEnd = myElemVTK2ObjIds.size();
 109   for(vtkIdType i = 0; i < iEnd; i++)
 110     if(myElemVTK2ObjIds[i] == theID)
 111       return i;
 112
 113   return -1;
 114 }
 115
 116 vtkIdType
 117 SALOME_ExtractGeometry
 118 ::GetNodeVTKId(vtkIdType theID)
 119 {
 120   if(!myStoreMapping || myIsDoneShallowCopy)
 121     return theID;
 122
 123   vtkIdType iEnd = myNodeVTK2ObjIds.size();
 124   for(vtkIdType i = 0; i < iEnd; i++)
 125     if(myNodeVTK2ObjIds[i] == theID)
 126       return i;
 127
 128   return -1;
 129 }
 130
 131
 132 //----------------------------------------------------------------------------
 133 vtkIdType
 134 SALOME_ExtractGeometry
 135 ::GetElemObjId(vtkIdType theVtkID)
 136 {
 137   if(!myStoreMapping || myIsDoneShallowCopy)
 138     return theVtkID;
 139
 140   if(theVtkID < myElemVTK2ObjIds.size())
 141     return myElemVTK2ObjIds[theVtkID];
 142
 143   return -1;
 144 }
 145
 146
 147 vtkIdType
 148 SALOME_ExtractGeometry
 149 ::GetNodeObjId(vtkIdType theVtkID)
 150 {
 151   if(!myStoreMapping || myIsDoneShallowCopy)
 152     return theVtkID;
 153
 154   if(theVtkID < myNodeVTK2ObjIds.size())
 155     return myNodeVTK2ObjIds[theVtkID];
 156
 157   return -1;
 158 }
 159
 160
 161 //----------------------------------------------------------------------------
 162 int
 163 SALOME_ExtractGeometry
 164 ::RequestData(vtkInformation *request,
 165               vtkInformationVector **inputVector,
 166               vtkInformationVector *outputVector)
 167 {
 168   // get the info objects
 169   vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
 170   vtkInformation *outInfo = outputVector->GetInformationObject(0);
 171
 172   // get the input and ouptut
 173   vtkDataSet *input = vtkDataSet::SafeDownCast(
 174     inInfo->Get(vtkDataObject::DATA_OBJECT()));
 175   vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
 176     outInfo->Get(vtkDataObject::DATA_OBJECT()));
 177
 178   myElemVTK2ObjIds.clear();
 179   myNodeVTK2ObjIds.clear();
 180   //
 181   myIsDoneShallowCopy = !this->ImplicitFunction;
 182
 183   if(!myIsDoneShallowCopy && myImplicitBoolean.GetPointer())
 184     if(vtkImplicitFunctionCollection* aFunction = myImplicitBoolean->GetFunction())
 185       myIsDoneShallowCopy = aFunction->GetNumberOfItems() == 0;
 186
 187   if(myIsDoneShallowCopy){
 188     output->ShallowCopy(input);
 189     return 1;
 190   }
 191
 192   return RequestData2(request,inputVector,outputVector);
 193 }
 194
 195
 196 //----------------------------------------------------------------------------
 197 int
 198 SALOME_ExtractGeometry
 199 ::RequestData2(vtkInformation *vtkNotUsed(request),
 200                vtkInformationVector **inputVector,
 201                vtkInformationVector *outputVector)
 202 {
 203   // get the info objects
 204   vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
 205   vtkInformation *outInfo = outputVector->GetInformationObject(0);
 206
 207   // get the input and ouptut
 208   vtkDataSet *input = vtkDataSet::SafeDownCast(
 209     inInfo->Get(vtkDataObject::DATA_OBJECT()));
 210   vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
 211     outInfo->Get(vtkDataObject::DATA_OBJECT()));
 212
 213   vtkIdType ptId, numPts, numCells, i, cellId, newCellId, newId, *pointMap;
 214   vtkIdList *cellPts;
 215   vtkCell *cell;
 216   int numCellPts;
 217   vtkFloatingPointType *x;
 218   vtkFloatingPointType multiplier;
 219   vtkPoints *newPts;
 220   vtkIdList *newCellPts;
 221   vtkPointData *pd = input->GetPointData();
 222   vtkCellData *cd = input->GetCellData();
 223   vtkPointData *outputPD = output->GetPointData();
 224   vtkCellData *outputCD = output->GetCellData();
 225   int npts;
 226   numCells = input->GetNumberOfCells();
 227   numPts = input->GetNumberOfPoints();
 228
 229   if ( ! this->ImplicitFunction )
 230     {
 231     vtkErrorMacro(<<"No implicit function specified");
 232     return 0;
 233     }
 234
 235   newCellPts = vtkIdList::New();
 236   newCellPts->Allocate(VTK_CELL_SIZE);
 237
 238   if ( this->ExtractInside )
 239     {
 240     multiplier = 1.0;
 241     }
 242   else
 243     {
 244     multiplier = -1.0;
 245     }
 246
 247   // Loop over all points determining whether they are inside the
 248   // implicit function. Copy the points and point data if they are.
 249   //
 250   pointMap = new vtkIdType[numPts]; // maps old point ids into new
 251   for (i=0; i < numPts; i++)
 252     {
 253     pointMap[i] = -1;
 254     }
 255
 256   output->Allocate(numCells/4); //allocate storage for geometry/topology
 257   newPts = vtkPoints::New();
 258   newPts->Allocate(numPts/4,numPts);
 259   outputPD->CopyAllocate(pd);
 260   outputCD->CopyAllocate(cd);
 261   vtkFloatArray *newScalars = NULL;
 262
 263   if(myStoreMapping){
 264     myElemVTK2ObjIds.reserve(numCells);
 265     myNodeVTK2ObjIds.reserve(numPts);
 266   }
 267
 268   if ( ! this->ExtractBoundaryCells )
 269     {
 270     for ( ptId=0; ptId < numPts; ptId++ )
 271       {
 272       x = input->GetPoint(ptId);
 273       if ( (this->ImplicitFunction->FunctionValue(x)*multiplier) < 0.0 )
 274         {
 275         newId = newPts->InsertNextPoint(x);
 276         pointMap[ptId] = newId;
 277         if(myStoreMapping)
 278           myNodeVTK2ObjIds.push_back(ptId);
 279         outputPD->CopyData(pd,ptId,newId);
 280         }
 281       }
 282     }
 283   else
 284     {
 285     // To extract boundary cells, we have to create supplemental information
 286     if ( this->ExtractBoundaryCells )
 287       {
 288       vtkFloatingPointType val;
 289       newScalars = vtkFloatArray::New();
 290       newScalars->SetNumberOfValues(numPts);
 291
 292       for (ptId=0; ptId < numPts; ptId++ )
 293         {
 294         x = input->GetPoint(ptId);
 295         val = this->ImplicitFunction->FunctionValue(x) * multiplier;
 296         newScalars->SetValue(ptId, val);
 297         if ( val < 0.0 )
 298           {
 299           newId = newPts->InsertNextPoint(x);
 300           pointMap[ptId] = newId;
 301           if(myStoreMapping)
 302             myNodeVTK2ObjIds.push_back(ptId);
 303           outputPD->CopyData(pd,ptId,newId);
 304           }
 305         }
 306       }
 307     }
 308
 309   // Now loop over all cells to see whether they are inside implicit
 310   // function (or on boundary if ExtractBoundaryCells is on).
 311   //
 312   for (cellId=0; cellId < numCells; cellId++)
 313     {
 314     cell = input->GetCell(cellId);
 315     cellPts = cell->GetPointIds();
 316     numCellPts = cell->GetNumberOfPoints();
 317
 318     newCellPts->Reset();
 319     if ( ! this->ExtractBoundaryCells ) //requires less work
 320       {
 321       for ( npts=0, i=0; i < numCellPts; i++, npts++)
 322         {
 323         ptId = cellPts->GetId(i);
 324         if ( pointMap[ptId] < 0 )
 325           {
 326           break; //this cell won't be inserted
 327           }
 328         else
 329           {
 330           newCellPts->InsertId(i,pointMap[ptId]);
 331           }
 332         }
 333       } //if don't want to extract boundary cells
 334
 335     else //want boundary cells
 336       {
 337       for ( npts=0, i=0; i < numCellPts; i++ )
 338         {
 339         ptId = cellPts->GetId(i);
 340         if ( newScalars->GetValue(ptId) <= 0.0 )
 341           {
 342           npts++;
 343           }
 344         }
 345       if ( npts > 0 )
 346         {
 347         for ( i=0; i < numCellPts; i++ )
 348           {
 349           ptId = cellPts->GetId(i);
 350           if ( pointMap[ptId] < 0 )
 351             {
 352             x = input->GetPoint(ptId);
 353             newId = newPts->InsertNextPoint(x);
 354             pointMap[ptId] = newId;
 355             if(myStoreMapping)
 356               myNodeVTK2ObjIds.push_back(ptId);
 357             outputPD->CopyData(pd,ptId,newId);
 358             }
 359           newCellPts->InsertId(i,pointMap[ptId]);
 360           }
 361         }//a boundary or interior cell
 362       }//if mapping boundary cells
 363
 364     if ( npts >= numCellPts || (this->ExtractBoundaryCells && npts > 0) )
 365       {
 366       newCellId = output->InsertNextCell(cell->GetCellType(),newCellPts);
 367       if(myStoreMapping)
 368         myElemVTK2ObjIds.push_back(cellId);
 369       outputCD->CopyData(cd,cellId,newCellId);
 370       }
 371     }//for all cells
 372
 373   // Update ourselves and release memory
 374   //
 375   delete [] pointMap;
 376   newCellPts->Delete();
 377   output->SetPoints(newPts);
 378   newPts->Delete();
 379
 380   if ( this->ExtractBoundaryCells )
 381     {
 382     newScalars->Delete();
 383     }
 384
 385   output->Squeeze();
 386
 387   return 1;
 388 }