diff options
Diffstat (limited to 'Build/source/utils/asymptote/base/three_tube.asy')
-rw-r--r-- | Build/source/utils/asymptote/base/three_tube.asy | 223 |
1 files changed, 134 insertions, 89 deletions
diff --git a/Build/source/utils/asymptote/base/three_tube.asy b/Build/source/utils/asymptote/base/three_tube.asy index c79ae3f619f..2eb5e4cf716 100644 --- a/Build/source/utils/asymptote/base/three_tube.asy +++ b/Build/source/utils/asymptote/base/three_tube.asy @@ -1,5 +1,6 @@ -void render(path3 s, real granularity=linegranularity, void f(path3, real)) +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) { @@ -65,8 +66,8 @@ rmf[] rmf(path3 g, real[] t) return R; } -surface bispline(real[][] z, real[][] p, real[][] q, real[][] r, - real[] x, real[] y, bool[][] cond={}) +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 @@ -87,62 +88,49 @@ surface bispline(real[][] z, real[][] p, real[][] q, real[][] r, } } - surface s=surface(count); - s.index=new int[n][m]; - int k=-1; + real[][][] s=new real[count][][]; + int k=0; for(int i=0; i < n; ++i) { + int ip=i+1; bool[] condi=all ? null : cond[i]; real xi=x[i]; real[] zi=z[i]; - real[] zp=z[i+1]; + real[] zp=z[ip]; real[] ri=r[i]; - real[] rp=r[i+1]; + real[] rp=r[ip]; real[] pi=p[i]; - real[] pp=p[i+1]; + real[] pp=p[ip]; real[] qi=q[i]; - real[] qp=q[i+1]; - real xp=x[i+1]; + real[] qp=q[ip]; + real xp=x[ip]; real hx=(xp-xi)/3; - int[] indexi=s.index[i]; for(int j=0; j < m; ++j) { real yj=y[j]; - real yp=y[j+1]; + int jp=j+1; + real yp=y[jp]; if(all || condi[j]) { - triple[][] P=array(4,array(4,O)); real hy=(yp-yj)/3; real hxy=hx*hy; - // x and y directions - for(int k=0; k < 4; ++k) { - P[0][k] += xi*X; - P[k][0] += yj*Y; - P[1][k] += (xp+2*xi)/3*X; - P[k][1] += (yp+2*yj)/3*Y; - P[2][k] += (2*xp+xi)/3*X; - P[k][2] += (2*yp+yj)/3*Y; - P[3][k] += xp*X; - P[k][3] += yp*Y; - } - // z: value - P[0][0] += zi[j]*Z; - P[3][0] += zp[j]*Z; - P[0][3] += zi[j+1]*Z; - P[3][3] += zp[j+1]*Z; - // z: first derivative - P[1][0] += (P[0][0].z+hx*pi[j])*Z; - P[1][3] += (P[0][3].z+hx*pi[j+1])*Z; - P[2][0] += (P[3][0].z-hx*pp[j])*Z; - P[2][3] += (P[3][3].z-hx*pp[j+1])*Z; - P[0][1] += (P[0][0].z+hy*qi[j])*Z; - P[3][1] += (P[3][0].z+hy*qp[j])*Z; - P[0][2] += (P[0][3].z-hy*qi[j+1])*Z; - P[3][2] += (P[3][3].z-hy*qp[j+1])*Z; - // z: second derivative - P[1][1] += (P[0][1].z+P[1][0].z-P[0][0].z+hxy*ri[j])*Z; - P[1][2] += (P[0][2].z+P[1][3].z-P[0][3].z-hxy*ri[j+1])*Z; - P[2][1] += (P[2][0].z+P[3][1].z-P[3][0].z-hxy*rp[j])*Z; - P[2][2] += (P[2][3].z+P[3][2].z-P[3][3].z+hxy*rp[j+1])*Z; - s.s[++k]=patch(P); - indexi[j]=k; + 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; } } } @@ -150,11 +138,11 @@ surface bispline(real[][] z, real[][] p, real[][] q, real[][] r, return s; } -// return the surface described by a real matrix f, interpolated with +// return the surface values described by a real matrix f, interpolated with // xsplinetype and ysplinetype. -surface surface(real[][] f, real[] x, real[] y, - splinetype xsplinetype=null, splinetype ysplinetype=xsplinetype, - bool[][] cond={}) +real[][][] bispline(real[][] f, real[] x, real[] y, + splinetype xsplinetype=null, + splinetype ysplinetype=xsplinetype, bool[][] cond={}) { real epsilon=sqrtEpsilon*norm(y); if(xsplinetype == null) @@ -176,9 +164,7 @@ surface surface(real[][] f, real[] x, real[] y, real[][] p=transpose(tp); for(int i=0; i < n; ++i) r[i]=clamped(d1[i],d2[i])(y,p[i]); - surface s=bispline(f,p,q,r,x,y,cond); - if(xsplinetype == periodic) s.ucyclic(true); - if(ysplinetype == periodic) s.vcyclic(true); + real[][][] s=bispline0(f,p,q,r,x,y,cond); return s; } @@ -248,11 +234,11 @@ surface surface(triple f(pair z), real[] u, real[] v, vperiodic(fz) ? periodic : notaknot}; } else if(vsplinetype.length != 3) abort("vsplinetype must have length 3"); - surface sx=surface(fx,ipt,jpt,usplinetype[0],vsplinetype[0],active); - surface sy=surface(fy,ipt,jpt,usplinetype[1],vsplinetype[1],active); - surface sz=surface(fz,ipt,jpt,usplinetype[2],vsplinetype[2],active); + 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.s.length); + surface s=surface(sx.length); s.index=new int[nu][nv]; int k=-1; for(int i=0; i < nu; ++i) { @@ -261,12 +247,20 @@ surface surface(triple f(pair z), real[] u, real[] v, indexi[j]=++k; } - for(int k=0; k < sx.s.length; ++k) { + for(int k=0; k < sx.length; ++k) { triple[][] Q=new triple[4][]; - for(int i=0; i < 4 ; ++i) - Q[i]=sequence(new triple(int j) { - return (sx.s[k].P[i][j].z,sy.s[k].P[i][j].z,sz.s[k].P[i][j].z); - },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); } @@ -282,8 +276,8 @@ surface surface(triple f(pair z), real[] u, real[] v, 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), + 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), @@ -294,12 +288,15 @@ path3 interp(path3 a, path3 b, real t) struct tube { surface s; - path3 center; + path3 center; // tube axis - void operator init(path3 p, real width, int sectors=4, - real granularity=linegranularity) { - sectors += sectors % 2; // Must be even. - int h=quotient(sectors,2); + 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 tube(path3, path3)=null) { real r=0.5*width; void generate(path3 p) { @@ -308,31 +305,41 @@ struct tube for(int i=0; i < n; ++i) { triple v=point(p,i); triple u=point(p,i+1)-v; - s.append(shift(v)*align(unit(u))*scale(r,r,abs(u))*unitcylinder); + 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),granularity, + 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)]; - return point(G,t.x)+r*(R.r*cos(t.y)-R.s*sin(t.y)); + 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]); } - real[] v=uniform(0,2pi,sectors); + 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[] Periodic={periodic,periodic,periodic}; + static splinetype[] Circular={circular,circular,circular}; if(T.length > 0) { - surface S=surface(f,sequence(T.length),v,Monotonic,Periodic); + surface S=surface(f,sequence(T.length),v,Monotonic,Circular); s.append(S); // Compute center of tube: @@ -341,19 +348,52 @@ struct tube triple[] pre=new triple[n+1]; triple[] point=new triple[n+1]; triple[] post=new triple[n+1]; + int[] index=S.index[0]; - pre[0]=point[0]=0.5*(S.s[index[0]].P[0][0]+S.s[index[h]].P[0][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 [][] P=S.s[index[0]].P; - triple [][] Q=S.s[index[h]].P; - post[i]=0.5*(P[1][0]+Q[1][0]); - pre[i+1]=0.5*(P[2][0]+Q[2][0]); - point[i+1]=0.5*(P[3][0]+Q[3][0]); + 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]; - post[n]=0.5*(S.s[index[0]].P[3][0]+S.s[index[h]].P[3][0]); - center=center&path3(pre,point,post,array(n+1,false),T.cyclic); + 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(tube != 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]; + tube(Center,path3(pre,point,post,b,T.cyclic)); + } } } } @@ -363,11 +403,16 @@ struct tube 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)); - s.append(shift(point(p,i))*t*align(dir(p,i,-1))*unithemisphere); - begin=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)); } } |