244 if (levelSet ==
nullptr) {
245 Logger::getInstance()
246 .addWarning(
"No level set was passed to FromSurfaceMesh.")
250 if (mesh ==
nullptr) {
251 Logger::getInstance()
252 .addWarning(
"No mesh was passed to FromSurfaceMesh.")
257 const auto &grid = levelSet->getGrid();
261 bool removeBoundaries[
D];
262 for (
unsigned i = 0; i <
D; ++i) {
263 if (!removeBoundaryTriangles[i] && grid.isBoundaryPeriodic(i)) {
264 removeBoundaries[i] =
false;
266 removeBoundaries[i] =
true;
270 std::vector<std::pair<viennahrle::Index<D>,
T>> points2;
274 typedef std::vector<std::pair<viennahrle::Index<D>, std::pair<T, T>>>
277 T gridDelta = grid.getGridDelta();
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;
286 std::vector<std::array<unsigned, D>> &elements =
287 mesh->template getElements<D>();
289 for (
unsigned currentElement = 0; currentElement < elements.size();
291 VectorType<T, D> nodes[
D];
292 VectorType<T, D> center;
293 std::fill(center.begin(), center.end(), 0);
295 std::bitset<2 * D> flags;
297 bool removeElement =
false;
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];
309 if (nodes[q][dim] > gridMin[dim])
311 if (nodes[q][dim] < gridMax[dim])
312 flags.reset(dim +
D);
314 center[dim] += nodes[q][dim];
316 if (removeBoundaries[dim] && (flags[dim] || flags[dim +
D])) {
317 removeElement =
true;
325 for (
int q = 0; q <
D; q++)
326 center[q] /=
static_cast<T>(
D);
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]);
341 viennahrle::Index<D> minIndex, maxIndex;
342 for (
int q = 0; q <
D; q++) {
344 static_cast<hrleIndexType
>(std::ceil(minNode[q] / gridDelta));
346 static_cast<hrleIndexType
>(std::floor(maxNode[q] / gridDelta));
351 for (
unsigned i = 0; i <
D; ++i) {
352 if (removeBoundaries[i]) {
353 minIndex[i] = std::max(minIndex[i], grid.getMinIndex(i));
355 maxIndex[i] = std::min(maxIndex[i], grid.getMaxGridPoint(i));
360 for (
int z = 0; z <
D; z++) {
361 VectorType<hrleIndexType,
D - 1> min_bb, max_bb;
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];
368 box bb(min_bb, max_bb);
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];
381 for (
int k = 1; k <
D; k++)
382 p[(k + z) %
D] = grid.gridPositionOfGlobalIndex(
383 (k + z) %
D, it_b[(k + z) %
D]);
386 int intersection_status = calculateGridlineTriangleIntersection(
387 p, nodes, z, intersection);
389 if (intersection_status != 0) {
391 assert(intersection >= minNode[z] &&
392 "Intersection should be inside Domain!");
393 assert(intersection <= maxNode[z] &&
394 "Intersection should be inside Domain!");
398 if (removeBoundaries[z] &&
399 intersection > gridDelta * (grid.getMaxIndex(z)))
401 if (removeBoundaries[z] &&
402 intersection < grid.getMinLocalCoordinate(z))
406 grid.globalCoordinate2LocalIndex(z, intersection);
408 auto floor =
static_cast<hrleIndexType
>(
409 std::floor(intersection2 - distanceEps));
410 auto ceil =
static_cast<hrleIndexType
>(
411 std::ceil(intersection2 + distanceEps));
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));
418 if (!removeBoundaries[z]) {
419 floor = grid.globalIndex2LocalIndex(z, floor);
420 ceil = grid.globalIndex2LocalIndex(z, ceil);
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];
429 for (it_b[z] = floor; it_b[z] <= ceil; ++(it_b[z])) {
431 T RealDistance = it_b[z] - intersection2;
433 T SignDistance = RealDistance;
435 SignDistance -= signEps * t[z];
437 if (intersection_status < 0) {
438 RealDistance = -RealDistance;
439 SignDistance = -SignDistance;
442 if (RealDistance < 0.)
443 SignDistance = -SignDistance;
445 RealDistance *= (1. - distanceEps * 1e-3);
446 if (RealDistance > 1.)
448 if (RealDistance < -1.)
451 if (RealDistance == 0.)
455 std::make_pair(grid.globalIndices2LocalIndices(it_b),
456 std::make_pair(SignDistance, RealDistance)));
463 std::sort(points.begin(), points.end());
466 auto it_points = points.begin();
467 while (it_points != points.end()) {
468 viennahrle::Index<D> tmp = it_points->first;
470 std::make_pair(it_points->first, it_points->second.second));
474 }
while ((it_points != points.end()) && (it_points->first == tmp));
478 levelSet->insertPoints(points2);
479 levelSet->getDomain().segment();
480 levelSet->finalize(2);