diff options
author | Karl Berry <karl@freefriends.org> | 2016-03-14 22:52:10 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2016-03-14 22:52:10 +0000 |
commit | 355bf72b81c6931539dd365430e78ebc10bcd0a2 (patch) | |
tree | ff9e855fb9eeb6cbc286d58e6cfe07bdb9cdf4ba /Build/source/utils/asymptote/base | |
parent | 1cdee1080e0ad8f88624629e00039d9fe89c1b54 (diff) |
asy 2.37 sources
git-svn-id: svn://tug.org/texlive/trunk@40025 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/utils/asymptote/base')
-rw-r--r-- | Build/source/utils/asymptote/base/bezulate.asy | 2 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/embed.asy | 11 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/math.asy | 12 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/plain_strings.asy | 21 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/smoothcontour3.asy | 684 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/three_surface.asy | 144 |
6 files changed, 515 insertions, 359 deletions
diff --git a/Build/source/utils/asymptote/base/bezulate.asy b/Build/source/utils/asymptote/base/bezulate.asy index d0b8958a4da..38d35718d4a 100644 --- a/Build/source/utils/asymptote/base/bezulate.asy +++ b/Build/source/utils/asymptote/base/bezulate.asy @@ -2,7 +2,7 @@ private real fuzz=sqrtEpsilon; real duplicateFuzz=1e-3; // Work around font errors. -real maxrefinements=7; +real maxrefinements=10; private real[][] intersections(pair a, pair b, path p) { diff --git a/Build/source/utils/asymptote/base/embed.asy b/Build/source/utils/asymptote/base/embed.asy index 588ad79c283..88495f12128 100644 --- a/Build/source/utils/asymptote/base/embed.asy +++ b/Build/source/utils/asymptote/base/embed.asy @@ -2,11 +2,6 @@ if(latex()) { usepackage("hyperref"); texpreamble("\hypersetup{"+settings.hyperrefOptions+"}"); usepackage("media9","bigfiles"); - texpreamble("\newif\ifplaybutton"); - texpreamble("\count255=\the\catcode`\@\makeatletter% -\@ifpackagelater{media9}{2013/11/15}{}{\playbuttontrue}% -\catcode`\@=\the\count255 -%"); } // For documentation of the options see @@ -19,11 +14,7 @@ string embedplayer(string name, string text="", string options="", if(width != 0) options += ",width="+(string) (width/pt)+"pt"; if(height != 0) options += ",height="+(string) (height/pt)+"pt"; return "% -\ifplaybutton% -\includemedia["+options+"]{"+text+"}{"+name+"}% -\else% -\includemedia[noplaybutton,"+options+"]{"+text+"}{"+name+"}% -\fi"; +\includemedia[noplaybutton,"+options+"]{"+text+"}{"+name+"}"; } // Embed media in pdf file diff --git a/Build/source/utils/asymptote/base/math.asy b/Build/source/utils/asymptote/base/math.asy index 8cfd61f4344..9d10e2b7794 100644 --- a/Build/source/utils/asymptote/base/math.asy +++ b/Build/source/utils/asymptote/base/math.asy @@ -437,3 +437,15 @@ real[] leastsquares(real[][] A, real[] b, bool warn=true) "a matrix with linearly dependent columns."); return solution; } + +// Namespace +struct rootfinder_settings { + static real roottolerance = 1e-4; +} + +real findroot(real f(real), real a, real b, + real tolerance=rootfinder_settings.roottolerance, + real fa=f(a), real fb=f(b)) +{ + return _findroot(f,a,b,tolerance,fa,fb); +} diff --git a/Build/source/utils/asymptote/base/plain_strings.asy b/Build/source/utils/asymptote/base/plain_strings.asy index a5c72b2035d..a88e3863416 100644 --- a/Build/source/utils/asymptote/base/plain_strings.asy +++ b/Build/source/utils/asymptote/base/plain_strings.asy @@ -228,4 +228,25 @@ string phantom(string s) return settings.tex != "none" ? "\phantom{"+s+"}" : ""; } +string[] spinner=new string[] {'|','/','-','\\'}; +spinner.cyclic=true; + +void progress(bool3 init=default) +{ + static int count=-1; + static int lastseconds=-1; + if(init == true) { + lastseconds=0; + write(stdout,' ',flush); + } else + if(init == default) { + int seconds=seconds(); + if(seconds > lastseconds) { + lastseconds=seconds; + write(stdout,'\b'+spinner[++count],flush); + } + } else + write(stdout,'\b',flush); +} + restricted int ocgindex=0; diff --git a/Build/source/utils/asymptote/base/smoothcontour3.asy b/Build/source/utils/asymptote/base/smoothcontour3.asy index 9b34fc7c098..4fe26654283 100644 --- a/Build/source/utils/asymptote/base/smoothcontour3.asy +++ b/Build/source/utils/asymptote/base/smoothcontour3.asy @@ -64,7 +64,7 @@ function[] bernstein = new function[] {B03, B13, B23, B33}; // the specified boundary path. However, the entries in the array // remain intact. patch patchwithnormals(path3 external, triple[] u0normals, triple[] u1normals, - triple[] v0normals, triple[] v1normals) + triple[] v0normals, triple[] v1normals) { assert(cyclic(external)); assert(length(external) == 4); @@ -213,8 +213,8 @@ patch patchwithnormals(path3 external, triple[] u0normals, triple[] u1normals, real[] solution = leastsquares(matrix, rightvector, warn=false); if (solution.length == 0) { // if the matrix was singular - write("Warning: unable to solve matrix for specifying normals " - + "on bezier patch. Using standard method."); + write("Warning: unable to solve matrix for specifying edge normals " + + "on bezier patch. Using coons patch."); return patch(external); } @@ -222,22 +222,95 @@ patch patchwithnormals(path3 external, triple[] u0normals, triple[] u1normals, for (int j = 1; j <= 2; ++j) { int position = 3 * (2 * (i-1) + (j-1)); controlpoints[i][j] = (solution[position], - solution[position+1], - solution[position+2]); + solution[position+1], + solution[position+2]); } } return patch(controlpoints); } -// A wrapper for the previous function when the normal direction +// This function attempts to produce a Bezier triangle +// with the specified boundary path and normal directions at the +// edge midpoints. The bezier triangle should be normal to +// n1 at point(external, 0.5), +// normal to n2 at point(external, 1.5), and +// normal to n3 at point(external, 2.5). +// The actual normal (as computed by the patch.normal() function) +// may be parallel to the specified normal, antiparallel, or +// even zero. +// +// A small amount of deviation is allowed in order to stabilize +// the algorithm (by keeping the mixed partials at the corners from +// growing too large). +patch trianglewithnormals(path3 external, triple n1, + triple n2, triple n3) { + assert(cyclic(external)); + assert(length(external) == 3); + // Use the formal symbols a3, a2b, abc, etc. to denote the control points, + // following the Wikipedia article on Bezier triangles. + triple a3 = point(external, 0), a2b = postcontrol(external, 0), + ab2 = precontrol(external, 1), b3 = point(external, 1), + b2c = postcontrol(external, 1), bc2 = precontrol(external, 2), + c3 = point(external, 2), ac2 = postcontrol(external, 2), + a2c = precontrol(external, 0); + + // Use orthogonal projection to ensure that the normal vectors are + // actually normal to the boundary path. + triple tangent = dir(external, 0.5); + n1 -= dot(n1,tangent)*tangent; + n1 = unit(n1); + + tangent = dir(external, 1.5); + n2 -= dot(n2,tangent)*tangent; + n2 = unit(n2); + + tangent = dir(external, 2.5); + n3 -= dot(n3,tangent)*tangent; + n3 = unit(n3); + + real wild = 2 * wildnessweight; + real[][] matrix = { {n1.x, n1.y, n1.z}, + {n2.x, n2.y, n2.z}, + {n3.x, n3.y, n3.z}, + { wild, 0, 0}, + { 0, wild, 0}, + { 0, 0, wild} }; + real[] rightvector = + { dot(n1, (a3 + 3a2b + 3ab2 + b3 - 2a2c - 2b2c)) / 4, + dot(n2, (b3 + 3b2c + 3bc2 + c3 - 2ab2 - 2ac2)) / 4, + dot(n3, (c3 + 3ac2 + 3a2c + a3 - 2bc2 - 2a2b)) / 4 }; + + // The inner control point that minimizes the sum of squares of + // the mixed partials on the corners. + triple tameinnercontrol = + ((a2b + a2c - a3) + (ab2 + b2c - b3) + (ac2 + bc2 - c3)) / 3; + rightvector.append(wild * new real[] + {tameinnercontrol.x, tameinnercontrol.y, tameinnercontrol.z}); + real[] solution = leastsquares(matrix, rightvector, warn=false); + if (solution.length == 0) { // if the matrix was singular + write("Warning: unable to solve matrix for specifying edge normals " + + "on bezier triangle. Using coons triangle."); + return patch(external); + } + triple innercontrol = (solution[0], solution[1], solution[2]); + return patch(external, innercontrol); +} + +// A wrapper for the previous functions when the normal direction // is given as a function of direction. The wrapper can also // accommodate cyclic boundary paths of between one and four // segments, although the results are best by far when there -// are four segments. +// are three or four segments. patch patchwithnormals(path3 external, triple normalat(triple)) { assert(cyclic(external)); assert(1 <= length(external) && length(external) <= 4); + if (length(external) == 3) { + triple n1 = normalat(point(external, 0.5)); + triple n2 = normalat(point(external, 1.5)); + triple n3 = normalat(point(external, 2.5)); + return trianglewithnormals(external, n1, n2, n3); + } while (length(external) < 4) external = external -- cycle; triple[] u0normals = new triple[3]; triple[] u1normals = new triple[3]; @@ -253,82 +326,6 @@ patch patchwithnormals(path3 external, triple normalat(triple)) { } /***********************************************/ -/*********** ROOT-FINDER UTILITY ***************/ -/***********************************************/ - -// Namespace -struct rootfinder_settings { - static real roottolerance = 1e-4; -} - -// Find a root for the specified continuous (but not -// necessarily differentiable) function. Whatever -// value t is returned, it is guaranteed that either -// t is within tolerance of a sign change, or -// abs(f(t)) <= 0.1 tolerance. -// An error is thrown if fa and fb are both positive -// or both negative. -// -// In the current implementation, binary search is interleaved -// with a modified version of linear interpolation. -real findroot(real f(real), real a, real b, - real tolerance = rootfinder_settings.roottolerance, - real fa = f(a), real fb = f(b)) -{ - if (fa == 0) return a; - if (fb == 0) return b; - real g(real); - if (fa < 0) { - assert(fb > 0); - g = f; - } else { - assert(fb < 0); - fa = -fa; - fb = -fb; - g = new real(real t) { return -f(t); }; - } - - real t = a; - real ft = fa; - - while (b - a > tolerance && abs(ft) > 0.1*tolerance) { - t = a + (b - a) / 2; - ft = g(t); - if (ft == 0) return t; - else if (ft > 0) { - b = t; - fb = ft; - } else if (ft < 0) { - a = t; - fa = ft; - } - - // linear interpolation - t = a - (b - a) / (fb - fa) * fa; - - // If the interpolated value is close to one edge of - // the interval, move it farther away from the edge in - // an effort to catch the root in the middle. - if (t - a < (b-a)/8) t = a + 2*(t-a); - else if (b - t < (b-a)/8) t = b - 2*(b-t); - - assert(t >= a && t <= b); - - ft = g(t); - if (ft == 0) return t; - else if (ft > 0) { - b = t; - fb = ft; - } else if (ft < 0) { - a = t; - fa = ft; - } - - } - return a - (b - a) / (fb - fa) * fa; -} - -/***********************************************/ /********* DUAL CUBE GRAPH UTILITY *************/ /***********************************************/ @@ -404,8 +401,8 @@ struct int_to_intset { void add(int key, int value) { for (int i = 0; i < keys.length; ++i) { if (keys[i] == key) { - values[i].add(value); - return; + values[i].add(value); + return; } } keys.push(key); @@ -499,9 +496,9 @@ int[] makecircle(edge[] edges) { if (adjacentvertices.length != 2) return null; for (int v : adjacentvertices) { if (v != lastvertex) { - lastvertex = currentvertex; - currentvertex = v; - break; + lastvertex = currentvertex; + currentvertex = v; + break; } } } while (currentvertex != startvertex); @@ -535,9 +532,10 @@ string operator cast(positionedvector vv) { // The angle, in degrees, between two vectors. real angledegrees(triple a, triple b) { - real lengthprod = max(abs(a) * abs(b), abs(dot(a,b))); + real dotprod = dot(a,b); + real lengthprod = max(abs(a) * abs(b), abs(dotprod)); if (lengthprod == 0) return 0; - return aCos(dot(a,b) / lengthprod); + return aCos(dotprod / lengthprod); } // A path (single curved segment) between two points. At each point @@ -560,10 +558,10 @@ path3 pathbetween(positionedvector v1, positionedvector v2) { // the span of v1 and v2. If v1 and v2 are dependent, returns an empty array // (if warn==false) or throws an error (if warn==true). real[] projecttospan_findcoeffs(triple toproject, triple v1, triple v2, - bool warn=false) { + bool warn=false) { real[][] matrix = {{v1.x, v2.x}, - {v1.y, v2.y}, - {v1.z, v2.z}}; + {v1.y, v2.y}, + {v1.z, v2.z}}; real[] desiredanswer = {toproject.x, toproject.y, toproject.z}; return leastsquares(matrix, desiredanswer, warn=warn); } @@ -573,7 +571,7 @@ real[] projecttospan_findcoeffs(triple toproject, triple v1, triple v2, // a >= mincoeff and b >= mincoeff. If v1 and v2 are linearly dependent, // return a random (positive) linear combination. triple projecttospan(triple toproject, triple v1, triple v2, - real mincoeff = 0.05) { + real mincoeff = 0.05) { real[] coeffs = projecttospan_findcoeffs(toproject, v1, v2, warn=false); real a, b; if (coeffs.length == 0) { @@ -639,25 +637,25 @@ path3 bisector(path3 edgecycle, int[] savevertices) { int opposite = i + 3; triple vec = unit(point(edgecycle, opposite) - point(edgecycle, i)); real[] coeffsbegin = projecttospan_findcoeffs(vec, forwarddirections[i], - backwarddirections[i]); + backwarddirections[i]); if (coeffsbegin.length == 0) continue; coeffsbegin[0] = max(coeffsbegin[0], mincoeff); coeffsbegin[1] = max(coeffsbegin[1], mincoeff); real[] coeffsend = projecttospan_findcoeffs(-vec, forwarddirections[opposite], - backwarddirections[opposite]); + backwarddirections[opposite]); if (coeffsend.length == 0) continue; coeffsend[0] = max(coeffsend[0], mincoeff); coeffsend[1] = max(coeffsend[1], mincoeff); real goodness = angles[i] * angles[opposite] * coeffsbegin[0] * coeffsend[0] - * coeffsbegin[1] * coeffsend[1]; + * coeffsbegin[1] * coeffsend[1]; if (goodness > maxgoodness) { maxgoodness = goodness; directionout = coeffsbegin[0] * forwarddirections[i] + - coeffsbegin[1] * backwarddirections[i]; + coeffsbegin[1] * backwarddirections[i]; directionin = -(coeffsend[0] * forwarddirections[opposite] + - coeffsend[1] * backwarddirections[opposite]); + coeffsend[1] * backwarddirections[opposite]); chosenindex = i; } } @@ -676,7 +674,7 @@ path3 bisector(path3 edgecycle, int[] savevertices) { // A path between two specified points (with specified normals) that lies // within a specified face of a rectangular solid. path3 pathinface(positionedvector v1, positionedvector v2, - triple facenorm, triple edge1normout, triple edge2normout) + triple facenorm, triple edge1normout, triple edge2normout) { triple dir1 = cross(v1.direction, facenorm); real dotprod = dot(dir1, edge1normout); @@ -707,15 +705,16 @@ triple normalout(int face) { // A path between two specified points (with specified normals) that lies // within a specified face of a rectangular solid. path3 pathinface(positionedvector v1, positionedvector v2, - int face, int edge1face, int edge2face) { + int face, int edge1face, int edge2face) { return pathinface(v1, v2, normalout(face), normalout(edge1face), - normalout(edge2face)); + normalout(edge2face)); } /***********************************************/ /******** DRAWING IMPLICIT SURFACES ************/ /***********************************************/ +// DEPRECATED // Quadrilateralization: // Produce a surface (array of *nondegenerate* Bezier patches) with a // specified three-segment boundary. The surface should approximate the @@ -725,11 +724,12 @@ path3 pathinface(positionedvector v1, positionedvector v2, // specified rectangular region, returns a length-zero array. // // Dividing a triangle into smaller quadrilaterals this way is opposite -// the usual trend in mathematics. However, the pathwithnormals algorithm -// does a poor job of choosing a good surface when the boundary path does +// the usual trend in mathematics. However, *before the introduction of bezier +// triangles,* the pathwithnormals algorithm +// did a poor job of choosing a good surface when the boundary path did // not consist of four positive-length segments. patch[] triangletoquads(path3 external, real f(triple), triple grad(triple), - triple a, triple b) { + triple a, triple b) { static real epsilon = 1e-3; assert(length(external) == 3); assert(cyclic(external)); @@ -805,6 +805,65 @@ patch[] triangletoquads(path3 external, real f(triple), triple grad(triple), patchwithnormals(quad2, grad)}; } +// Attempts to fill the path external (which should by a cyclic path consisting of +// three segments) with bezier triangle(s). Returns an empty array if it fails. +// +// In more detail: A single bezier triangle is computed using trianglewithnormals. The normals of +// the resulting triangle at the midpoint of each edge are computed. If any of these normals +// is in the negative f direction, the external triangle is subdivided into four external triangles +// and the same procedure is applied to each. If one or more of them has an incorrectly oriented +// edge normal, the function gives up and returns an empty array. +// +// Thus, the returned array consists of 0, 1, or 4 bezier triangles; no other array lengths +// are possible. +// +// This function assumes that the path orientation is consistent with f (and its gradient) +// -- i.e., that +// at a corner, (tangent in) x (tangent out) is in the positive f direction. +patch[] maketriangle(path3 external, real f(triple), + triple grad(triple), bool allowsubdivide = true) { + assert(cyclic(external)); + assert(length(external) == 3); + triple m1 = point(external, 0.5); + triple n1 = unit(grad(m1)); + triple m2 = point(external, 1.5); + triple n2 = unit(grad(m2)); + triple m3 = point(external, 2.5); + triple n3 = unit(grad(m3)); + patch beziertriangle = trianglewithnormals(external, n1, n2, n3); + if (dot(n1, beziertriangle.normal(0.5, 0)) >= 0 && + dot(n2, beziertriangle.normal(0.5, 0.5)) >= 0 && + dot(n3, beziertriangle.normal(0, 0.5)) >= 0) + return new patch[] {beziertriangle}; + + if (!allowsubdivide) return new patch[0]; + + positionedvector m1 = positionedvector(m1, n1); + positionedvector m2 = positionedvector(m2, n2); + positionedvector m3 = positionedvector(m3, n3); + path3 p12 = pathbetween(m1, m2); + path3 p23 = pathbetween(m2, m3); + path3 p31 = pathbetween(m3, m1); + patch[] triangles = maketriangle(p12 & p23 & p31 & cycle, f, grad=grad, + allowsubdivide=false); + if (triangles.length < 1) return new patch[0]; + + triangles.append(maketriangle(subpath(external, -0.5, 0.5) & reverse(p31) & cycle, + f, grad=grad, allowsubdivide=false)); + if (triangles.length < 2) return new patch[0]; + + triangles.append(maketriangle(subpath(external, 0.5, 1.5) & reverse(p12) & cycle, + f, grad=grad, allowsubdivide=false)); + if (triangles.length < 3) return new patch[0]; + + triangles.append(maketriangle(subpath(external, 1.5, 2.5) & reverse(p23) & cycle, + f, grad=grad, allowsubdivide=false)); + if (triangles.length < 4) return new patch[0]; + + return triangles; +} + + // Returns true if the point is "nonsingular" (in the sense that the magnitude // of the gradient is not too small) AND very close to the zero locus of f // (assuming f is locally linear). @@ -829,14 +888,14 @@ bool checkptincube(triple pt, triple a, triple b) { if (zmin > zmax) { real t = zmax; zmax=zmin; zmin=t; } return ((xmin <= pt.x) && (pt.x <= xmax) && - (ymin <= pt.y) && (pt.y <= ymax) && - (zmin <= pt.z) && (pt.z <= zmax)); + (ymin <= pt.y) && (pt.y <= ymax) && + (zmin <= pt.z) && (pt.z <= zmax)); } // A convenience function for combining the previous two tests. bool checkpt(triple testpt, real f(triple), triple grad(triple), - triple a, triple b) { + triple a, triple b) { return checkptincube(testpt, a, b) && check_fpt_zero(testpt, f, grad); } @@ -848,8 +907,8 @@ bool checkpt(triple testpt, real f(triple), triple grad(triple), // array, which merely indicates that the boundary cycle is too small // to be worth filling in. patch[] quadpatches(path3 edgecycle, positionedvector[] corners, - real f(triple), triple grad(triple), - triple a, triple b) { + real f(triple), triple grad(triple), + triple a, triple b, bool usetriangles) { assert(corners.cyclic); // The tolerance for considering two points "essentially identical." @@ -862,8 +921,8 @@ patch[] quadpatches(path3 edgecycle, positionedvector[] corners, if (corners.length == 2) return new patch[0]; corners.delete(i); edgecycle = subpath(edgecycle, 0, i) - & subpath(edgecycle, i+1, length(edgecycle)) - & cycle; + & subpath(edgecycle, i+1, length(edgecycle)) + & cycle; --i; assert(length(edgecycle) == corners.length); } @@ -875,9 +934,9 @@ patch[] quadpatches(path3 edgecycle, positionedvector[] corners, if (corners.length == 2) { // If the area is too small, just ignore it; otherwise, subdivide. real area0 = abs(cross(-dir(edgecycle, 0, sign=-1, normalize=false), - dir(edgecycle, 0, sign=1, normalize=false))); + dir(edgecycle, 0, sign=1, normalize=false))); real area1 = abs(cross(-dir(edgecycle, 1, sign=-1, normalize=false), - dir(edgecycle, 1, sign=1, normalize=false))); + dir(edgecycle, 1, sign=1, normalize=false))); if (area0 < areatolerance && area1 < areatolerance) return new patch[0]; else return null; } @@ -885,13 +944,14 @@ patch[] quadpatches(path3 edgecycle, positionedvector[] corners, for (int i = 0; i < length(edgecycle); ++i) { if (angledegrees(dir(edgecycle,i,sign=1), - dir(edgecycle,i+1,sign=-1)) > 80) { + dir(edgecycle,i+1,sign=-1)) > 80) { return null; } } if (length(edgecycle) == 3) { - patch[] toreturn = triangletoquads(edgecycle, f, grad, a, b); + patch[] toreturn = usetriangles ? maketriangle(edgecycle, f, grad) + : triangletoquads(edgecycle, f, grad, a, b); if (toreturn.length == 0) return null; else return toreturn; } @@ -912,8 +972,9 @@ patch[] quadpatches(path3 edgecycle, positionedvector[] corners, & reverse(middleguide) & cycle; if (length(edgecycle) == 5) { path3 secondpatch = middleguide - & subpath(edgecycle, bisectorindices[1], 5+bisectorindices[0]) & cycle; - toreturn = triangletoquads(secondpatch, f, grad, a, b); + & subpath(edgecycle, bisectorindices[1], 5+bisectorindices[0]) & cycle; + toreturn = usetriangles ? maketriangle(secondpatch, f, grad) + : triangletoquads(secondpatch, f, grad, a, b); if (toreturn.length == 0) return null; toreturn.push(patchwithnormals(firstpatch, grad)); } else { @@ -922,7 +983,7 @@ patch[] quadpatches(path3 edgecycle, positionedvector[] corners, & subpath(edgecycle, bisectorindices[1], 6+bisectorindices[0]) & cycle; toreturn = new patch[] {patchwithnormals(firstpatch, grad), - patchwithnormals(secondpatch, grad)}; + patchwithnormals(secondpatch, grad)}; } return toreturn; } @@ -933,8 +994,8 @@ vectorfunction nGrad(real f(triple)) { static real epsilon = 1e-3; return new triple(triple v) { return ( (f(v + epsilon*X) - f(v - epsilon*X)) / (2 epsilon), - (f(v + epsilon*Y) - f(v - epsilon*Y)) / (2 epsilon), - (f(v + epsilon*Z) - f(v - epsilon*Z)) / (2 epsilon) ); + (f(v + epsilon*Y) - f(v - epsilon*Y)) / (2 epsilon), + (f(v + epsilon*Z) - f(v - epsilon*Z)) / (2 epsilon) ); }; } @@ -953,18 +1014,18 @@ triple operator cast(evaluatedpoint p) { return p.pt; } // Compute the values of a function at every vertex of an nx by ny by nz // array of rectangular solids. evaluatedpoint[][][] make3dgrid(triple a, triple b, int nx, int ny, int nz, - real f(triple), bool allowzero = false) + real f(triple), bool allowzero = false) { evaluatedpoint[][][] toreturn = new evaluatedpoint[nx+1][ny+1][nz+1]; for (int i = 0; i <= nx; ++i) { for (int j = 0; j <= ny; ++j) { for (int k = 0; k <= nz; ++k) { - triple pt = (interp(a.x, b.x, i/nx), - interp(a.y, b.y, j/ny), - interp(a.z, b.z, k/nz)); - real value = f(pt); - if (value == 0 && !allowzero) value = 1e-5; - toreturn[i][j][k] = evaluatedpoint(pt, value); + triple pt = (interp(a.x, b.x, i/nx), + interp(a.y, b.y, j/ny), + interp(a.z, b.z, k/nz)); + real value = f(pt); + if (value == 0 && !allowzero) value = 1e-5; + toreturn[i][j][k] = evaluatedpoint(pt, value); } } } @@ -983,8 +1044,8 @@ T[][] slice(T[][] a, int start1, int end1, int start2, int end2) { return toreturn; } T[][][] slice(T[][][] a, int start1, int end1, - int start2, int end2, - int start3, int end3) { + int start2, int end2, + int start3, int end3) { T[][][] toreturn = new T[end1-start1][][]; for (int i = start1; i < end1; ++i) { toreturn[i-start1] = slice(a[i], start2, end2, start3, end3); @@ -1000,8 +1061,8 @@ T[][] slice(T[][] a, int start1, int end1, int start2, int end2) { return toreturn; } T[][][] slice(T[][][] a, int start1, int end1, - int start2, int end2, - int start3, int end3) { + int start2, int end2, + int start3, int end3) { T[][][] toreturn = new T[end1-start1][][]; for (int i = start1; i < end1; ++i) { toreturn[i-start1] = slice(a[i], start2, end2, start3, end3); @@ -1021,73 +1082,74 @@ struct gridwithzeros { triple grad(triple); real f(triple); int maxdepth; + bool usetriangles; // Populate the edges with zeros that have a sign change and are not already // populated. void fillzeros() { for (int j = 0; j < ny+1; ++j) { for (int k = 0; k < nz+1; ++k) { - real y = corners[0][j][k].pt.y; - real z = corners[0][j][k].pt.z; - real f_along_x(real t) { return f((t, y, z)); } - for (int i = 0; i < nx; ++i) { - if (xdirzeros[i][j][k] != null) continue; - evaluatedpoint start = corners[i][j][k]; - evaluatedpoint end = corners[i+1][j][k]; - if ((start.value > 0 && end.value > 0) || (start.value < 0 && end.value < 0)) - xdirzeros[i][j][k] = null; - else { - triple root = (0,y,z); - root += X * findroot(f_along_x, start.pt.x, end.pt.x, - fa=start.value, fb=end.value); - triple normal = grad(root); - xdirzeros[i][j][k] = positionedvector(root, normal); - } - } + real y = corners[0][j][k].pt.y; + real z = corners[0][j][k].pt.z; + real f_along_x(real t) { return f((t, y, z)); } + for (int i = 0; i < nx; ++i) { + if (xdirzeros[i][j][k] != null) continue; + evaluatedpoint start = corners[i][j][k]; + evaluatedpoint end = corners[i+1][j][k]; + if ((start.value > 0 && end.value > 0) || (start.value < 0 && end.value < 0)) + xdirzeros[i][j][k] = null; + else { + triple root = (0,y,z); + root += X * findroot(f_along_x, start.pt.x, end.pt.x, + fa=start.value, fb=end.value); + triple normal = grad(root); + xdirzeros[i][j][k] = positionedvector(root, normal); + } + } } } for (int i = 0; i < nx+1; ++i) { for (int k = 0; k < nz+1; ++k) { - real x = corners[i][0][k].pt.x; - real z = corners[i][0][k].pt.z; - real f_along_y(real t) { return f((x, t, z)); } - for (int j = 0; j < ny; ++j) { - if (ydirzeros[i][j][k] != null) continue; - evaluatedpoint start = corners[i][j][k]; - evaluatedpoint end = corners[i][j+1][k]; - if ((start.value > 0 && end.value > 0) || (start.value < 0 && end.value < 0)) - ydirzeros[i][j][k] = null; - else { - triple root = (x,0,z); - root += Y * findroot(f_along_y, start.pt.y, end.pt.y, - fa=start.value, fb=end.value); - triple normal = grad(root); - ydirzeros[i][j][k] = positionedvector(root, normal); - } - } + real x = corners[i][0][k].pt.x; + real z = corners[i][0][k].pt.z; + real f_along_y(real t) { return f((x, t, z)); } + for (int j = 0; j < ny; ++j) { + if (ydirzeros[i][j][k] != null) continue; + evaluatedpoint start = corners[i][j][k]; + evaluatedpoint end = corners[i][j+1][k]; + if ((start.value > 0 && end.value > 0) || (start.value < 0 && end.value < 0)) + ydirzeros[i][j][k] = null; + else { + triple root = (x,0,z); + root += Y * findroot(f_along_y, start.pt.y, end.pt.y, + fa=start.value, fb=end.value); + triple normal = grad(root); + ydirzeros[i][j][k] = positionedvector(root, normal); + } + } } } for (int i = 0; i < nx+1; ++i) { for (int j = 0; j < ny+1; ++j) { - real x = corners[i][j][0].pt.x; - real y = corners[i][j][0].pt.y; - real f_along_z(real t) { return f((x, y, t)); } - for (int k = 0; k < nz; ++k) { - if (zdirzeros[i][j][k] != null) continue; - evaluatedpoint start = corners[i][j][k]; - evaluatedpoint end = corners[i][j][k+1]; - if ((start.value > 0 && end.value > 0) || (start.value < 0 && end.value < 0)) - zdirzeros[i][j][k] = null; - else { - triple root = (x,y,0); - root += Z * findroot(f_along_z, start.pt.z, end.pt.z, - fa=start.value, fb=end.value); - triple normal = grad(root); - zdirzeros[i][j][k] = positionedvector(root, normal); - } - } + real x = corners[i][j][0].pt.x; + real y = corners[i][j][0].pt.y; + real f_along_z(real t) { return f((x, y, t)); } + for (int k = 0; k < nz; ++k) { + if (zdirzeros[i][j][k] != null) continue; + evaluatedpoint start = corners[i][j][k]; + evaluatedpoint end = corners[i][j][k+1]; + if ((start.value > 0 && end.value > 0) || (start.value < 0 && end.value < 0)) + zdirzeros[i][j][k] = null; + else { + triple root = (x,y,0); + root += Z * findroot(f_along_z, start.pt.z, end.pt.z, + fa=start.value, fb=end.value); + triple normal = grad(root); + zdirzeros[i][j][k] = positionedvector(root, normal); + } + } } } } @@ -1097,14 +1159,15 @@ struct gridwithzeros { // maximum subdivision depth. When a cube at maxdepth cannot be resolved to // patches, it is left empty. void operator init(int nx, int ny, int nz, - real f(triple), triple a, triple b, - int maxdepth = 6) { + real f(triple), triple a, triple b, + int maxdepth = 6, bool usetriangles) { this.nx = nx; this.ny = ny; this.nz = nz; grad = nGrad(f); this.f = f; this.maxdepth = maxdepth; + this.usetriangles = usetriangles; corners = make3dgrid(a, b, nx, ny, nz, f); xdirzeros = new positionedvector[nx][ny+1][nz+1]; ydirzeros = new positionedvector[nx+1][ny][nz+1]; @@ -1112,11 +1175,11 @@ struct gridwithzeros { for (int i = 0; i <= nx; ++i) { for (int j = 0; j <= ny; ++j) { - for (int k = 0; k <= nz; ++k) { - if (i < nx) xdirzeros[i][j][k] = null; - if (j < ny) ydirzeros[i][j][k] = null; - if (k < nz) zdirzeros[i][j][k] = null; - } + for (int k = 0; k <= nz; ++k) { + if (i < nx) xdirzeros[i][j][k] = null; + if (j < ny) ydirzeros[i][j][k] = null; + if (k < nz) zdirzeros[i][j][k] = null; + } } } @@ -1144,18 +1207,18 @@ struct gridwithzeros { corners = new evaluatedpoint[nx+1][ny+1][nz+1]; for (int i = 0; i <= nx; ++i) { for (int j = 0; j <= ny; ++j) { - for (int k = 0; k <= nz; ++k) { - if (i % 2 == 0 && j % 2 == 0 && k % 2 == 0) { - corners[i][j][k] = oldcorners[quotient(i,2)][quotient(j,2)][quotient(k,2)]; - } else { - triple pt = (interp(a.x, b.x, i/nx), - interp(a.y, b.y, j/ny), - interp(a.z, b.z, k/nz)); - real value = f(pt); - if (value == 0) value = 1e-5; - corners[i][j][k] = evaluatedpoint(pt, value); - } - } + for (int k = 0; k <= nz; ++k) { + if (i % 2 == 0 && j % 2 == 0 && k % 2 == 0) { + corners[i][j][k] = oldcorners[quotient(i,2)][quotient(j,2)][quotient(k,2)]; + } else { + triple pt = (interp(a.x, b.x, i/nx), + interp(a.y, b.y, j/ny), + interp(a.z, b.z, k/nz)); + real value = f(pt); + if (value == 0) value = 1e-5; + corners[i][j][k] = evaluatedpoint(pt, value); + } + } } } @@ -1163,23 +1226,23 @@ struct gridwithzeros { xdirzeros = new positionedvector[nx][ny+1][nz+1]; for (int i = 0; i < nx; ++i) { for (int j = 0; j < ny + 1; ++j) { - for (int k = 0; k < nz + 1; ++k) { - if (j % 2 != 0 || k % 2 != 0) { - xdirzeros[i][j][k] = null; - } else { - positionedvector zero = oldxdir[quotient(i,2)][quotient(j,2)][quotient(k,2)]; - if (zero == null) { - xdirzeros[i][j][k] = null; - continue; - } - real x = zero.position.x; - if (x > interp(a.x, b.x, i/nx) && x < interp(a.x, b.x, (i+1)/nx)) { - xdirzeros[i][j][k] = zero; - } else { - xdirzeros[i][j][k] = null; - } - } - } + for (int k = 0; k < nz + 1; ++k) { + if (j % 2 != 0 || k % 2 != 0) { + xdirzeros[i][j][k] = null; + } else { + positionedvector zero = oldxdir[quotient(i,2)][quotient(j,2)][quotient(k,2)]; + if (zero == null) { + xdirzeros[i][j][k] = null; + continue; + } + real x = zero.position.x; + if (x > interp(a.x, b.x, i/nx) && x < interp(a.x, b.x, (i+1)/nx)) { + xdirzeros[i][j][k] = zero; + } else { + xdirzeros[i][j][k] = null; + } + } + } } } @@ -1187,23 +1250,23 @@ struct gridwithzeros { ydirzeros = new positionedvector[nx+1][ny][nz+1]; for (int i = 0; i < nx+1; ++i) { for (int j = 0; j < ny; ++j) { - for (int k = 0; k < nz + 1; ++k) { - if (i % 2 != 0 || k % 2 != 0) { - ydirzeros[i][j][k] = null; - } else { - positionedvector zero = oldydir[quotient(i,2)][quotient(j,2)][quotient(k,2)]; - if (zero == null) { - ydirzeros[i][j][k] = null; - continue; - } - real y = zero.position.y; - if (y > interp(a.y, b.y, j/ny) && y < interp(a.y, b.y, (j+1)/ny)) { - ydirzeros[i][j][k] = zero; - } else { - ydirzeros[i][j][k] = null; - } - } - } + for (int k = 0; k < nz + 1; ++k) { + if (i % 2 != 0 || k % 2 != 0) { + ydirzeros[i][j][k] = null; + } else { + positionedvector zero = oldydir[quotient(i,2)][quotient(j,2)][quotient(k,2)]; + if (zero == null) { + ydirzeros[i][j][k] = null; + continue; + } + real y = zero.position.y; + if (y > interp(a.y, b.y, j/ny) && y < interp(a.y, b.y, (j+1)/ny)) { + ydirzeros[i][j][k] = zero; + } else { + ydirzeros[i][j][k] = null; + } + } + } } } @@ -1211,23 +1274,23 @@ struct gridwithzeros { zdirzeros = new positionedvector[nx+1][ny+1][nz]; for (int i = 0; i < nx + 1; ++i) { for (int j = 0; j < ny + 1; ++j) { - for (int k = 0; k < nz; ++k) { - if (i % 2 != 0 || j % 2 != 0) { - zdirzeros[i][j][k] = null; - } else { - positionedvector zero = oldzdir[quotient(i,2)][quotient(j,2)][quotient(k,2)]; - if (zero == null) { - zdirzeros[i][j][k] = null; - continue; - } - real z = zero.position.z; - if (z > interp(a.z, b.z, k/nz) && z < interp(a.z, b.z, (k+1)/nz)) { - zdirzeros[i][j][k] = zero; - } else { - zdirzeros[i][j][k] = null; - } - } - } + for (int k = 0; k < nz; ++k) { + if (i % 2 != 0 || j % 2 != 0) { + zdirzeros[i][j][k] = null; + } else { + positionedvector zero = oldzdir[quotient(i,2)][quotient(j,2)][quotient(k,2)]; + if (zero == null) { + zdirzeros[i][j][k] = null; + continue; + } + real z = zero.position.z; + if (z > interp(a.z, b.z, k/nz) && z < interp(a.z, b.z, (k+1)/nz)) { + zdirzeros[i][j][k] = zero; + } else { + zdirzeros[i][j][k] = null; + } + } + } } } @@ -1254,14 +1317,14 @@ struct gridwithzeros { void pushifnonnull(positionedvector v) { if (v != null) { - zeroedges.push(edge(currentface, nextface)); - zeros.push(v); + zeroedges.push(edge(currentface, nextface)); + zeros.push(v); } } positionedvector findzero(int face1, int face2) { edge e = edge(face1, face2); for (int i = 0; i < zeroedges.length; ++i) { - if (zeroedges[i] == e) return zeros[i]; + if (zeroedges[i] == e) return zeros[i]; } return null; } @@ -1303,7 +1366,7 @@ struct gridwithzeros { patch[] subdividecube() { if (!subdivide()) { - return new patch[0]; + return new patch[0]; } return draw(reportactive); } @@ -1324,19 +1387,31 @@ struct gridwithzeros { path3 edgecycle; for (int i = 0; i < faceorder.length; ++i) { path3 currentpath = pathinface(patchcorners[i], patchcorners[i+1], - faceorder[i+1], faceorder[i], - faceorder[i+2]); + faceorder[i+1], faceorder[i], + faceorder[i+2]); triple testpoint = point(currentpath, 0.5); if (!checkpt(testpoint, f, grad, corners[0][0][0], corners[1][1][1])) { - return subdividecube(); + return subdividecube(); } edgecycle = edgecycle & currentpath; } edgecycle = edgecycle & cycle; + + { // Ensure the outward normals are pointing in the same direction as the gradient. + triple tangentin = patchcorners[0].position - precontrol(edgecycle, 0); + triple tangentout = postcontrol(edgecycle, 0) - patchcorners[0].position; + triple normal = cross(tangentin, tangentout); + if (dot(normal, patchcorners[0].direction) < 0) { + edgecycle = reverse(edgecycle); + patchcorners = patchcorners[-sequence(patchcorners.length)]; + patchcorners.cyclic = true; + } + } + patch[] toreturn = quadpatches(edgecycle, patchcorners, f, grad, - corners[0][0][0], corners[1][1][1]); + corners[0][0][0], corners[1][1][1], usetriangles); if (alias(toreturn, null)) return subdividecube(); return toreturn; } @@ -1351,6 +1426,7 @@ struct gridwithzeros { cube.ny = 1; cube.nz = 1; cube.maxdepth = maxdepth; + cube.usetriangles = usetriangles; cube.corners = slice(corners,i,i+2,j,j+2,k,k+2); cube.xdirzeros = slice(xdirzeros,i,i+1,j,j+2,k,k+2); cube.ydirzeros = slice(ydirzeros,i,i+2,j,j+1,k,k+2); @@ -1372,6 +1448,7 @@ struct gridwithzeros { // grid will subdivide all the way to maxdepth if necessary to find points // on the surface. draw = new patch[](bool[] reportactive = null) { + if (alias(reportactive, null)) progress(true); // A list of all the patches not already drawn but known // to contain part of the surface. This "queue" is // actually implemented as stack for simplicity, since @@ -1382,49 +1459,49 @@ struct gridwithzeros { bool[][][] enqueued = new bool[nx][ny][nz]; for (int i = 0; i < enqueued.length; ++i) { for (int j = 0; j < enqueued[i].length; ++j) { - for (int k = 0; k < enqueued[i][j].length; ++k) { - enqueued[i][j][k] = false; - } + for (int k = 0; k < enqueued[i][j].length; ++k) { + enqueued[i][j][k] = false; + } } } void enqueue(int i, int j, int k) { if (i >= 0 && i < nx - && j >= 0 && j < ny - && k >= 0 && k < nz - && !enqueued[i][j][k]) { - queue.push((i,j,k)); - enqueued[i][j][k] = true; + && j >= 0 && j < ny + && k >= 0 && k < nz + && !enqueued[i][j][k]) { + queue.push((i,j,k)); + enqueued[i][j][k] = true; } if (!alias(reportactive, null)) { - if (i < 0) reportactive[XLOW] = true; - if (i >= nx) reportactive[XHIGH] = true; - if (j < 0) reportactive[YLOW] = true; - if (j >= ny) reportactive[YHIGH] = true; - if (k < 0) reportactive[ZLOW] = true; - if (k >= nz) reportactive[ZHIGH] = true; + if (i < 0) reportactive[XLOW] = true; + if (i >= nx) reportactive[XHIGH] = true; + if (j < 0) reportactive[YLOW] = true; + if (j >= ny) reportactive[YHIGH] = true; + if (k < 0) reportactive[ZLOW] = true; + if (k >= nz) reportactive[ZHIGH] = true; } } for (int i = 0; i < nx+1; ++i) { for (int j = 0; j < ny+1; ++j) { - for (int k = 0; k < nz+1; ++k) { - if (i < nx && xdirzeros[i][j][k] != null) { - for (int jj = j-1; jj <= j; ++jj) - for (int kk = k-1; kk <= k; ++kk) - enqueue(i, jj, kk); - } - if (j < ny && ydirzeros[i][j][k] != null) { - for (int ii = i-1; ii <= i; ++ii) - for (int kk = k-1; kk <= k; ++kk) - enqueue(ii, j, kk); - } - if (k < nz && zdirzeros[i][j][k] != null) { - for (int ii = i-1; ii <= i; ++ii) - for (int jj = j-1; jj <= j; ++jj) - enqueue(ii, jj, k); - } - } + for (int k = 0; k < nz+1; ++k) { + if (i < nx && xdirzeros[i][j][k] != null) { + for (int jj = j-1; jj <= j; ++jj) + for (int kk = k-1; kk <= k; ++kk) + enqueue(i, jj, kk); + } + if (j < ny && ydirzeros[i][j][k] != null) { + for (int ii = i-1; ii <= i; ++ii) + for (int kk = k-1; kk <= k; ++kk) + enqueue(ii, j, kk); + } + if (k < nz && zdirzeros[i][j][k] != null) { + for (int ii = i-1; ii <= i; ++ii) + for (int jj = j-1; jj <= j; ++jj) + enqueue(ii, jj, k); + } + } } } @@ -1448,9 +1525,9 @@ struct gridwithzeros { if (reportface[ZLOW]) enqueue(i,j,k-1); if (reportface[ZHIGH]) enqueue(i,j,k+1); surface.append(toappend); - if (settings.verbose > 1 && alias(reportactive, null)) write(stdout, '.'); + if (alias(reportactive, null)) progress(); } - if (settings.verbose > 1 && alias(reportactive, null)) write(stdout, '\n'); + if (alias(reportactive, null)) progress(false); return surface; }; } @@ -1478,22 +1555,25 @@ struct gridwithzeros { // maxdepth - the maximum depth to which the algorithm will subdivide in // an effort to find patches that closely approximate the true surface. surface implicitsurface(real f(triple) = null, real ff(real,real,real) = null, - triple a, triple b, - int n = nmesh, - bool keyword overlapedges = false, - int keyword nx=n, int keyword ny=n, - int keyword nz=n, - int keyword maxdepth = 8) { + triple a, triple b, + int n = nmesh, + bool keyword overlapedges = false, + int keyword nx=n, int keyword ny=n, + int keyword nz=n, + int keyword maxdepth = 8, + bool keyword usetriangles=true) { if (f == null && ff == null) abort("implicitsurface called without specifying a function."); if (f != null && ff != null) abort("Only specify one function when calling implicitsurface."); if (f == null) f = new real(triple w) { return ff(w.x, w.y, w.z); }; - gridwithzeros grid = gridwithzeros(nx, ny, nz, f, a, b, maxdepth=maxdepth); + gridwithzeros grid = gridwithzeros(nx, ny, nz, f, a, b, maxdepth=maxdepth, + usetriangles=usetriangles); patch[] patches = grid.draw(); if (overlapedges) { for (int i = 0; i < patches.length; ++i) { - triple center = patches[i].point(1/2,1/2); + triple center = (patches[i].triangular ? + patches[i].point(1/3, 1/3) : patches[i].point(1/2,1/2)); patches[i] = shift(center) * scale3(1.01) * shift(-center) * patches[i]; } } diff --git a/Build/source/utils/asymptote/base/three_surface.asy b/Build/source/utils/asymptote/base/three_surface.asy index 175ef562d80..c6d439e9fb0 100644 --- a/Build/source/utils/asymptote/base/three_surface.asy +++ b/Build/source/utils/asymptote/base/three_surface.asy @@ -68,7 +68,9 @@ struct patch { } triple Bu(int j, real u) {return bezier(P[0][j],P[1][j],P[2][j],P[3][j],u);} - triple BuP(int j, real u) {return bezierP(P[0][j],P[1][j],P[2][j],P[3][j],u);} + triple BuP(int j, real u) { + return bezierP(P[0][j],P[1][j],P[2][j],P[3][j],u); + } triple BuPP(int j, real u) { return bezierPP(P[0][j],P[1][j],P[2][j],P[3][j],u); } @@ -82,7 +84,9 @@ struct patch { } triple Bv(int i, real v) {return bezier(P[i][0],P[i][1],P[i][2],P[i][3],v);} - triple BvP(int i, real v) {return bezierP(P[i][0],P[i][1],P[i][2],P[i][3],v);} + triple BvP(int i, real v) { + return bezierP(P[i][0],P[i][1],P[i][2],P[i][3],v); + } triple BvPP(int i, real v) { return bezierPP(P[i][0],P[i][1],P[i][2],P[i][3],v); } @@ -99,13 +103,6 @@ struct patch { return bezier(Bu(0,u),Bu(1,u),Bu(2,u),Bu(3,u),v); } - triple pointtriangular(real u, real v) { - real w=1-u-v; - return w^2*(w*P[0][0]+3*(u*P[1][0]+v*P[1][1]))+ - u^2*(u*P[3][0]+3*(w*P[2][0]+v*P[3][1]))+ - 6*u*v*w*P[2][1]+v^2*(v*P[3][3]+3*(w*P[2][2]+u*P[3][2])); - } - // compute normal vectors for degenerate cases private triple normal0(real u, real v, real epsilon) { triple n=0.5*(cross(bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v), @@ -167,46 +164,92 @@ struct patch { return abs(n) > epsilon ? n : normal0(0,1,epsilon); } - triple normal00triangular() { - triple n=9*cross(P[1][0]-P[0][0],P[1][1]-P[0][0]); - real epsilon=fuzz*change2(P); - return abs(n) > epsilon ? n : normal0(0,0,epsilon); + triple pointtriangular(real u, real v) { + real w=1-u-v; + return w^2*(w*P[0][0]+3*(u*P[1][0]+v*P[1][1]))+ + u^2*(3*(w*P[2][0]+v*P[3][1])+u*P[3][0])+ + 6*u*v*w*P[2][1]+v^2*(3*(w*P[2][2]+u*P[3][2])+v*P[3][3]); } - triple normal10triangular() { - triple n=9*cross(P[3][0]-P[2][0],P[3][1]-P[2][0]); - real epsilon=fuzz*change2(P); - return abs(n) > epsilon ? n : normal0(1,0,epsilon); + triple bu(real u, real v) { + // Compute one-third of the directional derivative of a Bezier triangle + // in the u direction at (u,v). + real w=1-u-v; + return -w^2*P[0][0]+w*(w-2*u)*P[1][0]-2*w*v*P[1][1]+u*(2*w-u)*P[2][0]+ + 2*v*(w-u)*P[2][1]-v^2*P[2][2]+u^2*P[3][0]+2*u*v*P[3][1]+v^2*P[3][2]; } - triple normal01triangular() { - triple n=9*cross(P[3][2]-P[2][2],P[3][3]-P[2][2]); - real epsilon=fuzz*change2(P); - return abs(n) > epsilon ? n : normal0(0,1,epsilon); + triple buu(real u, real v) { + // Compute one-sixth of the second directional derivative of a Bezier + // triangle in the u direction at (u,v). + real w=1-u-v; + return w*P[0][0]+(u-2*w)*P[1][0]+v*P[1][1]+(w-2*u)*P[2][0]-2*v*P[2][1]+ + u*P[3][0]+v*P[3][1]; } - // Compute one-third of the directional derivative of a Bezier triangle in the u - // direction at point (u,v). - private triple bu(real u, real v) { + triple buuu() { + // Compute one-sixth of the third directional derivative of a Bezier + // triangle in the u direction at (u,v). + return -P[0][0]+3*P[1][0]-3*P[2][0]+P[3][0]; + } + + triple bv(real u, real v) { + // Compute one-third of the directional derivative of a Bezier triangle + // in the v direction at (u,v). real w=1-u-v; - return u*(w*2-u)*P[2][0]+2*v*(w-u)*P[2][1]+w*(w-2*u)*P[1][0]+ - u*(u*P[3][0]+2*v*P[3][1])+v*v*P[3][2]-w*(2*v*P[1][1]+w*P[0][0])- - v*v*P[2][2]; + return -w^2*P[0][0]-2*u*w*P[1][0]+w*(w-2*v)*P[1][1]-u^2*P[2][0]+ + 2*u*(w-v)*P[2][1]+v*(2*w-v)*P[2][2]+u*u*P[3][1]+2*u*v*P[3][2]+ + v^2*P[3][3]; } - // Compute one-third of the directional derivative of a Bezier triangle in the v - // direction at point (u,v). - private triple bv(real u, real v) { + triple bvv(real u, real v) { + // Compute one-sixth of the second directional derivative of a Bezier + // triangle in the v direction at (u,v). real w=1-u-v; - return u*2*(w-v)*P[2][1]+v*(2*w-v)*P[2][2]+w*(w-2*v)*P[1][1]+ - u*(u*P[3][1]+2*v*P[3][2])+v*v*P[3][3]-w*(2*u*P[1][0]+w*P[0][0])- - u*u*P[2][0]; + return w*P[0][0]+u*P[1][0]+(v-2*w)*P[1][1]-2*u*P[2][1]+(w-2*v)*P[2][2]+ + u*P[3][2]+v*P[3][3]; + } + + triple bvvv() { + // Compute one-sixth of the third directional derivative of a Bezier + // triangle in the v direction at (u,v). + return -P[0][0]+3*P[1][1]-3*P[2][2]+P[3][3]; + } + + // compute normal vectors for a degenerate Bezier triangle + private triple normaltriangular0(real u, real v, real epsilon) { + triple n=9*(cross(buu(u,v),bv(u,v))+ + cross(bu(u,v),bvv(u,v))); + return abs(n) > epsilon ? n : + 9*cross(buu(u,v),bvv(u,v))+ + 3*(cross(buuu(),bv(u,v))+cross(bu(u,v),bvvv())+ + cross(buuu(),bvv(u,v))+cross(buu(u,v),bvvv()))+ + cross(buuu(),bvvv()); } // Compute the normal of a Bezier triangle at (u,v) triple normaltriangular(real u, real v) { - // TODO: handle degeneracy - return 9*cross(bu(u,v),bv(u,v)); + triple n=9*cross(bu(u,v),bv(u,v)); + real epsilon=fuzz*change2(P); + return (abs(n) > epsilon) ? n : normal0(u,v,epsilon); + } + + triple normal00triangular() { + triple n=9*cross(P[1][0]-P[0][0],P[1][1]-P[0][0]); + real epsilon=fuzz*change2(P); + return abs(n) > epsilon ? n : normaltriangular0(0,0,epsilon); + } + + triple normal10triangular() { + triple n=9*cross(P[3][0]-P[2][0],P[3][1]-P[2][0]); + real epsilon=fuzz*change2(P); + return abs(n) > epsilon ? n : normaltriangular0(1,0,epsilon); + } + + triple normal01triangular() { + triple n=9*cross(P[3][2]-P[2][2],P[3][3]-P[2][2]); + real epsilon=fuzz*change2(P); + return abs(n) > epsilon ? n : normaltriangular0(0,1,epsilon); } pen[] colors(material m, light light=currentlight) { @@ -255,6 +298,8 @@ struct patch { normal10=normal10triangular; normal01=normal01triangular; colors=colorstriangular; + uequals=new path3(real u) {return nullpath3;}; + vequals=new path3(real u) {return nullpath3;}; } } @@ -448,7 +493,7 @@ patch reverse(patch s) assert(!s.triangular); patch S; S.P=transpose(s.P); - if(s.colors.length > 0) + if(s.colors.length > 0) S.colors=new pen[] {s.colors[0],s.colors[3],s.colors[2],s.colors[1]}; S.straight=s.straight; S.planar=s.planar; @@ -817,15 +862,15 @@ struct surface { return ucyclic() ? g&cycle : g; } - // A constructor for a possibly nonconvex simple cyclic path in a given plane. + // A constructor for a possibly nonconvex simple cyclic path in a given + // plane. void operator init(path p, triple plane(pair)=XYplane) { - bool straight=piecewisestraight(p); for(path g : regularize(p)) { if(length(g) == 3) { path3 G=path3(g,plane); - s.push(patch(G,coons3(G))); + s.push(patch(G,coons3(G),planar=true)); } else - s.push(patch(coons(g),plane,straight)); + s.push(patch(coons(g),plane,piecewisestraight(g))); } } @@ -1490,9 +1535,14 @@ void tensorshade(transform t=identity(), frame f, patch s, material m, light light=currentlight, projection P) { - if(s.triangular) s=tensor(s); + pen[] p; + if(s.triangular) { + p=s.colorstriangular(m,light); + p.push(p[0]); + s=tensor(s); + } else p=s.colors(m,light); tensorshade(f,box(t*s.min(P),t*s.max(P)),m.diffuse(), - s.colors(m,light),t*project(s.external(),P,1),t*project(s.internal(),P)); + p,t*project(s.external(),P,1),t*project(s.internal(),P)); } restricted pen[] nullpens={nullpen}; @@ -1536,18 +1586,19 @@ void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1, for(int p=depth.length-1; p >= 0; --p) { real[] a=depth[p]; int k=round(a[1]); + patch S=s.s[k]; pen meshpen=meshpen[k]; - if(!invisible(meshpen)) { + if(!invisible(meshpen) && !S.triangular) { if(group) begingroup3(f,meshname(name),render); meshpen=modifiers+meshpen; real step=nu == 0 ? 0 : 1/nu; for(int i=0; i <= nu; ++i) - draw(f,s.s[k].uequals(i*step),meshpen,meshlight,partname(i,render), + draw(f,S.uequals(i*step),meshpen,meshlight,partname(i,render), render); step=nv == 0 ? 0 : 1/nv; for(int j=0; j <= nv; ++j) - draw(f,s.s[k].vequals(j*step),meshpen,meshlight,partname(j,render), + draw(f,S.vequals(j*step),meshpen,meshlight,partname(j,render), render); if(group) endgroup3(f); @@ -1629,8 +1680,9 @@ void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1, pen modifiers; if(is3D()) modifiers=thin()+squarecap; for(int k=0; k < s.s.length; ++k) { + patch S=s.s[k]; pen meshpen=meshpen[k]; - if(!invisible(meshpen)) { + if(!invisible(meshpen) && !S.triangular) { meshpen=modifiers+meshpen; real step=nu == 0 ? 0 : 1/nu; for(int i=0; i <= nu; ++i) |