Inconsistent mesh generation from BoundaryDiscretizeRegion and ToBoundaryMesh at creases

The name of the pictureThe name of the pictureThe name of the pictureClash Royale CLAN TAG#URR8PPP











up vote
5
down vote

favorite
1












I'm trying to use BoundaryDiscretizeRegion and/or ToBoundaryMesh to generate polygonal models for 3d printing, but I'm having trouble getting consistently high-quality results. Here's a simple example:



r = 1;
t = 0.2;
h = 0.3;
shape1 = ImplicitRegion[(z < t && Norm[x, y] <= r) ||
r - t <= Norm[x, y] <= r, x, -r, r, y, -r, r, z, 0, 1];
shell = ImplicitRegion[
 r - t <= Norm[x, y] <= r, x, -r, r, y, -r, r, z, t, 1];
base = ImplicitRegion[
 Norm[x, y] <= r, x, -r, r, y, -r, r, z, 0, t];
hole = ImplicitRegion[
 Norm[x, y] <= h, x, -r, r, y, -r, r, z, 0, t];
shape2 = RegionDifference[base, hole];
shape3 = RegionUnion[shell, shape2];
opts = Sequence[SphericalRegion -> True, ImageSize -> Medium,
 Axes -> True, Boxed -> True];
Map[BoundaryDiscretizeRegion[#, opts] &, shape1, shape2, shape3]


inconsistent mesh quality



BoundaryMeshRegion has generated nice sharp edges without artifacts between the surfaces in most places, except at the inside corner between the base and shell in shape1 and shape3, and around the rim of the hole in shape3. Surprisingly, the rim of the same hole in shape2 is sharp with no artifacts.



How can I get nice sharp corners in those places as well? Note that techniques specific to these particular cylindrical shapes may be of some interest, but I want something that will generalize to more complex shapes.



Some things I have tried:



  • Decreasing MaxCellMeasure. This makes the artifacts smaller, but does not eliminate them, and cannot make them small enough to be unnoticeable before the kernel is unable to process the request.

  • Adjusting AccuracyGoal and PrecisionGoal. These seemed to have no effect.

  • The BoundaryGroupingThreshold option of ToBoundaryMesh sounded promising, but it didn't seem to do the trick, although the documentation is a bit vague so I'm not sure I was using it correctly.

  • Generating meshes for each of the solid pieces (e.g. shell, base) individually then using RegionUnion to join them. RegionUnion was not able to join the regions.

  • I tried a completely different approach of generating individual regions for each of the surfaces - outside of shell, inside of shell, rim, top of base, bottom of base - individually and then joining them. The resulting meshes had holes, incorrectly oriented faces, overlaps,and other problems that required cleanup, which wasn't even possible in some cases. It's also a very tedious approach.

EDIT: following added as a follow-up to @kh40tika suggestion to express result as a boolean combination of convex regions, and to discretize before doing the boolean operations. It's an interesting approach, but I'm having trouble applying it to the general case that I'm interested in.



My first step was to generalize Cylinder to an equivalent ImplicitRegion:



r = 1;
t = 0.2;
rdiff = RegionDifference;
rplus = RegionUnion;
(*vertcyl[rad_,height_,offset_]:=Cylinder[0,0,offset,0,0,offset+height,rad]*)
vertcyl[rad_, height_, offset_] := ImplicitRegion[
 Norm[x, y] <= rad,
 x, -rad, rad, y, -rad, rad, z, offset, offset + height]
shape1outer = vertcyl[r, 1, 0];
shape1inner = vertcyl[r - t, 1, t];

shape1outer, shape1inner = DiscretizeRegion /@ shape1outer, shape1inner;

shape1 = rdiff[shape1outer, shape1inner];


However that fails with this message:



BoundaryMeshRegion::bsuncl: The boundary surface is not closed because the edges Line[2802,2799,2240,1811,2328,2315,781,778,1591,1589,1734,1736,1237,1159,1196,1183,1179,1166,3874,3873,985,978,3898,3896,1966,1967,2789,2798,1964,1965,214,10,386,387,2769,2767,1589,1579,380,378,4734,4732,733,734,218,3690,416,786,2778,2779,<<1>>,5656,5658,10,11,2044,398,1733,1554,5658,5659,930,931,453,1309,5644,27,3260,3648,398,399,1806,919,4732,4733,3,1,411,725,687,669,989,988,5456,5458,2798,2797,3654,3651,2320,2321,794,795,50,36,1726,1729,706,724,<<825>>] only come from a single face.


So it appears that DiscretizeRegion may have some special knowledge of Cylinder? Is there a way to make your approach work using something more general than Cylinder?



Another problem I see is that in the general case I don't think it will be possible to use convex regions to generate the result. Here's an example of a more general case I'm interested in:



r[z_] := 0.6 + 0.4 z^2
t = 0.2;
h = 0.3; 
shape1 = ImplicitRegion[
 (z < t && Norm[x, y] <= r[z]) || r[z] - t <= Norm[x, y] <= r[z],
 x, -1, 1, y, -1, 1, z, 0, 1];
shell = ImplicitRegion[
 r[z] - t <= Norm[x, y] <= r[z],
 x, -1, 1, y, -1, 1, z, t, 1];
base = ImplicitRegion[
 Norm[x, y] <= r[z],
 x, -1, 1, y, -1, 1, z, 0, t];
hole = ImplicitRegion[
 Norm[x, y] <= h,
 x, -1, 1, y, -1, 1, z, 0, t];
shape2 = RegionDifference[base, hole];
shape3 = RegionUnion[shell, shape2];
opts = Sequence[SphericalRegion -> True, ImageSize -> Medium, Axes -> True, Boxed -> True];
Map[BoundaryDiscretizeRegion[#, opts] &, shape1, shape2, shape3]


enter image description here



I don't think it's possible to represent this as a boolean combination of convex regions, is it?



Here's another example. This one follows the proposed rules (convex region, discretize before doing boolean operations), but the generated mesh is not what I'm looking for:



DiscretizeRegion[
 ImplicitRegion[
 Norm[x, y] <= 1 && z <= 1 + 0.1 x + 0.1 y,
 x, -1, 1, y, -1, 1, z, 0, 1]]


bad mesh for convex region




mesh 3d-printing




share|improve this question






















  • Concerning the BoundaryGroupingThreshold option: this is to attribute markers to groups of contiguous boundary segments/faces. The threshold does not change the form of a boundary; it changes what ToBoundayMesh considers as segments/faces that point in the same direction. I tried to clarify this in the documentation.
    – user21
    Sep 3 at 5:43










  • Thanks @user21, that helps clarify. Do you have any insight to the more general problem?
    – Bruce Lucas
    Sep 3 at 12:54










  • I'd report this to the support.
    – user21
    Sep 4 at 6:14










  • Thanks, I will do that. I suspect it's a limitation of the algorithm, although I think it should be possible.
    – Bruce Lucas
    Sep 5 at 0:41














up vote
5
down vote

favorite
1












I'm trying to use BoundaryDiscretizeRegion and/or ToBoundaryMesh to generate polygonal models for 3d printing, but I'm having trouble getting consistently high-quality results. Here's a simple example:



r = 1;
t = 0.2;
h = 0.3;
shape1 = ImplicitRegion[(z < t && Norm[x, y] <= r) ||
r - t <= Norm[x, y] <= r, x, -r, r, y, -r, r, z, 0, 1];
shell = ImplicitRegion[
 r - t <= Norm[x, y] <= r, x, -r, r, y, -r, r, z, t, 1];
base = ImplicitRegion[
 Norm[x, y] <= r, x, -r, r, y, -r, r, z, 0, t];
hole = ImplicitRegion[
 Norm[x, y] <= h, x, -r, r, y, -r, r, z, 0, t];
shape2 = RegionDifference[base, hole];
shape3 = RegionUnion[shell, shape2];
opts = Sequence[SphericalRegion -> True, ImageSize -> Medium,
 Axes -> True, Boxed -> True];
Map[BoundaryDiscretizeRegion[#, opts] &, shape1, shape2, shape3]


inconsistent mesh quality



BoundaryMeshRegion has generated nice sharp edges without artifacts between the surfaces in most places, except at the inside corner between the base and shell in shape1 and shape3, and around the rim of the hole in shape3. Surprisingly, the rim of the same hole in shape2 is sharp with no artifacts.



How can I get nice sharp corners in those places as well? Note that techniques specific to these particular cylindrical shapes may be of some interest, but I want something that will generalize to more complex shapes.



Some things I have tried:



  • Decreasing MaxCellMeasure. This makes the artifacts smaller, but does not eliminate them, and cannot make them small enough to be unnoticeable before the kernel is unable to process the request.

  • Adjusting AccuracyGoal and PrecisionGoal. These seemed to have no effect.

  • The BoundaryGroupingThreshold option of ToBoundaryMesh sounded promising, but it didn't seem to do the trick, although the documentation is a bit vague so I'm not sure I was using it correctly.

  • Generating meshes for each of the solid pieces (e.g. shell, base) individually then using RegionUnion to join them. RegionUnion was not able to join the regions.

  • I tried a completely different approach of generating individual regions for each of the surfaces - outside of shell, inside of shell, rim, top of base, bottom of base - individually and then joining them. The resulting meshes had holes, incorrectly oriented faces, overlaps,and other problems that required cleanup, which wasn't even possible in some cases. It's also a very tedious approach.

EDIT: following added as a follow-up to @kh40tika suggestion to express result as a boolean combination of convex regions, and to discretize before doing the boolean operations. It's an interesting approach, but I'm having trouble applying it to the general case that I'm interested in.



My first step was to generalize Cylinder to an equivalent ImplicitRegion:



r = 1;
t = 0.2;
rdiff = RegionDifference;
rplus = RegionUnion;
(*vertcyl[rad_,height_,offset_]:=Cylinder[0,0,offset,0,0,offset+height,rad]*)
vertcyl[rad_, height_, offset_] := ImplicitRegion[
 Norm[x, y] <= rad,
 x, -rad, rad, y, -rad, rad, z, offset, offset + height]
shape1outer = vertcyl[r, 1, 0];
shape1inner = vertcyl[r - t, 1, t];

shape1outer, shape1inner = DiscretizeRegion /@ shape1outer, shape1inner;

shape1 = rdiff[shape1outer, shape1inner];


However that fails with this message:



BoundaryMeshRegion::bsuncl: The boundary surface is not closed because the edges Line[2802,2799,2240,1811,2328,2315,781,778,1591,1589,1734,1736,1237,1159,1196,1183,1179,1166,3874,3873,985,978,3898,3896,1966,1967,2789,2798,1964,1965,214,10,386,387,2769,2767,1589,1579,380,378,4734,4732,733,734,218,3690,416,786,2778,2779,<<1>>,5656,5658,10,11,2044,398,1733,1554,5658,5659,930,931,453,1309,5644,27,3260,3648,398,399,1806,919,4732,4733,3,1,411,725,687,669,989,988,5456,5458,2798,2797,3654,3651,2320,2321,794,795,50,36,1726,1729,706,724,<<825>>] only come from a single face.


So it appears that DiscretizeRegion may have some special knowledge of Cylinder? Is there a way to make your approach work using something more general than Cylinder?



Another problem I see is that in the general case I don't think it will be possible to use convex regions to generate the result. Here's an example of a more general case I'm interested in:



r[z_] := 0.6 + 0.4 z^2
t = 0.2;
h = 0.3; 
shape1 = ImplicitRegion[
 (z < t && Norm[x, y] <= r[z]) || r[z] - t <= Norm[x, y] <= r[z],
 x, -1, 1, y, -1, 1, z, 0, 1];
shell = ImplicitRegion[
 r[z] - t <= Norm[x, y] <= r[z],
 x, -1, 1, y, -1, 1, z, t, 1];
base = ImplicitRegion[
 Norm[x, y] <= r[z],
 x, -1, 1, y, -1, 1, z, 0, t];
hole = ImplicitRegion[
 Norm[x, y] <= h,
 x, -1, 1, y, -1, 1, z, 0, t];
shape2 = RegionDifference[base, hole];
shape3 = RegionUnion[shell, shape2];
opts = Sequence[SphericalRegion -> True, ImageSize -> Medium, Axes -> True, Boxed -> True];
Map[BoundaryDiscretizeRegion[#, opts] &, shape1, shape2, shape3]


enter image description here



I don't think it's possible to represent this as a boolean combination of convex regions, is it?



Here's another example. This one follows the proposed rules (convex region, discretize before doing boolean operations), but the generated mesh is not what I'm looking for:



DiscretizeRegion[
 ImplicitRegion[
 Norm[x, y] <= 1 && z <= 1 + 0.1 x + 0.1 y,
 x, -1, 1, y, -1, 1, z, 0, 1]]


bad mesh for convex region




mesh 3d-printing




share|improve this question






















  • Concerning the BoundaryGroupingThreshold option: this is to attribute markers to groups of contiguous boundary segments/faces. The threshold does not change the form of a boundary; it changes what ToBoundayMesh considers as segments/faces that point in the same direction. I tried to clarify this in the documentation.
    – user21
    Sep 3 at 5:43










  • Thanks @user21, that helps clarify. Do you have any insight to the more general problem?
    – Bruce Lucas
    Sep 3 at 12:54










  • I'd report this to the support.
    – user21
    Sep 4 at 6:14










  • Thanks, I will do that. I suspect it's a limitation of the algorithm, although I think it should be possible.
    – Bruce Lucas
    Sep 5 at 0:41












up vote
5
down vote

favorite
1









up vote
5
down vote

favorite
1






1





I'm trying to use BoundaryDiscretizeRegion and/or ToBoundaryMesh to generate polygonal models for 3d printing, but I'm having trouble getting consistently high-quality results. Here's a simple example:



r = 1;
t = 0.2;
h = 0.3;
shape1 = ImplicitRegion[(z < t && Norm[x, y] <= r) ||
r - t <= Norm[x, y] <= r, x, -r, r, y, -r, r, z, 0, 1];
shell = ImplicitRegion[
 r - t <= Norm[x, y] <= r, x, -r, r, y, -r, r, z, t, 1];
base = ImplicitRegion[
 Norm[x, y] <= r, x, -r, r, y, -r, r, z, 0, t];
hole = ImplicitRegion[
 Norm[x, y] <= h, x, -r, r, y, -r, r, z, 0, t];
shape2 = RegionDifference[base, hole];
shape3 = RegionUnion[shell, shape2];
opts = Sequence[SphericalRegion -> True, ImageSize -> Medium,
 Axes -> True, Boxed -> True];
Map[BoundaryDiscretizeRegion[#, opts] &, shape1, shape2, shape3]


inconsistent mesh quality



BoundaryMeshRegion has generated nice sharp edges without artifacts between the surfaces in most places, except at the inside corner between the base and shell in shape1 and shape3, and around the rim of the hole in shape3. Surprisingly, the rim of the same hole in shape2 is sharp with no artifacts.



How can I get nice sharp corners in those places as well? Note that techniques specific to these particular cylindrical shapes may be of some interest, but I want something that will generalize to more complex shapes.



Some things I have tried:



  • Decreasing MaxCellMeasure. This makes the artifacts smaller, but does not eliminate them, and cannot make them small enough to be unnoticeable before the kernel is unable to process the request.

  • Adjusting AccuracyGoal and PrecisionGoal. These seemed to have no effect.

  • The BoundaryGroupingThreshold option of ToBoundaryMesh sounded promising, but it didn't seem to do the trick, although the documentation is a bit vague so I'm not sure I was using it correctly.

  • Generating meshes for each of the solid pieces (e.g. shell, base) individually then using RegionUnion to join them. RegionUnion was not able to join the regions.

  • I tried a completely different approach of generating individual regions for each of the surfaces - outside of shell, inside of shell, rim, top of base, bottom of base - individually and then joining them. The resulting meshes had holes, incorrectly oriented faces, overlaps,and other problems that required cleanup, which wasn't even possible in some cases. It's also a very tedious approach.

EDIT: following added as a follow-up to @kh40tika suggestion to express result as a boolean combination of convex regions, and to discretize before doing the boolean operations. It's an interesting approach, but I'm having trouble applying it to the general case that I'm interested in.



My first step was to generalize Cylinder to an equivalent ImplicitRegion:



r = 1;
t = 0.2;
rdiff = RegionDifference;
rplus = RegionUnion;
(*vertcyl[rad_,height_,offset_]:=Cylinder[0,0,offset,0,0,offset+height,rad]*)
vertcyl[rad_, height_, offset_] := ImplicitRegion[
 Norm[x, y] <= rad,
 x, -rad, rad, y, -rad, rad, z, offset, offset + height]
shape1outer = vertcyl[r, 1, 0];
shape1inner = vertcyl[r - t, 1, t];

shape1outer, shape1inner = DiscretizeRegion /@ shape1outer, shape1inner;

shape1 = rdiff[shape1outer, shape1inner];


However that fails with this message:



BoundaryMeshRegion::bsuncl: The boundary surface is not closed because the edges Line[2802,2799,2240,1811,2328,2315,781,778,1591,1589,1734,1736,1237,1159,1196,1183,1179,1166,3874,3873,985,978,3898,3896,1966,1967,2789,2798,1964,1965,214,10,386,387,2769,2767,1589,1579,380,378,4734,4732,733,734,218,3690,416,786,2778,2779,<<1>>,5656,5658,10,11,2044,398,1733,1554,5658,5659,930,931,453,1309,5644,27,3260,3648,398,399,1806,919,4732,4733,3,1,411,725,687,669,989,988,5456,5458,2798,2797,3654,3651,2320,2321,794,795,50,36,1726,1729,706,724,<<825>>] only come from a single face.


So it appears that DiscretizeRegion may have some special knowledge of Cylinder? Is there a way to make your approach work using something more general than Cylinder?



Another problem I see is that in the general case I don't think it will be possible to use convex regions to generate the result. Here's an example of a more general case I'm interested in:



r[z_] := 0.6 + 0.4 z^2
t = 0.2;
h = 0.3; 
shape1 = ImplicitRegion[
 (z < t && Norm[x, y] <= r[z]) || r[z] - t <= Norm[x, y] <= r[z],
 x, -1, 1, y, -1, 1, z, 0, 1];
shell = ImplicitRegion[
 r[z] - t <= Norm[x, y] <= r[z],
 x, -1, 1, y, -1, 1, z, t, 1];
base = ImplicitRegion[
 Norm[x, y] <= r[z],
 x, -1, 1, y, -1, 1, z, 0, t];
hole = ImplicitRegion[
 Norm[x, y] <= h,
 x, -1, 1, y, -1, 1, z, 0, t];
shape2 = RegionDifference[base, hole];
shape3 = RegionUnion[shell, shape2];
opts = Sequence[SphericalRegion -> True, ImageSize -> Medium, Axes -> True, Boxed -> True];
Map[BoundaryDiscretizeRegion[#, opts] &, shape1, shape2, shape3]


enter image description here



I don't think it's possible to represent this as a boolean combination of convex regions, is it?



Here's another example. This one follows the proposed rules (convex region, discretize before doing boolean operations), but the generated mesh is not what I'm looking for:



DiscretizeRegion[
 ImplicitRegion[
 Norm[x, y] <= 1 && z <= 1 + 0.1 x + 0.1 y,
 x, -1, 1, y, -1, 1, z, 0, 1]]


bad mesh for convex region




mesh 3d-printing




share|improve this question














I'm trying to use BoundaryDiscretizeRegion and/or ToBoundaryMesh to generate polygonal models for 3d printing, but I'm having trouble getting consistently high-quality results. Here's a simple example:



r = 1;
t = 0.2;
h = 0.3;
shape1 = ImplicitRegion[(z < t && Norm[x, y] <= r) ||
r - t <= Norm[x, y] <= r, x, -r, r, y, -r, r, z, 0, 1];
shell = ImplicitRegion[
 r - t <= Norm[x, y] <= r, x, -r, r, y, -r, r, z, t, 1];
base = ImplicitRegion[
 Norm[x, y] <= r, x, -r, r, y, -r, r, z, 0, t];
hole = ImplicitRegion[
 Norm[x, y] <= h, x, -r, r, y, -r, r, z, 0, t];
shape2 = RegionDifference[base, hole];
shape3 = RegionUnion[shell, shape2];
opts = Sequence[SphericalRegion -> True, ImageSize -> Medium,
 Axes -> True, Boxed -> True];
Map[BoundaryDiscretizeRegion[#, opts] &, shape1, shape2, shape3]


inconsistent mesh quality



BoundaryMeshRegion has generated nice sharp edges without artifacts between the surfaces in most places, except at the inside corner between the base and shell in shape1 and shape3, and around the rim of the hole in shape3. Surprisingly, the rim of the same hole in shape2 is sharp with no artifacts.



How can I get nice sharp corners in those places as well? Note that techniques specific to these particular cylindrical shapes may be of some interest, but I want something that will generalize to more complex shapes.



Some things I have tried:



  • Decreasing MaxCellMeasure. This makes the artifacts smaller, but does not eliminate them, and cannot make them small enough to be unnoticeable before the kernel is unable to process the request.

  • Adjusting AccuracyGoal and PrecisionGoal. These seemed to have no effect.

  • The BoundaryGroupingThreshold option of ToBoundaryMesh sounded promising, but it didn't seem to do the trick, although the documentation is a bit vague so I'm not sure I was using it correctly.

  • Generating meshes for each of the solid pieces (e.g. shell, base) individually then using RegionUnion to join them. RegionUnion was not able to join the regions.

  • I tried a completely different approach of generating individual regions for each of the surfaces - outside of shell, inside of shell, rim, top of base, bottom of base - individually and then joining them. The resulting meshes had holes, incorrectly oriented faces, overlaps,and other problems that required cleanup, which wasn't even possible in some cases. It's also a very tedious approach.

EDIT: following added as a follow-up to @kh40tika suggestion to express result as a boolean combination of convex regions, and to discretize before doing the boolean operations. It's an interesting approach, but I'm having trouble applying it to the general case that I'm interested in.



My first step was to generalize Cylinder to an equivalent ImplicitRegion:



r = 1;
t = 0.2;
rdiff = RegionDifference;
rplus = RegionUnion;
(*vertcyl[rad_,height_,offset_]:=Cylinder[0,0,offset,0,0,offset+height,rad]*)
vertcyl[rad_, height_, offset_] := ImplicitRegion[
 Norm[x, y] <= rad,
 x, -rad, rad, y, -rad, rad, z, offset, offset + height]
shape1outer = vertcyl[r, 1, 0];
shape1inner = vertcyl[r - t, 1, t];

shape1outer, shape1inner = DiscretizeRegion /@ shape1outer, shape1inner;

shape1 = rdiff[shape1outer, shape1inner];


However that fails with this message:



BoundaryMeshRegion::bsuncl: The boundary surface is not closed because the edges Line[2802,2799,2240,1811,2328,2315,781,778,1591,1589,1734,1736,1237,1159,1196,1183,1179,1166,3874,3873,985,978,3898,3896,1966,1967,2789,2798,1964,1965,214,10,386,387,2769,2767,1589,1579,380,378,4734,4732,733,734,218,3690,416,786,2778,2779,<<1>>,5656,5658,10,11,2044,398,1733,1554,5658,5659,930,931,453,1309,5644,27,3260,3648,398,399,1806,919,4732,4733,3,1,411,725,687,669,989,988,5456,5458,2798,2797,3654,3651,2320,2321,794,795,50,36,1726,1729,706,724,<<825>>] only come from a single face.


So it appears that DiscretizeRegion may have some special knowledge of Cylinder? Is there a way to make your approach work using something more general than Cylinder?



Another problem I see is that in the general case I don't think it will be possible to use convex regions to generate the result. Here's an example of a more general case I'm interested in:



r[z_] := 0.6 + 0.4 z^2
t = 0.2;
h = 0.3; 
shape1 = ImplicitRegion[
 (z < t && Norm[x, y] <= r[z]) || r[z] - t <= Norm[x, y] <= r[z],
 x, -1, 1, y, -1, 1, z, 0, 1];
shell = ImplicitRegion[
 r[z] - t <= Norm[x, y] <= r[z],
 x, -1, 1, y, -1, 1, z, t, 1];
base = ImplicitRegion[
 Norm[x, y] <= r[z],
 x, -1, 1, y, -1, 1, z, 0, t];
hole = ImplicitRegion[
 Norm[x, y] <= h,
 x, -1, 1, y, -1, 1, z, 0, t];
shape2 = RegionDifference[base, hole];
shape3 = RegionUnion[shell, shape2];
opts = Sequence[SphericalRegion -> True, ImageSize -> Medium, Axes -> True, Boxed -> True];
Map[BoundaryDiscretizeRegion[#, opts] &, shape1, shape2, shape3]


enter image description here



I don't think it's possible to represent this as a boolean combination of convex regions, is it?



Here's another example. This one follows the proposed rules (convex region, discretize before doing boolean operations), but the generated mesh is not what I'm looking for:



DiscretizeRegion[
 ImplicitRegion[
 Norm[x, y] <= 1 && z <= 1 + 0.1 x + 0.1 y,
 x, -1, 1, y, -1, 1, z, 0, 1]]


bad mesh for convex region





mesh 3d-printing





share|improve this question













share|improve this question




share|improve this question








edited Sep 3 at 17:23

























asked Sep 3 at 1:48









Bruce Lucas

1714




1714











  • Concerning the BoundaryGroupingThreshold option: this is to attribute markers to groups of contiguous boundary segments/faces. The threshold does not change the form of a boundary; it changes what ToBoundayMesh considers as segments/faces that point in the same direction. I tried to clarify this in the documentation.
    – user21
    Sep 3 at 5:43










  • Thanks @user21, that helps clarify. Do you have any insight to the more general problem?
    – Bruce Lucas
    Sep 3 at 12:54










  • I'd report this to the support.
    – user21
    Sep 4 at 6:14










  • Thanks, I will do that. I suspect it's a limitation of the algorithm, although I think it should be possible.
    – Bruce Lucas
    Sep 5 at 0:41
















  • Concerning the BoundaryGroupingThreshold option: this is to attribute markers to groups of contiguous boundary segments/faces. The threshold does not change the form of a boundary; it changes what ToBoundayMesh considers as segments/faces that point in the same direction. I tried to clarify this in the documentation.
    – user21
    Sep 3 at 5:43










  • Thanks @user21, that helps clarify. Do you have any insight to the more general problem?
    – Bruce Lucas
    Sep 3 at 12:54










  • I'd report this to the support.
    – user21
    Sep 4 at 6:14










  • Thanks, I will do that. I suspect it's a limitation of the algorithm, although I think it should be possible.
    – Bruce Lucas
    Sep 5 at 0:41















Concerning the BoundaryGroupingThreshold option: this is to attribute markers to groups of contiguous boundary segments/faces. The threshold does not change the form of a boundary; it changes what ToBoundayMesh considers as segments/faces that point in the same direction. I tried to clarify this in the documentation.
– user21
Sep 3 at 5:43




Concerning the BoundaryGroupingThreshold option: this is to attribute markers to groups of contiguous boundary segments/faces. The threshold does not change the form of a boundary; it changes what ToBoundayMesh considers as segments/faces that point in the same direction. I tried to clarify this in the documentation.
– user21
Sep 3 at 5:43












Thanks @user21, that helps clarify. Do you have any insight to the more general problem?
– Bruce Lucas
Sep 3 at 12:54




Thanks @user21, that helps clarify. Do you have any insight to the more general problem?
– Bruce Lucas
Sep 3 at 12:54












I'd report this to the support.
– user21
Sep 4 at 6:14




I'd report this to the support.
– user21
Sep 4 at 6:14












Thanks, I will do that. I suspect it's a limitation of the algorithm, although I think it should be possible.
– Bruce Lucas
Sep 5 at 0:41




Thanks, I will do that. I suspect it's a limitation of the algorithm, although I think it should be possible.
– Bruce Lucas
Sep 5 at 0:41










2 Answers
2






active

oldest

votes

















up vote
4
down vote













Try this:



r = 1;
t = 0.2;
h = 0.3;
rdiff = RegionDifference;
rplus = RegionUnion;
vertcyl[rad_, height_, offset_] := 
 Cylinder[0, 0, offset, 0, 0, offset + height, rad]
shape1outer = vertcyl[r, 1, 0];
shape1inner = vertcyl[r - t, 1, t];
shellouter = vertcyl[r, 1 - t, t];
shellinner = vertcyl[r - t, 2, 0];
base = vertcyl[r, t, 0];
hole = vertcyl[h, t, 0];

shape1outer, shape1inner, shellouter, shellinner, base, hole = 
 DiscretizeRegion /@ shape1outer, shape1inner, shellouter, 
 shellinner, base, hole;

shape1 = rdiff[shape1outer, shape1inner];
shape2 = rdiff[base, hole];
shell = rdiff[shellouter, shellinner];
finalshape = rdiff[shape1, hole];

finalshape // RegionPlot3D


final result



Summary:



  • Always use convex region for CSG(union, subtract etc.) operations.


  • Turn convex region into polygonal mesh before CSG.


  • Make sure no faces overlap each other during CSG.






share|improve this answer






















  • Thanks for the suggestion, it's an interesting approach. However I'm having trouble applying it to the more general cases I'm interested in; I've edited my question to clarify. (Also, to clarify, you said to turn the region into a polygonal mesh, but DiscretizeRegion generates a polyhedral mesh, doesn't it? The only practical consequence is extra processing time to generate the polyhedra that aren't needed for the final result; Export to STL handles this by just exporting the polyhedral boundary.)
    – Bruce Lucas
    Sep 3 at 12:57

















up vote
2
down vote













It appears Mathematica isn't currently capable of computing good surface meshes from implicit regions, so I worked some more on constructing surface meshes directly. Here's how the example above can be constructed (building some more general machinery along the way).



First we define a primitive for generating surface patches:



(* emit a quadrilateral face as two triangles
 so degenerate triangles can be culled *) 
face[a_, b_, c_, d_] := Polygon[a, c, d], Polygon[a, b, c];

(* construct a surface patch by varying u,v from umin,vmin
 to umax,vmax with nu,nv steps. f is an expression
 over u and v that gives the x,y,z location of the patch at u,v *) 
patch[f_, u_, umin_, umax_, nu_, v_, vmin_, vmax_, nv_] := 
 Quiet[MeshRegion[
 (* coordinates *)
 Flatten[Table[
 f /. u -> umin + (umax - umin) (i - 1) / nu, 
 v -> vmin + (vmax - vmin) (j - 1) / nv,
 i, nu + 1,
 j, nv + 1],
 1],
 (* faces specified using indexes into the coordinate list *) 
 Flatten[Table[
 face[
 (i - 1) (nv + 1) + j,
 i (nv + 1) + j,
 i (nv + 1) + (j + 1),
 (i - 1) (nv + 1) + (j + 1)],
 i, nu,
 j, nv],
 2]
 ], MeshRegion::dgcellr (* quietly cull degenerate triangles *)];

(* example *)
patch[u + v, u - v, u v, u, 0, 1, 10, v, 0, 1, 10]


example patch



Now we build up some shapes as surface patches:



(* tube: cylindrical surface of possibly varying radius from p going in the z direction *)
ztube[p_, r_, nθ_, zmin_, zmax_, nz_] := patch[
 r Cos[θ], r Sin[θ], z + p ,
 θ, 0, 2 π, nθ,
 z, zmin, zmax, nz];

(* ring in xy plane with center at p, inner radius r0, outer radius r1 *)
zring[p_, r0_, r1_, nθ_] := patch[
 r Cos[θ], r Sin[θ], 0 + p,
 θ, 0, 2 π, nθ,
 r, r0, r1, 1];

(* cylindrical surface composed of a tube and two disks (rings with inner radius 0) *)
zcyl[p_, r_, h_, nθ_, nz_] := RegionUnion[
 ztube[p, r, nθ, 0, h, nz],
 zring[p, 0, r /. z -> 0, nθ],
 zring[p + 0, 0, h, 0, r /. z -> h, nθ]];

(* example *)
zcyl[0, 0, 0, 1, 1, 20, 1]


example cylinder



We can use these primitives to build up a more complex shape consisting of the union of 5 separate surfaces:



r[z_] := 0.6 + 0.4 z^2; (* varying radius *)
t := 0.1; (* thickness *)
nθ := 100; (* number of radial segments *)
nz := 10; (* number of segments in z direction *)
shape = RegionUnion[
 ztube[0, 0, 0, r[z], nθ, 0, 1, nz], (* outer side *)
 zring[0, 0, 0, 0, r[0], nθ], (* outer bottom *)
 ztube[0, 0, 0, r[z] - t, nθ, t, 1, nz], (* inner side *)
 zring[0, 0, t, r[t] - t, 0, nθ], (* inner bottom *)
 zring[0, 0, 1, r[1] - t, r[1], nθ] (* rim *)
 ], MeshCells[shape, 2] // Length


more complex shape



Note that in order for this union of surfaces to describe a proper solid boundary, e.g. for 3d printing, the edges of adjacent faces must coincide. For example, we have to use the same number of radial segments nθ for all the surfaces.



The above shapes are surfaces that happen to bound a solid region, but Mathematica doesn't know that. If we try to perform a region boolean operation on these shapes Mathematica doesn't know how to do it because they are 2-d surfaces embedded in 3-d:



hole = zcyl[0, 0, 0, 0.1, t, 20, 1];
RegionDifference[shape, hole]


can't do it



However, provided the surface is a proper closed surface, we can tell Mathematica to treat it as a surface that bounds a solid, and then we can do boolean region operations on the solids:



(* turn a region into the solid bounded by that region *)
solid[region_] := BoundaryMeshRegion[MeshCoordinates[region], MeshCells[region, 2]];
shapeSolid = solid[shape];
holeSolid = solid[hole];
withHole = RegionDifference[shapeSolid, holeSolid],
 MeshCells[withHole, 2] // Length


csg



However the drawback of doing boolean operations over solids described using mesh boundaries is an explosion in the number of faces - computing this single region difference has turned the 4400 faces of the original shape into 22,788 faces. I suspect Mathematica is using BSP-based algorithms for this, and this explosion of faces seems to be a characteristic of systems that use BSP algorithms, like OpenSCAD and OpenJSCAD. It's better to build up the shapes from explicit surfaces if possible.






share|improve this answer








New contributor




Bruce Lucas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.

















  • Great (+1). Since you do have a symbolic representation you could try this to down sample the boundary; I can not guarantee that this will work but it might be a worth a shot. If you do try it, let me know if it works.
    – user21
    Sep 5 at 6:10










Your Answer




StackExchange.ifUsing("editor", function ()
return StackExchange.using("mathjaxEditing", function ()
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix)
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
);
);
, "mathjax-editing");

StackExchange.ready(function()
var channelOptions =
tags: "".split(" "),
id: "387"
;
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function()
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled)
StackExchange.using("snippets", function()
createEditor();
);

else
createEditor();

);

function createEditor()
StackExchange.prepareEditor(
heartbeatType: 'answer',
convertImagesToLinks: false,
noModals: false,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
bindNavPrevention: true,
postfix: "",
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
);



);













 

draft saved


draft discarded


















StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f181130%2finconsistent-mesh-generation-from-boundarydiscretizeregion-and-toboundarymesh-at%23new-answer', 'question_page');

);

Post as a guest

















2 Answers
2






active

oldest

votes








2 Answers
2






active

oldest

votes









active

oldest

votes






active

oldest

votes








up vote
4
down vote













Try this:



r = 1;
t = 0.2;
h = 0.3;
rdiff = RegionDifference;
rplus = RegionUnion;
vertcyl[rad_, height_, offset_] := 
 Cylinder[0, 0, offset, 0, 0, offset + height, rad]
shape1outer = vertcyl[r, 1, 0];
shape1inner = vertcyl[r - t, 1, t];
shellouter = vertcyl[r, 1 - t, t];
shellinner = vertcyl[r - t, 2, 0];
base = vertcyl[r, t, 0];
hole = vertcyl[h, t, 0];

shape1outer, shape1inner, shellouter, shellinner, base, hole = 
 DiscretizeRegion /@ shape1outer, shape1inner, shellouter, 
 shellinner, base, hole;

shape1 = rdiff[shape1outer, shape1inner];
shape2 = rdiff[base, hole];
shell = rdiff[shellouter, shellinner];
finalshape = rdiff[shape1, hole];

finalshape // RegionPlot3D


final result



Summary:



  • Always use convex region for CSG(union, subtract etc.) operations.


  • Turn convex region into polygonal mesh before CSG.


  • Make sure no faces overlap each other during CSG.






share|improve this answer






















  • Thanks for the suggestion, it's an interesting approach. However I'm having trouble applying it to the more general cases I'm interested in; I've edited my question to clarify. (Also, to clarify, you said to turn the region into a polygonal mesh, but DiscretizeRegion generates a polyhedral mesh, doesn't it? The only practical consequence is extra processing time to generate the polyhedra that aren't needed for the final result; Export to STL handles this by just exporting the polyhedral boundary.)
    – Bruce Lucas
    Sep 3 at 12:57














up vote
4
down vote













Try this:



r = 1;
t = 0.2;
h = 0.3;
rdiff = RegionDifference;
rplus = RegionUnion;
vertcyl[rad_, height_, offset_] := 
 Cylinder[0, 0, offset, 0, 0, offset + height, rad]
shape1outer = vertcyl[r, 1, 0];
shape1inner = vertcyl[r - t, 1, t];
shellouter = vertcyl[r, 1 - t, t];
shellinner = vertcyl[r - t, 2, 0];
base = vertcyl[r, t, 0];
hole = vertcyl[h, t, 0];

shape1outer, shape1inner, shellouter, shellinner, base, hole = 
 DiscretizeRegion /@ shape1outer, shape1inner, shellouter, 
 shellinner, base, hole;

shape1 = rdiff[shape1outer, shape1inner];
shape2 = rdiff[base, hole];
shell = rdiff[shellouter, shellinner];
finalshape = rdiff[shape1, hole];

finalshape // RegionPlot3D


final result



Summary:



  • Always use convex region for CSG(union, subtract etc.) operations.


  • Turn convex region into polygonal mesh before CSG.


  • Make sure no faces overlap each other during CSG.






share|improve this answer






















  • Thanks for the suggestion, it's an interesting approach. However I'm having trouble applying it to the more general cases I'm interested in; I've edited my question to clarify. (Also, to clarify, you said to turn the region into a polygonal mesh, but DiscretizeRegion generates a polyhedral mesh, doesn't it? The only practical consequence is extra processing time to generate the polyhedra that aren't needed for the final result; Export to STL handles this by just exporting the polyhedral boundary.)
    – Bruce Lucas
    Sep 3 at 12:57












up vote
4
down vote










up vote
4
down vote









Try this:



r = 1;
t = 0.2;
h = 0.3;
rdiff = RegionDifference;
rplus = RegionUnion;
vertcyl[rad_, height_, offset_] := 
 Cylinder[0, 0, offset, 0, 0, offset + height, rad]
shape1outer = vertcyl[r, 1, 0];
shape1inner = vertcyl[r - t, 1, t];
shellouter = vertcyl[r, 1 - t, t];
shellinner = vertcyl[r - t, 2, 0];
base = vertcyl[r, t, 0];
hole = vertcyl[h, t, 0];

shape1outer, shape1inner, shellouter, shellinner, base, hole = 
 DiscretizeRegion /@ shape1outer, shape1inner, shellouter, 
 shellinner, base, hole;

shape1 = rdiff[shape1outer, shape1inner];
shape2 = rdiff[base, hole];
shell = rdiff[shellouter, shellinner];
finalshape = rdiff[shape1, hole];

finalshape // RegionPlot3D


final result



Summary:



  • Always use convex region for CSG(union, subtract etc.) operations.


  • Turn convex region into polygonal mesh before CSG.


  • Make sure no faces overlap each other during CSG.






share|improve this answer














Try this:



r = 1;
t = 0.2;
h = 0.3;
rdiff = RegionDifference;
rplus = RegionUnion;
vertcyl[rad_, height_, offset_] := 
 Cylinder[0, 0, offset, 0, 0, offset + height, rad]
shape1outer = vertcyl[r, 1, 0];
shape1inner = vertcyl[r - t, 1, t];
shellouter = vertcyl[r, 1 - t, t];
shellinner = vertcyl[r - t, 2, 0];
base = vertcyl[r, t, 0];
hole = vertcyl[h, t, 0];

shape1outer, shape1inner, shellouter, shellinner, base, hole = 
 DiscretizeRegion /@ shape1outer, shape1inner, shellouter, 
 shellinner, base, hole;

shape1 = rdiff[shape1outer, shape1inner];
shape2 = rdiff[base, hole];
shell = rdiff[shellouter, shellinner];
finalshape = rdiff[shape1, hole];

finalshape // RegionPlot3D


final result



Summary:



  • Always use convex region for CSG(union, subtract etc.) operations.


  • Turn convex region into polygonal mesh before CSG.


  • Make sure no faces overlap each other during CSG.







share|improve this answer














share|improve this answer



share|improve this answer








edited Sep 3 at 3:15

























answered Sep 3 at 3:05









kh40tika

22129




22129











  • Thanks for the suggestion, it's an interesting approach. However I'm having trouble applying it to the more general cases I'm interested in; I've edited my question to clarify. (Also, to clarify, you said to turn the region into a polygonal mesh, but DiscretizeRegion generates a polyhedral mesh, doesn't it? The only practical consequence is extra processing time to generate the polyhedra that aren't needed for the final result; Export to STL handles this by just exporting the polyhedral boundary.)
    – Bruce Lucas
    Sep 3 at 12:57
















  • Thanks for the suggestion, it's an interesting approach. However I'm having trouble applying it to the more general cases I'm interested in; I've edited my question to clarify. (Also, to clarify, you said to turn the region into a polygonal mesh, but DiscretizeRegion generates a polyhedral mesh, doesn't it? The only practical consequence is extra processing time to generate the polyhedra that aren't needed for the final result; Export to STL handles this by just exporting the polyhedral boundary.)
    – Bruce Lucas
    Sep 3 at 12:57















Thanks for the suggestion, it's an interesting approach. However I'm having trouble applying it to the more general cases I'm interested in; I've edited my question to clarify. (Also, to clarify, you said to turn the region into a polygonal mesh, but DiscretizeRegion generates a polyhedral mesh, doesn't it? The only practical consequence is extra processing time to generate the polyhedra that aren't needed for the final result; Export to STL handles this by just exporting the polyhedral boundary.)
– Bruce Lucas
Sep 3 at 12:57




Thanks for the suggestion, it's an interesting approach. However I'm having trouble applying it to the more general cases I'm interested in; I've edited my question to clarify. (Also, to clarify, you said to turn the region into a polygonal mesh, but DiscretizeRegion generates a polyhedral mesh, doesn't it? The only practical consequence is extra processing time to generate the polyhedra that aren't needed for the final result; Export to STL handles this by just exporting the polyhedral boundary.)
– Bruce Lucas
Sep 3 at 12:57










up vote
2
down vote













It appears Mathematica isn't currently capable of computing good surface meshes from implicit regions, so I worked some more on constructing surface meshes directly. Here's how the example above can be constructed (building some more general machinery along the way).



First we define a primitive for generating surface patches:



(* emit a quadrilateral face as two triangles
 so degenerate triangles can be culled *) 
face[a_, b_, c_, d_] := Polygon[a, c, d], Polygon[a, b, c];

(* construct a surface patch by varying u,v from umin,vmin
 to umax,vmax with nu,nv steps. f is an expression
 over u and v that gives the x,y,z location of the patch at u,v *) 
patch[f_, u_, umin_, umax_, nu_, v_, vmin_, vmax_, nv_] := 
 Quiet[MeshRegion[
 (* coordinates *)
 Flatten[Table[
 f /. u -> umin + (umax - umin) (i - 1) / nu, 
 v -> vmin + (vmax - vmin) (j - 1) / nv,
 i, nu + 1,
 j, nv + 1],
 1],
 (* faces specified using indexes into the coordinate list *) 
 Flatten[Table[
 face[
 (i - 1) (nv + 1) + j,
 i (nv + 1) + j,
 i (nv + 1) + (j + 1),
 (i - 1) (nv + 1) + (j + 1)],
 i, nu,
 j, nv],
 2]
 ], MeshRegion::dgcellr (* quietly cull degenerate triangles *)];

(* example *)
patch[u + v, u - v, u v, u, 0, 1, 10, v, 0, 1, 10]


example patch



Now we build up some shapes as surface patches:



(* tube: cylindrical surface of possibly varying radius from p going in the z direction *)
ztube[p_, r_, nθ_, zmin_, zmax_, nz_] := patch[
 r Cos[θ], r Sin[θ], z + p ,
 θ, 0, 2 π, nθ,
 z, zmin, zmax, nz];

(* ring in xy plane with center at p, inner radius r0, outer radius r1 *)
zring[p_, r0_, r1_, nθ_] := patch[
 r Cos[θ], r Sin[θ], 0 + p,
 θ, 0, 2 π, nθ,
 r, r0, r1, 1];

(* cylindrical surface composed of a tube and two disks (rings with inner radius 0) *)
zcyl[p_, r_, h_, nθ_, nz_] := RegionUnion[
 ztube[p, r, nθ, 0, h, nz],
 zring[p, 0, r /. z -> 0, nθ],
 zring[p + 0, 0, h, 0, r /. z -> h, nθ]];

(* example *)
zcyl[0, 0, 0, 1, 1, 20, 1]


example cylinder



We can use these primitives to build up a more complex shape consisting of the union of 5 separate surfaces:



r[z_] := 0.6 + 0.4 z^2; (* varying radius *)
t := 0.1; (* thickness *)
nθ := 100; (* number of radial segments *)
nz := 10; (* number of segments in z direction *)
shape = RegionUnion[
 ztube[0, 0, 0, r[z], nθ, 0, 1, nz], (* outer side *)
 zring[0, 0, 0, 0, r[0], nθ], (* outer bottom *)
 ztube[0, 0, 0, r[z] - t, nθ, t, 1, nz], (* inner side *)
 zring[0, 0, t, r[t] - t, 0, nθ], (* inner bottom *)
 zring[0, 0, 1, r[1] - t, r[1], nθ] (* rim *)
 ], MeshCells[shape, 2] // Length


more complex shape



Note that in order for this union of surfaces to describe a proper solid boundary, e.g. for 3d printing, the edges of adjacent faces must coincide. For example, we have to use the same number of radial segments nθ for all the surfaces.



The above shapes are surfaces that happen to bound a solid region, but Mathematica doesn't know that. If we try to perform a region boolean operation on these shapes Mathematica doesn't know how to do it because they are 2-d surfaces embedded in 3-d:



hole = zcyl[0, 0, 0, 0.1, t, 20, 1];
RegionDifference[shape, hole]


can't do it



However, provided the surface is a proper closed surface, we can tell Mathematica to treat it as a surface that bounds a solid, and then we can do boolean region operations on the solids:



(* turn a region into the solid bounded by that region *)
solid[region_] := BoundaryMeshRegion[MeshCoordinates[region], MeshCells[region, 2]];
shapeSolid = solid[shape];
holeSolid = solid[hole];
withHole = RegionDifference[shapeSolid, holeSolid],
 MeshCells[withHole, 2] // Length


csg



However the drawback of doing boolean operations over solids described using mesh boundaries is an explosion in the number of faces - computing this single region difference has turned the 4400 faces of the original shape into 22,788 faces. I suspect Mathematica is using BSP-based algorithms for this, and this explosion of faces seems to be a characteristic of systems that use BSP algorithms, like OpenSCAD and OpenJSCAD. It's better to build up the shapes from explicit surfaces if possible.






share|improve this answer








New contributor




Bruce Lucas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.

















  • Great (+1). Since you do have a symbolic representation you could try this to down sample the boundary; I can not guarantee that this will work but it might be a worth a shot. If you do try it, let me know if it works.
    – user21
    Sep 5 at 6:10














up vote
2
down vote













It appears Mathematica isn't currently capable of computing good surface meshes from implicit regions, so I worked some more on constructing surface meshes directly. Here's how the example above can be constructed (building some more general machinery along the way).



First we define a primitive for generating surface patches:



(* emit a quadrilateral face as two triangles
 so degenerate triangles can be culled *) 
face[a_, b_, c_, d_] := Polygon[a, c, d], Polygon[a, b, c];

(* construct a surface patch by varying u,v from umin,vmin
 to umax,vmax with nu,nv steps. f is an expression
 over u and v that gives the x,y,z location of the patch at u,v *) 
patch[f_, u_, umin_, umax_, nu_, v_, vmin_, vmax_, nv_] := 
 Quiet[MeshRegion[
 (* coordinates *)
 Flatten[Table[
 f /. u -> umin + (umax - umin) (i - 1) / nu, 
 v -> vmin + (vmax - vmin) (j - 1) / nv,
 i, nu + 1,
 j, nv + 1],
 1],
 (* faces specified using indexes into the coordinate list *) 
 Flatten[Table[
 face[
 (i - 1) (nv + 1) + j,
 i (nv + 1) + j,
 i (nv + 1) + (j + 1),
 (i - 1) (nv + 1) + (j + 1)],
 i, nu,
 j, nv],
 2]
 ], MeshRegion::dgcellr (* quietly cull degenerate triangles *)];

(* example *)
patch[u + v, u - v, u v, u, 0, 1, 10, v, 0, 1, 10]


example patch



Now we build up some shapes as surface patches:



(* tube: cylindrical surface of possibly varying radius from p going in the z direction *)
ztube[p_, r_, nθ_, zmin_, zmax_, nz_] := patch[
 r Cos[θ], r Sin[θ], z + p ,
 θ, 0, 2 π, nθ,
 z, zmin, zmax, nz];

(* ring in xy plane with center at p, inner radius r0, outer radius r1 *)
zring[p_, r0_, r1_, nθ_] := patch[
 r Cos[θ], r Sin[θ], 0 + p,
 θ, 0, 2 π, nθ,
 r, r0, r1, 1];

(* cylindrical surface composed of a tube and two disks (rings with inner radius 0) *)
zcyl[p_, r_, h_, nθ_, nz_] := RegionUnion[
 ztube[p, r, nθ, 0, h, nz],
 zring[p, 0, r /. z -> 0, nθ],
 zring[p + 0, 0, h, 0, r /. z -> h, nθ]];

(* example *)
zcyl[0, 0, 0, 1, 1, 20, 1]


example cylinder



We can use these primitives to build up a more complex shape consisting of the union of 5 separate surfaces:



r[z_] := 0.6 + 0.4 z^2; (* varying radius *)
t := 0.1; (* thickness *)
nθ := 100; (* number of radial segments *)
nz := 10; (* number of segments in z direction *)
shape = RegionUnion[
 ztube[0, 0, 0, r[z], nθ, 0, 1, nz], (* outer side *)
 zring[0, 0, 0, 0, r[0], nθ], (* outer bottom *)
 ztube[0, 0, 0, r[z] - t, nθ, t, 1, nz], (* inner side *)
 zring[0, 0, t, r[t] - t, 0, nθ], (* inner bottom *)
 zring[0, 0, 1, r[1] - t, r[1], nθ] (* rim *)
 ], MeshCells[shape, 2] // Length


more complex shape



Note that in order for this union of surfaces to describe a proper solid boundary, e.g. for 3d printing, the edges of adjacent faces must coincide. For example, we have to use the same number of radial segments nθ for all the surfaces.



The above shapes are surfaces that happen to bound a solid region, but Mathematica doesn't know that. If we try to perform a region boolean operation on these shapes Mathematica doesn't know how to do it because they are 2-d surfaces embedded in 3-d:



hole = zcyl[0, 0, 0, 0.1, t, 20, 1];
RegionDifference[shape, hole]


can't do it



However, provided the surface is a proper closed surface, we can tell Mathematica to treat it as a surface that bounds a solid, and then we can do boolean region operations on the solids:



(* turn a region into the solid bounded by that region *)
solid[region_] := BoundaryMeshRegion[MeshCoordinates[region], MeshCells[region, 2]];
shapeSolid = solid[shape];
holeSolid = solid[hole];
withHole = RegionDifference[shapeSolid, holeSolid],
 MeshCells[withHole, 2] // Length


csg



However the drawback of doing boolean operations over solids described using mesh boundaries is an explosion in the number of faces - computing this single region difference has turned the 4400 faces of the original shape into 22,788 faces. I suspect Mathematica is using BSP-based algorithms for this, and this explosion of faces seems to be a characteristic of systems that use BSP algorithms, like OpenSCAD and OpenJSCAD. It's better to build up the shapes from explicit surfaces if possible.






share|improve this answer








New contributor




Bruce Lucas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.

















  • Great (+1). Since you do have a symbolic representation you could try this to down sample the boundary; I can not guarantee that this will work but it might be a worth a shot. If you do try it, let me know if it works.
    – user21
    Sep 5 at 6:10












up vote
2
down vote










up vote
2
down vote









It appears Mathematica isn't currently capable of computing good surface meshes from implicit regions, so I worked some more on constructing surface meshes directly. Here's how the example above can be constructed (building some more general machinery along the way).



First we define a primitive for generating surface patches:



(* emit a quadrilateral face as two triangles
 so degenerate triangles can be culled *) 
face[a_, b_, c_, d_] := Polygon[a, c, d], Polygon[a, b, c];

(* construct a surface patch by varying u,v from umin,vmin
 to umax,vmax with nu,nv steps. f is an expression
 over u and v that gives the x,y,z location of the patch at u,v *) 
patch[f_, u_, umin_, umax_, nu_, v_, vmin_, vmax_, nv_] := 
 Quiet[MeshRegion[
 (* coordinates *)
 Flatten[Table[
 f /. u -> umin + (umax - umin) (i - 1) / nu, 
 v -> vmin + (vmax - vmin) (j - 1) / nv,
 i, nu + 1,
 j, nv + 1],
 1],
 (* faces specified using indexes into the coordinate list *) 
 Flatten[Table[
 face[
 (i - 1) (nv + 1) + j,
 i (nv + 1) + j,
 i (nv + 1) + (j + 1),
 (i - 1) (nv + 1) + (j + 1)],
 i, nu,
 j, nv],
 2]
 ], MeshRegion::dgcellr (* quietly cull degenerate triangles *)];

(* example *)
patch[u + v, u - v, u v, u, 0, 1, 10, v, 0, 1, 10]


example patch



Now we build up some shapes as surface patches:



(* tube: cylindrical surface of possibly varying radius from p going in the z direction *)
ztube[p_, r_, nθ_, zmin_, zmax_, nz_] := patch[
 r Cos[θ], r Sin[θ], z + p ,
 θ, 0, 2 π, nθ,
 z, zmin, zmax, nz];

(* ring in xy plane with center at p, inner radius r0, outer radius r1 *)
zring[p_, r0_, r1_, nθ_] := patch[
 r Cos[θ], r Sin[θ], 0 + p,
 θ, 0, 2 π, nθ,
 r, r0, r1, 1];

(* cylindrical surface composed of a tube and two disks (rings with inner radius 0) *)
zcyl[p_, r_, h_, nθ_, nz_] := RegionUnion[
 ztube[p, r, nθ, 0, h, nz],
 zring[p, 0, r /. z -> 0, nθ],
 zring[p + 0, 0, h, 0, r /. z -> h, nθ]];

(* example *)
zcyl[0, 0, 0, 1, 1, 20, 1]


example cylinder



We can use these primitives to build up a more complex shape consisting of the union of 5 separate surfaces:



r[z_] := 0.6 + 0.4 z^2; (* varying radius *)
t := 0.1; (* thickness *)
nθ := 100; (* number of radial segments *)
nz := 10; (* number of segments in z direction *)
shape = RegionUnion[
 ztube[0, 0, 0, r[z], nθ, 0, 1, nz], (* outer side *)
 zring[0, 0, 0, 0, r[0], nθ], (* outer bottom *)
 ztube[0, 0, 0, r[z] - t, nθ, t, 1, nz], (* inner side *)
 zring[0, 0, t, r[t] - t, 0, nθ], (* inner bottom *)
 zring[0, 0, 1, r[1] - t, r[1], nθ] (* rim *)
 ], MeshCells[shape, 2] // Length


more complex shape



Note that in order for this union of surfaces to describe a proper solid boundary, e.g. for 3d printing, the edges of adjacent faces must coincide. For example, we have to use the same number of radial segments nθ for all the surfaces.



The above shapes are surfaces that happen to bound a solid region, but Mathematica doesn't know that. If we try to perform a region boolean operation on these shapes Mathematica doesn't know how to do it because they are 2-d surfaces embedded in 3-d:



hole = zcyl[0, 0, 0, 0.1, t, 20, 1];
RegionDifference[shape, hole]


can't do it



However, provided the surface is a proper closed surface, we can tell Mathematica to treat it as a surface that bounds a solid, and then we can do boolean region operations on the solids:



(* turn a region into the solid bounded by that region *)
solid[region_] := BoundaryMeshRegion[MeshCoordinates[region], MeshCells[region, 2]];
shapeSolid = solid[shape];
holeSolid = solid[hole];
withHole = RegionDifference[shapeSolid, holeSolid],
 MeshCells[withHole, 2] // Length


csg



However the drawback of doing boolean operations over solids described using mesh boundaries is an explosion in the number of faces - computing this single region difference has turned the 4400 faces of the original shape into 22,788 faces. I suspect Mathematica is using BSP-based algorithms for this, and this explosion of faces seems to be a characteristic of systems that use BSP algorithms, like OpenSCAD and OpenJSCAD. It's better to build up the shapes from explicit surfaces if possible.






share|improve this answer








New contributor




Bruce Lucas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.









It appears Mathematica isn't currently capable of computing good surface meshes from implicit regions, so I worked some more on constructing surface meshes directly. Here's how the example above can be constructed (building some more general machinery along the way).



First we define a primitive for generating surface patches:



(* emit a quadrilateral face as two triangles
 so degenerate triangles can be culled *) 
face[a_, b_, c_, d_] := Polygon[a, c, d], Polygon[a, b, c];

(* construct a surface patch by varying u,v from umin,vmin
 to umax,vmax with nu,nv steps. f is an expression
 over u and v that gives the x,y,z location of the patch at u,v *) 
patch[f_, u_, umin_, umax_, nu_, v_, vmin_, vmax_, nv_] := 
 Quiet[MeshRegion[
 (* coordinates *)
 Flatten[Table[
 f /. u -> umin + (umax - umin) (i - 1) / nu, 
 v -> vmin + (vmax - vmin) (j - 1) / nv,
 i, nu + 1,
 j, nv + 1],
 1],
 (* faces specified using indexes into the coordinate list *) 
 Flatten[Table[
 face[
 (i - 1) (nv + 1) + j,
 i (nv + 1) + j,
 i (nv + 1) + (j + 1),
 (i - 1) (nv + 1) + (j + 1)],
 i, nu,
 j, nv],
 2]
 ], MeshRegion::dgcellr (* quietly cull degenerate triangles *)];

(* example *)
patch[u + v, u - v, u v, u, 0, 1, 10, v, 0, 1, 10]


example patch



Now we build up some shapes as surface patches:



(* tube: cylindrical surface of possibly varying radius from p going in the z direction *)
ztube[p_, r_, nθ_, zmin_, zmax_, nz_] := patch[
 r Cos[θ], r Sin[θ], z + p ,
 θ, 0, 2 π, nθ,
 z, zmin, zmax, nz];

(* ring in xy plane with center at p, inner radius r0, outer radius r1 *)
zring[p_, r0_, r1_, nθ_] := patch[
 r Cos[θ], r Sin[θ], 0 + p,
 θ, 0, 2 π, nθ,
 r, r0, r1, 1];

(* cylindrical surface composed of a tube and two disks (rings with inner radius 0) *)
zcyl[p_, r_, h_, nθ_, nz_] := RegionUnion[
 ztube[p, r, nθ, 0, h, nz],
 zring[p, 0, r /. z -> 0, nθ],
 zring[p + 0, 0, h, 0, r /. z -> h, nθ]];

(* example *)
zcyl[0, 0, 0, 1, 1, 20, 1]


example cylinder



We can use these primitives to build up a more complex shape consisting of the union of 5 separate surfaces:



r[z_] := 0.6 + 0.4 z^2; (* varying radius *)
t := 0.1; (* thickness *)
nθ := 100; (* number of radial segments *)
nz := 10; (* number of segments in z direction *)
shape = RegionUnion[
 ztube[0, 0, 0, r[z], nθ, 0, 1, nz], (* outer side *)
 zring[0, 0, 0, 0, r[0], nθ], (* outer bottom *)
 ztube[0, 0, 0, r[z] - t, nθ, t, 1, nz], (* inner side *)
 zring[0, 0, t, r[t] - t, 0, nθ], (* inner bottom *)
 zring[0, 0, 1, r[1] - t, r[1], nθ] (* rim *)
 ], MeshCells[shape, 2] // Length


more complex shape



Note that in order for this union of surfaces to describe a proper solid boundary, e.g. for 3d printing, the edges of adjacent faces must coincide. For example, we have to use the same number of radial segments nθ for all the surfaces.



The above shapes are surfaces that happen to bound a solid region, but Mathematica doesn't know that. If we try to perform a region boolean operation on these shapes Mathematica doesn't know how to do it because they are 2-d surfaces embedded in 3-d:



hole = zcyl[0, 0, 0, 0.1, t, 20, 1];
RegionDifference[shape, hole]


can't do it



However, provided the surface is a proper closed surface, we can tell Mathematica to treat it as a surface that bounds a solid, and then we can do boolean region operations on the solids:



(* turn a region into the solid bounded by that region *)
solid[region_] := BoundaryMeshRegion[MeshCoordinates[region], MeshCells[region, 2]];
shapeSolid = solid[shape];
holeSolid = solid[hole];
withHole = RegionDifference[shapeSolid, holeSolid],
 MeshCells[withHole, 2] // Length


csg



However the drawback of doing boolean operations over solids described using mesh boundaries is an explosion in the number of faces - computing this single region difference has turned the 4400 faces of the original shape into 22,788 faces. I suspect Mathematica is using BSP-based algorithms for this, and this explosion of faces seems to be a characteristic of systems that use BSP algorithms, like OpenSCAD and OpenJSCAD. It's better to build up the shapes from explicit surfaces if possible.







share|improve this answer








New contributor




Bruce Lucas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.









share|improve this answer



share|improve this answer






New contributor




Bruce Lucas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.









answered Sep 5 at 1:34









Bruce Lucas

1714




1714




New contributor




Bruce Lucas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.





New contributor





Bruce Lucas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.






Bruce Lucas is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.











  • Great (+1). Since you do have a symbolic representation you could try this to down sample the boundary; I can not guarantee that this will work but it might be a worth a shot. If you do try it, let me know if it works.
    – user21
    Sep 5 at 6:10
















  • Great (+1). Since you do have a symbolic representation you could try this to down sample the boundary; I can not guarantee that this will work but it might be a worth a shot. If you do try it, let me know if it works.
    – user21
    Sep 5 at 6:10















Great (+1). Since you do have a symbolic representation you could try this to down sample the boundary; I can not guarantee that this will work but it might be a worth a shot. If you do try it, let me know if it works.
– user21
Sep 5 at 6:10




Great (+1). Since you do have a symbolic representation you could try this to down sample the boundary; I can not guarantee that this will work but it might be a worth a shot. If you do try it, let me know if it works.
– user21
Sep 5 at 6:10

















 

draft saved


draft discarded















































 


draft saved


draft discarded














StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f181130%2finconsistent-mesh-generation-from-boundarydiscretizeregion-and-toboundarymesh-at%23new-answer', 'question_page');

);

Post as a guest
































































Comments

Post a Comment

Popular posts from this blog

What does second last employer means? [closed]

Image
Clash Royale CLAN TAG #URR8PPP .everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty margin-bottom:0; up vote -3 down vote favorite Say, for example I have worked in the following companies, Company A => 2010 - 2011 Company B => 2011 - 2012 Company C => 2012 - Till Date Company D => Have offer Now assume Company D is asking to submit documents from my second last employer. So from the list of companies given, which company would be qualified as second last employer? software-industry employer share | improve this question asked Aug 10 '14 at 5:16 Rajaprabhu Aravindasamy 559 1 6 13 closed as off-topic by Jim G., jmort253 ♦ Aug 11 '14 at 1:58 This question appears to be off-topic. The users who voted to close gave this specific reason: "Questions seeking advice on company-specific regulations, agreements, or policies should be directed to your manager or HR department. Questions that address...
Read more

Installing NextGIS Connect into QGIS 3?

Image
Clash Royale CLAN TAG #URR8PPP .everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty margin-bottom:0; up vote 1 down vote favorite Just acquired free version of NEXTGIS and trying to install the plugin for QGIS. It tells me it was not designed for version 3 of QGIS. How do I set up the QGIS-NEXTGIS interface then? qgis qgis-plugins nextgis share | improve this question edited 3 hours ago PolyGeo ♦ 51.9k 17 77 233 asked 4 hours ago IanGS 16 4 New contributor IanGS is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct. add a comment  |  up vote 1 down vote favorite Just acquired free version of NEXTGIS and trying to install the plugin for QGIS. It tells me it was not designed for version 3 of QGIS. How do I set up the QGIS-NEXTGIS interface then? qgis qgis-plugin...
Read more

One-line joke

A one-liner is a joke that is delivered in a single line. A good one-liner is said to be pithy - concise and meaningful. [1] Comedians and actors use this comedic method as part of their act, e.g. Jimmy Carr, Tommy Cooper, Rodney Dangerfield, Norm Macdonald, Ken Dodd, Stewart Francis, Zach Galifianakis, Mitch Hedberg, Anthony Jeselnik, Milton Jones, Shappi Khorsandi, Jay London, Mark Linn-Baker, Demetri Martin, Groucho Marx, Dan Mintz, Emo Philips, Tim Vine, Steven Wright and Henny Youngman and Arnold Schwarzenegger. [ citation needed ] Many fictional characters are also known to deliver one-liners, including James Bond, who usually includes pithy and laconic quips after disposing of a villain. [ citation needed ] Examples "A baby seal walks into a club." "A dyslexic man walks into a bra." (Anonymous) "There are three types of people, those who can count and those who can't." "Venison's dear, isn't it?" (Jimmy Carr) "An escala...
Read more