17
$\begingroup$

Context

I am interested in converting surfaces into solids (so I can make a 3D mesh out of them using ToElementMesh)

Say I have the following cool surface

f[t_] := With[{s = 3 t/2}, {(2 + Cos[s]) Cos[t], (2 + Cos[s]) Sin[t], Sin[s]} - {2, 0, 0}]
v1[t_] := Cross[f'[t], {0, 0, 1}] // Normalize
v2[t_] := Cross[f'[t], v1[t]] // Normalize
g[t_, θ_] := f[t] + (Cos[θ] v1[t] + Sin[θ] v2[t])/2
gr = ParametricPlot3D[Evaluate @ g[t, θ], {t, 0, 4 Pi}, {θ, 0, 2 Pi}, 
                      Mesh -> None, MaxRecursion -> 4, Boxed -> False, Axes -> False];

which I pinched from this answer.

Mathematica graphics

Question

How would like to make a solid out of what is in it?

I have a hint that there is a simple answer to this problem in the context of the new computational geometry tools of Mathematica 10 but I have not found it.

In some sense I am after the opposite of BoundaryMesh

$\endgroup$
3
  • $\begingroup$ Try DiscretizeGraphics $\endgroup$
    – user21
    Jun 11, 2015 at 16:23
  • $\begingroup$ In[491]:= gr // DiscretizeGraphics Out[491]= EmptyRegion[3]; so its a bit disapointing $\endgroup$
    – chris
    Jun 11, 2015 at 16:23
  • 3
    $\begingroup$ You could convert it to a ParametricRegion and then use ToElementMesh - some examples are in the documentation, see also the tutorial to element mesh generation. $\endgroup$
    – user21
    Jun 11, 2015 at 16:33

3 Answers 3

17
$\begingroup$

Method 1: Construct mesh elements manually

We can triangulate a periodic quad-lattice on the surface:

n = {180, 20};  (* number of points in each direction *)
pts = Table[
   g[4. Pi/n[[1]] t, 2. Pi/ n[[2]] θ], {t, n[[1]]}, {θ, n[[2]]}];

idcs = {{{1, 2, 4}, {1, 4, 3}}, {{1, 2, 3}, {2, 4, 3}}};  (* for a diamond pattern *)
tri = 1 +
  Array[Function[quad, quad[[#]] & /@ idcs[[Mod[+##, 2, 1]]]][Tuples[
        {Mod[#1 + {-1, 0}, Dimensions[pts][[1]]],
         Mod[#2 + {-1, 0}, Dimensions[pts][[2]]]}].{Length[First@pts], 1}] &,
       Most@Dimensions[pts]];

Needs["NDSolve`FEM`"].
bmesh = ToBoundaryMesh[
  "Coordinates" -> Flatten[pts, 1],
  "BoundaryElements" -> {TriangleElement[Flatten[tri, 2]]}
  ]

emesh = ToElementMesh[bmesh]
(*
  ElementMesh[{{-4.86396, 1.5}, {-3.32664, 3.32664}, {-1.5, 1.5}},
   {TetrahedronElement["<" 21544 ">"]}]
*)

MeshRegion@emesh

Mathematica graphics

Notes:

Tuple generates the indices of the quadrilaterals in each row (#1) & column (#2`), wrapping around at the end of the domain to close up the tube:

Tuples[{Mod[#1 + {-1, 0}, Dimensions[pts][[1]]],
        Mod[#2 + {-1, 0}, Dimensions[pts][[2]]]}]

The dot product with {Length[First@pts], 1} converts a {row, column} pair to the index of the corresponding point in pts.

We triangulate the quadrilaterals by alternating which diagonal is used. The variable idcs contains two lists of two triangles, representing both ways. The integers themselves are indices of the quadrilateral (given by Tuples[..].{Length[..], 1} above).

idcs = {{{1, 2, 4}, {1, 4, 3}},   (* 1-4 diagonal *)
        {{1, 2, 3}, {2, 4, 3}}};  (* 2-3 diagonal *)

Method 2: Mesh the domain and apply parametrization

This takes more advantage of FEM/ElementMesh capabilities. First mesh the domain. Then apply g to map domain mesh coordinates onto to the surface. Use these coordinates and the domain mesh to construct a boundary mesh on the surface. Finally, use ToElementMesh to construct a mesh of the solid.

 Needs["NDSolve`FEM`"]

tscale = 4; θscale = 0.5;                (* scale roughly proportional to speeds *)
dom = ToElementMesh[FullRegion[2], {{0, tscale}, {0, θscale}},         (* domain *)
   MaxCellMeasure -> {"Area" -> 0.001}];
coords = g[4 Pi #1/tscale, 2 Pi #2/θscale] & @@@ dom["Coordinates"];   (* apply g *)
bmesh2 = ToBoundaryMesh[
   "Coordinates" -> coords,
   "BoundaryElements" -> dom["MeshElements"]
   ];
emesh2 = ToElementMesh@ bmesh2
(*
  ElementMesh[{{-4.86407, 1.5}, {-3.32362, 3.32362}, {-1.49991, 1.49991}},
   {TetrahedronElement["<" 5581 ">"]}]
*)

MeshRegion@ emesh2

Mathematica graphics

Note: I thought I would have to glue the boundaries together by hand, but ToBoundaryMesh seems to handle it for me. :D

$\endgroup$
11
  • $\begingroup$ Thanks! I like method 2 better. $\endgroup$
    – chris
    Jun 12, 2015 at 7:02
  • 1
    $\begingroup$ @MichaelE2...learned such a lot from this +1, of course...thanks :) $\endgroup$
    – ubpdqn
    Jun 12, 2015 at 7:24
  • $\begingroup$ @chris You're welcome. I like method 2 better, too. But I came up with it second and after some time had passed. I was afraid I'd have to remap coordinates & indices to glue the boundaries together. $\endgroup$
    – Michael E2
    Jun 12, 2015 at 14:36
  • $\begingroup$ Is there a way to achieve the same thing if there's no known equation of the surface? $\endgroup$
    – RunnyKine
    Jul 10, 2015 at 21:43
  • $\begingroup$ @RunnyKine If the surface is not given by equations, how is it given? (If it's already a valid boundary mesh, then you're in business, aren't you? If it's self-intersecting, then you're in trouble. If it's not closed, then perhaps you could stick on a patch. The ParametricPlot3D has the last problem and maybe the first.) $\endgroup$
    – Michael E2
    Jul 11, 2015 at 0:34
8
$\begingroup$

Here is an attempt to use MichaelE2's method 2 but only using built-in functions with no need to load the FEM package.

tscale = 4; θscale = 0.5;
domain = DiscretizeRegion[FullRegion[2], {{0, tscale}, {0, θscale}}, 
                          MaxCellMeasure -> {"Area" -> 0.0005}]
coords = g[4 Pi #1/tscale, 2 Pi #2/θscale] & @@@ MeshCoordinates[domain];  
         (* This is the same g defined in the question *)
mr = MeshRegion[coords, MeshCells[domain, 2]];

Here is where it gets tricky with built-in functions. I first convert the MeshRegion mr to a Graphics3D object using Show then discretize the boundary and finally use TriangulateMesh to get a solid. BoundaryMeshRegion seems to fail here, hence, the workaround.

tm = TriangulateMesh @ BoundaryDiscretizeGraphics @ Show @ mr

Mathematica graphics

MeshCellCount[tm, 3] > 0    (* check that there are 3D cells *)

True

$\endgroup$
7
  • 1
    $\begingroup$ Nice! FEM has become a bit like black magic between undocumented functions, failing ones etc… :-) $\endgroup$
    – chris
    Jul 21, 2015 at 6:36
  • $\begingroup$ @user21 Why does BoundaryMeshRegion[coords, MeshCells[domain, 2]] complain about a boundary not being closed? $\endgroup$
    – RunnyKine
    Jul 21, 2015 at 11:02
  • $\begingroup$ @RunnyKine This is the message I get : BoundaryMeshRegion::bsuncl: "The boundary surface is not closed because the edges Line[{{719,1},{97,719},{1,97}}] only come from a single face." $\endgroup$ Apr 3, 2016 at 14:11
  • $\begingroup$ @AntonAntonov. I'm guessing you're using a much newer version of Mathematica. A lot of things that worked with these region functions in previous versions are now broken. If you check this site you'll see I've found quite a few. $\endgroup$
    – RunnyKine
    Apr 5, 2016 at 4:06
  • $\begingroup$ @RunnyKine Yes, I am using the newest version. Thank you for your clarification. I like that you programmed Michael E2's approach with built-in functions. $\endgroup$ Apr 5, 2016 at 4:31
6
$\begingroup$

For a closed surface such as this one, a slight modification of the function MakeTriangleMesh[] in this answer can be used:

MakeTriangleMesh[vl_List, {closedu : (True | False) : False,
                           closedv : (True | False) : False}, opts___] := 
       Module[{dims = Most[Dimensions[vl]], v = vl, idx},
              idx = Partition[Range[Times @@ dims], Last[dims]];
              If[TrueQ[closedu],
                 v = Most[v]; idx = Append[Most[idx], First[idx]]];
              If[TrueQ[closedv],
                 v = Most /@ v; idx = Composition[Append[#, First[#]] &, Most] /@ idx];
              BoundaryMeshRegion[Apply[Join, vl], 
                                 Triangle[Flatten[Apply[{Append[Reverse[#1], Last[#2]], 
                                                         Prepend[#2, First[#1]]} &, 
                                 Partition[idx, {2, 2}, {1, 1}], {2}], 2]], opts]] /;
       ArrayQ[vl, 3]

With[{m = 101, n = 31},
     MakeTriangleMesh[N[Table[Evaluate @ g[t, θ],
                              {t, 0, 4 π, 4 π/(m - 1)}, {θ, 0, 2 π, 2 π/(n - 1)}]],
                      {True, True}]]

meshed tube

Compute the volume:

Volume[%]
   24.674785447032324
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.