summaryrefslogtreecommitdiff
path: root/graphics/asymptote/base/graph3.asy
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2020-03-04 03:02:32 +0000
committerNorbert Preining <norbert@preining.info>2020-03-04 03:02:32 +0000
commit19d25b8009801aa98ea2f46b45c37c257f990491 (patch)
treec0da70fbab26feff40e87df289ae797917df8419 /graphics/asymptote/base/graph3.asy
parentb78f5a6705512314d6fe42fd5144c5f8e5d484ce (diff)
CTAN sync 202003040302
Diffstat (limited to 'graphics/asymptote/base/graph3.asy')
-rw-r--r--graphics/asymptote/base/graph3.asy220
1 files changed, 218 insertions, 2 deletions
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,