TTK
Loading...
Searching...
No Matches
OneSkeleton.cpp
Go to the documentation of this file.
1#include <OneSkeleton.h>
2#include <boost/container/small_vector.hpp>
3
4using namespace ttk;
5
7 setDebugMsgPrefix("OneSkeleton");
8}
9
10// 2D cells (triangles)
12 const std::vector<std::array<SimplexId, 2>> &edgeList,
13 const FlatJaggedArray &edgeStars,
14 const CellArray &cellArray,
15 FlatJaggedArray &edgeLinks) const {
16
17#ifndef TTK_ENABLE_KAMIKAZE
18 if(edgeList.empty())
19 return -1;
20 if(edgeStars.size() != edgeList.size())
21 return -2;
22#endif
23
24 Timer t;
25
26 const SimplexId edgeNumber = edgeStars.size();
27 std::vector<SimplexId> offsets(edgeNumber + 1);
28 // one vertex per star
29 std::vector<SimplexId> links(edgeStars.dataSize());
30
32 "Building edge links", 0, 0, threadNumber_, ttk::debug::LineMode::REPLACE);
33
34#ifdef TTK_ENABLE_OPENMP
35#pragma omp parallel for num_threads(threadNumber_)
36#endif // TTK_ENABLE_OPENMP
37 for(SimplexId i = 0; i < edgeNumber; i++) {
38 // copy the edgeStars offsets array
39 offsets[i] = edgeStars.offset(i);
40
41 for(SimplexId j = 0; j < edgeStars.size(i); j++) {
42 // for each cell/triangle in edge i's star, get the opposite vertex
43 for(int k = 0; k < 3; k++) {
44 const auto v = cellArray.getCellVertex(edgeStars.get(i, j), k);
45 if(v != edgeList[i][0] && v != edgeList[i][1]) {
46 // edge i does not contain vertex i
47 links[offsets[i] + j] = v;
48 break;
49 }
50 }
51 }
52 }
53
54 // don't forget the last offset
55 offsets[edgeNumber] = edgeStars.offset(edgeNumber);
56
57 edgeLinks.setData(std::move(links), std::move(offsets));
58
59 printMsg("Built " + std::to_string(edgeNumber) + " edge links", 1,
61
62 return 0;
63}
64
65// 3D cells (tetrahedron)
67 const std::vector<std::array<SimplexId, 2>> &edgeList,
68 const FlatJaggedArray &edgeStars,
69 const std::vector<std::array<SimplexId, 6>> &cellEdges,
70 FlatJaggedArray &edgeLinks) const {
71
72#ifndef TTK_ENABLE_KAMIKAZE
73 if(edgeList.empty())
74 return -1;
75 if((edgeStars.empty()) || (edgeStars.size() != edgeList.size()))
76 return -2;
77 if(cellEdges.empty())
78 return -3;
79#endif
80
81 Timer t;
82
83 const SimplexId edgeNumber = edgeStars.size();
84 std::vector<SimplexId> offsets(edgeNumber + 1);
85 // one edge per star
86 std::vector<SimplexId> links(edgeStars.dataSize());
87
89 "Building edge links", 0, 0, threadNumber_, debug::LineMode::REPLACE);
90
91#ifdef TTK_ENABLE_OPENMP
92#pragma omp parallel for num_threads(threadNumber_)
93#endif // TTK_ENABLE_OPENMP
94 for(SimplexId i = 0; i < edgeNumber; i++) {
95 // copy the edgeStars offsets array
96 offsets[i] = edgeStars.offset(i);
97
98 // current edge vertices
99 const auto &e = edgeList[i];
100 for(SimplexId j = 0; j < edgeStars.size(i); j++) {
101 const auto c = edgeStars.get(i, j);
102 for(size_t k = 0; k < cellEdges[c].size(); k++) {
103 // cell edge id
104 const auto ceid = cellEdges[c][k];
105 // cell edge vertices
106 const auto &ce = edgeList[ceid];
107
108 if(ce[0] != e[0] && ce[0] != e[1] && ce[1] != e[0] && ce[1] != e[1]) {
109 // ce and e have no vertex in common
110 links[offsets[i] + j] = ceid;
111 break;
112 }
113 }
114 }
115 }
116
117 // don't forget the last offset
118 offsets[edgeNumber] = edgeStars.offset(edgeNumber);
119
120 edgeLinks.setData(std::move(links), std::move(offsets));
121
122 printMsg("Built " + std::to_string(edgeNumber) + " edge links", 1,
124
125 return 0;
126}
127
128using edgeType = std::array<SimplexId, 2>;
129
130template <std::size_t n>
131std::array<edgeType, n> getLocalEdges(const CellArray &ttkNotUsed(cellArray),
132 const SimplexId ttkNotUsed(cid)) {
133 return {};
134}
135
136template <>
137std::array<edgeType, 1> getLocalEdges(const CellArray &cellArray,
138 const SimplexId cid) {
139 return {
140 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 0)),
141 static_cast<SimplexId>(cellArray.getCellVertex(cid, 1))},
142 };
143}
144
145template <>
146std::array<edgeType, 3> getLocalEdges(const CellArray &cellArray,
147 const SimplexId cid) {
148 // triangle case: {0-1}, {0-2}, {1-2}
149 return {
150 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 0)),
151 static_cast<SimplexId>(cellArray.getCellVertex(cid, 1))},
152
153 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 0)),
154 static_cast<SimplexId>(cellArray.getCellVertex(cid, 2))},
155
156 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 1)),
157 static_cast<SimplexId>(cellArray.getCellVertex(cid, 2))},
158 };
159}
160
161template <>
162std::array<edgeType, 4> getLocalEdges(const CellArray &cellArray,
163 const SimplexId cid) {
164 // quad case: {0-1}, {1-2}, {2-3}, {3-0}
165 return {
166 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 0)),
167 static_cast<SimplexId>(cellArray.getCellVertex(cid, 1))},
168 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 1)),
169 static_cast<SimplexId>(cellArray.getCellVertex(cid, 2))},
170 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 2)),
171 static_cast<SimplexId>(cellArray.getCellVertex(cid, 3))},
172 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 3)),
173 static_cast<SimplexId>(cellArray.getCellVertex(cid, 0))},
174 };
175}
176
177template <>
178std::array<edgeType, 6> getLocalEdges(const CellArray &cellArray,
179 const SimplexId cid) {
180 // tet case: {0-1}, {0-2}, {0-3}, {1-2}, {1-3}, {2-3}
181 return {
182 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 0)),
183 static_cast<SimplexId>(cellArray.getCellVertex(cid, 1))},
184 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 0)),
185 static_cast<SimplexId>(cellArray.getCellVertex(cid, 2))},
186 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 0)),
187 static_cast<SimplexId>(cellArray.getCellVertex(cid, 3))},
188 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 1)),
189 static_cast<SimplexId>(cellArray.getCellVertex(cid, 2))},
190 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 1)),
191 static_cast<SimplexId>(cellArray.getCellVertex(cid, 3))},
192 edgeType{static_cast<SimplexId>(cellArray.getCellVertex(cid, 2)),
193 static_cast<SimplexId>(cellArray.getCellVertex(cid, 3))},
194 };
195}
196
197template <std::size_t n>
199 const SimplexId &vertexNumber,
200 const CellArray &cellArray,
201 std::vector<std::array<SimplexId, 2>> &edgeList,
202 FlatJaggedArray &edgeStars,
203 std::vector<std::array<SimplexId, n>> &cellEdgeList) const {
204
205 Timer t;
206
207 // check parameters consistency (we need n to be consistent with the
208 // dimensionality of the mesh)
209 const size_t dim = cellArray.getCellVertexNumber(0) - 1;
210 // n == 4 -> non simplicial mesh (quadrangulation), skip check
211 if(n != dim * (dim + 1) / 2 && n != 4) {
212 this->printErr("Wrong template parameter (" + std::to_string(n)
213 + " edges per " + std::to_string(dim) + "D cell)");
214 this->printErr("Cannot build edge list");
215 return -1;
216 }
217
218 printMsg("Building edges", 0, 0, 1, ttk::debug::LineMode::REPLACE);
219
220 const SimplexId cellNumber = cellArray.getNbCells();
221 cellEdgeList.resize(cellNumber);
222
223 struct EdgeData {
224 // the id of the edge higher vertex
225 SimplexId highVert{};
226 // the edge id
227 SimplexId id{};
228 EdgeData(SimplexId hv, SimplexId i) : highVert{hv}, id{i} {
229 }
230 };
231
232 using boost::container::small_vector;
233 // for each vertex, a vector of EdgeData
234 std::vector<small_vector<EdgeData, 8>> edgeTable(vertexNumber);
235
236 const int timeBuckets = std::min<ttk::SimplexId>(10, cellNumber);
237 SimplexId edgeCount{};
238
239 for(SimplexId cid = 0; cid < cellNumber; cid++) {
240
241 // id of edge in cell
242 SimplexId ecid{};
243 const auto localEdges{getLocalEdges<n>(cellArray, cid)};
244
245 for(const auto &le : localEdges) {
246 // edge processing
247 SimplexId v0 = le[0];
248 SimplexId v1 = le[1];
249 if(v0 > v1) {
250 std::swap(v0, v1);
251 }
252 auto &vec = edgeTable[v0];
253 const auto pos
254 = std::find_if(vec.begin(), vec.end(),
255 [&](const EdgeData &a) { return a.highVert == v1; });
256 if(pos == vec.end()) {
257 // not found in edgeTable: new edge
258 vec.emplace_back(EdgeData{v1, edgeCount});
259 cellEdgeList[cid][ecid] = edgeCount;
260 edgeCount++;
261 } else {
262 // found an existing edge
263 cellEdgeList[cid][ecid] = pos->id;
264 }
265 ecid++;
266 }
267 if(debugLevel_ >= (int)(debug::Priority::INFO)) {
268 if(!(cid % ((cellNumber) / timeBuckets)))
269 printMsg("Building edges", (cid / (float)cellNumber),
271 }
272 }
273
274 // allocate & fill edgeList in parallel
275 edgeList.resize(edgeCount);
276
277#ifdef TTK_ENABLE_OPENMP
278#pragma omp parallel for num_threads(this->threadNumber_)
279#endif // TTK_ENABLE_OPENMP
280 for(SimplexId i = 0; i < vertexNumber; ++i) {
281 const auto &etable = edgeTable[i];
282 for(const auto &data : etable) {
283 edgeList[data.id] = {i, data.highVert};
284 }
285 }
286
287 // return cellEdgeList to get edgeStars
288 std::vector<SimplexId> offsets(edgeCount + 1);
289 // number of cells processed per edge
290 std::vector<SimplexId> starIds(edgeCount);
291
292 // store number of cells per edge
293 for(const auto &ce : cellEdgeList) {
294 for(const auto eid : ce) {
295 offsets[eid + 1]++;
296 }
297 }
298
299 // compute partial sum of number of cells per edge
300 for(size_t i = 1; i < offsets.size(); ++i) {
301 offsets[i] += offsets[i - 1];
302 }
303
304 // allocate flat edge stars vector
305 std::vector<SimplexId> edgeSt(offsets.back());
306
307 // fill flat neighbors vector using offsets and neighbors count vectors
308 for(size_t i = 0; i < cellEdgeList.size(); ++i) {
309 const auto &ce{cellEdgeList[i]};
310 for(const auto eid : ce) {
311 edgeSt[offsets[eid] + starIds[eid]] = i;
312 starIds[eid]++;
313 }
314 }
315
316 // fill FlatJaggedArray struct
317 edgeStars.setData(std::move(edgeSt), std::move(offsets));
318
319 printMsg(
320 "Built " + std::to_string(edgeCount) + " edges", 1, t.getElapsedTime(), 1);
321
322 // ethaneDiolMedium.vtu, 70Mtets, hal9000 (12coresHT)
323 // 1 thread: 10.4979 s
324 // 24 threads: 12.3994 s [not efficient in parallel]
325
326 return 0;
327}
328
329// explicit template instantiation for 1D cells (edges)
330template int OneSkeleton::buildEdgeList<1>(
331 const SimplexId &vertexNumber,
332 const CellArray &cellArray,
333 std::vector<std::array<SimplexId, 2>> &edgeList,
334 FlatJaggedArray &edgeStars,
335 std::vector<std::array<SimplexId, 1>> &cellEdgeList) const;
336
337// explicit template instantiation for 2D cells (triangles)
338template int OneSkeleton::buildEdgeList<3>(
339 const SimplexId &vertexNumber,
340 const CellArray &cellArray,
341 std::vector<std::array<SimplexId, 2>> &edgeList,
342 FlatJaggedArray &edgeStars,
343 std::vector<std::array<SimplexId, 3>> &cellEdgeList) const;
344
345// explicit template instantiation for 2D cells (quads)
346template int OneSkeleton::buildEdgeList<4>(
347 const SimplexId &vertexNumber,
348 const CellArray &cellArray,
349 std::vector<std::array<SimplexId, 2>> &edgeList,
350 FlatJaggedArray &edgeStars,
351 std::vector<std::array<SimplexId, 4>> &cellEdgeList) const;
352
353// explicit template instantiation for 3D cells (tetrathedron)
354template int OneSkeleton::buildEdgeList<6>(
355 const SimplexId &vertexNumber,
356 const CellArray &cellArray,
357 std::vector<std::array<SimplexId, 2>> &edgeList,
358 FlatJaggedArray &edgeStars,
359 std::vector<std::array<SimplexId, 6>> &cellEdgeList) const;
#define ttkNotUsed(x)
Mark function/method parameters that are not used in the function body at all.
Definition: BaseClass.h:47
std::array< SimplexId, 2 > edgeType
std::array< edgeType, n > getLocalEdges(const CellArray &ttkNotUsed(cellArray), const SimplexId ttkNotUsed(cid))
int threadNumber_
Definition: BaseClass.h:95
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
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>>
size_t dataSize() const
Returns the size of the data_ member.
SimplexId get(SimplexId id, SimplexId local) const
Returns the data inside the sub-vectors.
void setData(std::vector< SimplexId > &&data, std::vector< SimplexId > &&offsets)
Set internal data from pre-existing vectors.
SimplexId size(SimplexId id) const
Get the size of a particular sub-vector.
SimplexId offset(SimplexId id) const
Get the offset of a particular sub-vector.
bool empty() const
If the underlying buffers are empty.
int buildEdgeLinks(const std::vector< std::array< SimplexId, 2 > > &edgeList, const FlatJaggedArray &edgeStars, const CellArray &cellArray, FlatJaggedArray &edgeLinks) const
Definition: OneSkeleton.cpp:11
int buildEdgeList(const SimplexId &vertexNumber, const CellArray &cellArray, std::vector< std::array< SimplexId, 2 > > &edgeList, FlatJaggedArray &edgeStars, std::vector< std::array< SimplexId, n > > &cellEdgeList) const
double getElapsedTime()
Definition: Timer.h:15
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22