MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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