From 19d25b8009801aa98ea2f46b45c37c257f990491 Mon Sep 17 00:00:00 2001 From: Norbert Preining Date: Wed, 4 Mar 2020 03:02:32 +0000 Subject: CTAN sync 202003040302 --- graphics/asymptote/base/graph3.asy | 220 ++++++++++++++++++++++++++++++++++++- 1 file changed, 218 insertions(+), 2 deletions(-) (limited to 'graphics/asymptote/base/graph3.asy') diff --git a/graphics/asymptote/base/graph3.asy b/graphics/asymptote/base/graph3.asy index ed09830974..1fe63490c9 100644 --- a/graphics/asymptote/base/graph3.asy +++ b/graphics/asymptote/base/graph3.asy @@ -1547,6 +1547,33 @@ path3[] segment(triple[] v, bool[] cond, interpolate3 join=operator --) segment.length); } +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]; + for(int j=0; j < m; ++j) { + real norm=0; + for(int i=0; i < n; ++i) + norm=max(norm,abs(a[i][j])); + real epsilon=sqrtEpsilon*norm; + 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; + for(int i=0; i < n; ++i) { + real[] ai=a[i]; + real epsilon=sqrtEpsilon*norm(ai); + if(abs(ai[0]-ai[m]) > epsilon) return false; + } + return true; +} + bool uperiodic(triple[][] a) { int n=a.length; if(n == 0) return false; @@ -1614,7 +1641,7 @@ surface surface(triple[][] f, bool[][] cond={}) if(uperiodic(f)) s.ucyclic(true); if(vperiodic(f)) s.vcyclic(true); } - + return s; } @@ -1694,10 +1721,114 @@ surface bispline(real[][] z, real[][] p, real[][] q, real[][] r, } } } - + + return s; +} + +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]; + bool[] condp=cond[i+1]; + for(int j=0; j < m; ++j) + if(all || (condi[j] && condi[j+1] && condp[j] && condp[j+1])) + ++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]; + bool[] condp=all ? null : cond[i+1]; + for(int j=0; j < m; ++j) { + if(all || (condi[j] && condi[j+1] && condp[j] && condp[j+1])) { + 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); +} + // return the surface described by a real matrix f, interpolated with // xsplinetype and ysplinetype. surface surface(real[][] f, real[] x, real[] y, @@ -1802,6 +1933,91 @@ surface surface(triple f(pair z), pair a, pair b, int nu=nmesh, int nv=nu, return surface(v,active); } +// 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]); + if(!all) activei[j]=cond(z); + triple f=f(z); + 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; +} + // return the surface described by a parametric function f over box(a,b), // interpolated with usplinetype and vsplinetype. surface surface(triple f(pair z), pair a, pair b, int nu=nmesh, int nv=nu, -- cgit v1.2.3