moab
MixedCellConnectivity.h
Go to the documentation of this file.
00001 #ifndef __smoab_MixedCellConnectivity_h
00002 #define __smoab_MixedCellConnectivity_h
00003 
00004 #include "vtkCellType.h"
00005 #include <algorithm>
00006 
00007 namespace
00008 {
00009 
00010 template<int N> struct QuadratricOrdering{};
00011 
00012 template<> struct QuadratricOrdering<VTK_QUADRATIC_WEDGE>
00013 {
00014   static const int NUM_VERTS = 15;
00015   void reorder(vtkIdType* connectivity) const
00016   {
00017     std::swap_ranges(connectivity+9,connectivity+12,connectivity+12);
00018   }
00019 };
00020 
00021 template<> struct QuadratricOrdering<VTK_TRIQUADRATIC_HEXAHEDRON>
00022 {
00023   static const int NUM_VERTS = 27;
00024   void reorder(vtkIdType* connectivity) const
00025   {
00026     std::swap_ranges(connectivity+12,connectivity+16,connectivity+16);
00027 
00028     //move 20 to 22
00029     //move 22 to 23
00030     //move 23 to 20
00031 
00032     //swap 20 with 22
00033     std::swap(connectivity[20],connectivity[23]);
00034 
00035     //swap 22 with 23
00036     std::swap(connectivity[22],connectivity[23]);
00037   }
00038 };
00039 
00040 template<typename QuadraticOrdering>
00041 void FixQuadraticIdOrdering(vtkIdType* connectivity, vtkIdType numCells,
00042                             QuadraticOrdering& ordering)
00043 {
00044   //skip the first index that holds the length of the cells
00045   //if we skip it once here, and than properly increment it makes the code
00046   //far easier
00047   connectivity+=1;
00048   for(vtkIdType i=0; i < numCells; ++i)
00049     {
00050     ordering.reorder(connectivity);
00051     connectivity += ordering.NUM_VERTS + 1;
00052     }
00053 }
00054 }
00055 
00056 namespace smoab
00057 {
00058 
00059 class MixedCellConnectivity
00060 {
00061 public:
00062   MixedCellConnectivity(smoab::Range const& cells, moab::Interface* moab):
00063     Connectivity(),
00064     UniquePoints(),
00065     Info()
00066     {
00067     int count = 0;
00068     const std::size_t cellSize=cells.size();
00069     while(count != cellSize)
00070       {
00071       EntityHandle* connectivity;
00072       int numVerts=0, iterationCount=0;
00073       //use the highly efficent calls, since we know that are of the same dimension
00074       moab->connect_iterate(cells.begin()+count,
00075                             cells.end(),
00076                             connectivity,
00077                             numVerts,
00078                             iterationCount);
00079       //if we didn't read anything, break!
00080       if(iterationCount == 0)
00081         {
00082         break;
00083         }
00084 
00085       //identify the cell type that we currently have,
00086       //store that along with the connectivity in a temp storage vector
00087       const moab::EntityType type = moab->type_from_handle(*cells.begin()+count);
00088 
00089       //while all these cells are contiously of the same type,
00090       //quadric hexs in vtk have 20 points, but moab has 21 so we
00091       //need to store this difference
00092       int numVTKVerts = numVerts;
00093       int vtkCellType = smoab::vtkCellType(type,numVTKVerts);
00094 
00095       RunLengthInfo info = { vtkCellType, numVerts, (numVerts-numVTKVerts), iterationCount };
00096       this->Info.push_back(info);
00097       this->Connectivity.push_back(connectivity);
00098 
00099       count += iterationCount;
00100       }
00101     }
00102 
00103   //----------------------------------------------------------------------------
00104   void compactIds(vtkIdType& numCells, vtkIdType& connectivityLength)
00105     {
00106     //converts all the ids to be ordered starting at zero, and also
00107     //keeping the orginal logical ordering. Stores the result of this
00108     //operation in the unstrucutred grid that is passed in
00109 
00110     //lets determine the total length of the connectivity
00111     connectivityLength = 0;
00112     numCells = 0;
00113     for(InfoConstIterator i = this->Info.begin();
00114         i != this->Info.end();
00115         ++i)
00116       {
00117       connectivityLength += (*i).numCells * (*i).numVerts;
00118       numCells += (*i).numCells;
00119       }
00120 
00121     this->UniquePoints.reserve(connectivityLength);
00122 
00123     this->copyConnectivity(this->UniquePoints);
00124     std::sort(this->UniquePoints.begin(),this->UniquePoints.end());
00125 
00126     typedef std::vector<EntityHandle>::iterator EntityIterator;
00127     EntityIterator newEnd = std::unique(this->UniquePoints.begin(),
00128                                         this->UniquePoints.end());
00129 
00130     const std::size_t newSize = std::distance(this->UniquePoints.begin(),newEnd);
00131     this->UniquePoints.resize(newSize);
00132     }
00133 
00134   //----------------------------------------------------------------------------
00135   void moabPoints(smoab::Range& range) const
00136     {
00137     //from the documentation a reverse iterator is the fastest way
00138     //to insert into a range.
00139     std::copy(this->UniquePoints.rbegin(),
00140               this->UniquePoints.rend(),
00141               moab::range_inserter(range));
00142     }
00143 
00144   //----------------------------------------------------------------------------
00145   //copy the connectivity from the moab held arrays to the user input vector
00146   void copyConnectivity(std::vector<EntityHandle>& output) const
00147     {
00148     //walk the info to find the length of each sub connectivity array,
00149     //and insert them into the vector, ordering is implied by the order
00150     //the connecitivy sub array are added to this class
00151     ConnConstIterator c = this->Connectivity.begin();
00152     for(InfoConstIterator i = this->Info.begin();
00153         i != this->Info.end();
00154         ++i,++c)
00155       {
00156       //remember our Connectivity is a vector of pointers whose
00157       //length is held in the info vector.
00158       const int numUnusedPoints = (*i).numUnusedVerts;
00159       if(numUnusedPoints==0)
00160         {
00161         const int connLength = (*i).numCells * (*i).numVerts;
00162         std::copy(*c,*c+connLength,std::back_inserter(output));
00163         }
00164       else
00165         {
00166         //we have cell connectivity that we need to skip,
00167         //so we have to manual copy each cells connectivity
00168         const int size = (*i).numCells;
00169         const int numPoints = (*i).numVerts;
00170         for(int i=0; i < size; ++i)
00171           {
00172           std::copy(*c,*c+numPoints,std::back_inserter(output));
00173           }
00174         c+=numPoints + (*i).numUnusedVerts;
00175         }
00176 
00177       }
00178     }
00179 
00180   //copy the information from this contianer to a vtk cell array, and
00181   //related lookup information
00182   void copyToVtkCellInfo(vtkIdType* cellArray,
00183                          vtkIdType* cellLocations,
00184                          unsigned char* cellTypes) const
00185     {
00186     vtkIdType currentVtkConnectivityIndex = 0;
00187     ConnConstIterator c = this->Connectivity.begin();
00188     for(InfoConstIterator i = this->Info.begin();
00189         i != this->Info.end();
00190         ++i, ++c)
00191       {
00192       //for this group of the same cell type we need to fill the cellTypes
00193       const int numCells = (*i).numCells;
00194       const int numVerts = (*i).numVerts;
00195 
00196       std::fill_n(cellTypes,
00197                   numCells,
00198                   static_cast<unsigned char>((*i).type));
00199 
00200       //for each cell in this collection that have the same type
00201       //grab the raw array now, so we can properly increment for each vert in each cell
00202       EntityHandle* moabConnectivity = *c;
00203       for(int j=0;j < numCells; ++j)
00204         {
00205         cellLocations[j]= currentVtkConnectivityIndex;
00206 
00207         //cell arrays start and end are different, since we
00208         //have to account for element that states the length of each cell
00209         cellArray[0]=numVerts;
00210 
00211         for(int k=0; k < numVerts; ++k, ++moabConnectivity )
00212           {
00213           //this is going to be a root of some failures when we start
00214           //reading really large datasets under 32bit.
00215 
00216 
00217           //fyi, don't use a range ds for unique points, distance
00218           //function is horribly slow they need to override it
00219           EntityConstIterator result = std::lower_bound(
00220                                          this->UniquePoints.begin(),
00221                                          this->UniquePoints.end(),
00222                                          *moabConnectivity);
00223           std::size_t newId = std::distance(this->UniquePoints.begin(),
00224                                             result);
00225           cellArray[k+1] = static_cast<vtkIdType>(newId);
00226           }
00227 
00228         //skip any extra unused points, which is currnetly only
00229         //the extra center point in moab quadratic hex
00230         moabConnectivity+=(*i).numUnusedVerts;
00231 
00232         currentVtkConnectivityIndex += numVerts+1;
00233         cellArray += numVerts+1;
00234         }
00235 
00236       //For Tri-Quadratic-Hex and Quadratric-Wedge Moab and VTK
00237       //Differ on the order of the edge ids. For wedge we need to swap
00238       //indices 9,10,11 with 12,13,14 for each cell. For Hex we sawp
00239       //12,13,14,15 with 16,17,18,19
00240       int vtkCellType = (*i).type;
00241       vtkIdType* connectivity = cellArray - (numCells * (numVerts+1));
00242       if(vtkCellType == VTK_TRIQUADRATIC_HEXAHEDRON)
00243         {
00244         ::QuadratricOrdering<VTK_TRIQUADRATIC_HEXAHEDRON> newOrdering;
00245         ::FixQuadraticIdOrdering(connectivity, numCells, newOrdering);
00246         }
00247       else if(vtkCellType == VTK_QUADRATIC_WEDGE)
00248         {
00249         ::QuadratricOrdering<VTK_QUADRATIC_WEDGE> newOrdering;
00250         ::FixQuadraticIdOrdering(connectivity, numCells, newOrdering);
00251         }
00252 
00253       cellLocations += numCells;
00254       cellTypes += numCells;
00255       }
00256 
00257     }
00258 
00259 private:
00260   std::vector<EntityHandle*> Connectivity;
00261   std::vector<EntityHandle> UniquePoints;
00262 
00263   struct RunLengthInfo{ int type; int numVerts; int numUnusedVerts; int numCells; };
00264   std::vector<RunLengthInfo> Info;
00265 
00266   typedef std::vector<EntityHandle>::const_iterator EntityConstIterator;
00267   typedef std::vector<EntityHandle*>::const_iterator ConnConstIterator;
00268   typedef std::vector<RunLengthInfo>::const_iterator InfoConstIterator;
00269 };
00270 }
00271 #endif // __smoab_MixedCellConnectivity_h
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines