WolframResearch / OpenCascadeLink

Open source package for OpenCascadeLink, which is bundled with Wolfram Language products as of version 12.1
Other
31 stars 8 forks source link

OpenCascadeShapeUnion/Intersection/Difference return wrong result in some cases. #5

Closed ruebenko closed 2 years ago

ruebenko commented 3 years ago

This returns a wrong result:

    Needs["OpenCascadeLink`"]
    Needs["NDSolve`FEM`"]
    shape = OpenCascadeShape[CapsuleShape[{{0, -2, 0}, {0, 0, 1}}, 1]];
    r = ImplicitRegion[x^4 + y^4 + z^4 - 1 == 0, {x, y, z}];
    bmesh = ToBoundaryMesh[r, {{-2, 2}, {-2, 2}, {-2, 2}}];
    shape2 = OpenCascadeShape[bmesh];
    union = OpenCascadeShapeUnion[shape2, shape];
    OpenCascadeShapeSurfaceMeshToBoundaryMesh[union]["Wireframe"]

but these

    r = ImplicitRegion[x^4 + y^4 + z^4 - 1 <= 0, {x, y, z}];
    r = ImplicitRegion[x^4 + y^4 + z^4 - 1 < 0, {x, y, z}];
    r = ImplicitRegion[x^4 + y^4 + z^4 - 1 > 0, {x, y, z}];
    r = ImplicitRegion[x^4 + y^4 + z^4 - 1 >= 0, {x, y, z}];

return the correct result. For the broken case either use

    r = ImplicitRegion[x^4 + y^4 + z^4 - 1 <= 0, {x, y, z}];

Or this workaround can be used:

mp = MeshPrimitives[DiscretizeRegion[r], 2];
shape2 = OpenCascadeShapeSewing[
   OpenCascadeShape[Polygon[#]] & /@ (Reverse /@ mp[[All, 1]])];
ralphrmartin commented 3 years ago

I dont think your assertion that r = ImplicitRegion[x^4 + y^4 + z^4 - 1 > 0, {x, y, z}]; r = ImplicitRegion[x^4 + y^4 + z^4 - 1 >= 0, {x, y, z}]; work correctly is true.

Unless I have made a mistake, using r = ImplicitRegion[x^4 + y^4 + z^4 - 1 > 0, {x, y, z}]; gives me exactly the same shape result as r = ImplicitRegion[x^4 + y^4 + z^4 - 1 < 0, {x, y, z}];

However, as the solid is on the opposite side in these two cases, I would NOT expect the results to be the same. In the second case the implicit region should be a cube with a hole inside it the same shape as the first case. When unioned with the capsule, we should see the flat sides of the cube.

ralphrmartin commented 3 years ago

Adding to what I said, the problem appears to be a proper notion of solidity.

When you specify "<", "<=" or ">", ">=", it is obvious which side of the surface you consider to be solid, and hence which way the normal should face. (I presume OpenCascade has a well defined convention that normals always point out of (or into) the solid). However, when you create the region "==0" this is a 2D surface, and it is quite arbitrary which side the normal should face. The only person who can tell is the user, in general, so the user needs to be able to specify. The same is true for any 2D surface the user creates using Wolfram Language.

ruebenko commented 3 years ago

I get the same result for <= and < and the same result for > and >=. Unfortunately, github will not let me upload images. Please check again. (Perhaps you need to revert the previous "fix"?)

Conning the == case, I think this is true if the resulting surface is open, but in this case it is closed and it is possible to figure this out and have the proper normal set. While strictly speaking for the == case a closed surface could have arbitrary normal (then this would not be a bug) there is a strong argument for doing the 'right thing" in that the volume enclosed is the inside.

ruebenko commented 3 years ago

I'd like to point you to Mathematica Stackexchange that site is much better suited to ask questions and isolating issues than this github issues thing. I am also over there and will see what you possibly post.

ralphrmartin commented 3 years ago

My apologies. You are right that both <= and < and > and >= do work correctly. It seems that previously, my kernel stopped unexpectedly, and I looked at a stale result.

I completely agree with you that if the surface is simple and closed, there is a "natural" orientation for the == case. However, there does need to be a way for the user to control what happens in the open surface case. Furthermore, there are more tricky cases too, like a product surface comprising one sphere inside another, or self-intersecting surfaces, where what is correct is less obvious (winding numbers may need considering, etc).

Sorry, but I do not agree about StackExchange as it has no way of finding out when an issue has been resolved and closed.

ruebenko commented 2 years ago

This is fixed with the latest update. I have also added a new function OpenCascdeShapeFix that does shape healing and, among other things, can fix the orientation of wrongly oriented surfaces.