ViennaLS
Loading...
Searching...
No Matches
lsFromSurfaceMesh.hpp
Go to the documentation of this file.
1#pragma once
2
4
5#include <hrleTypes.hpp>
6
7#include <lsDomain.hpp>
8#include <lsMesh.hpp>
9
10namespace viennals {
11
12using namespace viennacore;
13
15template <class T, int D> class FromSurfaceMesh {
16 using hrleIndexType = viennahrle::IndexType;
17
19 class box {
20 VectorType<hrleIndexType, D - 1> xMin, xMax;
21
22 public:
25 box(const VectorType<hrleIndexType, D - 1> &idx0,
26 const VectorType<hrleIndexType, D - 1> &idx1) {
27 xMin = Min(idx0, idx1);
28 xMax = Max(idx0, idx1);
29 }
30
32 box(const box &box) {
33 xMin = box.xMin;
34 xMax = box.xMax;
35 }
36
38 explicit box(const VectorType<hrleIndexType, D - 1> &idx) {
39 xMin = idx;
40 xMax = idx;
41 }
42
43 box &operator=(const box &b) = default;
44
46 bool is_empty() const {
47 for (int i = 0; i < D - 1; i++)
48 if (xMax[i] < xMin[i])
49 return true;
50 return false;
51 }
52
54 const VectorType<hrleIndexType, D - 1> &min() const { return xMin; }
56 const VectorType<hrleIndexType, D - 1> &max() const { return xMax; }
57
59 class iterator {
60 VectorType<hrleIndexType, D - 1> pos;
61 const box &b;
62
63 public:
64 explicit iterator(const box &bx) : pos(bx.min()), b(bx) {}
65
67 int i;
68 for (i = 0; i < D - 2; i++) {
69 ++pos[i];
70 if (pos[i] <= b.xMax[i]) {
71 break;
72 } else {
73 pos[i] = b.xMin[i];
74 }
75 }
76 if (i == D - 2)
77 ++pos[i];
78 return *this;
79 }
80
82 iterator tmp = *this;
83 ++(*this);
84 return tmp;
85 }
86
87 bool is_finished() const { return (pos[D - 2] > b.xMax[D - 2]); }
88
89 const VectorType<hrleIndexType, D - 1> &operator*() const { return pos; }
90 };
91 };
92
93 SmartPointer<Domain<T, D>> levelSet = SmartPointer<Domain<T, D>>::New();
94 SmartPointer<Mesh<T>> mesh = SmartPointer<Mesh<T>>::New();
95 // bool removeBoundaryTriangles = true;
96 std::array<bool, 3> removeBoundaryTriangles{true, true, true};
97 T boundaryEps = 1e-5;
98 T distanceEps = 1e-4;
99 T signEps = 1e-6;
100
107 int calculateGridlineTriangleIntersection(VectorType<T, 2> Point,
108 const VectorType<T, 2> *c, int dir,
109 T &intersection) {
110
111 bool inside_pos = true;
112 bool inside_neg = true;
113
114 T A[2];
115
116 const int dirA = (dir + 1) % 2;
117
118 for (int k = 0; k < 2; k++) {
119 A[k] = c[(k + 1) % 2][dirA] - Point[dirA];
120 if (k == dir)
121 A[k] = -A[k];
122 if (A[k] < T(0))
123 inside_pos = false;
124 if (A[k] > T(0))
125 inside_neg = false;
126 }
127
128 if ((!inside_pos && inside_neg) || (inside_pos && !inside_neg)) {
129 T sum = A[0] + A[1];
130 int k = 0;
131 for (int i = 1; i < 2; ++i) {
132 if (inside_pos) {
133 if (A[i] > A[k])
134 k = i;
135 } else {
136 if (A[i] < A[k])
137 k = i;
138 }
139 }
140 intersection = c[k][dir] +
141 (c[(k + 1) % 2][dir] - c[k][dir]) * (A[(k + 1) % 2] / sum);
142 return (inside_pos) ? 1 : -1;
143 } else {
144 return 0;
145 }
146 }
147
154 int calculateGridlineTriangleIntersection(VectorType<T, 3> Point,
155 const VectorType<T, 3> *c, int dir,
156 T &intersection) {
157
158 bool inside_pos = true;
159 bool inside_neg = true;
160
161 T A[3];
162
163 const int dirA = (dir + 1) % 3;
164 const int dirB = (dir + 2) % 3;
165
166 for (int k = 0; k < 3; k++) {
167 bool swapped =
168 (c[(k + 1) % 3] <
169 c[(k + 2) % 3]); // necessary to guarantee anti commutativity
170 const VectorType<T, 3> &v1 = (swapped) ? c[(k + 2) % 3] : c[(k + 1) % 3];
171 const VectorType<T, 3> &v2 = (swapped) ? c[(k + 1) % 3] : c[(k + 2) % 3];
172
173 A[k] = (v1[dirA] - Point[dirA]) * (v2[dirB] - Point[dirB]) -
174 (v2[dirA] - Point[dirA]) * (v1[dirB] - Point[dirB]);
175
176 if (swapped)
177 A[k] = -A[k];
178
179 if (A[k] < T(0))
180 inside_pos = false;
181 if (A[k] > T(0))
182 inside_neg = false;
183 }
184
185 if ((!inside_pos && inside_neg) || (inside_pos && !inside_neg)) {
186 T sum = A[0] + A[1] + A[2];
187 int k = 0;
188 for (int i = 1; i < 3; ++i) {
189 if (inside_pos) {
190 if (A[i] > A[k])
191 k = i;
192 } else {
193 if (A[i] < A[k])
194 k = i;
195 }
196 }
197 intersection =
198 c[k][dir] +
199 (c[(k + 1) % 3][dir] - c[k][dir]) * (A[(k + 1) % 3] / sum) +
200 (c[(k + 2) % 3][dir] - c[k][dir]) * (A[(k + 2) % 3] / sum);
201 return (inside_pos) ? 1 : -1;
202 } else {
203 return 0;
204 }
205 }
206
207public:
208 FromSurfaceMesh() = default;
209
210 FromSurfaceMesh(SmartPointer<Domain<T, D>> passedLevelSet,
211 SmartPointer<Mesh<T>> passedMesh,
212 bool passedRemoveBoundaryTriangles = true)
213 : levelSet(passedLevelSet), mesh(passedMesh),
214 removeBoundaryTriangles{passedRemoveBoundaryTriangles,
215 passedRemoveBoundaryTriangles,
216 passedRemoveBoundaryTriangles} {}
217
218 void setLevelSet(SmartPointer<Domain<T, D>> passedLevelSet) {
219 levelSet = passedLevelSet;
220 }
221
222 void setMesh(SmartPointer<Mesh<T>> passedMesh) { mesh = passedMesh; }
223
227 void setRemoveBoundaryTriangles(bool passedRemoveBoundaryTriangles) {
228 removeBoundaryTriangles.fill(passedRemoveBoundaryTriangles);
229 }
230
235 template <std::size_t N>
237 std::array<bool, N> passedRemoveBoundaryTriangles) {
238 for (unsigned i = 0; i < D && i < N; ++i) {
239 removeBoundaryTriangles[i] = passedRemoveBoundaryTriangles[i];
240 }
241 }
242
243 void apply() {
244 if (levelSet == nullptr) {
245 Logger::getInstance()
246 .addWarning("No level set was passed to FromSurfaceMesh.")
247 .print();
248 return;
249 }
250 if (mesh == nullptr) {
251 Logger::getInstance()
252 .addWarning("No mesh was passed to FromSurfaceMesh.")
253 .print();
254 return;
255 }
256
257 const auto &grid = levelSet->getGrid();
258
259 // caluclate which directions should apply removeBoundaryTriangles
260 // it only makes sense to keep boundary triangles for periodic BNC
261 bool removeBoundaries[D];
262 for (unsigned i = 0; i < D; ++i) {
263 if (!removeBoundaryTriangles[i] && grid.isBoundaryPeriodic(i)) {
264 removeBoundaries[i] = false;
265 } else {
266 removeBoundaries[i] = true;
267 }
268 }
269
270 std::vector<std::pair<viennahrle::Index<D>, T>> points2;
271
272 // setup list of grid points with distances to surface elements
273 {
274 typedef std::vector<std::pair<viennahrle::Index<D>, std::pair<T, T>>>
275 point_vector;
276 point_vector points;
277 T gridDelta = grid.getGridDelta();
278
279 VectorType<T, D> gridMin, gridMax;
280 for (unsigned i = 0; i < D; ++i) {
281 gridMin[i] = grid.getMinIndex(i) * gridDelta;
282 gridMax[i] = grid.getMaxIndex(i) * gridDelta;
283 }
284
285 // for each surface element do
286 std::vector<std::array<unsigned, D>> &elements =
287 mesh->template getElements<D>();
288
289 for (unsigned currentElement = 0; currentElement < elements.size();
290 currentElement++) {
291 VectorType<T, D> nodes[D]; // nodes of element
292 VectorType<T, D> center; // center point of triangle
293 std::fill(center.begin(), center.end(), 0);
294
295 std::bitset<2 * D> flags;
296 flags.set();
297 bool removeElement = false;
298
299 for (int dim = 0; dim < D; dim++) {
300 for (int q = 0; q < D; q++) {
301 nodes[q][dim] = mesh->nodes[elements[currentElement][q]][dim];
302 if (std::abs(nodes[q][dim] - gridMin[dim]) <
303 boundaryEps * gridDelta)
304 nodes[q][dim] = gridMin[dim];
305 if (std::abs(nodes[q][dim] - gridMax[dim]) <
306 boundaryEps * gridDelta)
307 nodes[q][dim] = gridMax[dim];
308
309 if (nodes[q][dim] > gridMin[dim])
310 flags.reset(dim);
311 if (nodes[q][dim] < gridMax[dim])
312 flags.reset(dim + D);
313
314 center[dim] += nodes[q][dim]; // center point calculation
315 }
316 if (removeBoundaries[dim] && (flags[dim] || flags[dim + D])) {
317 removeElement = true;
318 }
319 }
320
321 if (removeElement)
322 continue;
323
324 // center point calculation
325 for (int q = 0; q < D; q++)
326 center[q] /= static_cast<T>(D);
327
328 // VectorType<T,D> normal=NormalVector(c); //normal vector
329 // calculation
330
331 // determine the minimum and maximum nodes of the element, based on
332 // their coordinates
333 VectorType<T, D> minNode = nodes[0];
334 VectorType<T, D> maxNode = nodes[0];
335 for (int i = 1; i < D; ++i) {
336 minNode = Min(minNode, nodes[i]);
337 maxNode = Max(maxNode, nodes[i]);
338 }
339
340 // find indices which describe the triangle
341 viennahrle::Index<D> minIndex, maxIndex;
342 for (int q = 0; q < D; q++) {
343 minIndex[q] =
344 static_cast<hrleIndexType>(std::ceil(minNode[q] / gridDelta));
345 maxIndex[q] =
346 static_cast<hrleIndexType>(std::floor(maxNode[q] / gridDelta));
347 }
348
349 // Only iterate over indices which are actually inside the domain
350 // if boundary conditions should be ignored
351 for (unsigned i = 0; i < D; ++i) {
352 if (removeBoundaries[i]) {
353 minIndex[i] = std::max(minIndex[i], grid.getMinIndex(i));
354 // for Periodic BNC, MaxGridPoint is MaxIndex-1
355 maxIndex[i] = std::min(maxIndex[i], grid.getMaxGridPoint(i));
356 }
357 }
358
359 // each cartesian direction
360 for (int z = 0; z < D; z++) {
361 VectorType<hrleIndexType, D - 1> min_bb, max_bb;
362
363 for (int h = 0; h < D - 1; ++h) {
364 min_bb[h] = minIndex[(z + h + 1) % D];
365 max_bb[h] = maxIndex[(z + h + 1) % D];
366 }
367
368 box bb(min_bb, max_bb);
369
370 if (bb.is_empty())
371 continue;
372
373 for (typename box::iterator it_bb(bb); !it_bb.is_finished();
374 ++it_bb) {
375
376 viennahrle::Index<D> it_b;
377 for (int h = 0; h < D - 1; ++h)
378 it_b[(z + h + 1) % D] = (*it_bb)[h];
379
380 VectorType<T, D> p;
381 for (int k = 1; k < D; k++)
382 p[(k + z) % D] = grid.gridPositionOfGlobalIndex(
383 (k + z) % D, it_b[(k + z) % D]);
384
385 T intersection;
386 int intersection_status = calculateGridlineTriangleIntersection(
387 p, nodes, z, intersection);
388
389 if (intersection_status != 0) {
390 // if there is an intersection
391 assert(intersection >= minNode[z] &&
392 "Intersection should be inside Domain!");
393 assert(intersection <= maxNode[z] &&
394 "Intersection should be inside Domain!");
395 // intersection = std::max(intersection, minNode[z]);
396 // intersection = std::min(intersection, maxNode[z]);
397
398 if (removeBoundaries[z] &&
399 intersection > gridDelta * (grid.getMaxIndex(z)))
400 continue;
401 if (removeBoundaries[z] &&
402 intersection < grid.getMinLocalCoordinate(z))
403 continue;
404
405 T intersection2 =
406 grid.globalCoordinate2LocalIndex(z, intersection);
407
408 auto floor = static_cast<hrleIndexType>(
409 std::floor(intersection2 - distanceEps));
410 auto ceil = static_cast<hrleIndexType>(
411 std::ceil(intersection2 + distanceEps));
412
413 floor = std::max(floor, minIndex[z] - 1);
414 ceil = std::min(ceil, maxIndex[z] + 1);
415 floor = std::max(floor, grid.getMinIndex(z));
416 ceil = std::min(ceil, grid.getMaxIndex(z));
417
418 if (!removeBoundaries[z]) {
419 floor = grid.globalIndex2LocalIndex(z, floor);
420 ceil = grid.globalIndex2LocalIndex(z, ceil);
421 }
422
423 VectorType<T, D> t = center;
424 t[z] -= intersection;
425 for (int k = 1; k < D; k++)
426 t[(z + k) % D] -= p[(z + k) % D];
427 Normalize(t);
428
429 for (it_b[z] = floor; it_b[z] <= ceil; ++(it_b[z])) {
430
431 T RealDistance = it_b[z] - intersection2;
432
433 T SignDistance = RealDistance;
434
435 SignDistance -= signEps * t[z];
436
437 if (intersection_status < 0) {
438 RealDistance = -RealDistance;
439 SignDistance = -SignDistance;
440 }
441
442 if (RealDistance < 0.)
443 SignDistance = -SignDistance;
444
445 RealDistance *= (1. - distanceEps * 1e-3);
446 if (RealDistance > 1.)
447 RealDistance = 1.;
448 if (RealDistance < -1.)
449 RealDistance = -1.;
450
451 if (RealDistance == 0.)
452 RealDistance = 0.; // to avoid zeros with negative sign
453
454 points.push_back(
455 std::make_pair(grid.globalIndices2LocalIndices(it_b),
456 std::make_pair(SignDistance, RealDistance)));
457 }
458 }
459 }
460 }
461 }
462
463 std::sort(points.begin(), points.end()); // sort points lexicographically
464
465 // setup list of index/distance pairs for level set initialization
466 auto it_points = points.begin();
467 while (it_points != points.end()) {
468 viennahrle::Index<D> tmp = it_points->first;
469 points2.push_back(
470 std::make_pair(it_points->first, it_points->second.second));
471
472 do {
473 ++it_points;
474 } while ((it_points != points.end()) && (it_points->first == tmp));
475 }
476 }
477
478 levelSet->insertPoints(points2); // initialize level set function
479 levelSet->getDomain().segment(); // distribute points evenly across threads
480 levelSet->finalize(2);
481 }
482};
483
484// add all template specialisations for this class
485PRECOMPILE_PRECISION_DIMENSION(FromSurfaceMesh)
486
487} // namespace viennals
Class containing all information about the level set, including the dimensions of the domain,...
Definition lsDomain.hpp:27
Iterator over all grid points, contained by a box.
Definition lsFromSurfaceMesh.hpp:59
iterator & operator++()
Definition lsFromSurfaceMesh.hpp:66
iterator operator++(int)
Definition lsFromSurfaceMesh.hpp:81
const VectorType< hrleIndexType, D - 1 > & operator*() const
Definition lsFromSurfaceMesh.hpp:89
iterator(const box &bx)
Definition lsFromSurfaceMesh.hpp:64
bool is_finished() const
Definition lsFromSurfaceMesh.hpp:87
void setRemoveBoundaryTriangles(bool passedRemoveBoundaryTriangles)
Set whether all triangles outside the domain should be ignored (=true) or whether boundary conditions...
Definition lsFromSurfaceMesh.hpp:227
void setRemoveBoundaryTriangles(std::array< bool, N > passedRemoveBoundaryTriangles)
Set whether all triangles outside the domain should be ignored (=true) or whether boundary conditions...
Definition lsFromSurfaceMesh.hpp:236
FromSurfaceMesh(SmartPointer< Domain< T, D > > passedLevelSet, SmartPointer< Mesh< T > > passedMesh, bool passedRemoveBoundaryTriangles=true)
Definition lsFromSurfaceMesh.hpp:210
void setLevelSet(SmartPointer< Domain< T, D > > passedLevelSet)
Definition lsFromSurfaceMesh.hpp:218
void apply()
Definition lsFromSurfaceMesh.hpp:243
void setMesh(SmartPointer< Mesh< T > > passedMesh)
Definition lsFromSurfaceMesh.hpp:222
This class holds an explicit mesh, which is always given in 3 dimensions. If it describes a 2D mesh,...
Definition lsMesh.hpp:21
#define PRECOMPILE_PRECISION_DIMENSION(className)
Definition lsPreCompileMacros.hpp:24
Definition lsAdvect.hpp:36
constexpr int D
Definition pyWrap.cpp:70
double T
Definition pyWrap.cpp:68