TTK
Loading...
Searching...
No Matches
ThreeSkeleton.cpp
Go to the documentation of this file.
1#include <ThreeSkeleton.h>
2#include <boost/container/small_vector.hpp>
3
4using namespace ttk;
5
7 // TODO could use typeid for automatic value
8 setDebugMsgPrefix("ThreeSkeleton");
9}
10
12 const SimplexId &vertexNumber,
13 const CellArray &cellArray,
14 FlatJaggedArray &cellNeighbors,
15 FlatJaggedArray *triangleStars) const {
16
17 auto localTriangleStars = triangleStars;
18 FlatJaggedArray defaultTriangleStars{};
19 if(!localTriangleStars) {
20 localTriangleStars = &defaultTriangleStars;
21 }
22
23 if(localTriangleStars->empty()) {
24 TwoSkeleton twoSkeleton;
25 twoSkeleton.setThreadNumber(threadNumber_);
26 twoSkeleton.setDebugLevel(debugLevel_);
27 twoSkeleton.buildTriangleList(
28 vertexNumber, cellArray, nullptr, localTriangleStars);
29 }
30
31 Timer t;
32
33 printMsg("Building cell neighbors", 0, 0, 1, ttk::debug::LineMode::REPLACE);
34
35 const SimplexId cellNumber = cellArray.getNbCells();
36 const SimplexId triangleNumber = localTriangleStars->size();
37 std::vector<SimplexId> offsets(cellNumber + 1);
38 // number of neighbors processed per cell
39 std::vector<SimplexId> neighborsId(cellNumber);
40
41 for(SimplexId i = 0; i < triangleNumber; i++) {
42 if(localTriangleStars->size(i) == 2) {
43 // tetra cells in triangle i's star
44 const auto cs0 = localTriangleStars->get(i, 0);
45 const auto cs1 = localTriangleStars->get(i, 1);
46 offsets[cs0 + 1]++;
47 offsets[cs1 + 1]++;
48 }
49 }
50
51 // compute partial sum of number of neighbors per vertex
52 for(size_t i = 1; i < offsets.size(); ++i) {
53 offsets[i] += offsets[i - 1];
54 }
55
56 // allocate flat neighbors vector
57 std::vector<SimplexId> neighbors(offsets.back());
58
59 // fill flat neighbors vector using offsets and neighbors count vectors
60 for(SimplexId i = 0; i < triangleNumber; i++) {
61 if(localTriangleStars->size(i) == 2) {
62 // tetra cells in triangle i's star
63 const auto cs0 = localTriangleStars->get(i, 0);
64 const auto cs1 = localTriangleStars->get(i, 1);
65 neighbors[offsets[cs0] + neighborsId[cs0]] = cs1;
66 neighborsId[cs0]++;
67 neighbors[offsets[cs1] + neighborsId[cs1]] = cs0;
68 neighborsId[cs1]++;
69 }
70 }
71
72 // fill FlatJaggedArray struct
73 cellNeighbors.setData(std::move(neighbors), std::move(offsets));
74
75 printMsg("Built " + std::to_string(cellNumber) + " cell neighbors", 1,
76 t.getElapsedTime(), 1);
77
78 // ethaneDiol.vtu, 8.7Mtets, vger (4coresHT)
79 // 1 thread: 9.80 s
80 // 4 threads: 14.18157 s
81
82 // ethaneDiol.vtu, 8.7Mtets, richard (8coresHT)
83 // 1 thread: 9.6763 s
84 // 4 threads: 16.7489 s
85
86 // ethaneDiol.vtu, 8.7Mtets, vger (4coresHT), VTK implementation
87 // 1 thread: 14.9023 s
88 // 4 threads: 8.58048 s
89
90 // ethaneDiol.vtu, 8.7Mtets, richard (8coresHT), VTK implementation
91 // 1 thread: 15.1387 s
92 // 8 threads: 5.38286 s
93
94 return 0;
95}
96
98 const SimplexId &vertexNumber,
99 const CellArray &cellArray,
100 FlatJaggedArray &cellNeighbors,
101 FlatJaggedArray *vertexStars) const {
102
103 // TODO: ASSUME uniform mesh here!
104 if(cellArray.getNbCells() && cellArray.getCellVertexNumber(0) == 3) {
105 TwoSkeleton twoSkeleton;
106 twoSkeleton.setDebugLevel(debugLevel_);
107 twoSkeleton.setThreadNumber(threadNumber_);
108 return twoSkeleton.buildCellNeighborsFromVertices(
109 vertexNumber, cellArray, cellNeighbors, vertexStars);
110 }
111
112 if(cellArray.getNbCells() && cellArray.getCellVertexNumber(0) <= 2) {
113 // 1D
114 printErr("buildCellNeighborsFromVertices in 1D:");
115 printErr("Not implemented! TODO?!");
116 return -1;
117 }
118
119 auto localVertexStars = vertexStars;
120 FlatJaggedArray defaultVertexStars{};
121
122 if(!localVertexStars) {
123 localVertexStars = &defaultVertexStars;
124 }
125
126 if(localVertexStars->empty()) {
127
128 ZeroSkeleton zeroSkeleton;
129 zeroSkeleton.setThreadNumber(threadNumber_);
130 zeroSkeleton.setDebugLevel(debugLevel_);
131 zeroSkeleton.buildVertexStars(vertexNumber, cellArray, *localVertexStars);
132 }
133
134 Timer t;
135
136 printMsg("Building cell neighnors", 0, 0, threadNumber_,
138
139 const SimplexId cellNumber = cellArray.getNbCells();
140 using boost::container::small_vector;
141 // for each cell/tetra, a vector of neighbors
142 std::vector<small_vector<SimplexId, 4>> neighbors(cellNumber);
143
144#ifdef TTK_ENABLE_OPENMP
145#pragma omp parallel for num_threads(threadNumber_)
146#endif
147 for(SimplexId cid = 0; cid < cellNumber; cid++) {
148 const SimplexId nbVertCell = cellArray.getCellVertexNumber(cid);
149
150 // go triangle by triangle
151 for(SimplexId j = 0; j < nbVertCell; j++) {
152
153 SimplexId v0 = cellArray.getCellVertex(cid, j);
154 SimplexId v1 = cellArray.getCellVertex(cid, (j + 1) % nbVertCell);
155 SimplexId v2 = cellArray.getCellVertex(cid, (j + 2) % nbVertCell);
156
157 // perform an intersection of the 3 (sorted) star lists
158 SimplexId pos0 = 0, pos1 = 0, pos2 = 0;
159 SimplexId intersection = -1;
160
161 while(pos0 < localVertexStars->size(v0)
162 && pos1 < localVertexStars->size(v1)
163 && pos2 < localVertexStars->size(v2)) {
164
165 SimplexId biggest = localVertexStars->get(v0, pos0);
166 if(localVertexStars->get(v1, pos1) > biggest) {
167 biggest = localVertexStars->get(v1, pos1);
168 }
169 if(localVertexStars->get(v2, pos2) > biggest) {
170 biggest = localVertexStars->get(v2, pos2);
171 }
172
173 for(SimplexId l = pos0; l < localVertexStars->size(v0); l++) {
174 if(localVertexStars->get(v0, l) < biggest) {
175 pos0++;
176 } else {
177 break;
178 }
179 }
180 for(SimplexId l = pos1; l < localVertexStars->size(v1); l++) {
181 if(localVertexStars->get(v1, l) < biggest) {
182 pos1++;
183 } else {
184 break;
185 }
186 }
187 for(SimplexId l = pos2; l < localVertexStars->size(v2); l++) {
188 if(localVertexStars->get(v2, l) < biggest) {
189 pos2++;
190 } else {
191 break;
192 }
193 }
194
195 if(pos0 < localVertexStars->size(v0)
196 && pos1 < localVertexStars->size(v1)
197 && pos2 < localVertexStars->size(v2)) {
198
199 if((localVertexStars->get(v0, pos0)
200 == localVertexStars->get(v1, pos1))
201 && (localVertexStars->get(v0, pos0)
202 == localVertexStars->get(v2, pos2))) {
203
204 if(localVertexStars->get(v0, pos0) != cid) {
205 intersection = localVertexStars->get(v0, pos0);
206 break;
207 }
208
209 pos0++;
210 pos1++;
211 pos2++;
212 }
213 }
214 }
215
216 if(intersection != -1) {
217 neighbors[cid].emplace_back(intersection);
218 }
219 }
220 }
221
222 // convert to a FlatJaggedArray
223 cellNeighbors.fillFrom(neighbors);
224
225 printMsg("Built " + std::to_string(cellNumber) + " cell neighbors", 1,
227
228 // ethaneDiol.vtu, 8.7Mtets, richard (4coresHT)
229 // 1 thread: 9.39488 s
230 // 2 threads: 5.93132 s [faster than any other implementation]
231 // 4 threads: 3.3445 s [~x3, not too bad]
232
233 // ethaneDiol.vtu, 8.7Mtets, hal9000 (12coresHT)
234 // 1 thread: 4.4475 s [speedup on VTK ~2.5]
235 // 24 threads: 1.4488 s
236
237 // ethaneDiolMedium.vtu, 70Mtets, hal9000 (12coresHT)
238 // 1 thread: 29.4951 s [speedup on VTK ~3.3]
239 // 24 threads: 12.7058 s (parallel speed up ~3, speed up on VTK: 2)
240
241 // NOTE:
242 // when dealing with simple types such as int, it's better to have heap
243 // variable (instead of a vector<int> indexed by threadId).
244 // it is MUCH MUCH MUCH better
245 // especially with write access (cache misses)
246
247 return 0;
248}
int threadNumber_
Definition: BaseClass.h:95
virtual int setThreadNumber(const int threadNumber)
Definition: BaseClass.h:80
CellArray generic array of cells
LongSimplexId getCellVertex(const LongSimplexId cellId, const SimplexId localVertId) const
LongSimplexId getNbCells() const
Get the number of cells in the array.
SimplexId getCellVertexNumber(const LongSimplexId cellId) const
int debugLevel_
Definition: Debug.h:379
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
virtual int setDebugLevel(const int &debugLevel)
Definition: Debug.cpp:147
int printMsg(const std::string &msg, const debug::Priority &priority=debug::Priority::INFO, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cout) const
Definition: Debug.h:118
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:149
Replacement for std::vector<std::vector<SimplexId>>
void setData(std::vector< SimplexId > &&data, std::vector< SimplexId > &&offsets)
Set internal data from pre-existing vectors.
void fillFrom(const std::vector< T > &src, int threadNumber=1)
Fill buffers from a std::vector<std::vector<SimplexId>>
int buildCellNeighborsFromVertices(const SimplexId &vertexNumber, const CellArray &cellArray, FlatJaggedArray &cellNeighbors, FlatJaggedArray *vertexStars=nullptr) const
int buildCellNeighborsFromTriangles(const SimplexId &vertexNumber, const CellArray &cellArray, FlatJaggedArray &cellNeighbors, FlatJaggedArray *triangleStars=nullptr) const
double getElapsedTime()
Definition: Timer.h:15
TwoSkeleton processing package.
Definition: TwoSkeleton.h:24
int buildTriangleList(const SimplexId &vertexNumber, const CellArray &cellArray, std::vector< std::array< SimplexId, 3 > > *triangleList=nullptr, FlatJaggedArray *triangleStars=nullptr, std::vector< std::array< SimplexId, 4 > > *cellTriangleList=nullptr) const
int buildCellNeighborsFromVertices(const SimplexId &vertexNumber, const CellArray &cellArray, FlatJaggedArray &cellNeighbors, FlatJaggedArray *vertexStars=nullptr) const
Definition: TwoSkeleton.cpp:65
ZeroSkeleton processing package.
Definition: ZeroSkeleton.h:25
int buildVertexStars(const SimplexId &vertexNumber, const CellArray &cellArray, FlatJaggedArray &vertexStars) const
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22