Triangulation of surfaces

Ivan Cnop
Vrije universiteit Brussel
icnop@vub.ac.be

Problem

If we plot three dimensional surfaces, we may get bad results near discontinuities

f(x_, 
        y_) := -((10 x)/(Sqrt[(x^2) + (y^2)]))

<< Default3D`

Plot3D[f[x, y], {x, -3, 3}, {y, -3, 3}]

Increasing Plotpoints doesn't help:

Plot3D[f[x, y], {x, -3, 3}, {y, -3, 3}, 
      PlotPoints -> 31]

Plot3D[f[x, y], {x, -3, 3}, {y, -3, 3}, 
      PlotPoints -> 101]

Spinning

<< RealTime3D`

Plot3D[f[x, y], {x, -3, 3}, {y, -3, 3}, 
      PlotPoints -> 51]

The problem is that we only want more plotpoints near the discontinuity.

Solution: nonlinear grids

Basis: a  triangular grid

a = ( .5 + (Sqrt[3])/2  I)

pts = 
        Table[{ .2 k +  .2 a l}, {k, -20, 0}, {l, 0, 
            20}];

coords[z_] := {Re[z], Im[z]};

ListPlot[Map[coords, Flatten[pts]], 
      AspectRatio -> Automatic]

More dense around the origin

\(\(pts2 = Map[#\^3 &, pts];\)\)

ListPlot[Map[coords, Flatten[pts2]], 
      AspectRatio -> Automatic]

Selecting triangles

Flatten[Map[RotateLeft, pts]];

Flatten[RotateLeft[pts]];

tri = 
        Transpose[{Flatten[pts], Flatten[Map[RotateLeft, pts]], 
            Flatten[RotateLeft[pts]]}];

tri2 = 
        Transpose[{Flatten[pts2], Flatten[Map[RotateLeft, pts2]], 
            Flatten[RotateLeft[pts2]]}];

nearQ[{z1_, z2_, 
          z3_}] := (Abs[z1 - z2] + Abs[z3 - z2] + Abs[z1 - z3] < 
          20)

N[Select[tri, nearQ]];

N[Select[tri2, nearQ]];

polygon[{z1_, z2_, z3_}] := 
      Polygon[Map[coords, {z1, z2, z3}]]

Show[Graphics[Map[polygon, Select[tri2, nearQ]]]]

Using only half of the triangles will allow us to "look through" the surface

Computing values

Triples of points

f[z_] := f[Re[z], Im[z]]

f[{z1_, z2_, z3_}] := {f[z1], f[z2], f[z3]}

f[x_, y_] := -10 x /Sqrt[(x^2) + (y^2)]

Fixed viewpoint

<< Default3D`

polygon3D[{z1_, z2_, z3_}] := 
      Polygon[{{Re[z1], Im[z1], f[Re[z1], Im[z1]]}, {Re[z2], Im[z2], 
            f[Re[z2], Im[z2]]}, {Re[z3], Im[z3], 
            f[Re[z3], Im[z3]]}}]

Show[Graphics3D[Map[polygon3D, Select[N[tri2], nearQ]]], 
      ViewPoint -> {2, 3,  .5}, AmbientLight  ->  Hue[0]]

Spinning the result

<< RealTime3D`

Show[Graphics3D[Map[polygon3D, Select[N[tri2], nearQ]]]]

Other examples

f[x_, y_] := - 20 x  y/((x^2) + (y^2))

polygon3D[{z1_, z2_, z3_}] := 
      Polygon[{{Re[z1], Im[z1], f[Re[z1], Im[z1]]}, {Re[z2], Im[z2], 
            f[Re[z2], Im[z2]]}, {Re[z3], Im[z3], 
            f[Re[z3], Im[z3]]}}]

Show[Graphics3D[Map[polygon3D, Select[N[tri2], nearQ]]]]

Applications

If we have a numerical problem with strange behavior around a point, we can make grid points more dense around this point and obtain more reliable computations and graphics.