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