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