#528 causes invalid UpdateVert calls
#529
Comments
|
I ported it to C++ and reduced it to triangulator error: It looks something like this: Code for the sample: #include "samples.h"
namespace {
using namespace manifold;
using namespace glm;
constexpr float AtomicRadiusN2 = 0.65;
constexpr float BondPairN2 = 1.197;
constexpr float AtomicRadiusSi = 1.1;
constexpr float LatticeCellSizeSi = 5.4309;
constexpr float fccOffset = 0.25;
constexpr float AtomicRadiusC = 0.7;
constexpr float LatticeCellSizeC = 3.65;
constexpr float cellLenA = 2.464;
constexpr float cellLenB = cellLenA;
constexpr float cellLenC = 6.711;
constexpr float cellAngleA = 90;
constexpr float cellAngleB = cellAngleA;
constexpr float cellAngleC = 120;
constexpr float LayerSeperationC = 3.364;
Manifold bond(int fn = 16, vec3 p1 = {0, 0, 0}, vec3 p2 = {1, 1, 1},
float ar1 = 1.0, float ar2 = 2.0) {
float cyR = std::min(ar1, ar2) / 5.0;
float dist = length(p1 - p2);
vec3 cyC = (p1 + p2) / 2.0f;
float beta = degrees(acos((p1.z - p2.z) / dist));
float gamma = degrees(atan2(p1.y - p2.y, p1.x - p2.x));
vec3 rot = {0.0, beta, gamma};
return Manifold::Cylinder(dist, cyR, -1, fn, true)
.Rotate(rot.x, rot.y, rot.z)
.Translate(cyC);
}
Manifold bondPair(int fn = 16, float d = 0.0, float ar = 1.0) {
float axD = pow(d, 1.0 / 3.0);
vec3 p1 = {+axD, -axD, -axD};
vec3 p2 = {-axD, +axD, +axD};
Manifold sphere = Manifold::Sphere(ar, fn);
return sphere.Translate(p1) + sphere.Translate(p2) + bond(fn, p1, p2, ar, ar);
}
Manifold hexagonalClosePacked(int fn = 16, vec3 dst = {1.0, 1.0, 1.0},
float ar = 1.0) {
std::vector<Manifold> parts;
vec3 p1 = {0, 0, 0};
parts.push_back(Manifold::Sphere(ar, fn));
float baseAg = 30;
vec3 ag = {baseAg, baseAg + 120, baseAg + 240};
vec3 points[] = {{cosd(ag.x) * dst.x, sind(ag.x) * dst.x, 0},
{cosd(ag.y) * dst.y, sind(ag.y) * dst.y, 0},
{cosd(ag.z) * dst.z, sind(ag.z) * dst.z, 0}};
for (vec3 p2 : points) {
parts.push_back(Manifold::Sphere(ar, fn).Translate(p2));
parts.push_back(bond(fn, p1, p2, ar, ar));
}
return Manifold::BatchBoolean(parts, OpType::Add);
}
Manifold fccDiamond(int fn = 16, float ar = 1.0, float unitCell = 2.0,
float fccOffset = 0.25) {
std::vector<Manifold> parts;
float huc = unitCell / 2.0;
float od = fccOffset * unitCell;
vec3 interstitial[] = {
{+od, +od, +od}, {+od, -od, -od}, {-od, +od, -od}, {-od, -od, +od}};
vec3 corners[] = {{+huc, +huc, +huc},
{+huc, -huc, -huc},
{-huc, +huc, -huc},
{-huc, -huc, +huc}};
vec3 fcc[] = {{+huc, 0, 0}, {-huc, 0, 0}, {0, +huc, 0},
{0, -huc, 0}, {0, 0, +huc}, {0, 0, -huc}};
for (auto p : corners) parts.push_back(Manifold::Sphere(ar, fn).Translate(p));
for (auto p : fcc) parts.push_back(Manifold::Sphere(ar, fn).Translate(p));
for (auto p : interstitial)
parts.push_back(Manifold::Sphere(ar, fn).Translate(p));
vec3 bonds[][2] = {{interstitial[0], corners[0]}, {interstitial[0], fcc[0]},
{interstitial[0], fcc[2]}, {interstitial[0], fcc[4]},
{interstitial[1], corners[1]}, {interstitial[1], fcc[0]},
{interstitial[1], fcc[3]}, {interstitial[1], fcc[5]},
{interstitial[2], corners[2]}, {interstitial[2], fcc[1]},
{interstitial[2], fcc[2]}, {interstitial[2], fcc[5]},
{interstitial[3], corners[3]}, {interstitial[3], fcc[1]},
{interstitial[3], fcc[3]}, {interstitial[3], fcc[4]}};
for (auto b : bonds) parts.push_back(bond(fn, b[0], b[1], ar, ar));
return Manifold::BatchBoolean(parts, OpType::Add);
}
Manifold SiCell(int fn = 16, float x = 1.0, float y = 1.0, float z = 1.0) {
return fccDiamond(fn, AtomicRadiusSi, LatticeCellSizeSi, fccOffset)
.Translate({LatticeCellSizeSi * x, LatticeCellSizeSi * y,
LatticeCellSizeSi * z});
}
Manifold SiN2Cell(int fn = 16, float x = 1.0, float y = 1.0, float z = 1.0) {
float n2Offset = LatticeCellSizeSi / 8;
return bondPair(fn, BondPairN2, AtomicRadiusN2)
.Translate({LatticeCellSizeSi * x - n2Offset,
LatticeCellSizeSi * y + n2Offset,
LatticeCellSizeSi * z + n2Offset}) +
SiCell(fn, x, y, z);
}
Manifold GraphiteCell(int fn, vec3 xyz = {1.0, 1.0, 1.0}) {
vec3 loc = {(cellLenA * xyz.x * cosd(30) * 2),
((cellLenB * sind(30)) + cellLenC) * xyz.y, xyz.z};
return hexagonalClosePacked(fn, {cellLenA, cellLenB, cellLenC}, AtomicRadiusC)
.Translate(loc);
}
} // namespace
namespace manifold {
Manifold CondensedMatter(int fn) {
std::vector<Manifold> parts;
float siOffset = 3.0 * LatticeCellSizeSi / 8.0;
for (int x = -3; x <= 3; x++)
for (int y = -1; y <= 2; y++)
parts.push_back(
GraphiteCell(fn, {x + (y % 2 == 0 ? 0.0 : 0.5), y,
LayerSeperationC * 0.5 + LatticeCellSizeSi * 1.5})
.Translate({0, -siOffset, 0})
.Rotate(0, 0, 45));
float xyPlane[] = {-2, -1, 0, +1, +2};
for (float x : xyPlane)
for (float y : xyPlane) parts.push_back(SiN2Cell(fn, x, y, 1));
return Manifold::BatchBoolean(parts, OpType::Add);
}
} // namespace manifold |
|
Okay, so two different bugs? Looks like the more serious is that |
|
Yeah the segfault behavior is a very serious bug. I was not sure if the invalid triangulation can still preserve manifoldness so I did not make it a separate issue. |
|
I think I should add bound check for VecDH by default, it will slow down things a tiny bit but it is worth paying for the price. |
|
Yeah, probably an assert so it only compiles for debug mode, right? Because we're not going to be able to recover anyway. |
|
I can repro your triangulator error with |
|
I actually forgot how I got the segfault, I think it is related to compiler version (the undefined behavior may not cause segfault in some cases, but it is bad enough....). I do get address sanitizer errors. For the issue with Unique, it is because thrust Unique is buggy and can access invalid memory. That is the main reason I started switching the algorithms to pstl, but sadly OSX does not support pstl for now so we still default to thrust implementation in that case. |
|
Ah, I think I only got the access error when using openscad, the index was -1. Backtrace: |
|
Any chance you can write a test that'll fail this way on our CI? |
|
I think I mentioned this before, but this is another reason I'd like to switch SparseIndices from relying on sort and unique to just using a hash table. It should be faster, simpler, and probably allow us to remove these weird thrust calls entirely. |
Yeah I can try to do that tmr. I think this is related to how their spheres and cylinders are constructed, as they are using constructors in openscad instead of ours so there will be subtle differences. |
|
Ah yes, that makes sense. Now that we have a way to import 3D model files into our tests, hopefully you can just export the models from OpenSCAD and use them directly. |
|
while modifying the examples I found a bug: 582c85a originally this can evaluate to 0 when fn is not specified, and causes the triangulator to crash. |
|
No, I am unable to reproduce the problem in our tests... |
|
Hmm, that's troubling... Well, I'm working on the triangulator refactor now, though that'll likely take some time. I can't make much progress on the decimator until I have a repro - maybe we'll luck out and it'll be randomly fixed by another merge. |
|
OK I got a simpler reproduce, but it is non-deterministic... --- a/src/utilities/include/vec.h
+++ b/src/utilities/include/vec.h
@@ -59,12 +59,20 @@ class VecView {
operator VecView<const T>() const { return {ptr_, size_}; }
inline const T &operator[](int i) const {
- if (i < 0 || i >= size_) throw std::out_of_range("Vec out of range");
+ if (i < 0 || i >= size_) {
+ printf("i: %d, size: %d\n", i, size_);
+ __builtin_trap();
+ throw std::out_of_range("Vec out of range");
+ }
return ptr_[i];
}
inline T &operator[](int i) {
- if (i < 0 || i >= size_) throw std::out_of_range("Vec out of range");
+ if (i < 0 || i >= size_) {
+ printf("i: %d, size: %d\n", i, size_);
+ __builtin_trap();
+ throw std::out_of_range("Vec out of range");
+ }
return ptr_[i];
}The out of range are always -1 or -3, so I think we should probably add more check to make sure that the halfedge is not already removed... Also, I implemented multithreaded CollapseEdge, but the performance is not great: the performance is slightly better on my machine with 20 cores (50ms faster on average? over something that took 2400ms), and slower when the number of cores decreases. I double checked to make sure that my partition implementation should be pretty fast, so I guess the problem is just too hard to parallelize (at least for now). |
|
I ran I am seeing the test fail hard with address sanitizer enabled, but at |
|
Interesting, I guess the problem with concurrency issue is that it can be hard to reproduce... I actually tried to switch SparseIndices to a hash table, at least in Kernel12. The performance for Kernel12 is now much faster, but we first have to build the hash table and the build is pretty slow (using phmap, which is faster than std::unordered_map alreadh). There are several reasons that forbid using hash table everywhere:
|
|
Agreed we should use our existing parallel hashmap. Does it need a refactor? Despite all the code, what |
|
And for unique, it might be simpler to just iterate through the hashmap and add the |
|
Indeed, that should work. I was thinking about using pair as the key type but it is not really needed. |
In case this is helpful: |
|
Oh it seems that this is triggered by triangulation failure, fixing the triangulator should hopefully fix this... (reverting 8011fb6 fixes the issue) |
|
I sometimes got this as well: |
|
Thanks, that looks actionable! I'll take a look. |
|
I am now pretty sure that:
diff --git a/src/manifold/src/edge_op.cpp b/src/manifold/src/edge_op.cpp
index 4e3db36..9f69dff 100644
--- a/src/manifold/src/edge_op.cpp
+++ b/src/manifold/src/edge_op.cpp
@@ -158,7 +158,7 @@ void Manifold::Impl::SimplifyTopology() {
entries[i].index = i;
});
- sort(policy, entries.begin(), entries.end());
+ stable_sort(policy, entries.begin(), entries.end());
for (int i = 0; i < nbEdges - 1; ++i) {
if (entries[i].start == entries[i + 1].start &&
entries[i].end == entries[i + 1].end) {
diff --git a/src/manifold/src/sort.cpp b/src/manifold/src/sort.cpp
index 13deab4..007b305 100644
--- a/src/manifold/src/sort.cpp
+++ b/src/manifold/src/sort.cpp
@@ -261,6 +261,7 @@ void Manifold::Impl::Finish() {
if (halfedge_.size() == 0) return;
CompactProps();
+ if (halfedge_.size() % 6 != 0) __builtin_trap();
ASSERT(halfedge_.size() % 6 == 0, topologyErr,
"Not an even number of faces after sorting faces!");
diff --git a/src/manifold/src/edge_op.cpp b/src/manifold/src/edge_op.cpp
index 4e3db36..9161fe2 100644
--- a/src/manifold/src/edge_op.cpp
+++ b/src/manifold/src/edge_op.cpp
@@ -343,6 +351,8 @@ void Manifold::Impl::FormLoop(int current, int end) {
}
void Manifold::Impl::CollapseTri(const glm::ivec3& triEdge) {
+ if (halfedge_[triEdge[1]].pairedHalfedge == -1)
+ return;
int pair1 = halfedge_[triEdge[1]].pairedHalfedge;
int pair2 = halfedge_[triEdge[2]].pairedHalfedge;
halfedge_[pair1].pairedHalfedge = pair2;
@@ -355,6 +365,8 @@ void Manifold::Impl::CollapseTri(const glm::ivec3& triEdge) {
void Manifold::Impl::RemoveIfFolded(int edge) {
const glm::ivec3 tri0edge = TriOf(edge);
const glm::ivec3 tri1edge = TriOf(halfedge_[edge].pairedHalfedge);
+ if (halfedge_[tri0edge[1]].pairedHalfedge == -1)
+ return;
if (halfedge_[tri0edge[1]].endVert == halfedge_[tri1edge[1]].endVert) {
if (halfedge_[tri0edge[1]].pairedHalfedge == tri1edge[2]) {
if (halfedge_[tri0edge[2]].pairedHalfedge == tri1edge[1]) { |
|
Alternatively we can use sort, but only dedupe the smaller indices: diff --git a/src/manifold/src/edge_op.cpp b/src/manifold/src/edge_op.cpp
index 4e3db36..cae7f74 100644
--- a/src/manifold/src/edge_op.cpp
+++ b/src/manifold/src/edge_op.cpp
@@ -162,7 +162,8 @@ void Manifold::Impl::SimplifyTopology() {
for (int i = 0; i < nbEdges - 1; ++i) {
if (entries[i].start == entries[i + 1].start &&
entries[i].end == entries[i + 1].end) {
- DedupeEdge(entries[i].index);
+ DedupeEdge(std::min(entries[i].index, entries[i + 1].index));
+ entries[i + 1].index = std::max(entries[i].index, entries[i + 1].index);
numFlagged++;
}
}And I can confirm that 1, 2 and 3 are all independent from each other. Fixing one will not fix another. I feel like I am starting to use my brain now. |
|
I like your solution - let's get that into a PR, close this confusing thread and start some new issues for the remaining problems independently. |
|
And extra points if you can manage to add a test first that allows our CI to catch this problem(s) and verify a fix. |
For 2, one way would be to sort with descending ID in order to verify that ID ordering is indeed important. But I don't think it would be very helpful to put this to the CI. For 3, I don't think it is possible to reproduce in the CI (reliably and easily). The only way I can think of is to record a trace, but this can only let us reproduce a specific instance, and cannot prove that our fix indeed fixes the problem. |
|
Yeah, I was afraid of that, but figured I'd ask 😄 |
Not sure why, but it seems that when openscad is compiled with the current master after #528, https://gist.github.com/ochafik/70a6b15e982b7ccd5a79ff9afd99dbcf#file-condensed-matter-scad will break and access invalid memory in
UpdateVertcall:I tried compiling with #525 and it works fine. Maybe we should try to port the example to c++ and debug. We can use address sanitizer to try to debug this kind of issues.
There is still invalid access when parallelization is turned off.
The text was updated successfully, but these errors were encountered: