diff options
Diffstat (limited to 'Build/source/utils/asymptote/base')
-rw-r--r-- | Build/source/utils/asymptote/base/plain_prethree.asy | 15 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/rational.asy | 73 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/rationalSimplex.asy | 222 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/shaders/fragment.glsl | 85 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/shaders/vertex.glsl | 34 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/simplex.asy | 142 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/three.asy | 29 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/three_light.asy | 27 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/three_surface.asy | 38 |
9 files changed, 485 insertions, 180 deletions
diff --git a/Build/source/utils/asymptote/base/plain_prethree.asy b/Build/source/utils/asymptote/base/plain_prethree.asy index cfc0c3c7411..0d3921e3f0d 100644 --- a/Build/source/utils/asymptote/base/plain_prethree.asy +++ b/Build/source/utils/asymptote/base/plain_prethree.asy @@ -165,7 +165,6 @@ struct light { real[][] specular; pen background=nullpen; // Background color of the 3D canvas. real specularfactor; - bool viewport; // Are the lights specified (and fixed) in the viewport frame? triple[] position; // Only directional lights are currently implemented. transform3 T=identity(4); // Transform to apply to normal vectors. @@ -176,7 +175,7 @@ struct light { pen[] ambient=array(diffuse.length,black), pen[] specular=diffuse, pen background=nullpen, real specularfactor=1, - bool viewport=false, triple[] position) { + triple[] position) { int n=diffuse.length; assert(ambient.length == n && specular.length == n && position.length == n); @@ -192,21 +191,18 @@ struct light { this.position[i]=unit(position[i]); } this.specularfactor=specularfactor; - this.viewport=viewport; } void operator init(pen diffuse=white, pen ambient=black, pen specular=diffuse, - pen background=nullpen, real specularfactor=1, - bool viewport=false...triple[] position) { + pen background=nullpen, real specularfactor=1 ...triple[] position) { int n=position.length; operator init(array(n,diffuse),array(n,ambient),array(n,specular), - background,specularfactor,viewport,position); + background,specularfactor,position); } void operator init(pen diffuse=white, pen ambient=black, pen specular=diffuse, - pen background=nullpen, bool viewport=false, - real x, real y, real z) { - operator init(diffuse,ambient,specular,background,viewport,(x,y,z)); + pen background=nullpen, real x, real y, real z) { + operator init(diffuse,ambient,specular,background,(x,y,z)); } void operator init(explicit light light) { @@ -215,7 +211,6 @@ struct light { specular=copy(light.specular); background=light.background; specularfactor=light.specularfactor; - viewport=light.viewport; position=copy(light.position); } diff --git a/Build/source/utils/asymptote/base/rational.asy b/Build/source/utils/asymptote/base/rational.asy index a142063806d..d24da0d6645 100644 --- a/Build/source/utils/asymptote/base/rational.asy +++ b/Build/source/utils/asymptote/base/rational.asy @@ -133,6 +133,26 @@ bool operator ==(rational[][] r, rational[][] s) return all(sequence(new bool(int i) {return r[i] == s[i];},r.length)); } +bool[] operator <(rational[] r, rational s) +{ + return sequence(new bool(int i) {return r[i] < s;},r.length); +} + +bool[] operator >(rational[] r, rational s) +{ + return sequence(new bool(int i) {return r[i] > s;},r.length); +} + +bool[] operator <=(rational[] r, rational s) +{ + return sequence(new bool(int i) {return r[i] <= s;},r.length); +} + +bool[] operator >=(rational[] r, rational s) +{ + return sequence(new bool(int i) {return r[i] >= s;},r.length); +} + rational min(rational a, rational b) { return a <= b ? a : b; @@ -143,35 +163,47 @@ rational max(rational a, rational b) return a >= b ? a : b; } +string string(rational r) +{ + return r.q == 1 ? string(r.p) : string(r.p)+"/"+string(r.q); +} + +string texstring(rational r) +{ + if(r.q == 1) return string(r.p); + string s; + if(r.p < 0) s="-"; + return s+"\frac{"+string(abs(r.p))+"}{"+string(r.q)+"}"; +} -void write(string s="", rational r, suffix suffix=endl) { - if(r.q == 1) - write(s+string(r.p),suffix); - else - write(s+string(r.p)+"/"+string(r.q),suffix); +void write(file fout=stdout, string s="", rational r, suffix suffix=none) +{ + write(fout,s+string(r),suffix); } -void write(string s="", rational[] a, suffix suffix=endl) { +void write(file fout=stdout, string s="", rational[] a, suffix suffix=none) +{ if(s != "") - write(s); + write(fout,s); for(int i=0; i < a.length; ++i) { - write(i,none); - write(':\t',a[i]); + write(fout,i,none); + write(fout,':\t',a[i],endl); } - write(suffix); + write(fout,suffix); } -void write(string s="", rational[][] a, suffix suffix=endl) { +void write(file fout=stdout, string s="", rational[][] a, suffix suffix=none) +{ if(s != "") - write(s); + write(fout,s); for(int i=0; i < a.length; ++i) { rational[] ai=a[i]; for(int j=0; j < ai.length; ++j) { - write(ai[j],tab); + write(fout,ai[j],tab); } - write(); + write(fout,endl); } - write(suffix); + write(fout,suffix); } bool rectangular(rational[][] m) @@ -203,9 +235,18 @@ rational max(rational[] a) rational abs(rational r) { - return rational(abs(r.p),r.q,false); + return rational(abs(r.p),r.q,false); +} + +rational[] operator -(rational[] r) +{ + return sequence(new rational(int i) {return -r[i];},r.length); } +rational[][] rationalidentity(int n) +{ + return sequence(new rational[](int i) {return sequence(new rational(int j) {return j == i ? 1 : 0;},n);},n); +} /* rational r=rational(1,3)+rational(1,4); diff --git a/Build/source/utils/asymptote/base/rationalSimplex.asy b/Build/source/utils/asymptote/base/rationalSimplex.asy index fd7128b8025..5175b87ba74 100644 --- a/Build/source/utils/asymptote/base/rationalSimplex.asy +++ b/Build/source/utils/asymptote/base/rationalSimplex.asy @@ -1,6 +1,29 @@ // Rational simplex solver written by John C. Bowman and Pouria Ramazi, 2018. import rational; +void simplexTableau(rational[][] E, int[] Bindices, int I=-1, int J=-1) {} +void simplexPhase2() {} + +void simplexWrite(rational[][] E, int[] Bindicies, int, int) +{ + int m=E.length-1; + int n=E[0].length-1; + + write(E[m][n],tab); + for(int j=0; j < n; ++j) + write(E[m][j],tab); + write(); + + for(int i=0; i < m; ++i) { + write(E[i][n],tab); + for(int j=0; j < n; ++j) { + write(E[i][j],tab); + } + write(); + } + write(); +}; + struct simplex { static int OPTIMAL=0; static int UNBOUNDED=1; @@ -49,7 +72,7 @@ struct simplex { if(Em[J] < 0) break; if(J == N) - return 0; + break; int I=-1; rational M; @@ -65,28 +88,74 @@ struct simplex { rational e=E[i][J]; if(e > 0) { rational v=E[i][N]/e; - if(v <= M) {M=v; I=i;} + if(v < M) {M=v; I=i;} // Bland's rule: choose smallest argmin } } if(I == -1) return UNBOUNDED; // Can only happen in Phase 2. + simplexTableau(E,Bindices,I,J); + + // Generate new tableau Bindices[I]=J; + rowreduce(E,N,I,J); + } + return OPTIMAL; + } + + int iterateDual(rational[][] E, int N, int[] Bindices) { + while(true) { + // Find first negative entry in right (basic variable) column + rational[] Em=E[m]; + int I; + for(I=0; I < m; ++I) { + if(E[I][N] < 0) break; + } + + if(I == m) + break; + + int J=-1; + rational M; + for(int j=0; j < N; ++j) { + rational e=E[I][j]; + if(e < 0) { + M=-E[m][j]/e; + J=j; + break; + } + } + for(int j=J+1; j < N; ++j) { + rational e=E[I][j]; + if(e < 0) { + rational v=-E[m][j]/e; + if(v < M) {M=v; J=j;} // Bland's rule: choose smallest argmin + } + } + if(J == -1) + return INFEASIBLE; // Can only happen in Phase 2. + + simplexTableau(E,Bindices,I,J); // Generate new tableau + Bindices[I]=J; rowreduce(E,N,I,J); } - return 0; + return OPTIMAL; } // Try to find a solution x to Ax=b that minimizes the cost c^T x, - // where A is an m x n matrix, x is a vector of length n, b is a - // vector of length m, and c is a vector of length n. + // where A is an m x n matrix, x is a vector of n non-negative numbers, + // b is a vector of length m, and c is a vector of length n. + // Can set phase1=false if the last m columns of A form the identity matrix. void operator init(rational[] c, rational[][] A, rational[] b, - bool phase1=true) { - // Phase 1 + bool phase1=true, bool dual=false) { + if(dual) phase1=false; + // Phase 1 m=A.length; + if(m == 0) {case=INFEASIBLE; return;} n=A[0].length; + if(n == 0) {case=INFEASIBLE; return;} int N=phase1 ? n+m : n; rational[][] E=new rational[m+1][N+1]; @@ -98,7 +167,7 @@ struct simplex { for(int i=0; i < m; ++i) { rational[] Ai=A[i]; rational[] Ei=E[i]; - if(b[i] >= 0) { + if(b[i] >= 0 || dual) { for(int j=0; j < n; ++j) { rational Aij=Ai[j]; Ei[j]=Aij; @@ -126,7 +195,7 @@ struct simplex { rational sum=0; for(int i=0; i < m; ++i) { - rational B=abs(b[i]); + rational B=dual ? b[i] : abs(b[i]); E[i][N]=B; sum -= B; } @@ -136,27 +205,45 @@ struct simplex { for(int j=0; j < m; ++j) Em[n+j]=0; - int[] Bindices=sequence(new int(int x){return x;},m)+n; + int[] Bindices; if(phase1) { + Bindices=sequence(new int(int x){return x;},m)+n; iterate(E,N,Bindices); if(Em[J] != 0) { + simplexTableau(E,Bindices); case=INFEASIBLE; return; } - } - + } else Bindices=sequence(new int(int x){return x;},m)+n-m; + + rational[] cB=phase1 ? new rational[m] : c[n-m:n]; rational[][] D=phase1 ? new rational[m+1][n+1] : E; - rational[] Dm=D[m]; - rational[] cb=phase1 ? new rational[m] : c[n-m:n]; if(phase1) { + bool output=true; + // Drive artificial variables out of basis. + for(int i=0; i < m; ++i) { + int k=Bindices[i]; + if(k >= n) { + rational[] Ei=E[i]; + int j; + for(j=0; j < n; ++j) + if(Ei[j] != 0) break; + if(j == n) continue; + output=false; + simplexTableau(E,Bindices,i,j); + Bindices[i]=j; + rowreduce(E,n,i,j); + } + } + if(output) simplexTableau(E,Bindices); int ip=0; // reduced i for(int i=0; i < m; ++i) { int k=Bindices[i]; if(k >= n) continue; Bindices[ip]=k; - cb[ip]=c[k]; + cB[ip]=c[k]; rational[] Dip=D[ip]; rational[] Ei=E[i]; for(int j=0; j < n; ++j) @@ -171,27 +258,33 @@ struct simplex { Dip[j]=Em[j]; Dip[n]=Em[N]; - m=ip; - - for(int j=0; j < n; ++j) { - rational sum=0; - for(int k=0; k < m; ++k) - sum += cb[k]*D[k][j]; - Dm[j]=c[j]-sum; + if(m > ip) { + Bindices.delete(ip,m-1); + D.delete(ip,m-1); + m=ip; } + if(!output) simplexTableau(D,Bindices); + } - // Done with Phase 1 + rational[] Dm=D[m]; + for(int j=0; j < n; ++j) { + rational sum=0; + for(int k=0; k < m; ++k) + sum += cB[k]*D[k][j]; + Dm[j]=c[j]-sum; } - + rational sum=0; for(int k=0; k < m; ++k) - sum += cb[k]*D[k][n]; + sum += cB[k]*D[k][n]; Dm[n]=-sum; - if(iterate(D,n,Bindices) == UNBOUNDED) { - case=UNBOUNDED; - return; - } + simplexPhase2(); + + case=(dual ? iterateDual : iterate)(D,n,Bindices); + simplexTableau(D,Bindices); + if(case != OPTIMAL) + return; for(int j=0; j < n; ++j) x[j]=0; @@ -200,15 +293,16 @@ struct simplex { x[Bindices[k]]=D[k][n]; cost=-Dm[n]; - case=OPTIMAL; } // Try to find a solution x to sgn(Ax-b)=sgn(s) that minimizes the cost - // c^T x, where A is an m x n matrix, x is a vector of length n, b is a - // vector of length m, and c is a vector of length n. + // c^T x, where A is an m x n matrix, x is a vector of n non-negative + // numbers, b is a vector of length m, and c is a vector of length n. void operator init(rational[] c, rational[][] A, int[] s, rational[] b) { int m=A.length; + if(m == 0) {case=INFEASIBLE; return;} int n=A[0].length; + if(n == 0) {case=INFEASIBLE; return;} int count=0; for(int i=0; i < m; ++i) @@ -226,6 +320,9 @@ struct simplex { int k=0; + bool phase1=false; + bool dual=count == m && all(c >= 0); + for(int i=0; i < m; ++i) { rational[] ai=a[i]; for(int j=0; j < k; ++j) @@ -234,47 +331,32 @@ struct simplex { ai[n+k]=-s[i]; for(int j=k+1; j < count; ++j) ai[n+j]=0; - if(s[i] != 0) ++k; + int si=s[i]; + if(si == 0) phase1=true; + else { + ++k; + rational bi=b[i]; + if(bi == 0) { + if(si == 1) { + s[i]=-1; + for(int j=0; j < n+count; ++j) + ai[j]=-ai[j]; + } + } else if(si*bi > 0) { + if(dual && si == 1) { + b[i]=-bi; + s[i]=-1; + for(int j=0; j < n+count; ++j) + ai[j]=-ai[j]; + } else + phase1=true; + } + } } - bool phase1=!all(s == -1); - operator init(concat(c,array(count,rational(0))),a,b,phase1); + operator init(concat(c,array(count,rational(0))),a,b,phase1,dual); - if(case == OPTIMAL) + if(case == OPTIMAL && count > 0) x.delete(n,n+count-1); } } - -/* -simplex S=simplex(new rational[] {4,1,1}, - new rational[][] {{2,1,2},{3,3,1}}, - new rational[] {4,3}); - -simplex S=simplex(new rational[] {2,6,1,1}, - new rational[][] {{1,2,0,1},{1,2,1,1},{1,3,-1,2},{1,1,1,0}}, - new rational[] {6,7,7,5}); -simplex S=simplex(new rational[] {-10,-12,-12,0,0,0}, - new rational[][] {{1,2,2,1,0,0}, - {2,1,2,0,1,0}, - {2,2,1,0,0,1}}, - new rational[] {20,20,20}); - -simplex S=simplex(new rational[] {-10,-12,-12}, - new rational[][] {{1,2,2}, - {2,1,2}, - {2,2,1}}, - new int[] {0,0,-1}, - new rational[] {20,20,20}); - -simplex S=simplex(new rational[] {1,1,1,0}, - new rational[][] {{1,2,3,0}, - {-1,2,6,0}, - {0,4,9,0}, - {0,0,3,1}}, - new rational[] {3,2,5,1}); - -write(); -write("case:",S.case); -write("x:",S.x); -write("Cost=",S.cost); -*/ diff --git a/Build/source/utils/asymptote/base/shaders/fragment.glsl b/Build/source/utils/asymptote/base/shaders/fragment.glsl new file mode 100644 index 00000000000..b6a85597451 --- /dev/null +++ b/Build/source/utils/asymptote/base/shaders/fragment.glsl @@ -0,0 +1,85 @@ +struct Material
+{
+ vec4 diffuse,ambient,emissive,specular;
+ float shininess;
+};
+
+struct Light
+{
+ vec4 direction;
+ vec4 diffuse,ambient,specular;
+};
+
+uniform int nlights;
+
+uniform Light lights[Nlights];
+
+uniform MaterialBuffer {
+ Material Materials[Nmaterials];
+};
+
+in vec3 Normal;
+
+#ifdef EXPLICIT_COLOR
+in vec4 Color;
+#endif
+flat in int materialIndex;
+
+out vec4 outColor;
+
+void main()
+{
+ vec4 Diffuse;
+ vec4 Ambient;
+ vec4 Emissive;
+ vec4 Specular;
+ float Shininess;
+
+#ifdef EXPLICIT_COLOR
+ if(materialIndex < 0) {
+ int index=-materialIndex-1;
+ Material m=Materials[index];
+ Diffuse=Color;
+ Ambient=Color;
+ Emissive=vec4(0.0,0.0,0.0,1.0);
+ Specular=m.specular;
+ Shininess=m.shininess;
+ } else {
+ Material m=Materials[materialIndex];
+ Diffuse=m.diffuse;
+ Ambient=m.ambient;
+ Emissive=m.emissive;
+ Specular=m.specular;
+ Shininess=m.shininess;
+ }
+#else
+ Material m=Materials[materialIndex];
+ Diffuse=m.diffuse;
+ Ambient=m.ambient;
+ Emissive=m.emissive;
+ Specular=m.specular;
+ Shininess=m.shininess;
+#endif
+ // Phong-Blinn model
+ if(nlights > 0) {
+ vec3 diffuse=vec3(0,0,0);
+ vec3 specular=vec3(0,0,0);
+ vec3 ambient=vec3(0,0,0);
+ vec3 Z=vec3(0,0,1);
+
+ for(int i=0; i < nlights; ++i) {
+ vec3 L=normalize(lights[i].direction.xyz);
+ diffuse += lights[i].diffuse.rgb*abs(dot(Normal,L));
+ ambient += lights[i].ambient.rgb;
+ specular += pow(abs(dot(Normal,normalize(L+Z))),Shininess)*
+ lights[i].specular.rgb;
+ }
+
+ vec3 color=diffuse*Diffuse.rgb+
+ ambient*Ambient.rgb+
+ specular*Specular.rgb+
+ Emissive.rgb;
+ outColor=vec4(color,Diffuse[3]);
+ } else
+ outColor=Diffuse;
+}
diff --git a/Build/source/utils/asymptote/base/shaders/vertex.glsl b/Build/source/utils/asymptote/base/shaders/vertex.glsl new file mode 100644 index 00000000000..9e3feda6988 --- /dev/null +++ b/Build/source/utils/asymptote/base/shaders/vertex.glsl @@ -0,0 +1,34 @@ +in vec3 position;
+in vec3 normal;
+
+#ifdef EXPLICIT_COLOR
+in uint color;
+#endif
+
+in int material;
+
+uniform mat4 projViewMat;
+uniform mat4 viewMat;
+uniform mat4 normMat;
+
+out vec3 ViewPosition;
+out vec3 Normal;
+
+#ifdef EXPLICIT_COLOR
+out vec4 Color;
+#endif
+
+flat out int materialIndex;
+
+void main()
+{
+ gl_Position=projViewMat*vec4(position,1.0);
+ ViewPosition=(viewMat*vec4(position,1.0)).xyz;
+ Normal=normalize((normMat*vec4(normal,0)).xyz);
+
+#ifdef EXPLICIT_COLOR
+ Color=unpackUnorm4x8(color);
+#endif
+
+ materialIndex=material;
+}
diff --git a/Build/source/utils/asymptote/base/simplex.asy b/Build/source/utils/asymptote/base/simplex.asy index 16a495a082e..040621dd904 100644 --- a/Build/source/utils/asymptote/base/simplex.asy +++ b/Build/source/utils/asymptote/base/simplex.asy @@ -1,4 +1,4 @@ -// General simplex solver written by John C. Bowman and Pouria Ramazi, 2018. +// Real simplex solver written by John C. Bowman and Pouria Ramazi, 2018. struct simplex { static int OPTIMAL=0; @@ -49,7 +49,7 @@ struct simplex { if(Em[J] < 0) break; if(J == N) - return 0; + break; int I=-1; real M; @@ -65,28 +65,69 @@ struct simplex { real e=E[i][J]; if(e > epsilonA) { real v=E[i][N]/e; - if(v <= M) {M=v; I=i;} + if(v < M) {M=v; I=i;} // Bland's rule: choose smallest argmin } } if(I == -1) return UNBOUNDED; // Can only happen in Phase 2. + // Generate new tableau Bindices[I]=J; + rowreduce(E,N,I,J); + } + return OPTIMAL; + } + + int iterateDual(real[][] E, int N, int[] Bindices) { + while(true) { + // Find first negative entry in right (basic variable) column + real[] Em=E[m]; + int I; + for(I=0; I < m; ++I) { + if(E[I][N] < 0) break; + } + + if(I == m) + break; + + int J=-1; + real M; + for(int j=0; j < N; ++j) { + real e=E[I][j]; + if(e < epsilonA) { + M=-E[m][j]/e; + J=j; + break; + } + } + for(int j=J+1; j < N; ++j) { + real e=E[I][j]; + if(e < epsilonA) { + real v=-E[m][j]/e; + if(v < M) {M=v; J=j;} // Bland's rule: choose smallest argmin + } + } + if(J == -1) + return INFEASIBLE; // Can only happen in Phase 2. // Generate new tableau + Bindices[I]=J; rowreduce(E,N,I,J); } - return 0; + return OPTIMAL; } // Try to find a solution x to Ax=b that minimizes the cost c^T x, // where A is an m x n matrix, x is a vector of n non-negative numbers, // b is a vector of length m, and c is a vector of length n. - void operator init(real[] c, real[][] A, real[] b, bool phase1=true) { + // Can set phase1=false if the last m columns of A form the identity matrix. + void operator init(real[] c, real[][] A, real[] b, bool phase1=true, + bool dual=false) { + if(dual) phase1=false; static real epsilon=sqrt(realEpsilon); epsilonA=epsilon*norm(A); - // Phase 1 + // Phase 1 m=A.length; if(m == 0) {case=INFEASIBLE; return;} n=A[0].length; @@ -102,7 +143,7 @@ struct simplex { for(int i=0; i < m; ++i) { real[] Ai=A[i]; real[] Ei=E[i]; - if(b[i] >= 0) { + if(b[i] >= 0 || dual) { for(int j=0; j < n; ++j) { real Aij=Ai[j]; Ei[j]=Aij; @@ -130,7 +171,7 @@ struct simplex { real sum=0; for(int i=0; i < m; ++i) { - real B=abs(b[i]); + real B=dual ? b[i] : abs(b[i]); E[i][N]=B; sum -= B; } @@ -140,27 +181,40 @@ struct simplex { for(int j=0; j < m; ++j) Em[n+j]=0.0; - int[] Bindices=sequence(new int(int x){return x;},m)+n; + int[] Bindices; if(phase1) { + Bindices=sequence(new int(int x){return x;},m)+n; iterate(E,N,Bindices); if(abs(Em[J]) > epsilonA) { case=INFEASIBLE; return; } - } + } else Bindices=sequence(new int(int x){return x;},m)+n-m; + real[] cB=phase1 ? new real[m] : c[n-m:n]; real[][] D=phase1 ? new real[m+1][n+1] : E; - real[] Dm=D[m]; - real[] cb=phase1 ? new real[m] : c[n-m:n]; if(phase1) { + // Drive artificial variables out of basis. + for(int i=0; i < m; ++i) { + int k=Bindices[i]; + if(k >= n) { + real[] Ei=E[i]; + int j; + for(j=0; j < n; ++j) + if(Ei[j] != 0) break; + if(j == n) continue; + Bindices[i]=j; + rowreduce(E,n,i,j); + } + } int ip=0; // reduced i for(int i=0; i < m; ++i) { int k=Bindices[i]; if(k >= n) continue; Bindices[ip]=k; - cb[ip]=c[k]; + cB[ip]=c[k]; real[] Dip=D[ip]; real[] Ei=E[i]; for(int j=0; j < n; ++j) @@ -175,27 +229,29 @@ struct simplex { Dip[j]=Em[j]; Dip[n]=Em[N]; - m=ip; - - for(int j=0; j < n; ++j) { - real sum=0; - for(int k=0; k < m; ++k) - sum += cb[k]*D[k][j]; - Dm[j]=c[j]-sum; + if(m > ip) { + Bindices.delete(ip,m-1); + D.delete(ip,m-1); + m=ip; } + } - // Done with Phase 1 + real[] Dm=D[m]; + for(int j=0; j < n; ++j) { + real sum=0; + for(int k=0; k < m; ++k) + sum += cB[k]*D[k][j]; + Dm[j]=c[j]-sum; } - + real sum=0; for(int k=0; k < m; ++k) - sum += cb[k]*D[k][n]; + sum += cB[k]*D[k][n]; Dm[n]=-sum; - if(iterate(D,n,Bindices) == UNBOUNDED) { - case=UNBOUNDED; - return; - } + case=(dual ? iterateDual : iterate)(D,n,Bindices); + if(case != OPTIMAL) + return; for(int j=0; j < n; ++j) x[j]=0; @@ -204,7 +260,6 @@ struct simplex { x[Bindices[k]]=D[k][n]; cost=-Dm[n]; - case=OPTIMAL; } // Try to find a solution x to sgn(Ax-b)=sgn(s) that minimizes the cost @@ -232,6 +287,9 @@ struct simplex { int k=0; + bool phase1=false; + bool dual=count == m && all(c >= 0); + for(int i=0; i < m; ++i) { real[] ai=a[i]; for(int j=0; j < k; ++j) @@ -240,14 +298,32 @@ struct simplex { ai[n+k]=-s[i]; for(int j=k+1; j < count; ++j) ai[n+j]=0; - if(s[i] != 0) ++k; + int si=s[i]; + if(si == 0) phase1=true; + else { + ++k; + real bi=b[i]; + if(bi == 0) { + if(si == 1) { + s[i]=-1; + for(int j=0; j < n+count; ++j) + ai[j]=-ai[j]; + } + } else if(si*bi > 0) { + if(dual && si == 1) { + b[i]=-bi; + s[i]=-1; + for(int j=0; j < n+count; ++j) + ai[j]=-ai[j]; + } else + phase1=true; + } + } } - // bool phase1=!all(s == -1); // TODO: Check - bool phase1=true; - operator init(concat(c,array(count,0.0)),a,b,phase1); + operator init(concat(c,array(count,0.0)),a,b,phase1,dual); - if(case == OPTIMAL) + if(case == OPTIMAL && count > 0) x.delete(n,n+count-1); } } diff --git a/Build/source/utils/asymptote/base/three.asy b/Build/source/utils/asymptote/base/three.asy index 742ceda1c04..8b9d2137fc9 100644 --- a/Build/source/utils/asymptote/base/three.asy +++ b/Build/source/utils/asymptote/base/three.asy @@ -2021,7 +2021,7 @@ path3 arc(triple c, triple v1, triple v2, triple normal=O, bool direction=CCW) } string invalidnormal="invalid normal vector"; - real fuzz=sqrtEpsilon*max(abs(v1),abs(v2)); + real fuzz=sqrtEpsilon; if(abs(v1.z) > fuzz || abs(v2.z) > fuzz) abort(invalidnormal); @@ -2558,23 +2558,6 @@ private string Format(transform3 t, string sep=" ") Format(t[0][3])+sep+Format(t[1][3])+sep+Format(t[2][3]); } -string lightscript(light light) { - string script="for(var i=scene.lights.count-1; i >= 0; i--) - scene.lights.removeByIndex(i);"+'\n\n'; - for(int i=0; i < light.position.length; ++i) { - string Li="L"+string(i); - real[] diffuse=light.diffuse[i]; - script += Li+"=scene.createLight();"+'\n'+ - Li+".direction.set("+format(-light.position[i],",")+");"+'\n'+ - Li+".color.set("+format((diffuse[0],diffuse[1],diffuse[2]),",")+");"+'\n'; - } - // Work around initialization bug in Adobe Reader 8.0: - return script +" -scene.lightScheme=scene.LIGHT_MODE_HEADLAMP; -scene.lightScheme=scene.LIGHT_MODE_FILE; -"; -} - void writeJavaScript(string name, string preamble, string script) { file out=output(name); @@ -2625,10 +2608,6 @@ string embed3D(string prefix, string label=prefix, string text=label, shipout3(prefix,f); string name=prefix+".js"; - // Adobe Reader doesn't appear to support user-specified viewport lights. - bool lightscript=light.on() && !light.viewport; - if(lightscript) - writeJavaScript(name,lightscript(light),script); if(!settings.inlinetex && !prconly()) file3.push(prefix+".prc"); @@ -2637,7 +2616,7 @@ string embed3D(string prefix, string label=prefix, string text=label, transform3 inv=inverse(flipxz*P.T.modelview); string options3="3Dlights="+ - (light.on() ? (light.viewport ? "Headlamp" : "File") : "None"); + (light.on() ? "Headlamp" : "None"); if(defaultembed3Doptions != "") options3 += ","+defaultembed3Doptions; if((settings.render < 0 || !settings.embed) && settings.auto3D) @@ -2653,8 +2632,6 @@ string embed3D(string prefix, string label=prefix, string text=label, if(options != "") options3 += ","+options; if(settings.inlinetex) prefix=jobname(prefix); - if(lightscript) - options3 += ",add3Djscript="+prefix+".js"; options3 += ",add3Djscript=asylabels.js"; return text == "" ? Embed(prefix+".prc","",options3,width,height) : @@ -2909,7 +2886,7 @@ object embed(string prefix=outprefix(), string label=prefix, P.zoom,m,M,P.viewportshift, tinv*inv*shift(0,0,zcenter),Light.background(),Light.position, Light.diffuse,Light.ambient,Light.specular, - Light.viewport,view && !preview); + view && !preview); if(!preview) return F; } diff --git a/Build/source/utils/asymptote/base/three_light.asy b/Build/source/utils/asymptote/base/three_light.asy index 6cbd9b6f0be..3bd747dd4ec 100644 --- a/Build/source/utils/asymptote/base/three_light.asy +++ b/Build/source/utils/asymptote/base/three_light.asy @@ -86,36 +86,41 @@ pen color(triple normal, material m, light light, transform3 T=light.T) { real[] Ambient=rgba(m.ambient()); real[] Specular=rgba(m.specular()); real[] p=rgba(m.emissive()); + real[] diffuse={0,0,0,0}; + real[] ambient={0,0,0,0}; + real[] specular={0,0,0,0}; for(int i=0; i < position.length; ++i) { - triple L=light.viewport ? position[i] : T*position[i]; - real Ldotn=max(dot(normal,L),0); - p += light.ambient[i]*Ambient+Ldotn*light.diffuse[i]*Diffuse; - // Apply specularfactor to partially compensate non-pixel-based rendering. - if(Ldotn > 0) // Phong-Blinn model of specular reflection - p += dot(normal,unit(L+Z))^s*light.specularfactor* - light.specular[i]*Specular; + triple L=position[i]; + real dotproduct=abs(dot(normal,L)); + diffuse += dotproduct*light.diffuse[i]; + ambient += light.ambient[i]; + dotproduct=abs(dot(normal,unit(L+Z))); + // Phong-Blinn model of specular reflection + specular += dotproduct^s*light.specular[i]; } + p += diffuse*Diffuse; + p += ambient*Ambient; + // Apply specularfactor to partially compensate non-pixel-based rendering. + p += specular*Specular*light.specularfactor; return rgb(p[0],p[1],p[2])+opacity(opacity(m.diffuse())); } light operator * (transform3 t, light light) { light light=light(light); - if(!light.viewport) light.position=shiftless(t)*light.position; return light; } light operator cast(triple v) {return light(v);} -light Viewport=light(ambient=gray(0.1),specularfactor=3,viewport=true, - (0.25,-0.25,1)); +light Viewport=light(ambient=gray(0.1),specularfactor=3,(0.25,-0.25,1)); light White=light(new pen[] {rgb(0.38,0.38,0.45),rgb(0.6,0.6,0.67), rgb(0.5,0.5,0.57)},specularfactor=3, new triple[] {(-2,-1.5,-0.5),(2,1.1,-2.5),(-0.5,0,2)}); light Headlamp=light(gray(0.8),ambient=gray(0.1),specular=gray(0.7), - specularfactor=3,viewport=true,dir(42,48)); + specularfactor=3,dir(42,48)); currentlight=Headlamp; diff --git a/Build/source/utils/asymptote/base/three_surface.asy b/Build/source/utils/asymptote/base/three_surface.asy index 6e10246c774..b3eebee45cc 100644 --- a/Build/source/utils/asymptote/base/three_surface.asy +++ b/Build/source/utils/asymptote/base/three_surface.asy @@ -1409,11 +1409,32 @@ void draw3D(frame f, int type=0, patch s, triple center=O, material m, PRCshininess,s.colors,interaction.type,prc); } +int computeNormals(triple[] v, int[][] vi, triple[] n, int[][] ni) +{ + triple lastnormal=O; + for(int i=0; i < vi.length; ++i) { + int[] vii=vi[i]; + int[] nii=ni[i]; + triple normal=normal(new triple[] {v[vii[0]],v[vii[1]],v[vii[2]]}); + if(normal != lastnormal || n.length == 0) { + n.push(normal); + lastnormal=normal; + } + nii[0]=nii[1]=nii[2]=n.length-1; + } + return ni.length; +} + // Draw triangles on a frame. void draw(frame f, triple[] v, int[][] vi, triple[] n={}, int[][] ni={}, material m=currentpen, pen[] p={}, int[][] pi={}, light light=currentlight) { + bool normals=ni.length > 0; + if(!normals) { + ni=new int[vi.length][3]; + normals=computeNormals(v,vi,n,ni) > 0; + } if(p.length > 0) m=mean(p); m=material(m,light); @@ -1428,23 +1449,12 @@ void draw(picture pic=currentpicture, triple[] v, int[][] vi, triple[] n={}, int[][] ni={}, material m=currentpen, pen[] p={}, int[][] pi={}, light light=currentlight) { - bool colors=pi.length > 0; bool normals=ni.length > 0; - if(!colors && !normals) { - n=new triple[]; + if(!normals) { ni=new int[vi.length][3]; - triple lastnormal=O; - for(int i=0; i < vi.length; ++i) { - int[] vii=vi[i]; - int[] nii=ni[i]; - triple normal=normal(new triple[] {v[vii[0]],v[vii[1]],v[vii[2]]}); - if(normal != lastnormal || n.length == 0) { - n.push(normal); - lastnormal=normal; - } - nii[0]=nii[1]=nii[2]=n.length-1; - } + normals=computeNormals(v,vi,n,ni) > 0; } + bool colors=pi.length > 0; pic.add(new void(frame f, transform3 t, picture pic, projection P) { triple[] v=t*v; |