moab
|
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