TikZ – Draw a Smooth Surface

Question

I would like to draw a smooth surface like https://mathematica.stackexchange.com/questions/20912/how-to-draw-a-higher-genus-surface but using Asymptote, TikZ, Metapost, PStricks, or other TeX friendly software. The following comes out pretty rough. If I increase nx it crashes (a frequent problem for me using Asymptote).

settings.outformat = "pdf";
settings.render = 16;
size(8cm);
real R=1.2;
real r=.5;
real a=sqrt(2);
real Rsqr=R^2;
real rsqr=r^2;
real asqr=a^2;
real diffsqr=Rsqr-rsqr;
real sqrdiffsqr=diffsqr^2;
real twsumsqr=2*(rsqr+Rsqr);
real twdiffsqr=2*diffsqr;

import contour3;
currentprojection=perspective(4,6,8);
real f(real x, real y, real z) 
{ 
    real ysqr=y^2;
    real zsqr=z^2;
    real p=(-r - R + x)^2;
    real q=(r + R + x)^2;
    return
        -asqr 
        + 
        (
            sqrdiffsqr 
            - twsumsqr*(p + ysqr) 
            + twdiffsqr*zsqr
            + (p+ysqr+zsqr)^2
        )
        *
        (
            sqrdiffsqr 
            - twsumsqr*(q + ysqr) 
            + twdiffsqr*zsqr
            + (q+ysqr+zsqr)^2
        );
}

material opaq=material(diffusepen=gray(0.5), emissivepen=gray(0.4), specularpen=black);

draw(
    surface(
        contour3(
            f, 
            (-2*(r + R), -(r + R),-r - a),
            (2*(r + R), (r + R), r + a), 
            nx=60)
    ), 
    surfacepen=opaq
);

Answer 1

Interesting question. Putting two torus is not a big deal as soon as you have some 3D routines and a lot of 3D solutions can do it. The remark of Charles Staats is relevant.

I tried to make a sufficiently smooth transition between two torus. To do this I take an (x,y,z) parametrization of the left part which is possible for x belonging to [-R-a,-R], compute the interval in which y lies and then obtain the z value. The modification consists to add a x dependant function null in -R and having a null derivative in -R-a. From a mathematical point of view the result is C^1 (I think)

Second Solution. I remember the interpolation in Tchebychev points versus equidistant points. Because for x fixed the curve is closed in some parts of an ellipse, I change the y parametrization so that the points are more well-balanced (it was the reason I put sqrt(y) in the first attempt). So the solution is to take sin(pi/2*y) and then it is possible to use the smooth-spline functionality of Asymptote. Please find the (always dirty and to improve) code

size(200);
import graph3;

currentprojection=perspective(5,4,4);

real R=3;
real a=1;

triple f(pair t) {
  triple z;
  z= ((R+a*cos(t.y))*cos(t.x),(R+a*cos(t.y))*sin(t.x),a*sin(t.y));
  return z;
}
surface s=surface(f,(0,0),(2pi,2pi),Spline);
surface ss=shift((2*R+2*a,0,0))*surface(f,(0,0),(2pi,2pi),Spline);

draw(s,green,render(compression=Low,merge=true));
draw(ss,green,render(compression=Low,merge=true));

triple f11 (pair z)
{
  real y,ty,x,yy;
  triple zz;
  if ((z.y>=0) & (z.y<=1))
    {
      x=z.x+1/2*(z.x+R)^2/a;
      y=sin(pi*(z.y)/2)*sqrt((R+a)^2-x*x);
      zz=(-z.x,y,sqrt(a^2-(R-sqrt((x)^2+y^2))^2));
    }
  if ((z.y<0) & (z.y>=-1))
    {
      x=z.x+1/2*(z.x+R)^2/a;
      y=-sin(pi*abs(z.y)/2)*sqrt((R+a)^2-x*x);
      zz=(-z.x,y,sqrt(a^2-(R-sqrt((x)^2+y^2))^2));
    }
  if ((z.y>1) && (z.y<=2))
    {
      x=z.x+1/2*(z.x+R)^2/a;
      yy=2-z.y;
      y=sin(pi*(yy)/2)*sqrt((R+a)^2-x*x);
      zz=(-z.x,y,-sqrt(a^2-(R-sqrt((x)^2+y^2))^2));
    }
  if ((z.y>1) && (z.y<=2))
    {
      zz=(-z.x,y,-sqrt(a^2-(R-sqrt((x)^2+y^2))^2));
    }
  if ((z.y>=-2) && (z.y<-1))
    {
      x=z.x+1/2*(z.x+R)^2/a;
      yy=-2-z.y;
      y=-sin(pi*(abs(yy))/2)*sqrt((R+a)^2-x*x);
      zz=(-z.x,y,-sqrt(a^2-(R-sqrt((x)^2+y^2))^2));
    }
  return zz;
}

triple f21(pair z)
{
  triple z=f11(z);
  return (-z.x,z.y,z.z);
}
int N=20;
surface st=surface(f11,(-R-a,-2),(-R,2),N,2*N,Spline);
surface sstb=shift((2*R+2*a,0,0))*surface(f21,(-R-a,-2),(-R,2),N,2*N,Spline);
surface pont1=surface(st,sstb);
draw(pont1,green, render(compression=Low,merge=true));

You can observe that the two torus are drawn and a 'bridge' is added. The code is faster and the result looks like

First solution. However I cannot use the surface(...,Spline) routine to have a true smooth surface. I think it is due in fact to some artefacts of the Bézier patch parametrization (De Boor reference). Perhaps that using some monotonic option to Spline is the solution but I have to admit that I do not have the time. Im y opinion, with the standard surface(....,100,100), the result is satisfactory. I am sure that it is possible to improve it.

Please find the (dirty and to clean) code

size(200);
import graph3;

currentprojection=perspective(5,4,4);

real R=3;
real a=1;

triple f(pair t) {
  triple z;
  z= ((R+a*cos(t.y))*cos(t.x),(R+a*cos(t.y))*sin(t.x),a*sin(t.y));
  return z;
}


bool active(pair pos) {return (R+a*cos(pos.y))*cos(pos.x)<=3.2;}
bool active2(pair pos) {return (R+a*cos(pos.y))*cos(pos.x)>=-3.2;}

surface s=surface(f,(0,0),(2pi,2pi),200,200, active);
surface ss=shift((2*R+2*a,0,0))*surface(f,(0,0),(2pi,2pi),200,200,active2);

draw(s,green,render(compression=Low,merge=true));
draw(ss,green,render(compression=Low,merge=true));

triple f11 (pair z)
{
  real y,ty,x;
  if (z.y>=0)
    {
      x=z.x+1/2*(z.x+R)^2/a;
      y=sqrt(z.y)*sqrt((R+a)^2-x*x);
    }
  else
    {
      x=z.x+1/2*(z.x+R)^2/a;
      y=-sqrt(abs(z.y))*sqrt((R+a)^2-x*x);
    }
  return (-z.x,y,sqrt(a^2-(R-sqrt((x)^2+y^2))^2));
}

triple f12(pair z)
{
  triple z=f11(z);
  return (z.x,z.y,-z.z);
}

triple f21(pair z)
{
  triple z=f11(z);
  return (-z.x,z.y,z.z);
}
triple f22(pair z)
{
  triple z=f11(z);
  return (-z.x,z.y,-z.z);
}
int N=100;
surface st=surface(f11,(-R-a,-1),(-R,1),N,N);
surface sst=surface(f12,(-R-a,-1),(-R,1),N,N);
surface sstb=shift((2*R+2*a,0,0))*surface(f21,(-R-a,-1),(-R,1),N,N);
surface sstbm=shift((2*R+2*a,0,0))*surface(f22,(-R-a,-1),(-R,1),N,N);
surface pont1=surface(st,sst,sstb,sstbm);
draw(pont1,green);//blue);

and the result