diff options
Diffstat (limited to 'Master/texmf-dist/asymptote/three_tube.asy')
-rw-r--r-- | Master/texmf-dist/asymptote/three_tube.asy | 417 |
1 files changed, 417 insertions, 0 deletions
diff --git a/Master/texmf-dist/asymptote/three_tube.asy b/Master/texmf-dist/asymptote/three_tube.asy new file mode 100644 index 00000000000..aa5688f081e --- /dev/null +++ b/Master/texmf-dist/asymptote/three_tube.asy @@ -0,0 +1,417 @@ +void render(path3 s, void f(path3, real), render render=defaultrender) +{ + real granularity=render.tubegranularity; + void Split(triple z0, triple c0, triple c1, triple z1, real t0=0, real t1=1, + real depth=mantissaBits) { + if(depth > 0) { + real S=straightness(z0,c0,c1,z1); + if(S > 0) { + --depth; + if(S > max(granularity*max(abs(z0),abs(c0),abs(c1),abs(z1)))) { + triple m0=0.5*(z0+c0); + triple m1=0.5*(c0+c1); + triple m2=0.5*(c1+z1); + triple m3=0.5*(m0+m1); + triple m4=0.5*(m1+m2); + triple m5=0.5*(m3+m4); + real tm=0.5*(t0+t1); + Split(z0,m0,m3,m5,t0,tm,depth); + Split(m5,m4,m2,z1,tm,t1,depth); + return; + } + } + } + f(z0..controls c0 and c1..z1,t0); + } + Split(point(s,0),postcontrol(s,0),precontrol(s,1),point(s,1)); +} + +struct rmf +{ + triple p,r,t,s; + void operator init(triple p, triple r, triple t) + { + this.p=p; + this.r=r; + this.t=t; + s=cross(t,r); + } +} + +// Rotation minimizing frame +// http://www.cs.hku.hk/research/techreps/document/TR-2007-07.pdf +rmf[] rmf(path3 g, real[] t) +{ + rmf[] R=new rmf[t.length]; + triple d=dir(g,0); + R[0]=rmf(point(g,0),perp(d),d); + for(int i=1; i < t.length; ++i) { + rmf Ri=R[i-1]; + real t=t[i]; + triple p=point(g,t); + triple v1=p-Ri.p; + if(v1 != O) { + triple r=Ri.r; + triple u1=unit(v1); + triple ti=Ri.t; + triple tp=ti-2*dot(u1,ti)*u1; + ti=dir(g,t); + triple rp=r-2*dot(u1,r)*u1; + triple u2=unit(ti-tp); + rp=rp-2*dot(u2,rp)*u2; + R[i]=rmf(p,unit(rp),unit(ti)); + } else + R[i]=R[i-1]; + } + return R; +} + +private real[][][] bispline0(real[][] z, real[][] p, real[][] q, real[][] r, + real[] x, real[] y, bool[][] cond={}) +{ // z[i][j] is the value at (x[i],y[j]) + // p and q are the first derivatives with respect to x and y, respectively + // r is the second derivative ddu/dxdy + int n=x.length-1; + int m=y.length-1; + + bool all=cond.length == 0; + + int count; + if(all) + count=n*m; + else { + count=0; + for(int i=0; i < n; ++i) { + bool[] condi=cond[i]; + for(int j=0; j < m; ++j) + if(condi[j]) ++count; + } + } + + real[][][] s=new real[count][][]; + int k=0; + for(int i=0; i < n; ++i) { + int ip=i+1; + real xi=x[i]; + real xp=x[ip]; + real hx=(xp-xi)/3; + real[] zi=z[i]; + real[] zp=z[ip]; + real[] ri=r[i]; + real[] rp=r[ip]; + real[] pi=p[i]; + real[] pp=p[ip]; + real[] qi=q[i]; + real[] qp=q[ip]; + bool[] condi=all ? null : cond[i]; + for(int j=0; j < m; ++j) { + if(all || condi[j]) { + real yj=y[j]; + int jp=j+1; + real yp=y[jp]; + real hy=(yp-yj)/3; + real hxy=hx*hy; + real zij=zi[j]; + real zip=zi[jp]; + real zpj=zp[j]; + real zpp=zp[jp]; + real pij=hx*pi[j]; + real ppj=hx*pp[j]; + real qip=hy*qi[jp]; + real qpp=hy*qp[jp]; + real zippip=zip+hx*pi[jp]; + real zppmppp=zpp-hx*pp[jp]; + real zijqij=zij+hy*qi[j]; + real zpjqpj=zpj+hy*qp[j]; + + s[k]=new real[][] {{zij,zijqij,zip-qip,zip}, + {zij+pij,zijqij+pij+hxy*ri[j], + zippip-qip-hxy*ri[jp],zippip}, + {zpj-ppj,zpjqpj-ppj-hxy*rp[j], + zppmppp-qpp+hxy*rp[jp],zppmppp}, + {zpj,zpjqpj,zpp-qpp,zpp}}; + ++k; + } + } + } + + return s; +} + +// return the surface values described by a real matrix f, interpolated with +// xsplinetype and ysplinetype. +real[][][] bispline(real[][] f, real[] x, real[] y, + splinetype xsplinetype=null, + splinetype ysplinetype=xsplinetype, bool[][] cond={}) +{ + real epsilon=sqrtEpsilon*norm(y); + if(xsplinetype == null) + xsplinetype=(abs(x[0]-x[x.length-1]) <= epsilon) ? periodic : notaknot; + if(ysplinetype == null) + ysplinetype=(abs(y[0]-y[y.length-1]) <= epsilon) ? periodic : notaknot; + int n=x.length; int m=y.length; + real[][] ft=transpose(f); + real[][] tp=new real[m][]; + for(int j=0; j < m; ++j) + tp[j]=xsplinetype(x,ft[j]); + real[][] q=new real[n][]; + for(int i=0; i < n; ++i) + q[i]=ysplinetype(y,f[i]); + real[][] qt=transpose(q); + real[] d1=xsplinetype(x,qt[0]); + real[] d2=xsplinetype(x,qt[m-1]); + real[][] r=new real[n][]; + real[][] p=transpose(tp); + for(int i=0; i < n; ++i) + r[i]=clamped(d1[i],d2[i])(y,p[i]); + return bispline0(f,p,q,r,x,y,cond); +} + +bool uperiodic(real[][] a) { + int n=a.length; + if(n == 0) return false; + int m=a[0].length; + real[] a0=a[0]; + real[] a1=a[n-1]; + real epsilon=sqrtEpsilon*norm(a); + for(int j=0; j < m; ++j) + if(abs(a0[j]-a1[j]) > epsilon) return false; + return true; +} +bool vperiodic(real[][] a) { + int n=a.length; + if(n == 0) return false; + int m=a[0].length-1; + real epsilon=sqrtEpsilon*norm(a); + for(int i=0; i < n; ++i) + if(abs(a[i][0]-a[i][m]) > epsilon) return false; + return true; +} + +// return the surface described by a parametric function f evaluated at u and v +// and interpolated with usplinetype and vsplinetype. +surface surface(triple f(pair z), real[] u, real[] v, + splinetype[] usplinetype, splinetype[] vsplinetype=Spline, + bool cond(pair z)=null) +{ + int nu=u.length-1; + int nv=v.length-1; + real[] ipt=sequence(u.length); + real[] jpt=sequence(v.length); + real[][] fx=new real[u.length][v.length]; + real[][] fy=new real[u.length][v.length]; + real[][] fz=new real[u.length][v.length]; + + bool[][] active; + bool all=cond == null; + if(!all) active=new bool[u.length][v.length]; + + for(int i=0; i <= nu; ++i) { + real ui=u[i]; + real[] fxi=fx[i]; + real[] fyi=fy[i]; + real[] fzi=fz[i]; + bool[] activei=all ? null : active[i]; + for(int j=0; j <= nv; ++j) { + pair z=(ui,v[j]); + triple f=(all || (activei[j]=cond(z))) ? f(z) : O; + fxi[j]=f.x; + fyi[j]=f.y; + fzi[j]=f.z; + } + } + + if(usplinetype.length == 0) { + usplinetype=new splinetype[] {uperiodic(fx) ? periodic : notaknot, + uperiodic(fy) ? periodic : notaknot, + uperiodic(fz) ? periodic : notaknot}; + } else if(usplinetype.length != 3) abort("usplinetype must have length 3"); + + if(vsplinetype.length == 0) { + vsplinetype=new splinetype[] {vperiodic(fx) ? periodic : notaknot, + vperiodic(fy) ? periodic : notaknot, + vperiodic(fz) ? periodic : notaknot}; + } else if(vsplinetype.length != 3) abort("vsplinetype must have length 3"); + + real[][][] sx=bispline(fx,ipt,jpt,usplinetype[0],vsplinetype[0],active); + real[][][] sy=bispline(fy,ipt,jpt,usplinetype[1],vsplinetype[1],active); + real[][][] sz=bispline(fz,ipt,jpt,usplinetype[2],vsplinetype[2],active); + + surface s=surface(sx.length); + s.index=new int[nu][nv]; + int k=-1; + for(int i=0; i < nu; ++i) { + int[] indexi=s.index[i]; + for(int j=0; j < nv; ++j) + indexi[j]=++k; + } + + for(int k=0; k < sx.length; ++k) { + triple[][] Q=new triple[4][]; + real[][] Px=sx[k]; + real[][] Py=sy[k]; + real[][] Pz=sz[k]; + for(int i=0; i < 4 ; ++i) { + real[] Pxi=Px[i]; + real[] Pyi=Py[i]; + real[] Pzi=Pz[i]; + Q[i]=new triple[] {(Pxi[0],Pyi[0],Pzi[0]), + (Pxi[1],Pyi[1],Pzi[1]), + (Pxi[2],Pyi[2],Pzi[2]), + (Pxi[3],Pyi[3],Pzi[3])}; + } + s.s[k]=patch(Q); + } + + if(usplinetype[0] == periodic && usplinetype[1] == periodic && + usplinetype[1] == periodic) s.ucyclic(true); + + if(vsplinetype[0] == periodic && vsplinetype[1] == periodic && + vsplinetype[1] == periodic) s.vcyclic(true); + + return s; +} + +path3 interp(path3 a, path3 b, real t) +{ + int n=size(a); + return path3(sequence(new triple(int i) { + return interp(precontrol(a,i),precontrol(b,i),t);},n), + sequence(new triple(int i) {return interp(point(a,i),point(b,i),t);},n), + sequence(new triple(int i) {return interp(postcontrol(a,i), + postcontrol(b,i),t);},n), + sequence(new bool(int i) {return straight(a,i) && straight(b,i);},n), + cyclic(a) && cyclic(b)); +} + +struct tube +{ + surface s; + path3 center; // tube axis + + void Null(transform3) {} + void Null(transform3, bool) {} + + void operator init(path3 p, real width, render render=defaultrender, + void cylinder(transform3)=Null, + void sphere(transform3, bool half)=Null, + void pipe(path3, path3)=null) { + real r=0.5*width; + + void generate(path3 p) { + int n=length(p); + if(piecewisestraight(p)) { + for(int i=0; i < n; ++i) { + triple v=point(p,i); + triple u=point(p,i+1)-v; + transform3 t=shift(v)*align(unit(u))*scale(r,r,abs(u)); + s.append(t*unitcylinder); + cylinder(t); + } + center=center&p; + } else { + real[] T; + path3 G; + for(int i=0; i < n; ++i) + render(subpath(p,i,i+1), + new void(path3 g, real s) { + G=G&g; + T.push(i+s); + },render); + T.push(n); + T.cyclic=cyclic(p); + rmf[] rmf=rmf(p,T); + triple f(pair t) { + rmf R=rmf[round(t.x)]; + int n=round(t.y); + static real[] x={1,0,-1,0}; + static real[] y={0,1,0,-1}; + return point(G,t.x)+r*(R.r*x[n]-R.s*y[n]); + } + + static real[] v={0,1,2,3,0}; + static real[] circular(real[] x, real[] y) { + static real a=8/3*(sqrt(2)-1); + return a*periodic(x,y); + } + + static splinetype[] Monotonic={monotonic,monotonic,monotonic}; + static splinetype[] Circular={circular,circular,circular}; + if(T.length > 0) { + surface S=surface(f,sequence(T.length),v,Monotonic,Circular); + s.append(S); + + // Compute center of tube: + int n=S.index.length; + if(T.cyclic) --n; + triple[] pre=new triple[n+1]; + triple[] point=new triple[n+1]; + triple[] post=new triple[n+1]; + + int[] index=S.index[0]; + triple Point; + for(int m=0; m < 4; ++m) + Point += S.s[index[m]].P[0][0]; + pre[0]=point[0]=0.25*Point; + + for(int i=0; i < n; ++i) { + index=S.index[i]; + triple Pre,Point,Post; + for(int m=0; m < 4; ++m) { + triple [][] P=S.s[index[m]].P; + Post += P[1][0]; + Pre += P[2][0]; + Point += P[3][0]; + } + post[i]=0.25*Post; + pre[i+1]=0.25*Pre; + point[i+1]=0.25*Point; + + } + + index=S.index[n-1]; + triple Post; + for(int m=0; m < 4; ++m) + Post += S.s[index[m]].P[3][0]; + post[n]=0.25*Post; + + bool[] b=array(n+1,false); + path3 Center=path3(pre,point,post,b,T.cyclic); + center=center&Center; + + if(pipe != null) { // Compute path along tube + triple[] pre=new triple[n+1]; + triple[] point=new triple[n+1]; + triple[] post=new triple[n+1]; + pre[0]=point[0]=S.s[S.index[0][0]].P[0][0]; + for(int i=0; i < n; ++i) { + triple [][] P=S.s[S.index[i][0]].P; + post[i]=P[1][0]; + pre[i+1]=P[2][0]; + point[i+1]=P[3][0]; + } + post[n]=S.s[S.index[n-1][0]].P[3][0]; + pipe(Center,path3(pre,point,post,b,T.cyclic)); + } + } + } + } + + transform3 t=scale3(r); + bool cyclic=cyclic(p); + int begin=0; + int n=length(p); + for(int i=cyclic ? 0 : 1; i < n; ++i) + if(abs(dir(p,i,1)-dir(p,i,-1)) > sqrtEpsilon) { + generate(subpath(p,begin,i)); + triple dir=dir(p,i,-1); + s.append(shift(point(p,i))*t*align(dir)* + (dir != O ? unithemisphere : unitsphere)); + int L=length(center); + sphere(shift(point(center,L))*t*align(dir(center,L,-1)), + half=straight(p,i-1) && straight(p,i)); + begin=i; + } + generate(subpath(p,begin,n)); + } +} |