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

Long meetings (6-7 hours a day): Being “babysat” by supervisor

Image
Clash Royale CLAN TAG #URR8PPP .everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty margin-bottom:0; up vote 31 down vote favorite 4 I hold a PhD in mechanical engineering with a 7 year background in self-guided, self-paced research including 2 years as a graduate level course instructor and a subject matter expert for various engineering projects in a university setting in the United States. For the last 10-11 months, I am working in a research industry in France where my task is to debug C++ spaghetti code written over a period of several years by my current supervisor. It is a team of 2: just I and my supervisor. Diurnally, I am posed with a "bug of the day" to get out of this C++ code. My supervisor generally gives me about 15-30 minutes to solve it and then accosts me about whether or not the task is completed and then rinses and repeats this until the end of the day. Off-late, my supervisor has been holding 6-7 hour long "meetin
Read more

Is the Concept of Multiple Fantasy Races Scientifically Flawed? [closed]

Image
Clash Royale CLAN TAG #URR8PPP up vote 2 down vote favorite Many fantasy worlds have multiple humanoid species (elves, orcs, humans, etc) living in the same world, in addition to that, they often have the ability to interbreed. And I can't help but question their believability... fantasy-races reproduction interspecies-relations share | improve this question edited Aug 9 at 14:46 Michael Kjörling ♦ 21.1k 10 85 161 asked Aug 9 at 13:13 Chlodio 118 6 closed as unclear what you're asking by MichaelK, Gryphon, John, Vincent, Frostfyre Aug 9 at 15:25 Please clarify your specific problem or add additional details to highlight exactly what you need. As it's currently written, it’s hard to tell exactly what you're asking. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question. 3 Stack Exchange exist to
Read more

Confectionery

Image
Prepared foods made with a lot of sugar or other cabohydrates "Sweetmeat" redirects here. For the racehorse, see Sweetmeat (horse). Not to be confused with Sweetbread. This Kransekake is a traditional Scandinavian baker's confection. Confectionery is the art of making confections , which are food items that are rich in sugar and carbohydrates. Exact definitions are difficult. [1] In general, though, confectionery is divided into two broad and somewhat overlapping categories, bakers' confections and sugar confections. [2] Bakers' confectionery, also called flour confections , includes principally sweet pastries, cakes, and similar baked goods. Sugar confectionery includes candies ( sweets in British English), candied nuts, chocolates, chewing gum, bubble gum, pastillage, and other confections that are made primarily of sugar. In some cases, chocolate confections (confections made of chocolate) are treated as a separate category, as are sugar-free versions of
Read more