Skip to content

Commit

Permalink
New node checking
Browse files Browse the repository at this point in the history
  • Loading branch information
tobre1 committed Jan 9, 2025
1 parent be5fa3e commit 09556b8
Show file tree
Hide file tree
Showing 3 changed files with 150 additions and 62 deletions.
117 changes: 55 additions & 62 deletions include/viennals/lsToSurfaceMeshRefined.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ using namespace viennacore;

template <class LsNT, class MeshNT = LsNT, int D = 3>
class ToSurfaceMeshRefined {

typedef Domain<LsNT, D> lsDomainType;
typedef typename Domain<LsNT, D>::DomainType hrleDomainType;
typedef KDTree<LsNT, std::array<LsNT, 3>> kdTreeType;
Expand All @@ -29,9 +30,34 @@ class ToSurfaceMeshRefined {
SmartPointer<kdTreeType> kdTree = nullptr;

const MeshNT epsilon;
const MeshNT minNodeDistanceFactor = 0.2;

bool checkNodeDistance = true;

struct compareNodes {
bool operator()(const Vec3D<MeshNT> &nodeA,
const Vec3D<MeshNT> &nodeB) const {
if (nodeA[0] < nodeB[0] - minNodeDistance)
return true;
if (nodeA[0] > nodeB[0] + minNodeDistance)
return false;

if (nodeA[1] < nodeB[1] - minNodeDistance)
return true;
if (nodeA[1] > nodeB[1] + minNodeDistance)
return false;

if (nodeA[2] < nodeB[2] - minNodeDistance)
return true;
if (nodeA[2] > nodeB[2] + minNodeDistance)
return false;

return false;
}
};

public:
static MeshNT minNodeDistance;

ToSurfaceMeshRefined(double eps = 1e-12) : epsilon(eps) {}

ToSurfaceMeshRefined(SmartPointer<lsDomainType> passedLevelSet,
Expand All @@ -56,9 +82,11 @@ class ToSurfaceMeshRefined {
kdTree = passedKdTree;
}

void setMinNodeDistanceFactor(MeshNT factor) {
minNodeDistanceFactor = factor;
}
// void setMinNodeDistanceFactor(MeshNT factor) {
// minNodeDistanceFactor = factor;
// }

void setCheckNodeDistance(bool check) { checkNodeDistance = check; }

void apply() {
if (levelSet == nullptr) {
Expand All @@ -76,7 +104,6 @@ class ToSurfaceMeshRefined {

mesh->clear();
const auto gridDelta = levelSet->getGrid().getGridDelta();
const MeshNT minNodeDistance = gridDelta * minNodeDistanceFactor;
mesh->minimumExtent = Vec3D<MeshNT>{std::numeric_limits<MeshNT>::max(),
std::numeric_limits<MeshNT>::max(),
std::numeric_limits<MeshNT>::max()};
Expand All @@ -100,13 +127,16 @@ class ToSurfaceMeshRefined {
nodeContainerType;

nodeContainerType nodes[D];
minNodeDistance = gridDelta * 0.01;
std::map<Vec3D<MeshNT>, unsigned, compareNodes> nodeCoordinates;

typename nodeContainerType::iterator nodeIt;
lsInternal::MarchingCubes marchingCubes;

std::vector<Vec3D<MeshNT>> triangleCenters;
std::vector<Vec3D<MeshNT>> normals;
const bool buildKdTreeFlag = kdTree != nullptr;
const bool checkNodeFlag = checkNodeDistance;

// iterate over all active surface points
for (hrleConstSparseCellIterator<hrleDomainType> cellIt(
Expand Down Expand Up @@ -156,7 +186,8 @@ class ToSurfaceMeshRefined {
nodeIt = nodes[dir].find(d);
if (nodeIt != nodes[dir].end()) {
nod_numbers[n] = nodeIt->second;
} else { // if node does not exist yet
} else {
// if node does not exist yet
// calculate coordinate of new node
Vec3D<MeshNT> cc{}; // initialise with zeros
for (int z = 0; z < D; z++) {
Expand All @@ -167,10 +198,10 @@ class ToSurfaceMeshRefined {
cellIt.getIndices(z) +
BitMaskToVector<D, hrleIndexType>(p0)[z]);
} else {
MeshNT d0, d1;

d0 = static_cast<MeshNT>(cellIt.getCorner(p0).getValue());
d1 = static_cast<MeshNT>(cellIt.getCorner(p1).getValue());
MeshNT d0 =
static_cast<MeshNT>(cellIt.getCorner(p0).getValue());
MeshNT d1 =
static_cast<MeshNT>(cellIt.getCorner(p1).getValue());

// calculate the surface-grid intersection point
if (d0 == -d1) { // includes case where d0=d1=0
Expand All @@ -190,14 +221,22 @@ class ToSurfaceMeshRefined {
cc[z] = gridDelta * cc[z];
}

int nodeIdx = checkIfNodeExists(cc, minNodeDistance);
int nodeIdx = -1;
if (checkNodeFlag) {
auto checkNode = nodeCoordinates.find(cc);
if (checkNode != nodeCoordinates.end()) {
nodeIdx = checkNode->second;
}
}
if (nodeIdx >= 0) {
// node exists or close node exists
nod_numbers[n] = nodeIdx;
} else {
// insert new node
nod_numbers[n] =
mesh->insertNextNode(cc); // insert new surface node
nodes[dir][d] = nod_numbers[n];
nodeCoordinates[cc] = nod_numbers[n];

for (int a = 0; a < D; a++) {
if (cc[a] < mesh->minimumExtent[a])
mesh->minimumExtent[a] = cc[a];
Expand All @@ -214,7 +253,7 @@ class ToSurfaceMeshRefined {
mesh->nodes[nod_numbers[2]]);
LsNT norm = std::sqrt(normal[0] * normal[0] + normal[1] * normal[1] +
normal[2] * normal[2]);
if (norm > gridDelta * gridDelta * 1e-4) {
if (norm > epsilon) {
mesh->insertNextElement(nod_numbers); // insert new surface element
for (int d = 0; d < D; d++) {
normal[d] /= norm;
Expand Down Expand Up @@ -249,55 +288,6 @@ class ToSurfaceMeshRefined {
}
}

int checkIfNodeExists(const Vec3D<MeshNT> &node,
const MeshNT minNodeDistance) {
const auto &nodes = mesh->getNodes();
const uint N = nodes.size();
const int maxThreads = omp_get_max_threads();
const int numThreads = maxThreads > N / (10 * maxThreads) ? 1 : maxThreads;
std::vector<int> threadLocal(numThreads, -1);
std::atomic<bool> foundPoint(false);
const uint numPointsPerThread = N / numThreads;

#pragma omp parallel shared(node, nodes, threadLocal) num_threads(numThreads)
{
int threadId = omp_get_thread_num();

uint i = threadId * numPointsPerThread;
const uint endPointIdx =
threadId == numThreads - 1 ? N : (threadId + 1) * numPointsPerThread;

while (i < endPointIdx && !foundPoint) {
if (nodeClose(node, nodes[i], minNodeDistance)) {
threadLocal[threadId] = i;
foundPoint = true;
}
i++;
}
}

int idx = -1;
for (size_t i = 0; i < threadLocal.size(); i++) {
if (threadLocal[i] >= 0) {
idx = threadLocal[i];
break;
}
}

return idx;
}

static bool nodeClose(const Vec3D<MeshNT> &nodeA, const Vec3D<MeshNT> &nodeB,
const MeshNT distance) {
const auto nodeDist = std::abs(nodeA[0] - nodeB[0]) +
std::abs(nodeA[1] - nodeB[1]) +
std::abs(nodeA[2] - nodeB[2]);
if (nodeDist < distance)
return true;

return false;
}

static bool inline triangleMisformed(
const std::array<unsigned, D> &nod_numbers) {
if constexpr (D == 3) {
Expand All @@ -321,4 +311,7 @@ class ToSurfaceMeshRefined {
}
};

template <class LsNT, class MeshNT, int D>
MeshNT ToSurfaceMeshRefined<LsNT, MeshNT, D>::minNodeDistance = 0.1;

} // namespace viennals
7 changes: 7 additions & 0 deletions tests/RefinedSurfaceMesh/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
project(RefinedSurfaceMesh LANGUAGES CXX)

add_executable(${PROJECT_NAME} "${PROJECT_NAME}.cpp")
target_link_libraries(${PROJECT_NAME} PRIVATE ViennaLS)

add_dependencies(ViennaLS_Tests ${PROJECT_NAME})
add_test(NAME ${PROJECT_NAME} COMMAND $<TARGET_FILE:${PROJECT_NAME}>)
88 changes: 88 additions & 0 deletions tests/RefinedSurfaceMesh/RefinedSurfaceMesh.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,88 @@
#include <chrono>
#include <iostream>

#include <lsBooleanOperation.hpp>
#include <lsDomain.hpp>
#include <lsExpand.hpp>
#include <lsFromSurfaceMesh.hpp>
#include <lsMakeGeometry.hpp>
#include <lsToSurfaceMeshRefined.hpp>
#include <lsToVoxelMesh.hpp>
#include <lsVTKWriter.hpp>

/**
Minimal example showing how to write different
meshes created by the algorithms lsToVoxelMesh and lsToSurfaceMesh.
\example Make3DSphere.cpp
*/

namespace ls = viennals;

int main() {

constexpr int D = 3;

omp_set_num_threads(1);

double gridDelta = 0.4;

auto sphere1 = ls::SmartPointer<ls::Domain<double, D>>::New(
gridDelta); //, boundaryCons);

auto sphere2 = ls::SmartPointer<ls::Domain<double, D>>::New(
gridDelta); //, boundaryCons);
double origin[3] = {5., 0., 0.};
double radius = 7.3;

ls::MakeGeometry<double, D>(
sphere1, ls::SmartPointer<ls::Sphere<double, D>>::New(origin, radius))
.apply();
origin[0] = -5.;
ls::MakeGeometry<double, D>(
sphere2, ls::SmartPointer<ls::Sphere<double, D>>::New(origin, radius))
.apply();
ls::BooleanOperation<double, D>(sphere1, sphere2,
ls::BooleanOperationEnum::UNION)
.apply();

std::cout << "Number of points: " << sphere1->getDomain().getNumberOfPoints()
<< std::endl;
auto mesh = ls::SmartPointer<ls::Mesh<>>::New();

ls::ToSurfaceMeshRefined<double, double, D>(sphere1, mesh).apply();
ls::VTKWriter<double>(mesh, "test-refined.vtp").apply();

std::cout << "Refined mesh written to test-refined.vtp" << std::endl;
std::cout << "Number of points: " << mesh->nodes.size() << std::endl;
std::cout << "Number of triangles: " << mesh->triangles.size() << std::endl;

auto &nodes = mesh->nodes;
auto minNodeDistance =
ls::ToSurfaceMeshRefined<double, double, D>::minNodeDistance;
for (unsigned i = 0; i < nodes.size(); ++i) {
for (unsigned j = i + 1; j < nodes.size(); ++j) {
double dist = 0.;
// manhattan distance
for (unsigned k = 0; k < D; ++k) {
dist += std::abs(nodes[i][k] - nodes[j][k]);
}
if (dist < minNodeDistance) {
std::cout << "Distance between nodes " << i << " and " << j
<< " is smaller than minNodeDistance: " << dist << std::endl;
}
}
}
std::cout << "Minimum node distance: " << minNodeDistance << std::endl;

ls::ToSurfaceMeshRefined<double, double, D> noCheck(sphere1, mesh);
noCheck.setCheckNodeDistance(false);
noCheck.apply();

std::cout << "Not refined mesh written to test-not-refined.vtp" << std::endl;
std::cout << "Number of points: " << mesh->nodes.size() << std::endl;
std::cout << "Number of triangles: " << mesh->triangles.size() << std::endl;

ls::VTKWriter<double>(mesh, "test-no-refined.vtp").apply();

return 0;
}

0 comments on commit 09556b8

Please sign in to comment.