(* Version dated 12/18/2000 7/2/97 *) (*************************************************************************** * de Casteljau's algorithm for rational triangular and rectangular surface patches (11/25/2002) **************************************************************************) (*************************************************************************** REMEMBER: To run programs involving rectangular surfaces, you need to input radecas.m All points must have a weight. If polynomial surface, weight = 1. **************************************************************************) (* Given a reference triangle (r, s, t), the barycentric coordinates *) (* of r are (1, 0, 0), of s are (0, 1, 0), and of t are (0, 0, 1) *) (* In terms of an affine frame, t is the origin, r is on the x-axis *) (* and s is on the y-axis. Wrt this frame, if p has barycentric coords *) (* (lambda, mu, nu), then the affine coordinates are (lambda, mu) *) (* Pattern of calls for the various rendering procedures *) (* if light = 1, light triangulation trianglight, *) (* else ltriang *) (* *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* sqrender[{net__}, mm_, a_, b_, c_, d_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* This is wrt the square [c, a] x [d, b] *) (* *) (* rendersq[{net__}, mm_, a_, b_, c_, d_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_,debug_] := *) (* Same as sqrender, but wrt [a, b] x [c, d], more intuitive order *) (* *) (* trender[{net__}, mm_, {newtrig_}, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* render1[{net__}, mm_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* render3[{net__}, mm_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* render4[{net__}, mm_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz,debug_] *) (* *) (* render5[{net__}, mm_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* render6[{net__}, mm_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* render34[{net__}, mm_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* render56[{net__}, mm_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* render3to6[{net__}, mm_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_,debug] *) (* *) (* lrender renders the surface by u-curves and v-curves *) (* wrt square [-1, 1] x [-1, 1] *) (* *) (* lrender[{net__},mm_,n_,ni_,lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* *) (* lfrender renders a closed surface by u-curves and v-curves *) (* *) (* lfrender[{net__},mm_,n_,ni_,lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* render[{net__}, mm_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* render3col[{net__}, mm_, n_, light_, fcnet_, *) (* lx_, mx_,ly_,my_,lz_,mz_,debug_] := *) (* *) (* new2net[{net__}, mm_, a_, b_, n_, *) (* lx_,mx_,ly_,my_,lz_,mz_, debug_] *) (* *) (* subdiv4[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* To render a square patch *) (* recrender[{net__}, p_, q_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* rfrac6[{net__}, mm_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_] *) (* *) (* sqfractal[{net__}, mm_, a_, b_, c_, d_, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_] *) (* *) (* tfractal[{net__}, mm_, {newtrig__}, n_, light_, fcnet_, *) (* lx_,mx_,ly_,my_,lz_,mz_] *) (* *) (* fractal4[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,light_] *) (* *) (* fractal4 is the same as tfractal for (r,s,t) *) (* *) (* fractal2[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,light_] *) (* *) (* fractal2 is the same as sqfractal for standard square *) (* *) (* fastdiv6 is the "spider-style" subdivision *) (* fastdiv6[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,light_] *) (* *) (* fastdiv4 is the "diamond-style" subdivision *) (* fastdiv4[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,light_] *) (* *) (* clocurv3D renders a closed curve on a surface *) (* specified by points a, b in the parameter plane *) (* clocurv3D[{net__},mm_,n_,a_,b_,lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* scurv3D renders a curve segment on a surface, specified by *) (* points a, b in the parameter plane *) (* scurv3D[{net__},mm_,n_,a_,b_,lx_,mx_,ly_,my_,lz_,mz_,debug_] *) (* *) (* Some useful reference triangles *) reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; reftrig = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; reftrig4 = {{1, -1, 1}, {-1, 1, 1}, {-1, -1, 3}}; reftrigb = {{1, -1}, {1, 1}, {-1, -1}} (* To turn on and off the simulated lighting *) (* light = 4 or 2 or > 10 corresponds to explicit colors, otherwise simulated illumination *) (* light = 4 corresponds to four colors, light = 2 to coloring the inside and outside *) (* light = 20 corresponds to one patch blue, one patch red *) lightkind[light_] := Block[ {res}, If[light === 4 || light === 2 || light > 10, res = Lighting -> False , res = Lighting -> True ]; res ]; (* To turn on and off coloring the inside and outside *) (* light = 2 corresponds to coloring the inside and the outside of the surface, red and green *) facecolors[light_] := Block[ {res}, If[light === 2, res = {FaceForm[RGBColor[1, 0, 0], RGBColor[0, 1, 0]]}, res = {Thickness[0.0006]} ]; res[[1]] ]; (* To build a full triangulation from a triangular net *) (* in terms of line segments *) triangul[{cnet__},m_] := Block[ {cc = {cnet}, row1, row2, n, i, j, lseg = {}, pt, edge}, ( row2 = {}; Do[ pt = cc[[j]]; row2 = Append[row2, pt], {j, 1, m + 1} ]; Do[n = m + 2 - i; cc = Drop[cc, n]; row1 = row2; row2 = {}; Do[ pt = cc[[j]]; row2 = Append[row2, pt], {j, 1, n - 1} ]; Do[ edge = {row1[[j]], row1[[j + 1]]}; lseg = Append[lseg, edge], {j, 1, n - 1} ]; Do[ edge = {row2[[j]], row1[[j]]}; lseg = Append[lseg, edge]; edge = {row2[[j]], row1[[j + 1]]}; lseg = Append[lseg, edge], {j, 1, n - 1} ], {i, 1, m} ]; lseg ) ]; (* To build a full triangulation from a triangular net *) (* in terms of solid nontransparent triangles *) soltriangul[{cnet__},m_] := Block[ {cc = {cnet}, row1, row2, n, i, j, lseg = {}, pt, triangle}, ( row2 = {}; Do[ pt = cc[[j]]; row2 = Append[row2, pt], {j, 1, m + 1} ]; Do[n = m + 2 - i; cc = Drop[cc, n]; row1 = row2; row2 = {}; Do[ pt = cc[[j]]; row2 = Append[row2, pt], {j, 1, n - 1} ]; Do[ triangle = {row1[[j]], row2[[j]], row1[[j + 1]]}; lseg = Append[lseg, triangle], {j, 1, n - 1} ]; Do[ triangle = {row2[[j]], row2[[j + 1]], row1[[j + 1]]}; lseg = Append[lseg, triangle], {j, 1, n - 2} ], {i, 1, m} ]; lseg ) ]; (* To build a partial triangulation from a triangular net *) (* omits intermediate horizontal segments *) ftriangul[{cnet__},m_] := Block[ {cc = {cnet}, row1, row2, n, i, j, lseg = {}, pt, edge}, (row2 = {}; Do[ pt = cc[[j]]; row2 = Append[row2, pt], {j, 1, m + 1} ]; Do[n = m + 2 - i; cc = Drop[cc, n]; row1 = row2; row2 = {}; Do[ pt = cc[[j]]; row2 = Append[row2, pt], {j, 1, n - 1} ]; Do[ If[n === m + 1, edge = {row1[[j]], row1[[j + 1]]}; lseg = Append[lseg, edge] ], {j, 1, n - 1} ]; Do[ edge = {row2[[j]], row1[[j]]}; lseg = Append[lseg, edge]; edge = {row2[[j]], row1[[j + 1]]}; lseg = Append[lseg, edge], {j, 1, n - 1} ], {i, 1, m} ]; lseg ) ]; (* Links the nodes of the left and right sides of a triangular net *) connects[{cnet__},m_] := Block[ {cc = {cnet}, row1, row2, n, j, lseg = {}, pt1, pt2, edge}, ( row1 = {}; row2 = {}; Do[n = m + 2 - j; pt1 = cc[[1]]; pt2 = cc[[n]]; row1 = Prepend[row1, pt1]; row2 = Prepend[row2, pt2]; If[j =!= m + 1, cc = Drop[cc,n]], {j, 1, m + 1} ]; Do[ edge = {row1[[j]], row1[[j + 1]]}; lseg = Append[lseg, edge], {j, 1, m} ]; Do[ edge = {row2[[j]], row2[[j + 1]]}; lseg = Append[lseg, edge], {j, 1, m} ]; lseg ) ]; (* Links the nodes of the three outer edges of a triangular net *) (* First collects the nodes on the base (s, t) (NOTE; s before t) in row3 *) (* Then collects the nodes on the sides (t, r) and (s, r) in row1 and row2 *) (* Returns the Join of the lists in the order (s, t), (t, r), (s, r) *) trilink[{cnet__},m_] := Block[ {cc = {cnet}, row1, row2, row3, n, j, lseg = {}, pt1, pt2, edge}, ( row1 = {}; row2 = {}; row3 = {}; Do[ pt1 = cc[[j]]; row3 = Append[row3, pt1], {j, 1, m + 1} ]; Do[n = m + 2 - j; pt1 = cc[[1]]; pt2 = cc[[n]]; row1 = Prepend[row1, pt1]; row2 = Prepend[row2, pt2]; If[j =!= m + 1, cc = Drop[cc,n]], {j, 1, m + 1} ]; Do[ edge = {row1[[j]], row1[[j + 1]]}; lseg = Append[lseg, edge], {j, 1, m} ]; Do[ edge = {row2[[j]], row2[[j + 1]]}; lseg = Append[lseg, edge], {j, 1, m} ]; Do[ edge = {row3[[j]], row3[[j + 1]]}; lseg = Append[lseg, edge], {j, 1, m} ]; lseg ) ]; (* To build a list of Line segments from a list of edges *) buildmesh[{cnet__}] := Block[ {cc = {cnet}, ll, edge, i, lseg = {}, res}, (ll = Length[cc]; Do[ edge = cc[[i]]; If[Length[edge] === 2, lseg = Append[lseg, Line[edge]], lseg = Append[lseg, {RGBColor[1,0,0], Point[edge]}] ], {i, 1, ll} ]; lseg ) ]; (* To build a list of Opaque triangles from a list of triples of nodes *) solbuildmesh[{cnet__}] := Block[ {cc = {cnet}, ll, triangle, i, lseg = {}, res}, (ll = Length[cc]; Do[ triangle = cc[[i]]; If[Length[triangle] === 3 && Length[triangle[[1]]] === 3, lseg = Append[lseg, Polygon[triangle]], lseg = Append[lseg, {RGBColor[1,0,0], Point[triangle]}] ], {i, 1, ll} ]; lseg ) ]; (* To test that a triangular net really has size (m + 1)(m + 2)/2, *) (* when the polar degree is m *) testm[cnet__, m_] := Block[ {cc = cnet, ll, s, i, stop, res}, (ll = Length[cc]; s = N[Sqrt[2 * ll]]; i = 1; stop = 1; res = 1; While[i <= s && stop === 1, If[(res === ll) && (i === m + 1), stop = 0, i = i + 1; res = res + i] ]; If[stop === 1, Print["*** Net Size: ", ll, " Inconsistent with Surface Degree: ", m, " ***"]]; stop ) ]; (* Computes the polar degree m associated with a net of sixe (m + 1)*(m + 2)/2 *) netsize[cnet__] := Block[ {cc = cnet, ll, s, i, stop, res}, (ll = Length[cc]; s = N[Sqrt[2 * ll]]; i = 1; stop = 1; res = 1; While[i <= s && stop === 1, If[(res === ll), stop = 0, i = i + 1; res = res + i] ]; res = i - 1; res ) ]; (* Builds a list of line segments for the union of lists of nets *) (* This version calls connects *) (* Also adds the three corners of *) (* each pyramidal net as list elements *) listriang[{cnet__},m_] := Block[ {cc = {cnet}, ll, i, lineseg = {}, edglist, linelis, edge, pt}, ( ll = Length[cc]; Do[ edgelist = cc[[i]]; linelis = connects[edgelist,m]; edge = Last[linelis]; pt = Last[edge]; linelis = Append[linelis,pt]; edge = First[linelis]; pt = First[edge]; linelis = Append[linelis,pt]; edge = linelis[[m]]; pt = Last[edge]; linelis = Append[linelis,pt]; lineseg = Join[lineseg, linelis], {i, 1, ll} ]; lineseg = buildmesh[lineseg]; lineseg ) ]; (* Builds a list of line segments for the union of lists of nets *) (* This version calls trilink *) (* Also adds the three corners of *) (* each pyramidal net as list elements *) (* trilink returns the lists in the order (s, t), (t, r), (s, r) *) listrig[{cnet__},m_] := Block[ {cc = {cnet}, ll, i, lineseg = {}, edglist, linelis, edge, pt, res}, ( ll = Length[cc]; Do[ edgelist = cc[[i]]; linelis = trilink[edgelist,m]; edge = Last[linelis]; pt = Last[edge]; (* this is r *) linelis = Append[linelis,pt]; edge = First[linelis]; pt = First[edge]; (* this is s *) linelis = Append[linelis,pt]; edge = linelis[[m]]; pt = Last[edge]; (* this is t *) linelis = Append[linelis,pt]; lineseg = Join[lineseg, linelis], {i, 1, ll} ]; lineseg = buildmesh[lineseg]; lineseg ) ]; (* Builds a list of line segments for the union of lists of nets *) (* This version calls trilink *) (* Also adds the three corners of *) (* each pyramidal net as list elements *) (* uses different colors on sides (blue, orange, red, green *) (* trilink returns the lists in the order (s, t), (t, r), (s, r) *) trianglight[{cnet__},oldm_] := Block[ {cc = {cnet}, ll, l, i, j, lineseg = {}, edglist, m, linelis, edge, pt}, ( ll = Length[cc]; l = ll/4; Do[ Do[ m = netsize[cc[[j]]]; edgelist= projnet[cc[[j]]]; linelis = trilink[edgelist,m]; edge = Last[linelis]; pt = Last[edge]; (* this is r *) linelis = Append[linelis,pt]; edge = First[linelis]; pt = First[edge]; (* this is s *) linelis = Append[linelis,pt]; edge = linelis[[m]]; pt = Last[edge]; (* this is t *) linelis = Append[linelis,pt]; linelis = buildmesh[linelis]; If[j === 1, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[j === 2, linelis = Prepend[linelis, RGBColor[1, 0.8, 0]]]; (* orange *) If[j === 3, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) If[j === 4, linelis = Prepend[linelis, RGBColor[0, 0.7, 0]]]; (* green *) lineseg = Join[lineseg, linelis], {j, 1, 4} ]; cc = Drop[cc, 4], {i, 1, l} ]; lineseg ) ]; (* Calls segltriang if switch === -1 (full triangulation by line segments), otherwise if switch === 4 triangulation by solid colored triangles, otherwise ltriangsolid (default lighting, or two colors). Only called when switch =!= 1 *) (* When switch === 3, omit vertex points *) (* When switch > 10, unique color, blue = 11, red = 12, green = 13, yellow = 14 *) ltriang[switch_,{cnet__},oldm_] := Block[ {cc = {cnet}, res}, ( If[switch === -1, res = segltriang[cc,oldm], If[switch === 4, res = ltriang4solid[cc,oldm], If[switch > 10, res = ltriangcolor[cc,oldm,switch], res = ltriangsolid[switch,cc,oldm] ] ] ]; res ) ]; (* Builds a list of line segments for the union of lists of nets *) (* This version calls triangul *) (* Also adds the three corners of *) (* each pyramidal net as list elements *) ltriang1[{cnet__},oldm_] := Block[ {cc = {cnet}, ll, m, i, lineseg = {}, edglist, linelis, edge, pt}, ( ll = Length[cc]; Do[ m = netsize[cc[[i]]]; edgelist= projnet[cc[[i]]]; linelis = triangul[edgelist,m]; edge = Last[linelis]; pt = First[edge]; linelis = Append[linelis,pt]; edge = First[linelis]; pt = First[edge]; linelis = Append[linelis,pt]; edge = linelis[[m]]; pt = Last[edge]; linelis = Append[linelis,pt]; lineseg = Join[lineseg, linelis], {i, 1, ll} ]; lineseg = buildmesh[lineseg]; lineseg ) ]; (* Builds a list of line segments for the union of lists of nets *) (* This version calls triangul *) (* Also adds the three corners of *) (* each pyramidal net as list elements *) (* uses different colors on four subpatches *) (* actually, patch 2 gets the green and patch 4 the orange! *) segltriang[{cnet__},oldm_] := Block[ {cc = {cnet}, ll, l, i, j, lineseg = {}, edglist, linelis, m, edge, pt}, ( ll = Length[cc]; l = ll/4; Do[ Do[ m = netsize[cc[[j]]]; edgelist= projnet[cc[[j]]]; linelis = triangul[edgelist,m]; edge = Last[linelis]; pt = First[edge]; linelis = Append[linelis,pt]; edge = First[linelis]; pt = First[edge]; linelis = Append[linelis,pt]; edge = linelis[[m]]; pt = Last[edge]; linelis = Append[linelis,pt]; linelis = buildmesh[linelis]; If[j === 1, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[j === 2, linelis = Prepend[linelis, RGBColor[1, 0.8, 0]]]; (* orange *) If[j === 3, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) If[j === 4, linelis = Prepend[linelis, RGBColor[0, 0.7, 0]]]; (* green *) lineseg = Join[lineseg, linelis], {j, 1, 4} ]; cc = Drop[cc, 4], {i, 1, l} ]; lineseg ) ]; (* Builds a list of solid triangles for the union of lists of nets *) (* This version calls soltriangul *) (* Also adds the three corners of *) (* each pyramidal net as list elements *) (* If switch === 3, do not construct corner points *) ltriangsolid[switch_,{cnet__},oldm_] := Block[ {cc = {cnet}, l, i, lineseg = {}, edglist, linelis, m, edge, pt}, ( l = Length[cc]; Do[ m = netsize[cc[[i]]]; edgelist= projnet[cc[[i]]]; linelis = soltriangul[edgelist,m]; If[switch =!= 3, edge = Last[linelis]; pt = edge[[2]]; (* top corner *) linelis = Append[linelis,pt] ]; If[switch =!= 3, edge = First[linelis]; pt = First[edge]; (* lower left corner *) linelis = Append[linelis,pt] ]; If[switch =!= 3, edge = linelis[[m]]; pt = Last[edge]; (* lower right corner *) linelis = Append[linelis,pt] ]; linelis = solbuildmesh[linelis]; lineseg = Join[lineseg, linelis], {i, 1, l} ]; lineseg ) ]; (* Builds a list of solid triangles for the union of lists of nets *) (* This version calls soltriangul *) (* Also adds the three corners of *) (* each pyramidal net as list elements *) (* uses different colors on four subpatches *) ltriang4solid[{cnet__},oldm_] := Block[ {cc = {cnet}, ll, l, i, j, lineseg = {}, edglist, linelis, m, edge, pt}, ( ll = Length[cc]; l = ll/4; Do[ Do[ m = netsize[cc[[j]]]; edgelist= projnet[cc[[j]]]; linelis = soltriangul[edgelist,m]; edge = Last[linelis]; pt = edge[[2]]; (* top corner *) linelis = Append[linelis,pt]; edge = First[linelis]; pt = First[edge]; (* lower left corner *) linelis = Append[linelis,pt]; edge = linelis[[m]]; pt = Last[edge]; (* lower right corner *) linelis = Append[linelis,pt]; linelis = solbuildmesh[linelis]; If[j === 1, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[j === 2, linelis = Prepend[linelis, RGBColor[1, 1, 0]]]; (* yellow *) If[j === 3, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) If[j === 4, linelis = Prepend[linelis, RGBColor[0, 1, 0]]]; (* green *) lineseg = Join[lineseg, linelis], {j, 1, 4} ]; cc = Drop[cc, 4], {i, 1, l} ]; lineseg ) ]; (* Builds a list of solid triangles for the union of lists of nets *) (* This version calls soltriangul *) (* assigns a unique color to all subpatches *) (* Either blue = 11, red = 12, green = 13, or yellow = 14 *) ltriangcolor[{cnet__},oldm_,col_] := Block[ {cc = {cnet}, ll, i, lineseg = {}, edglist, linelis, m, color}, ( color = col; ll = Length[cc]; Do[ m = netsize[cc[[i]]]; edgelist= projnet[cc[[i]]]; linelis = soltriangul[edgelist,m]; linelis = solbuildmesh[linelis]; If[color === 11, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[color === 12, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) If[color === 13, linelis = Prepend[linelis, RGBColor[0, 1, 0]]]; (* green *) If[color === 14, linelis = Prepend[linelis, RGBColor[1, 1, 0]]]; (* yellow *) lineseg = Join[lineseg, linelis], {i, 1, ll} ]; lineseg ) ]; (* Calls segfractriang if switch === -1 (full triangulation by line segments), and ltriangsolid otherwise (triangulation by solid colored triangles Only called when switch =!= 1 *) fractriang[switch_,{cnet__},oldm_] := Block[ {cc = {cnet}, res}, ( If[switch === -1, res = segfractriang[cc,oldm], If[switch === 4, res = ltriang3solid[cc,oldm], If[switch > 10, res = ltriangcolor[cc,oldm,switch], res = ltriangsolid[switch,cc,oldm] ] ] ]; res ) ]; (* Builds a list of line segments for the union of lists of nets *) (* This version calls triangul *) (* Also adds the three corners of *) (* each pyramidal net as list elements *) (* uses different colors on three subpatches *) (* this version is for fractal effect *) segfractriang[{cnet__},oldm_] := Block[ {cc = {cnet}, ll, l, i, j, lineseg = {}, edglist, linelis, m, edge, pt}, ( ll = Length[cc]; l = ll/3; Do[ Do[ m = netsize[cc[[j]]]; edgelist= projnet[cc[[j]]]; linelis = triangul[edgelist,m]; edge = Last[linelis]; pt = First[edge]; linelis = Append[linelis,pt]; edge = First[linelis]; pt = First[edge]; linelis = Append[linelis,pt]; edge = linelis[[m]]; pt = Last[edge]; linelis = Append[linelis,pt]; linelis = buildmesh[linelis]; If[j === 1, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[j === 2, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) If[j === 3, linelis = Prepend[linelis, RGBColor[0, 0.7, 0]]]; (* green *) lineseg = Join[lineseg, linelis], {j, 1, 3} ]; cc = Drop[cc, 3], {i, 1, l} ]; lineseg ) ]; (* Builds a list of solid triangles for the union of lists of nets *) (* This version for fractal effect calls soltriangul *) (* Also adds the three corners of *) (* each pyramidal net as list elements *) (* uses different colors on three subpatches *) ltriang3solid[{cnet__},oldm_] := Block[ {cc = {cnet}, ll, l, i, j, lineseg = {}, edglist, linelis, m, edge, pt}, ( ll = Length[cc]; l = ll/3; Do[ Do[ m = netsize[cc[[j]]]; edgelist= projnet[cc[[j]]]; linelis = soltriangul[edgelist,m]; edge = Last[linelis]; pt = edge[[2]]; (* top corner *) linelis = Append[linelis,pt]; edge = First[linelis]; pt = First[edge]; (* lower left corner *) linelis = Append[linelis,pt]; edge = linelis[[m]]; pt = Last[edge]; linelis = Append[linelis,pt]; linelis = solbuildmesh[linelis]; If[j === 1, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[j === 2, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) If[j === 3, linelis = Prepend[linelis, RGBColor[0, 1, 0]]]; (* green *) lineseg = Join[lineseg, linelis], {j, 1, 3} ]; cc = Drop[cc, 3], {i, 1, l} ]; lineseg ) ]; (* builds a list of line segments for the union of lists of nets *) (* This version calls ftriangul *) (* Also adds the three corners of *) (* each pyramidal net as list elements *) lftriang[{cnet__},m_] := Block[ {cc = {cnet}, ll, i, lineseg = {}, edglist, linelis, edge, pt}, ( ll = Length[cc]; Do[ edgelist = cc[[i]]; linelis = ftriangul[edgelist,m]; edge = Last[linelis]; pt = First[edge]; linelis = Append[linelis,pt]; edge = First[linelis]; pt = First[edge]; linelis = Append[linelis,pt]; edge = linelis[[m]]; pt = Last[edge]; linelis = Append[linelis,pt]; lineseg = Join[lineseg, linelis], {i, 1, ll} ]; lineseg = buildmesh[lineseg]; lineseg ) ]; (* Converts a triangular net into a triangular matrix *) convtomat[{cnet__},m_] := Block[ {cc = {cnet}, i, j, n, mat = {}, row}, ( Do[n = m + 2 - i; row = {}; Do[ pt = cc[[j]]; row = Append[row, pt], {j, 1, n} ]; cc = Drop[cc, n]; mat = Append[mat, row], {i, 1, m + 1} ]; mat ) ]; (* To transpose a triangular net wrt to left edge *) (* Converts (r, s, t) to (s, r, t) *) transnetj[{cnet__},m_] := Block[ {cc = {cnet}, i, j, n, aux1, aux2, row, pt, res}, ( aux1 = {}; res = {}; Do[ pt = cc[[j]]; aux1 = Append[aux1, {pt}], {j, 1, m + 1} ]; cc = Drop[cc, m + 1]; Do[n = m + 1 - i; Do[ pt = cc[[j]]; aux2 = Append[aux1[[j]], pt]; aux1 = ReplacePart[aux1, aux2, j], {j, 1, n} ]; cc = Drop[cc, n], {i, 1, m} ]; Do[ row = aux1[[j]]; res = Join[res, row], {j, 1, m + 1} ]; res ) ]; (* To rotate a triangular net *) (* Converts (r, s, t) to (s, t, r) *) transnetk[{cnet__},m_] := Block[ {cc = {cnet}, i, j, n, aux1, aux2, row, pt, res}, ( aux1 = {}; res = {}; Do[ pt = cc[[j]]; aux1 = Append[aux1, {pt}], {j, 1, m + 1} ]; cc = Drop[cc, m + 1]; Do[n = m + 1 - i; Do[ pt = cc[[j]]; aux2 = Prepend[aux1[[j]], pt]; aux1 = ReplacePart[aux1, aux2, j], {j, 1, n} ]; cc = Drop[cc, n], {i, 1, m} ]; Do[ row = aux1[[j]]; res = Join[res, row], {j, 1, m + 1} ]; res ) ]; (* de Casteljau for triangular patches, with subdivision into 3 subpatches *) (* Returns the nets art, ast, ars *) sdecas3[{net__}, mm_, {a__}] := Block[ {cc = {net}, dd, row, row0, row1, row2, barcor = {a}, net0, cctr = {}, ccts = {}, ccsr = {}, m, i, j, k, l, pt, res, art, ast, ars}, (m = mm; net0 = convtomat[cc,m]; (* Print[" netsize m in sdecas3 = ", m]; *) row0 = {}; row1 = {}; row2 = {}; Do[ row0 = Append[row0, net0[[i + 1, 1]]], {i, 0, m} ]; Do[ row1 = Append[row1, net0[[1, j + 1]]], {j, 0, m} ]; Do[ row2 = Append[row2, net0[[i + 1, m - i + 1]]], {i, 0, m} ]; cctr = row0; ccts = row1; ccsr = row2; Do[ dd = {}; Do[ row = {}; Do[ pt = barcor[[1]] * net0[[i + 2, j + 1]] + barcor[[2]] * net0[[i + 1, j + 2]] + barcor[[3]] * net0[[i + 1, j + 1]]; row = Append[row, pt], {j, 0, m - i - l} ]; dd = Join[dd, row], {i, 0, m - l} ]; If[m - l =!= 0, net0 = convtomat[dd,m - l]; row0 = {}; row1 = {}; row2 = {}; Do[ row0 = Append[row0, net0[[i + 1, 1]]], {i, 0, m - l} ]; Do[ row1 = Append[row1, net0[[1, j + 1]]], {j, 0, m - l} ]; Do[ row2 = Append[row2, net0[[i + 1, m - l - i + 1]]], {i, 0, m - l} ]; cctr = Join[cctr, row0]; ccts = Join[ccts, row1]; ccsr = Join[ccsr, row2], cctr = Join[cctr, dd]; ccts = Join[ccts, dd]; ccsr = Join[ccsr, dd] ], {l, 1, m} ]; ast = ccts; (* rat = transnetj[cctr, m]; *) (* rsa = transnetk[ccsr, m]; *) art = cctr; ars = ccsr; res = Join[{art},{ast}]; res = Join[res, {ars}]; res ) ]; (* de Casteljau for triangular patches, with subdivision into 3 subpatches *) (* Returns the layers of the pyramidal set of polar values *) decas3a[{net__}, mm_, {a__}] := Block[ {cc = {net}, dd, row, barcor = {a}, net0, m, i, j, k, l, pt, res}, (m = mm; res = {cc}; net0 = convtomat[cc,m]; Print[" netsize m in decas3a = ", m]; Do[ dd = {}; Do[ row = {}; Do[ pt = barcor[[1]] * net0[[i + 2, j + 1]] + barcor[[2]] * net0[[i + 1, j + 2]] + barcor[[3]] * net0[[i + 1, j + 1]]; row = Append[row, pt], {j, 0, m - i - l} ]; dd = Join[dd, row], {i, 0, m - l} ]; If[m - l =!= 0, net0 = convtomat[dd,m - l] ]; res = Join[res, {dd}], {l, 1, m} ]; res ) ]; (* de Casteljau for triangular patches, with subdivision into 3 subpatches *) (* this version returns the control net for the ref triangle (a, s, t) *) subdecas3ra[{net__}, mm_, {a__}] := Block[ {cc = {net}, dd, row, row0, row1, row2, barcor = {a}, net0, cctr = {}, ccts = {}, ccsr = {}, m, i, j, k, l, pt}, (m = mm; net0 = convtomat[cc,m]; row0 = {}; row1 = {}; row2 = {}; Do[ row1 = Append[row1, net0[[1, j + 1]]], {j, 0, m} ]; ccts = row1; Do[ dd = {}; Do[ row = {}; Do[ pt = barcor[[1]] * net0[[i + 2, j + 1]] + barcor[[2]] * net0[[i + 1, j + 2]] + barcor[[3]] * net0[[i + 1, j + 1]]; row = Append[row, pt], {j, 0, m - i - l} ]; dd = Join[dd, row], {i, 0, m - l} ]; If[m - l =!= 0, net0 = convtomat[dd,m - l]; row0 = {}; row1 = {}; row2 = {}; Do[ row1 = Append[row1, net0[[1, j + 1]]], {j, 0, m - l} ]; ccts = Join[ccts, row1], ccts = Join[ccts, dd] ], {l, 1, m} ]; ccts ) ]; (* de Casteljau for triangular patches, with subdivision into 3 subpatches *) (* this version returns the control net for the ref triangle (a, r, t) *) subdecas3sa[{net__}, mm_, {a__}] := Block[ {cc = {net}, dd, row, row0, row1, row2, barcor = {a}, net0, cctr = {}, ccts = {}, ccsr = {}, m, i, j, k, l, pt}, (m = mm; net0 = convtomat[cc,m]; row0 = {}; row1 = {}; row2 = {}; Do[ row0 = Append[row0, net0[[i + 1, 1]]], {i, 0, m} ]; cctr = row0; Do[ dd = {}; Do[ row = {}; Do[ pt = barcor[[1]] * net0[[i + 2, j + 1]] + barcor[[2]] * net0[[i + 1, j + 2]] + barcor[[3]] * net0[[i + 1, j + 1]]; row = Append[row, pt], {j, 0, m - i - l} ]; dd = Join[dd, row], {i, 0, m - l} ]; If[m - l =!= 0, net0 = convtomat[dd,m - l]; row0 = {}; row1 = {}; row2 = {}; Do[ row0 = Append[row0, net0[[i + 1, 1]]], {i, 0, m - l} ]; cctr = Join[cctr, row0], cctr = Join[cctr, dd] ], {l, 1, m} ]; cctr ) ]; (* de Casteljau for triangular patches, with subdivision into 3 subpatches *) (* this version returns the control net for the ref triangle (a, r, s) *) subdecas3ta[{net__}, mm_, {a__}] := Block[ {cc = {net}, dd, row, row0, row1, row2, barcor = {a}, net0, cctr = {}, ccts = {}, ccsr = {}, m, i, j, k, l, pt}, (m = mm; net0 = convtomat[cc,m]; row0 = {}; row1 = {}; row2 = {}; Do[ row2 = Append[row2, net0[[i + 1, m - i + 1]]], {i, 0, m} ]; ccsr = row2; Do[ dd = {}; Do[ row = {}; Do[ pt = barcor[[1]] * net0[[i + 2, j + 1]] + barcor[[2]] * net0[[i + 1, j + 2]] + barcor[[3]] * net0[[i + 1, j + 1]]; row = Append[row, pt], {j, 0, m - i - l} ]; dd = Join[dd, row], {i, 0, m - l} ]; If[m - l =!= 0, net0 = convtomat[dd,m - l]; row0 = {}; row1 = {}; row2 = {}; Do[ row2 = Append[row2, net0[[i + 1, m - l - i + 1]]], {i, 0, m - l} ]; ccsr = Join[ccsr, row2], ccsr = Join[ccsr, dd] ], {l, 1, m} ]; ccsr ) ]; (* polar value by de Casteljau for triangular patches *) poldecas3[{net__}, mm_, {a__}] := Block[ {cc = {net}, dd, barcor = {a}, net0, row, m, i, j, k, l, pt, res}, (m = mm; net0 = convtomat[cc,m]; Do[ dd = {}; Do[ row = {}; Do[ pt = barcor[[l,1]] * net0[[i + 2, j + 1]] + barcor[[l,2]] * net0[[i + 1, j + 2]] + barcor[[l,3]] * net0[[i + 1, j + 1]]; row = Append[row, pt], {j, 0, m - i - l} ]; dd = Join[dd, row], {i, 0, m - l} ]; If[m - l =!= 0, net0 = convtomat[dd,m - l], res = dd ], {l, 1, m} ]; res ) ]; (* Computes a 3 X 3 determinant *) (* Assuming a, b, c are column vectors *) ddet3[a_, b_, c_] := Block[ {a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, res}, ( a1 = a[[1]]; a2 = a[[2]]; a3 = a[[3]]; b1 = b[[1]]; b2 = b[[2]]; b3 = b[[3]]; c1 = c[[1]]; c2 = c[[2]]; c3 = c[[3]]; d1 = b2*c3 - b3*c2; d2 = b1*c3 - b3*c1; d3 = b1*c2 - b2*c1; res = a1*d1 - a2*d2 + a3*d3; res ) ]; (* Computes a 3 X 3 determinant *) (* Assuming a, b, c are row vectors *) det3[a_, b_, c_] := Block[ {a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, res}, ( a1 = a[[1]]; a2 = a[[2]]; a3 = a[[3]]; b1 = b[[1]]; b2 = b[[2]]; b3 = b[[3]]; c1 = c[[1]]; c2 = c[[2]]; c3 = c[[3]]; d1 = b1*c2 - b2*c1; d2 = a1*c2 - a2*c1; d3 = a1*b2 - a2*b1; res = a3*d1 - b3*d2 + c3*d3; res ) ]; (* Solves a 3 X 3 linear system *) (* Assuming a, b, c, d are column vectors *) solve3[a_, b_, c_, d_] := Block[ {a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, x1, x2, x3, dd, res}, (dd = det3[a, b, c]; If[dd === 0. || dd === 0, Print["*** Null Determinant ***"]; dd = 10^(-10)]; x1 = det3[d, b, c]/dd; x2 = det3[a, d, c]/dd; x3 = det3[a, b, d]/dd; res = {x1, x2, x3}; res ) ]; (* computes new barycentric coordinates (u1, u2, u3) of a point b *) (* with old barycentric coordinates (l2, m2, n2) w.r.t *) (* ref triangle (r, s, t), w.r.t new ref triangle (a, s, t) *) (* where a has coords (l1, m1, n1) w.r.t ref triangle (r, s, t *) (* WARNING: l1 must be NONZERO *) newbaryb[{a__}, {b__}] := Block[ {aa = {a}, bb = {b}, u1, u2, u3, l1, m1, n1, l2, m2, n2, res}, ( l1 = aa[[1]]; m1 = aa[[2]]; n1 = aa[[3]]; l2 = bb[[1]]; m2 = bb[[2]]; n2 = bb[[3]]; If[l1 === 0. || l1 === 0, Print["*** Error, l1 = 0"]; Return["*** Error, l1 = 0"]]; u1 = l2/l1; u2 = m2 - m1 * u1; u3 = n2 - n1 * u1; res = {u1, u2, u3}; res ) ]; (* computes new barycentric coordinates (v1, v2, v3) of a point c *) (* with old barycentric coordinates (l3, m3, n3) w.r.t *) (* ref triangle (r, s, t), w.r.t new ref triangle (a, b, t) *) (* where a has coords (l1, m1, n1) w.r.t ref triangle (r, s, t) *) (* and b has coords (l2, m2, n2) w.r.t ref triangle (r, s, t) *) (* WARNING: a and b cannot be collinear with t *) newbaryc[{a__}, {b__}, {c__}] := Block[ {aa = {a}, bb = {b}, cc = {c}, v1, v2, v3, l1, m1, n1, l2, m2, n2, l3, m3, n3, dd, res}, ( l1 = aa[[1]]; m1 = aa[[2]]; n1 = aa[[3]]; l2 = bb[[1]]; m2 = bb[[2]]; n2 = bb[[3]]; l3 = cc[[1]]; m3 = cc[[2]]; n3 = cc[[3]]; dd = l1 * m2 - l2 * m1; If[dd === 0. || dd === 0, Print["*** Error, dd = 0"]; Return["*** Error, dd = 0"]]; v1 = (l3 * m2 - l2 * m3)/dd; v2 = (l1 * m3 - l3 * m1)/dd; v3 = n3 - n1 * v1 - n2 * v2; res = {v1, v2, v3}; res ) ]; (* computes new control net, given barycentric coords for *) (* new reference triangle (a, b, c) w.r.t (r, s, t) *) (* The algorithm first computes the control net w.r.t *) (* (a, s, t) using subdecas3ra, then the control net w.r.t *) (* (a, b, t) using subdecas3sa and transnetj, and *) (* the control net w.r.t (a, b, c) using subdecas3ta and *) (* transnetk *) (* this version calls newcnet3old when l1 = 0 or d = 0 *) znewcnet3[{net__}, m_, {reftrig__}] := Block[ {cc = {net}, newtrig = {reftrig}, neta, netb, newnet, a, b, c, nb, nc, rr, ss, tt, args, i, j, k, l, res, l1, m1, l2, m2, d}, (newnet = {}; a = newtrig[[1]]; b = newtrig[[2]]; c = newtrig[[3]]; l1 = a[[1]]; m1 = a[[2]]; l2 = b[[1]]; m2 = b[[2]]; d = l1 * m2 - l2 * m1; (* Print[" d in newcnet3: ", d]; *) If[d === 0. || d === 0, Print["** In newcnet3, d = 0 **"]]; If[l1 === 0. || l1 === 0, Print["** In newcnet3, l1 = 0 **"]]; If[l1 === 0. || d === 0.|| l1 === 0 || d === 0, newnet = newcnet3old[cc, m, newtrig], neta = subdecas3ra[cc, m, a]; nb = newbaryb[a, b]; netb = subdecas3sa[neta, m, nb]; netb = transnetj[netb, m]; nc = newbaryc[a, b, c]; newnet = subdecas3ta[netb, m, nc]; newnet = transnetk[newnet, m] ]; newnet ) ]; (* computes new control net, given barycentric coords for *) (* new reference triangle (a, b, c) w.r.t (r, s, t) *) (* In the simplest case where a, s, t are not collinear and *) (* b is not on (a, t), *) (* the algorithm first computes the control net w.r.t *) (* (a, s, t) using subdecas3ra, then the control net w.r.t *) (* (a, b, t) using subdecas3sa and transnetj, and *) (* the control net w.r.t (a, b, c) using subdecas3ta and *) (* transnetk *) (* The other cases are also treated *) (* this version is faster than the previous one *) (* neta = (a, s, t), or (a, r, t), or (a, r, s) *) (* The function returns the net (a, b, c) *) newcnet3[{net__}, m_, {reftrig__}] := Block[ {cc = {net}, newtrig = {reftrig}, neta, netb, newnet, a, b, c, nb, nc, rr, ss, tt}, (newnet = {}; a = newtrig[[1]]; b = newtrig[[2]]; c = newtrig[[3]]; rr = {1, 0, 0}; ss = {0, 1, 0}; tt = {0, 0, 1}; If[collin[a, ss, tt] === 0, (* In this case, a is not on (s, t). Want (a, s, t) *) (* Print[" a NOT on (s, t) "]; *) neta = subdecas3ra[cc, m, a]; (* Now, the net is (a, s, t) *) newnet = fnewaux[neta, m, a, b, c, ss, tt], (* In this case, a is on (s, t) *) (* Print[" a IS on (s, t) "]; *) If[a === ss, (* In this case, a is on (s, t) and a = s. Want (a, r, t) *) (* Print[" a on (s, t) and a = ss"]; *) neta = transnetj[cc, m]; (* Now, the net is (a, r, t) *) newnet = fnewaux[neta, m, a, b, c, rr, tt], (* In this case, a is on (s, t) and a <> s. Want (a, r, s) *) (* Print[" a on (s, t) and a <> ss"]; *) neta = subdecas3ta[cc, m, a]; (* Now, the net is (a, r, s) *) newnet = fnewaux[neta, m, a, b, c, rr, ss] ] ]; newnet ) ]; (* This routine is used by newcnet3 *) (* Initially, neta = (a,s,t) *) fnewaux[{net__}, m_, a_, b_, c_, ss_, tt_] := Block[ {neta = {net}, netb, newnet, nb, nc}, ( If[collin[b, a, tt] === 0, (* In this case, b is not on (a, t). Want (a, b, t) *) (* Print[" b NOT on (a, t) "]; *) nb = solve3[a, ss, tt, b]; netb = subdecas3sa[neta, m, nb]; netb = transnetj[netb, m]; (* Now, the net is (a, b, t) *) nc = solve3[a, b, tt, c]; newnet = subdecas3ta[netb, m, nc]; newnet = transnetk[newnet, m], (* In this case, b is on (a, t). Want (a, b, s) *) (* Print[" b IS on (a, t) "]; *) neta = transnetk[neta, m]; neta = transnetk[neta, m]; (* Now, the net (a, s, t) is (t, a, s) *) nb = solve3[tt, a, ss, b]; netb = subdecas3ra[neta, m, nb]; netb = transnetj[netb, m]; (* Now, the net is (a, b, s) *) nc = solve3[a, b, ss, c]; newnet = subdecas3ta[netb, m, nc]; newnet = transnetk[newnet, m] ]; newnet ) ]; (* computes new control net, given barycentric coords for *) (* new reference triangle (a, b, c) *) (* slower version *) newcnet3old[{net__}, m_, {reftrig__}] := Block[ {cc = {net}, newtrig = {reftrig}, a, newnet, rr, ss, tt, args, i, j, k, l}, (newnet = {}; rr = newtrig[[1]]; ss = newtrig[[2]]; tt = newtrig[[3]]; Print[" In newcnet3old!"]; Do[ Do[args = {}; Do[ args = Append[args, rr], {l, 1, i} ]; Do[ args = Append[args, ss], {l, 1, j} ]; Do[ args = Append[args, tt], {l, 1, m - i - j} ]; a = poldecas3[cc,m,args]; newnet = Join[newnet,a], {j, 0, m - i} ], {i, 0, m} ]; newnet ) ]; (* a = (t + s)/2, b = (t + r)/2, c = (r + s)/2 *) (* Subdivides into four subpatches as follows: first compute the *) (* point F((t + s)/2) and the three subpatches art, ast, ars, *) (* discard the subpatch ast (2nd subpatch) then compute *) (* the point F((t + r)/2) on the subpatch art, and the three *) (* subpatches bat, brt, bar, discard the subpatch brt *) (* (2nd subpatch), and finally compute the point F((r + s)/2) on *) (* the subpatch ars and the three subpatches cas, crs, car, *) (* and discard the subpatch crs (second subpatch). *) (* The join of bat, bar, cas, car is returned *) (* diamond-style subdivision *) splitdecas4[{net__}, mm_] := Block[ {cc = {net}, newnet = {}, barcor1, barcor2, anet, art, ars, bat, bar, cas, car, net1, net2, m = mm}, (barcor1 = {1/2, 0, 1/2}; barcor2 = {0, 1/2, 1/2}; Print[" netsize m in splitdecas4 = ", m]; anet = sdecas3[cc, m, barcor2]; art = anet[[1]]; ars = anet[[3]]; anet = sdecas3[art, m, barcor2]; bat = anet[[1]]; bar = anet[[3]]; net1 = Join[{bat}, {bar}]; anet = sdecas3[ars, m, barcor2]; cas = anet[[1]]; car = anet[[3]]; net2 = Join[{cas}, {car}]; newnet = Join[net1,net2]; newnet ) ]; (* computes additional control points to have a nicer triangulation *) (* in the case of splitting into 4 patches using splitdecas4 *) split24[{net__}, mm_] := Block[ {cc = {net}, dd, net1, net2, net3, net4, netabc, netbrc, triangle1, triangle2, m, i, l, ll, surf}, (m = mm; triangle1 = {{1, 0, 0}, {0, 1, 0}, {1, 1, -1}}; triangle2 = {{1, 0, 0}, {2, 0, -1}, {1, 1, -1}}; l = Length[cc]; ll = l/4; dd = {}; Do[ net1 = cc[[1]]; net2 = cc[[2]]; net3 = cc[[3]]; net4 = cc[[4]]; netabc = newcnet3[net1, m, triangle1]; netbrc = newcnet3[net1, m, triangle2]; dd = Join[dd, {net1}]; dd = Join[dd, {netbrc}]; dd = Join[dd, {netabc}]; dd = Join[dd, {net3}]; cc = Drop[cc, 4], {i, 1, ll} ]; dd ) ]; (* g = (t + r + s)/3, a = (t + s)/2, b = (t + r)/2, c = (r + s)/2 *) (* Subdivides into six subpatches: first compute the *) (* point F((t + r + s)/3) and the three subpatches *) (* grt, gst, grs. Then compute the *) (* point F((t + r)/2) and the three subpatches bgt, brt, bgr, *) (* discard the subpatch brt (2nd subpatch) then compute *) (* the point F((t + s)/2) on the subpatch gst, and the three *) (* subpatches agt, ast, ags, discard the subpatch ast *) (* (2nd subpatch), and finally compute the point F((s + r)/2) on *) (* the subpatch grs and the three subpatches cgs, crs, cgr, *) (* and discard the subpatch crs (2nd subpatch). *) (* Returns the list bgt, bgr, agt, ags, cgs, cgr *) (* spider-web style subdivision *) splitdecas6[{net__}, mm_] := Block[ {cc = {net}, newnet = {}, barcor1, barcor2, g, anet, grt, gst, grs, bgt, bgr, agt, ags, cgs, cgr, n1, n2, n3, m = mm}, (barcor1 = {1/2, 0, 1/2}; barcor2 = {0, 1/2, 1/2}; g = {1/3, 1/3, 1/3}; Print[" netsize m in splitdecas6 = ", m]; anet = sdecas3[cc, m, g]; grt = anet[[1]]; gst = anet[[2]]; grs = anet[[3]]; anet = sdecas3[grt, m, barcor2]; bgt = anet[[1]]; bgr = anet[[3]]; n1 = Join[{bgt}, {bgr}]; anet = sdecas3[gst, m, barcor2]; agt = anet[[1]]; ags = anet[[3]]; n2 = Join[{agt}, {ags}]; anet = sdecas3[grs, m, barcor2]; cgs = anet[[1]]; cgr = anet[[3]]; n3 = Join[{cgs}, {cgr}]; newnet = Join[n1,n2]; newnet = Join[newnet,n3]; newnet ) ]; (* experimenting with subdivision *) (* this version seems correct, and yields a strange triangulation! *) splitd4[{net__}, mm_] := Block[ {cc = {net}, newnet = {}, barcor1, barcor2, barcor3, anet, ccts1, ccsr1, lnet2, lnet3, m = mm, res}, (barcor1 = {1/2, 0, 1/2}; barcor2 = {0, 1/2, 1/2}; barcor3 = {1/2, 1/2, 0}; newnet = sdecas3[cc, m, barcor2]; lnet2 = newnet[[1]]; lnet3 = newnet[[3]]; res = Join[{lnet2},{lnet3}]; res ) ]; (* Subdivides into four subpatches, by dividing the base triangle *) (* into four subtriangles using the middles of the original sides *) subdecas12[{net__}, mm_] := Block[ {cc = {net}, newnet = {}, barcor1, barcor2, barcor3, barcor4, anet, m = mm}, (barcor1 = {{1/2, 0, 1/2}, {0, 1/2, 1/2}, {0, 0, 1}}; barcor2 = {{1/2, 0, 1/2}, {1/2, 1/2, 0}, {1, 0, 0}}; barcor3 = {{1/2, 0, 1/2}, {0, 1/2, 1/2}, {1/2, 1/2, 0}}; barcor4 = {{1/2, 1/2, 0}, {0, 1/2, 1/2}, {0, 1, 0}}; anet = newcnet3[cc, m, barcor1]; newnet = Append[newnet, anet]; anet = newcnet3[cc, m, barcor2]; newnet = Append[newnet, anet]; anet = newcnet3[cc, m, barcor3]; newnet = Append[newnet, anet]; anet = newcnet3[cc, m, barcor4]; newnet = Append[newnet, anet]; newnet ) ]; (* Subdivides into four subpatches, by dividing the base triangle *) (* into four subtriangles using the middles of the original sides *) (* uses a tricky scheme involving 5 de Casteljau steps *) (* The triangles are abt, sca, bac, and crb *) nsubdecas45[{net__}, oldm_, debug_] := Block[ {art = {net}, newnet = {}, barcor1, barcor2, barcor3, bat, cba, rcb, sca, abt, bac, crb, m = oldm}, (barcor1 = {1, 1, -1}; barcor2 = {0, 1/2, 1/2}; barcor3 = {1, -1, 1}; bat = subdecas3sa[art, m, barcor2]; cba = subdecas3ta[bat, m, barcor1]; rcb = subdecas3ta[cba, m, barcor1]; sca = subdecas3sa[cba, m, barcor3]; abt = transnetj[bat, m]; bac = transnetk[cba, m]; crb = transnetj[rcb, m]; newnet = Join[{abt}, {sca}]; newnet = Join[newnet, {bac}]; newnet = Join[newnet, {crb}]; newnet ) ]; (* Subdivides into four subpatches, by dividing the base triangle *) (* into four subtriangles using the middles of the original sides *) (* uses a tricky scheme involving 5 de Casteljau steps *) (* basic version that does not try to resolve singularities *) mainsubdecas45[{net__}, oldm_, debug_] := Block[ {cc = {net}, newnet, m, art, a}, (a = {0, 1/2, 1/2}; m = netsize[cc]; art = subdecas3sa[cc, m, a]; If[debug === -20, Print["*** Left corner patch art: ", art]]; newnet = nsubdecas4[art, m, debug]; newnet ) ]; (* Subdivides into four subpatches, by dividing the base triangle *) (* into four subtriangles using the middles of the original sides *) (* uses a tricky scheme involving 4 de Casteljau steps *) (* basic version that does not try to resolve singularities *) (* The triangles are abt, sca, bac, and crb *) mainsubdecas4[{net__}, oldm_, debug_] := Block[ {cc = {net}, newnet, m, art, ars, barcor2}, (m = netsize[cc]; (* Print[" calling mainsubdecas4 with netsize = ", m]; *) barcor1 = {-1, 1, 1}; barcor2 = {0, 1/2, 1/2}; newnet = sdecas3[cc, m, barcor2]; art = newnet[[1]]; ars = newnet[[3]]; newnet = nsubdecas4[art, ars, m, debug]; newnet ) ]; (* Subdivides into four subpatches, by dividing the base triangle *) (* into four subtriangles using the middles of the original sides *) (* uses a tricky scheme involving 4 de Casteljau steps *) (* basic version that does not try to resolve singularities *) (* The triangles are abt, sca, bac, and crb *) nsubdecas4[{net1__}, {net2__}, oldm_, debug_] := Block[ {newnet, m, art = {net1}, ars = {net2}, bat, bar, cas, sca, cbr, cba, bac, crb, barcor1, barcor2}, (m = netsize[art]; (* Print[" calling mainsubdecas4 with netsize = ", m]; *) barcor1 = {-1, 1, 1}; barcor2 = {0, 1/2, 1/2}; If[debug === -20, Print["*** art: ",art]]; If[debug === -20, Print["*** ars: ",ars]]; newnet = sdecas3[art, m, barcor2]; bat = newnet[[1]]; bar = newnet[[3]]; abt = transnetj[bat, m]; If[debug === -20, Print["*** abt: ",abt]]; If[debug === -20, Print["*** bar: ",bar]]; cas = subdecas3sa[ars, m, barcor2]; sca = transnetk[cas, m]; sca = transnetk[sca, m]; newnet = sdecas3[bar, m, barcor1]; cbr = newnet[[1]]; cba = newnet[[3]]; bac = transnetk[cba, m]; crb = transnetj[cbr, m]; crb = transnetk[crb, m]; newnet = Join[{abt}, {sca}]; newnet = Join[newnet, {bac}]; newnet = Join[newnet, {crb}]; newnet ) ]; (* To test that a rectangular net really has size (p + 1)(q + 1) *) (* when the polar bidegree is (p, q) *) testpq[cnet__, p_, q_] := Block[ {cc = cnet, ll, stop}, (ll = Length[cc]; If[ll === (p + 1)*(q + 1), stop = 0, stop = 1]; If[stop === 1, Print["*** Net Size: ", ll, " Inconsistent with Surface bidegree: (", p, ", ",q, ") ***"]]; stop ) ]; (* computes a rectangular net of degree ((m - nm), m) from a triangular net of degree m *) (* by "blowing up" a degenerate corner. Simple version wrt standard frame (r, s, t) *) nsqblowup[{oldnet__}, mm_, nm_, debug_] := Block[ {net = {oldnet}, newnet, m = mm, rr, ss, temp, izargs, jzargs, pt, u, v, w, i, j, ii, permset, perm, setm, pi, kk, ll, xx, lowi}, (rr = 0; ss = 1; lowi = nm; Print["In sqblowup, m = ", m, ", nm = ", nm]; newnet = {}; setm = Table[ii, {ii, 1, m}]; If[debug === -20, Print["*** Left corner patch in blowup: ", net]]; Do[izargs = {}; Do[ izargs = Append[izargs, rr], {ii, 1, m - i} ]; Do[ izargs = Append[izargs, ss], {ii, m + 1 - i, m} ]; Do[ jzargs = {}; xx = 0; Do[ jzargs = Append[jzargs, rr], {ii, 1, m - j} ]; Do[ jzargs = Append[jzargs, ss], {ii, m + 1 - j, m} ]; If[debug === 20, Print[" izargs = ", izargs]]; If[debug === 20, Print[" jzargs = ", jzargs]]; permset = Permutations[setm]; ll = Length[permset]; Do[ perm = permset[[kk]]; args = {}; If[debug === 20, Print[" perm = ", perm]]; Do[ pi = perm[[ii]]; u = izargs[[ii]]*(1 - jzargs[[pi]]); v = izargs[[ii]]*jzargs[[pi]]; w = 1 - izargs[[ii]]; temp = {u, v, w}; args = Append[args, temp], {ii, 1, m} ]; If[debug === 20, Print[" args = ", args]]; pt = poldecas3[net, m, args]; xx = xx + pt, {kk, 1, ll} ]; pt = xx/ll; If[debug === 20, Print[" xx/ll: ", pt]]; newnet = Join[newnet, pt]; Print["In sqblowup, i = ", i, ", j = ", j], {j, 0, m} ], {i, lowi, m} ]; If[debug === 20, Print[" newnet = ", newnet]]; newnet ) ]; (* computes a rectangular net of degree ((m - nm), m) from a triangular net of degree m *) (* by "blowing up" a degenerate corner *) (* General version, blowing up a triangular patch based on the triangle *) (* ((r_1, r_2), (s_1, s_2), (t_1, t_2)), wrt the standard frame (r, s, t) *) sqblowup[{oldnet__}, {reftrig__}, mm_, nm_, debug_] := Block[ {net = {oldnet}, newnet, m = mm, rr, ss, ntrig = {reftrig}, temp, izargs, jzargs, pt, u, v, w, i, j, ii, permset, perm, setm, pi, kk, ll, xx, lowi, r1, r2, s1, s2, t1, t2}, (rr = 0; ss = 1; lowi = nm; Print["In sqblowup, m = ", m, ", nm = ", nm]; r1 = ntrig[[1, 1]]; r2 = ntrig[[1, 2]]; s1 = ntrig[[2, 1]]; s2 = ntrig[[2, 2]]; t1 = ntrig[[3, 1]]; t2 = ntrig[[3, 2]]; Print[" r1 = ", r1, ", r2 = ", r2, ", s1 = ", s1, ", s2 = ", s2, ", t1 = ", t1, ", t2 = ", t2]; newnet = {}; setm = Table[ii, {ii, 1, m}]; If[debug === -20, Print["*** Left corner patch in blowup: ", net]]; Do[izargs = {}; Do[ izargs = Append[izargs, rr], {ii, 1, m - i} ]; Do[ izargs = Append[izargs, ss], {ii, m + 1 - i, m} ]; Do[ jzargs = {}; xx = 0; Do[ jzargs = Append[jzargs, rr], {ii, 1, m - j} ]; Do[ jzargs = Append[jzargs, ss], {ii, m + 1 - j, m} ]; If[debug === 20, Print[" izargs = ", izargs]]; If[debug === 20, Print[" jzargs = ", jzargs]]; permset = Permutations[setm]; ll = Length[permset]; Do[ perm = permset[[kk]]; args = {}; If[debug === 20, Print[" perm = ", perm]]; Do[ pi = perm[[ii]]; u = (s1 - r1)*izargs[[ii]]*jzargs[[pi]] + (r1 - t1)*izargs[[ii]] + t1; v = (s2 - r2)*izargs[[ii]]*jzargs[[pi]] + (r2 - t2)*izargs[[ii]] + t2; w = 1 - u - v; temp = {u, v, w}; args = Append[args, temp], {ii, 1, m} ]; If[debug === 20, Print[" args = ", args]]; pt = poldecas3[net, m, args]; xx = xx + pt, {kk, 1, ll} ]; pt = xx/ll; If[debug === 20, Print[" xx/ll: ", pt]]; newnet = Join[newnet, pt]; Print["In sqblowup, i = ", i, ", j = ", j], {j, 0, m} ], {i, lowi, m} ]; If[debug === 20, Print[" newnet = ", newnet]]; newnet ) ]; (* converts a triangular net of degree m to a rectangular net of degree (m, m) *) (* wrt to the standard Cartesian frame ((0, 0), ((1, 0), (0, 1))) *) trtosq[{oldnet__}, mm_, debug_] := Block[ {net = {oldnet}, newnet, m = mm, rr, ss, temp, izargs, jzargs, pt, u, v, w, i, j, ii, permset, perm, setm, pi, kk, ll, xx}, (rr = 0; ss = 1; newnet = {}; setm = Table[ii, {ii, 1, m}]; If[debug === -20, Print["*** Input net: ", net]]; Do[izargs = {}; Do[ izargs = Append[izargs, rr], {ii, 1, m - i} ]; Do[ izargs = Append[izargs, ss], {ii, m + 1 - i, m} ]; Do[ jzargs = {}; xx = 0; Do[ jzargs = Append[jzargs, rr], {ii, 1, m - j} ]; Do[ jzargs = Append[jzargs, ss], {ii, m + 1 - j, m} ]; If[debug === 20, Print[" izargs = ", izargs]]; If[debug === 20, Print[" jzargs = ", jzargs]]; permset = Permutations[setm]; ll = Length[permset]; Do[ perm = permset[[kk]]; args = {}; If[debug === 20, Print[" perm = ", perm]]; Do[ pi = perm[[ii]]; u = izargs[[ii]]; v = jzargs[[pi]]; w = 1 - u - v; temp = {u, v, w}; args = Append[args, temp], {ii, 1, m} ]; If[debug === 20, Print[" args = ", args]]; pt = poldecas3[net, m, args]; xx = xx + pt, {kk, 1, ll} ]; pt = xx/ll; If[debug === 20, Print[" xx/ll: ", pt]]; newnet = Join[newnet, pt]; Print["In trtosq, i = ", i, " j = ", j], {j, 0, m} ], {i, 0, m} ]; If[debug === 20, Print[" newnet = ", newnet]]; newnet ) ]; (* to convert a triangular net of degree m to a rectangular net of degree (m, m) *) (* This program is used with recrender *) convtrsq[{oldnet__}, m_, debug_] := Block[ {net0 = {oldnet}, net, newnet, res1, res}, net = maptohat[net0]; Print[" net before conversion: ", net]; res1 = trtosq[net, m, debug]; res = piomeg[res1]; Print[" Square net: ", res]; res ]; (* to convert a triangular net of degree m to a rectangular net of *) (* degree ((m - nm), m) *) (* The rectangular net is obtained by blowing up of the triangular patch *) (* This program is used with recrender *) blowtrsq[{oldnet__}, m_, nm_, reftrig_, debug_] := Block[ {net0 = {oldnet}, net, newnet, res1, res}, net = maptohat[net0]; Print[" net before conversion: ", net]; res1 = sqblowup[net, reftrig, m, nm, debug]; res = piomeg[res1]; Print[" Square net: ", res]; res ]; (* to convert a triangular net of degree m to a rectangular net of degree (m, m) *) (* After taking a new reference triangle reftrig *) (* This program is used with recrender *) gconvtrsq[{oldnet__}, m_, reftrig_, debug_] := Block[ {net0 = {oldnet}, net, net2, newnet, res1, res}, net = maptohat[net0]; Print[" Net before conversion: ", net]; net2 = newcnet3[net, m, reftrig]; Print[" net2: ", net2]; res1 = trtosq[net2, m, debug]; res = piomeg[res1]; Print[" Square net: ", res]; res ]; (* converts a rectangular net of degree (p, q) to *) (* a triangular net of degree m = p + q, wrt the standard frame (r, s, t), where r = (1, 0, 0), s = (0, 1, 0), t = (0, 0, 1) *) rectotr[{oldnet__}, p_, q_, debug_] := Block[ {net = {oldnet}, i, j, newnet, Isets, I1, J1, ii, jj, l, m, setm, pt1, pt2, pt, lambda, mu, xx, u, v, rr, ss, tt, r1, s1, r2, s2}, rr = {1, 0}; ss = {0, 1}; tt = {0, 0}; r1 = 0; r2 = 0; s1 = 1; s2 = 1; newnet = {}; m = p + q; (* Print["p = ", p," q = ", q]; *) setm = Table[i, {i, m}]; Do[ Do[lambda = {}; mu = {}; u = {}; v = {}; Do[ lambda = Append[lambda, rr[[1]]]; mu = Append[mu, rr[[2]]], {ii, 1, i} ]; Do[ lambda = Append[lambda, ss[[1]]]; mu = Append[mu, ss[[2]]], {ii, 1, j} ]; Do[ lambda = Append[lambda, tt[[1]]]; mu = Append[mu, tt[[2]]], {ii, 1, m - i - j} ]; Isets = subsets[setm, p]; l = Length[Isets]; pt = 0; Do[ I1 = Isets[[ii]]; J1 = relcomp[setm, I1]; If[debug === 20, Print[" I1 = ", I1]]; If[debug === 20, Print[" J1 = ", J1]]; u = {}; v = {}; Do[ kk = I1[[jj]]; pt1 = lambda[[kk]]; u = Append[u, pt1], {jj, 1, p} ]; If[debug === 20, Print[" u = ", u]]; Do[ kk = J1[[jj]]; pt2 = mu[[kk]]; v = Append[v, pt2], {jj, 1, q} ]; If[debug === 20, Print[" v = ", v]]; xx = recpoldecas[net, p, q, r1, s1, r2, s2, u, v, debug]; pt = pt + xx, {ii, 1, l} ]; xx = pt/Binomial[m,p]; newnet = Join[newnet, {xx}]; Print["In rectotr, i = ", i, ", j = ", j], {j, 0, m - i} ], {i, 0, m} ]; If[debug === 20, Print[" newnet = ", newnet]]; newnet ]; (* to convert a triangular net of degree m to a rectangular net of degree m x m *) (* This program is used with recrender *) convrectr[{oldnet__}, p_, q_, debug_] := Block[ {net0 = {oldnet}, net, net2, newnet, res1, res}, net = maptohat[net0]; Print[" net before conversion: ", net]; res1 = rectotr[net, p, q, debug]; res = piomeg[res1]; Print[" Triangular net: ", res]; res ]; (* converts a rectangular net of degree (p, q) to *) (* a triangular net of degree m = p + q, wrt the frame (r, s, d) *) (* where r = (1, 0, 0), s = (0, 1, 0), d = (1, 1, -1) *) rectotr2[{oldnet__}, p_, q_, debug_] := Block[ {net = {oldnet}, i, j, newnet, Isets, I1, J1, ii, jj, l, m, setm, pt1, pt2, pt, lambda, mu, xx, u, v, rr, ss, tt, r1, s1, r2, s2}, rr = {1, 0}; ss = {0, 1}; tt = {1, 1}; r1 = 0; r2 = 0; s1 = 1; s2 = 1; newnet = {}; m = p + q; (* Print["p = ", p," q = ", q]; *) setm = Table[i, {i, m}]; Do[ Do[lambda = {}; mu = {}; u = {}; v = {}; Do[ lambda = Append[lambda, rr[[1]]]; mu = Append[mu, rr[[2]]], {ii, 1, i} ]; Do[ lambda = Append[lambda, ss[[1]]]; mu = Append[mu, ss[[2]]], {ii, 1, j} ]; Do[ lambda = Append[lambda, tt[[1]]]; mu = Append[mu, tt[[2]]], {ii, 1, m - i - j} ]; Isets = subsets[setm, p]; l = Length[Isets]; pt = 0; Do[ I1 = Isets[[ii]]; J1 = relcomp[setm, I1]; If[debug === 20, Print[" I1 = ", I1]]; If[debug === 20, Print[" J1 = ", J1]]; u = {}; v = {}; Do[ kk = I1[[jj]]; pt1 = lambda[[kk]]; u = Append[u, pt1], {jj, 1, p} ]; If[debug === 20, Print[" u = ", u]]; Do[ kk = J1[[jj]]; pt2 = mu[[kk]]; v = Append[v, pt2], {jj, 1, q} ]; If[debug === 20, Print[" v = ", v]]; xx = recpoldecas[net, p, q, r1, s1, r2, s2, u, v, debug]; pt = pt + xx, {ii, 1, l} ]; xx = pt/Binomial[m,p]; newnet = Join[newnet, {xx}]; Print["In rectotr, i = ", i, ", j = ", j], {j, 0, m - i} ], {i, 0, m} ]; If[debug === 20, Print[" newnet = ", newnet]]; newnet ]; (* Computation of a polar value for a rectangular net of degree (p, q) *) recpoldecas[{oldnet__}, p_, q_, r1_, s1_, r2_, s2_, u_, v_, debug_] := Block[ {net = {oldnet}, i, j, bistar, bstarj, pt, res}, bistar = {}; Do[ bstarj = {}; Do[ pt = net[[(q + 1)*i + j + 1]]; bstarj = Append[bstarj, pt], {j, 0, q} ]; If[debug === 20, Print[" bstarj: ", bstarj]]; pt = bdecas[bstarj, v, r2, s2]; bistar = Append[bistar, pt], {i, 0, p} ]; If[debug === 20, Print[" bistar: ", bistar]]; res = bdecas[bistar, u, r1, s1]; If[debug === 20, Print[" polar value: ", res]]; res ]; (* to test the conversion programs *) tconvert[{oldnet__}, m_, debug_] := Block[ {net0 = {oldnet}, newnet, res1, res2, res}, net = maptohat[net0]; Print[" Net before conversion: ", net]; res1 = trtosq[net, m, debug]; Print[" Square net: ", res1]; res2 = rectotr[res1, m, m, debug]; Print[" Triangular net: ", res2]; res = piomeg[res2]; res ]; (* to test the blowup programs *) testblow2[{oldnet__}, m_, nm_, debug_] := Block[ {net0 = {oldnet}, theta, theta1, theta2, lnet, newnet1a, newnet1b, newnet, res1, res2, res, trig}, net = maptohat[net0]; Print[" Net before blowing up: ", net]; trig = {{1, 0}, {0, 1}, {0, 0}}; res1 = sqblowup[net, trig, m, nm, debug]; Print[" Square net after blowing up: ", res1]; res2 = rectotr[res1, m - nm, m, debug]; res = piomeg[res2]; res ]; (* to test the blowup programs with theta1, etc, triangular net *) testblow[{oldnet__}, m_, nm_, debug_] := Block[ {net0 = {oldnet}, theta, theta1, theta2, lnet, newnet1a, newnet1b, newnet, res1, res2, res3, res4, res, trig, a}, a = {0, 1/2, 1/2}; lnet = net123[net0, m, debug]; theta = thetanet[lnet, m, debug]; theta1 = theta[[1]]; theta2 = theta[[2]]; net = subdecas3sa[theta1, m, a]; (* to get rta from art, i.e, put the singularity at the origin *) net = transnetk[net, m]; Print[" Net before blowing up: ", net]; trig = {{1, 0}, {0, 0}, {0, 1/2}}; res1 = sqblowup[theta1, trig, m, nm, debug]; Print[" rectangular res1 = ", res1]; res2 = rectotr[res1, m - nm, m, debug]; Print[" res2 = ", res2]; net = subdecas3ra[theta1, m, a]; net = transnetk[net, m]; Print[" Net before blowing up: ", net]; trig = {{1, 0}, {0, 1}, {0, 1/2}}; res3 = sqblowup[theta1, trig, m, nm, debug]; Print[" rectangular res3 = ", res3]; res4 = rectotr[res3, m - nm, m, debug]; Print[" res4 = ", res4]; res2 = piomeg[res2]; res4 = piomeg[res4]; res = Join[{res2}, {res4}]; res ]; (* to test the blowup programs with theta1, etc, rectangular net *) testblow3[{oldnet__}, m_, nm_, debug_] := Block[ {net0 = {oldnet}, theta, theta1, theta2, lnet, newnet1a, newnet1b, newnet, res1, res2, res3, res4, res, trig, a}, a = {0, 1/2, 1/2}; lnet = net123[net0, m, debug]; theta = thetanet[lnet, m, debug]; theta1 = theta[[1]]; theta2 = theta[[2]]; net = subdecas3sa[theta1, m, a]; (* to get rta from art, i.e, put the singularity at the origin *) net = transnetk[net, m]; Print[" Net before blowing up: ", net]; trig = {{1, 0}, {0, 0}, {0, 1/2}}; res1 = sqblowup[theta1, trig, m, nm, debug]; Print[" res1 = ", res1]; net = subdecas3ra[theta1, m, a]; net = transnetk[net, m]; Print[" Net before blowing up: ", net]; trig = {{1, 0}, {0, 1}, {0, 1/2}}; res2 = sqblowup[theta1, trig, m, nm, debug]; Print[" res2 = ", res2]; res1 = piomeg[res1]; res2 = piomeg[res2]; res = Join[{res1}, {res2}]; res ]; (* Subdivides into four subpatches, by dividing the base triangle *) (* into four subtriangles using the middles of the original sides *) (* uses a tricky scheme involving 5 de Casteljau steps *) (* also takes care of net with a degenerate upper corner *) (* by performing a blow up *) subdecas4and5[{net__}, oldm_, debug_] := Block[ {basenet = {net}, newnet1a, newnet1b, newnet1, newnet2a, newnet2b, newnet2, newnet, m, pair, i, j, a, anet, art, ars, nm, reftrig, trig}, (a = {0, 1/2, 1/2}; reftrig = {{1, 0, 0}, {0, 1, 0}, {1, 1, -1}}; m = netsize[basenet]; Print[" netsize m in subdecas4 = ", m]; art = subdecas3sa[basenet, m, a]; If[debug === -20, Print["*** Left corner patch art: ", art]]; pair = zerocorner[art, m, debug]; (* Print["*** pair in subdecas4: ", pair]; *) If[pair[[1]] === -1, newnet = nsubdecas45[art, m, debug], (* Exceptional case where top corner triangle is degenerate *) Print[" ** Degenerate corner in net art: ", art]; Print[" Attempting to repair the net by blowing up the corner."]; i = pair[[1]]; nm = m - i; (* to get rta from art, i.e, put the singularity at the origin *) anet = transnetk[art, m]; trig = {{1, 0}, {0, 0}, {0, 1/2}}; newnet1a = sqblowup[basenet, trig, m, nm, debug]; Print[" ** Blowing up to a (", i, " , ", m, ")-rectangular patch done"]; Print[" ** Now converting to a triangular patch of degre ", 2*m - nm]; newnet1a = rectotr[newnet1a, m - nm, m, debug]; Print[" ** Repair done, newnet1a : ", newnet1a]; newnet1b = newcnet3[newnet1a, 2*m - nm, reftrig]; Print[" ** Repair done, newnet1b: ", newnet1b]; newnet1 = Join[{newnet1a}, {newnet1b}]; (* Now, treats the other half ars *) ars = subdecas3ra[basenet, m, a]; pair = zerocorner[ars, m, debug]; If[pair[[1]] === -1, newnet2 = ars, If[debug === -20, Print["*** Right corner patch ars: ", ars]]; Print[" ** Degenerate corner in net ars: ", ars]; Print[" Attempting to repair the net by blowing up the corner."]; i = pair[[1]]; nm = m - i; (* to get rsa from ars, i.e, put the singularity at the origin *) anet = transnetk[ars, m]; trig = {{1, 0}, {0, 1}, {0, 1/2}}; newnet2a = sqblowup[basenet, trig, m, nm, debug]; Print[" ** Blowing up to a (", i, " , ", m, ")-rectangular patch done"]; Print[" ** Now converting to a triangular patch of degree ", 2*m - nm]; newnet2a = rectotr[newnet2a, m - nm, m, debug]; Print[" ** Repair done, newnet2a: ", newnet2a]; newnet2b = newcnet3[newnet2a, 2*m - nm, reftrig]; Print[" ** Repair done, newnet2b: ", newnet2b]; newnet2 = Join[{newnet2a}, {newnet2b}] ]; newnet = Join[newnet1, newnet2] ]; If[debug === -20, Print["Subdivided nets: ", newnet]]; newnet ) ]; (* Subdivides into four subpatches, by dividing the base triangle *) (* into four subtriangles using the middles of the original sides *) (* uses a tricky scheme involving 4 de Casteljau steps *) (* also takes care of net with a degenerate upper corner *) (* by performing a blow up *) subdecas4[{net__}, oldm_, debug_] := Block[ {basenet = {net}, newnet1a, newnet1b, newnet1, newnet2a, newnet2b, newnet2, newnet, m, pair, i, j, a, anet, art, ars, nm, reftrig, trig}, (a = {0, 1/2, 1/2}; reftrig = {{1, 0, 0}, {0, 1, 0}, {1, 1, -1}}; m = netsize[basenet]; Print[" netsize m in subdecas4 = ", m]; newnet = sdecas3[basenet, m, a]; art = newnet[[1]]; ars = newnet[[3]]; If[debug === -20, Print["*** Left corner patch art: ", art]]; pair = zerocorner[art, m, debug]; (* Print["*** pair in subdecas4: ", pair]; *) If[pair[[1]] === -1, newnet = nsubdecas4[art, ars, m, debug], (* Exceptional case where top corner triangle is degenerate *) Print[" ** Degenerate corner in net art: ", art]; Print[" Attempting to repair the net by blowing up the corner."]; i = pair[[1]]; nm = m - i; (* to get rta from art, i.e, put the singularity at the origin *) anet = transnetk[art, m]; trig = {{1, 0}, {0, 0}, {0, 1/2}}; newnet1a = sqblowup[basenet, trig, m, nm, debug]; Print[" ** Blowing up to a (", i, " , ", m, ")-rectangular patch done"]; Print[" ** Now converting to a triangular patch of degre ", 2*m - nm]; newnet1a = rectotr[newnet1a, m - nm, m, debug]; Print[" ** Repair done, newnet1a : ", newnet1a]; newnet1b = newcnet3[newnet1a, 2*m - nm, reftrig]; Print[" ** Repair done, newnet1b: ", newnet1b]; newnet1 = Join[{newnet1a}, {newnet1b}]; (* Now, treats the other half ars *) pair = zerocorner[ars, m, debug]; If[pair[[1]] === -1, newnet2 = ars, If[debug === -20, Print["*** Right corner patch ars: ", ars]]; Print[" ** Degenerate corner in net ars: ", ars]; Print[" Attempting to repair the net by blowing up the corner."]; i = pair[[1]]; nm = m - i; (* to get rsa from ars, i.e, put the singularity at the origin *) anet = transnetk[ars, m]; trig = {{1, 0}, {0, 1}, {0, 1/2}}; newnet2a = sqblowup[basenet, trig, m, nm, debug]; Print[" ** Blowing up to a (", i, " , ", m, ")-rectangular patch done"]; Print[" ** Now converting to a triangular patch of degree ", 2*m - nm]; newnet2a = rectotr[newnet2a, m - nm, m, debug]; Print[" ** Repair done, newnet2a: ", newnet2a]; newnet2b = newcnet3[newnet2a, 2*m - nm, reftrig]; Print[" ** Repair done, newnet2b: ", newnet2b]; newnet2 = Join[{newnet2a}, {newnet2b}] ]; newnet = Join[newnet1, newnet2] ]; If[debug === -20, Print["Subdivided nets: ", newnet]]; newnet ) ]; (* Tries to repair a zero corner when the top nonzero row is consistent, ie *) (* when all elements are equal *) fzerocorner[{net0__}, pair_, m_, debug_] := Block[ {oldnet = {net0}, dn, i, a, j, res}, ( If[Length[pair] === 3, i = pair[[2]]; a = pair[[3]]; dn = ((m - i)*(m - i + 1))/2; res = Drop[oldnet, -dn]; Do[ res = Append[res, a], {j, 1, dn} ]; Print["Repaired net in fzerocorner: ", res], res = oldnet ]; res ) ]; (* Detects whether a corner of a net contains zero vectors *) (* And if so, if all the elements on the top nonzero row are identical *) zerocorner[{net__}, oldm_, debug_] := Block[ {cc = {net}, anet, newnet = {}, m = oldm, pt, stop1, stop2, i, j, res, flag}, (anet = convtomat[cc, m]; (* checks it top corner is zero *) stop1 = 0; i = m; While[stop1 === 0 && i >= 0, j = 0; stop2 = 0; While[stop2 === 0 && j <= m - i, pt = anet[[i+1, j+1]]; If[nearzero[pt] === 1, j = j + 1, stop2 = 1] ]; If[stop2 === 0, i = i - 1, stop1 = 1] ]; (* checks it top nonzero row is consistent *) If[stop1 === 1 && i < m, flag = 0; j = 0; While[flag === 0 && j < m - i, If[anet[[i + 1, j + 1]] === anet[[i + 1, j + 2]], j = j + 1, flag = 1 ] ]; If[flag === 1, res = {i, j}, res = {-1, i, anet[[i + 1, 1]]}], res = {-1} ]; res ) ]; (* Fractal version a la Sierpinski Gasket *) (* Subdivides into four subpatches, by dividing the base triangle *) (* into four subtriangles using the middles of the original sides *) (* uses a tricky scheme involving 5 de Casteljau steps *) (* omits the central triangle, for fractal effect *) (* returns abt, sca, crb *) fracdecas4[{oldnet__}, oldm_] := Block[ {net = {oldnet}, newnet = {}, barcor1, barcor2, barcor3, art, bat, cba, rcb, sca, abt, crb, m = oldm}, (barcor1 = {1, 1, -1}; barcor2 = {0, 1/2, 1/2}; barcor3 = {1, -1, 1}; Print[" netsize m in fracdecas4 = ", m]; art = subdecas3sa[net, m, barcor2]; bat = subdecas3sa[art, m, barcor2]; cba = subdecas3ta[bat, m, barcor1]; rcb = subdecas3ta[cba, m, barcor1]; sca = subdecas3sa[cba, m, barcor3]; abt = transnetj[bat, m]; crb = transnetj[rcb, m]; newnet = Join[{abt}, {sca}]; (* omits central patch bac *) newnet = Join[newnet, {crb}]; newnet ) ]; (* performs a subdivision step on each net in a list *) (* using splitdecas4, i.e., splitting the reference triangle *) (* into two halves, and then again each half into two halves *) itersplit4[{net__}, m_] := Block[ {cnet = {net}, lnet = {}, l, i}, (l = Length[cnet]; Do[ lnet = Join[lnet, splitdecas4[cnet[[i]], m]] , {i, 1, l} ]; lnet ) ]; (* performs a subdivision step on each net in a list *) (* using splitdecas6, i.e., splitting into six patches *) (* First split into three triangles using the center of gravity *) (* and then split each subtriangle into two *) itersplit6[{net__}, m_] := Block[ {cnet = {net}, lnet = {}, l, i}, (l = Length[cnet]; Do[ lnet = Join[lnet, splitdecas6[cnet[[i]], m]] , {i, 1, l} ]; lnet ) ]; (* performs a subdivision step on each net in a list *) (* using sdecas3, i.e., splitting into three patches *) (* This version does NOT converge to a triangulation *) (* Since it does not make progress on the sides of *) (* reference triangles *) itersub3[{net__}, m_, {a__}] := Block[ {cnet = {net}, lnet = {}, aa = {a}, l, i}, (l = Length[cnet]; Do[ lnet = Join[lnet, sdecas3[cnet[[i]], m, aa]] , {i, 1, l} ]; lnet ) ]; (* performs a subdivision step on each net in a list *) (* using subdecas4, i.e., splitting the reference triangle *) (* into a regular pattern of four subtriangles *) (* This version calls subdecas4, (5 de Casteljau steps) *) itersub4[{net__}, m_, debug_] := Block[ {cnet = {net}, lnet = {}, l, i}, (l = Length[cnet]; Do[ lnet = Join[lnet, subdecas4[cnet[[i]], m, debug]] , {i, 1, l} ]; lnet ) ]; (* performs a subdivision step on each net in a list *) (* using fracdecas4, i.e., splitting into four patches *) (* and omitting the central triangle for fractal effect *) iterfrac4[{net__}, m_] := Block[ {cnet = {net}, lnet = {}, l, i}, (l = Length[cnet]; Do[ lnet = Join[lnet, fracdecas4[cnet[[i]], m]] , {i, 1, l} ]; lnet ) ]; (* To project a point in the hat space back onto the affine space *) (* version assigning a point far away when w near zero *) wprojpt[pt_] := Block[ {h, w}, w = Last[pt]; h = Drop[pt, -1]; If[Abs[w] > 10^(-20), If[Abs[w] <= 10^(-16), Print["*** Warning: Point with very small weight: ", pt," ***"]]; h = h/w, Print["*** Warning: Point at infinity!: ", pt, " ***"]; h = 10^(20) * h ]; h ]; (* To project a point in the hat space back onto the affine space *) (* version returning the origin when w near infinity *) projpt[pt_] := Block[ {h, w}, w = Last[pt]; h = Drop[pt, -1]; If[Abs[w] > 10^(-20), If[Abs[w] <= 10^(-16), Print["*** Warning: Point with very small weight: ", pt," ***"]]; h = h/w, Print["*** Warning: Point at infinity!: ", pt, " ***"]; h = 0 * h ]; h ]; (* To project down onto affine space, keeping the weights *) piomeg[{net__}] := Block[ {snet = {net}, newnet = {}, pt, newpt, j, l2, h, w}, l2 = Length[snet]; Do[ pt = snet[[j]]; w = Last[pt]; h = Drop[pt, -1]; If[w =!= 0 && w =!= 0., h = h/w]; newpt = Append[h, w]; newnet = Append[newnet,newpt], {j, 1, l2} ]; newnet ]; (* To project a net in the hat space back onto the affine space *) projnet[{net__}] := Block[ {snet = {net}, newnet = {}, pt, h, j, l2, flag}, l2 = Length[snet]; Do[ pt = snet[[j]]; h = projpt[pt]; newnet = Append[newnet,h], {j, 1, l2} ]; newnet ]; (* To compute the correct value of a degenerate point *) repairpt[{oldnet__},mm_] := Block[ {oldn = {oldnet}, pp, pol1, pol2, pta, ptb, curv1, lcurv}, pta = {1, 1, -1}; ptb = {-1, -1, 3}; (* Exceptional case where a point is near the zero vector *) pol1 = curvcpoly[oldn, mm, pta, ptb]; pol2 = negpoly[pol1]; curv1 = gsubdiv[pol2, -1, 1, -1, 1, 1]; lcurv = curv1[[1]]; pp = Last[lcurv]; pp ]; (* To project a list of nets in the hat space back onto the affine space *) projlis[{netlis__}] := Block[ {slis = {netlis}, newlis = {}, anet, newnet, j, l2}, l2 = Length[slis]; Do[ anet = slis[[j]]; newnet = projnet[anet]; newlis = Append[newlis,newnet], {j, 1, l2} ]; newlis ]; (* To map a net in the affine space into the hat space *) maptohat[{net__}] := Block[ {lnet = {net}, newnet = {}, pt, h, w, i, l1}, l1 = Length[lnet]; Do[ pt = lnet[[i]]; w = Last[pt]; h = Drop[pt, -1]; If[w =!= 0, h = w * h; pt = Append[h, w] ]; newnet = Append[newnet,pt], {i, 1, l1} ]; newnet ]; (* To strip the weights from a list of points *) strip[{net__}] := Block[ {snet = {net}, newnet = {}, j, l2, h}, l2 = Length[snet]; Do[ pt = snet[[j]]; h = Drop[pt, -1]; newnet = Append[newnet,h], {j, 1, l2} ]; newnet ]; (* performs n subdivision steps using itersub3, i.e., splits into 3 patches *) (* This version does NOT converge to a triangulation *) rsubdiv3[{net__}, m_, n_] := Block[ {newnet = {net}, i, a}, ( a = {1/3, 1/3, 1/3}; newnet = {newnet}; Do[ newnet = itersub3[newnet, m, a], {i, 1, n} ]; newnet ) ]; (* performs n subdivision steps using itersub4, i.e., *) (* splits the reference triangle into a regular pattern *) (* of four triangles *) rsubdiv4[{net__}, m_, n_, debug_] := Block[ {innet = {net}, i}, ( (* When Length[innet] === 2, we have two symmetric nets *) If[Length[innet] === 2, newnet = innet , newnet = {innet} ]; Do[ newnet = itersub4[newnet, m, debug], {i, 1, n} ]; newnet ) ]; (* performs n subdivision steps using iterfrac4, i.e., splits into 4 patches *) (* and does fractal effect by omitting central triangle *) rfracdiv4[{net__}, m_, n_] := Block[ {innet = {net}, i}, ( (* When Length[innet] === 2, we have two symmetric nets *) If[Length[innet] === 2, newnet = innet, newnet = {innet} ]; Do[ newnet = iterfrac4[newnet, m], {i, 1, n} ]; newnet ) ]; (* performs n subdivision steps using itersplit4 *) (* i.e., splits into 4 subpatches, by splitting the reference triangle *) (* into two halves, and then again each half into two halves *) subdivrec4[{net__}, m_, n_] := Block[ {newnet = {net}, i}, ( newnet = {newnet}; Do[ newnet = itersplit4[newnet, m], {i, 1, n} ]; newnet ) ]; (* performs n subdivision steps using itersplit4 *) (* i.e., splits into 4 subpatches, and finishes with split24 *) subdivr4[{net__}, m_, n_] := Block[ {res, newnet = {net}, i}, ( newnet = {newnet}; Do[ newnet = itersplit4[newnet, m], {i, 1, n} ]; res = split24[newnet, m]; res ) ]; (* performs n subdivision steps using itersplit6 *) (* i.e., splits into 6 subpatches *) subdivrec6[{net__}, m_, n_] := Block[ {res = {}, newnet = {net}, i}, ( res = {newnet}; Do[ res = itersplit6[res, m], {i, 1, n} ]; (* res = projlis[newnet]; *) res ) ]; (* performs n subdivision steps using first itersub3, *) (* i.e., splits into 3 patches, followed by itersub4, *) (* i.e., splits into 4 patches, according to regular pattern *) rsubdivmix[{net__}, m_, n_, debug_] := Block[ {newnet = {net}, i, a}, ( a = {1/3, 1/3, 1/3}; newnet = {newnet}; Do[ newnet = itersub3[newnet, m, a]; newnet = itersub4[newnet, m, debug], {i, 1, n} ]; newnet ) ]; (* performs n subdivision steps using first itersub4, *) (* i.e., splits into 4 patches, according to regular pattern *) (* i.e., splits into 4 patches, followed by itersub3, *) rsubdivmixb[{net__}, m_, n_, debug_] := Block[ {newnet = {net}, i, a}, ( a = {1/3, 1/3, 1/3}; newnet = {newnet}; Do[ newnet = itersub4[newnet, m, debug]; newnet = itersub3[newnet, m, a], {i, 1, n} ]; newnet ) ]; (* prepares triangular control net *) controlnet3D[{net__},m_] := Block[ {res, net2}, net2 = strip[{net}]; net2 = triangul[net2,m]; net2 = buildmesh[net2]; res = Prepend[net2, RGBColor[0,1,0]]; res ]; (* To display the result of n subdivision steps using rsubdiv3 *) (* i.e. splitting into 3 patches recursively *) subdiv3[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net2, surf}, net2 = maptohat[{net}]; surf = rsubdiv3[net2, mm, n]; If[debug === 20, Print["surf = ", surf]]; surf = fractriang[debug,surf,mm]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splits the reference triangle into a regular pattern *) (* of four triangles *) subdiv4[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net2, surf, stop, oldnet}, stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net2 = maptohat[{net}]; oldnet = net2; surf = rsubdiv4[net2, mm, n, debug]; surf = ltriang[debug,surf,mm]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting the reference triangle into 4 subtriangles *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) render1[{net__}, mm_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net2, surf, image, stop, oldnet}, stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net2 = maptohat[{net}]; oldnet = net2; surf = rsubdiv4[net2, mm, n, debug]; If[light === 1, surf = trianglight[surf,mm], surf = ltriang[light,surf,mm]]; If[fcnet === 1, image = {surf, controlnet3D[{net},mm]}, image = surf ]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using rfracdiv4 *) (* i.e. splitting the reference triangle into 4 subtriangles *) (* and performs fractal effect *) fractal4[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,light_] := Block[ {net2, surf}, net2 = maptohat[{net}]; surf = rfracdiv4[net2, mm, n]; surf = fractriang[light,surf,mm]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using splitdecas4 *) (* i.e. subdividing 3 times, getting four subpatches *) (* The triangulation is somewhat odd, lots of diamonds *) fastdiv4[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,light_] := Block[ {net2, surf}, net2 = maptohat[{net}]; surf = subdivrec4[net2, mm, n]; surf = ltriang[light,surf,mm]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using splitdecas6 *) (* i.e. subdividing 4 times, getting six subpatches *) (* spider-web style subdivision *) fastdiv6[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net2, surf,oldnet}, net2 = maptohat[{net}]; oldnet = net2; surf = subdivrec6[net2, mm, n]; surf = fractriang[debug,surf,mm]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[debug], DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using *) (* rsubdiv3 i.e. splitting into 3 patches, followed by rsubdiv4 *) (* i.e. splitting into 4 patches. *) mixsubdiv[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net2, surf}, net2 = maptohat[{net}]; surf = rsubdivmix[net2, mm, n, debug]; If[debug === 20, Print["surf = ", surf]]; surf = fractriang[debug,surf,mm]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, first using a new *) (* reference triangle *) subd4[{net__}, mm_, {newtrig__}, n_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net2, surf, oldnet}, net2 = maptohat[{net}]; oldnet = net2; net2 = newcnet3[net2,mm,{newtrig}]; Print["New Net: ", net2]; surf = ltriang[debug,rsubdiv4[net2, mm, n, debug],mm]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* flips the sign of control points in the hat space *) signflip[{pt__}] := Block[ {apt = {pt}, h, w, res}, w = Last[apt]; h = Drop[apt,-1]; If[w === 0, h = -h, w = -w]; res = Append[h, w]; res ]; (* Computes new nets 3, 4, 5, 6, involved in showing all the patches *) (* *) new4net[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,light_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, surf1, surf2, surf3, surf4, surf5, surf6, surf, oldnet, reftrig1, reftrig2, reftrig3, pt, i, j, k, debug}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; net0 = maptohat[{net}]; oldnet = net0; debug = 0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; anet = convtomat[net2,mm]; theta1 = {}; Do[ Do[ pt = anet[[j + 1, mm - i - j + 1]]; If[OddQ[i + j], pt = -pt]; theta1 = Append[theta1, pt], {j, 0, mm - i} ], {i, 0, mm} ]; net3 = newcnet3[net0,mm,reftrig3]; Print["New Net 3: ", net3]; anet = convtomat[net3,mm]; theta2 = {}; Do[ Do[ pt = anet[[i + 1, j + 1]]; If[OddQ[mm - i - j], pt = -pt]; theta2 = Append[theta2, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net theta1: ", theta1]; Print["New Net theta2: ", theta2]; rho1 = {}; Do[ Do[ pt = anet[[j + 1, mm - i - j + 1]]; If[OddQ[j], pt = -pt]; rho1 = Append[rho1, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net rho1: ", rho1]; anet = convtomat[net1,mm]; rho2 = {}; Do[ Do[ pt = anet[[mm - i - j + 1, i + 1]]; If[OddQ[mm - j], pt = -pt]; rho2 = Append[rho2, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net rho2: ", rho2]; surf3 = rsubdiv4[theta1, mm, n, debug]; surf3 = ltriang[light,surf3,mm]; Print["Third patch done! "]; surf4 = rsubdiv4[theta2, mm, n, debug]; surf4 = ltriang[light,surf4,mm]; Print["Fourth patch done! "]; surf5 = rsubdiv4[rho1, mm, n, debug]; surf5 = ltriang[light,surf5,mm]; Print["Fifth patch done! "]; surf6 = rsubdiv4[rho2, mm, n, debug]; surf6 = ltriang[light,surf6,mm]; Print["Sixth patch done! "]; surf = Join[surf3, surf4]; surf = Join[surf, surf5]; surf = Join[surf, surf6]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* Computes new nets 5 and 6 involved in showing all the patches *) (* *) new56net[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, surf1, surf2, surf3, surf4, surf5, surf6, surf, oldnet, reftrig1, reftrig2, reftrig3, pt, i, j, k, debug}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; net0 = maptohat[{net}]; oldnet = net0; debug = 0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; net3 = newcnet3[net0,mm,reftrig3]; Print["New Net 3: ", net3]; anet = convtomat[net3,mm]; rho1 = {}; Do[ Do[ pt = anet[[j + 1, mm - i - j + 1]]; If[OddQ[j], pt = -pt]; rho1 = Append[rho1, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net rho1: ", rho1]; anet = convtomat[net1,mm]; rho2 = {}; Do[ Do[ pt = anet[[mm - i - j + 1, i + 1]]; If[OddQ[mm - j], pt = -pt]; rho2 = Append[rho2, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net rho2: ", rho2]; surf5 = rsubdiv4[rho1, mm, n, debug]; surf5 = ltriang[light,surf5,mm]; Print["Fifth patch done! "]; surf6 = rsubdiv4[rho2, mm, n, debug]; surf6 = ltriang[light,surf6,mm]; Print["Sixth patch done! "]; surf = Join[surf5, surf6]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* Computes new nets 3 and 4 involved in showing all the patches *) (* *) new34net[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,light_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, surf1, surf2, surf3, surf4, surf5, surf6, surf, oldnet, reftrig1, reftrig2, reftrig3, pt, i, j, k, debug}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; debug = 0; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; anet = convtomat[net2,mm]; theta1 = {}; Do[ Do[ pt = anet[[j + 1, mm - i - j + 1]]; If[OddQ[i + j], pt = - pt]; theta1 = Append[theta1, pt], {j, 0, mm - i} ], {i, 0, mm} ]; net3 = newcnet3[net0,mm,reftrig3]; Print["New Net 3: ", net3]; anet = convtomat[net3,mm]; theta2 = {}; Do[ Do[ pt = anet[[i + 1, j + 1]]; If[OddQ[mm - i - j], pt = -pt]; theta2 = Append[theta2, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net theta1: ", theta1]; Print["New Net theta2: ", theta2]; surf3 = rsubdiv4[theta1, mm, n, debug]; surf3 = ltriang[light,surf3,mm]; Print["Third patch done! "]; surf4 = rsubdiv4[theta2, mm, n, debug]; surf4 = ltriang[light,surf4,mm]; Print["Fourth patch done! "]; surf = Join[surf3, surf4]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* Computes new net 4 involved in showing all the patches *) (* *) new4anet[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,light_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, surf1, surf2, surf3, surf4, surf5, surf6, surf, oldnet, reftrig1, reftrig2, reftrig3, pt, i, j, k, debug}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; net0 = maptohat[{net}]; oldnet = net0; debug = 0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; net3 = newcnet3[net0,mm,reftrig3]; Print["New Net 3: ", net3]; anet = convtomat[net3,mm]; theta2 = {}; Do[ Do[ pt = anet[[i + 1, j + 1]]; If[OddQ[mm - i - j], pt = -pt]; theta2 = Append[theta2, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net theta2: ", theta2]; surf4 = rsubdiv4[theta2, mm, n, debug]; surf = ltriang[light,surf4,mm]; Print["Fourth patch done! "]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* Computes new net 3 involved in showing all the patches *) (* *) new3anet[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,light_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, surf1, surf2, surf3, surf4, surf5, surf6, surf, oldnet, reftrig1, reftrig2, reftrig3, pt, i, j, k, debug}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; net0 = maptohat[{net}]; oldnet = net0; debug = 0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; anet = convtomat[net2,mm]; theta1 = {}; Do[ Do[ pt = anet[[j + 1, mm - i - j + 1]]; If[OddQ[i + j], pt = - pt]; theta1 = Append[theta1, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net theta1: ", theta1]; surf3 = rsubdiv4[theta1, mm, n, debug]; surf = ltriang[light,surf3,mm]; Print["Third patch done! "]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* Computes the six new nets involved in showing all the patches *) (* and displays the closed surface *) new6net[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,light_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, surf1, surf2, surf3, surf4, surf5, surf6, surf, oldnet, reftrig1, reftrig2, reftrig3, pt, i, j, k, debug}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; net0 = maptohat[{net}]; oldnet = net0; debug = 0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; anet = convtomat[net2,mm]; theta1 = {}; Do[ Do[ pt = anet[[j + 1, mm - i - j + 1]]; If[OddQ[i + j], pt = -pt]; theta1 = Append[theta1, pt], {j, 0, mm - i} ], {i, 0, mm} ]; net3 = newcnet3[net0,mm,reftrig3]; Print["New Net 3: ", net3]; anet = convtomat[net3,mm]; theta2 = {}; Do[ Do[ pt = anet[[i + 1, j + 1]]; If[OddQ[mm - i - j], pt = -pt]; theta2 = Append[theta2, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net theta1: ", theta1]; Print["New Net theta2: ", theta2]; rho1 = {}; Do[ Do[ pt = anet[[j + 1, mm - i - j + 1]]; If[OddQ[j], pt = -pt]; rho1 = Append[rho1, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net rho1: ", rho1]; anet = convtomat[net1,mm]; rho2 = {}; Do[ Do[ pt = anet[[mm - i - j + 1, i + 1]]; If[OddQ[mm - j], pt = -pt]; rho2 = Append[rho2, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net rho2: ", rho2]; surf1 = rsubdiv4[net1, mm, n, debug]; surf1 = ltriang[light,surf1,mm]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; surf2 = ltriang[light,surf2,mm]; Print["Second patch done! "]; surf3 = rsubdiv4[theta1, mm, n, debug]; surf3 = ltriang[light,surf3,mm]; Print["Third patch done! "]; surf4 = rsubdiv4[theta2, mm, n, debug]; surf4 = ltriang[light,surf4,mm]; Print["Fourth patch done! "]; surf5 = rsubdiv4[rho1, mm, n, debug]; surf5 = ltriang[light,surf5,mm]; Print["Fifth patch done! "]; surf6 = rsubdiv4[rho2, mm, n, debug]; surf6 = ltriang[light,surf6,mm]; Print["Sixth patch done! "]; surf = Join[surf1, surf2]; surf = Join[surf, surf3]; surf = Join[surf, surf4]; surf = Join[surf, surf5]; surf = Join[surf, surf6]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* Computes the six new nets involved in showing all the patches *) (* and displays the closed surface *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) render[{net__}, mm_, n_, light_, fcnet_, lx_, mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, oldnet, surf1, surf2, surf3, surf4, surf5, surf6, surf, image, stop, theta, rho, lnet, reftrig1, reftrig2, reftrig3, pt, i, j, k}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; lnet = net123[{net}, mm, debug]; net1 = lnet[[1]]; net2 = lnet[[2]]; theta = thetanet[lnet, mm, debug]; theta1 = theta[[1]]; theta2 = theta[[2]]; rho = rhonet[lnet, mm, debug]; rho1 = rho[[1]]; rho2 = rho[[2]]; surf1 = rsubdiv4[net1, mm, n, debug]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = ltriang[light,surf1,mm]]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = ltriang[light,surf2,mm]]; Print["Second patch done! "]; surf3 = rsubdiv4[theta1, mm, n, debug]; If[light === 1, surf3 = trianglight[surf3,mm], surf3 = ltriang[light,surf3,mm]]; Print["Third patch done! "]; surf4 = rsubdiv4[theta2, mm, n, debug]; If[light === 1, surf4 = trianglight[surf4,mm], surf4 = ltriang[light,surf4,mm]]; Print["Fourth patch done! "]; surf5 = rsubdiv4[rho1, mm, n, debug]; If[light === 1, surf5 = trianglight[surf5,mm], surf5 = ltriang[light,surf5,mm]]; Print["Fifth patch done! "]; surf6 = rsubdiv4[rho2, mm, n, debug]; If[light === 1, surf6 = trianglight[surf6,mm], surf6 = ltriang[light,surf6,mm]]; Print["Sixth patch done! "]; surf = Join[surf1, surf2]; surf = Join[surf, surf3]; surf = Join[surf, surf4]; surf = Join[surf, surf5]; surf = Join[surf, surf6]; Print["Ready to display! "]; If[fcnet === 1, image = {surf, controlnet3D[{net},mm]}, image = surf ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* Computes the six new nets involved in showing all the patches *) (* and displays the closed surface *) (* if light = 1, uses light triangulation trianglight, *) (* This version colors patches 1,2 in blue, 3,4 in red, 5,6 in green *) (* if fcnet = 1, displays control net, otherwise skip it *) render3col[{net__}, mm_, n_, light_, fcnet_, lx_, mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, oldnet, surf1, surf2, surf3, surf4, surf5, surf6, surf, image, stop, theta, rho, lnet, reftrig1, reftrig2, reftrig3, pt, i, j, k}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; lnet = net123[{net}, mm, debug]; net1 = lnet[[1]]; net2 = lnet[[2]]; theta = thetanet[lnet, mm, debug]; theta1 = theta[[1]]; theta2 = theta[[2]]; rho = rhonet[lnet, mm, debug]; rho1 = rho[[1]]; rho2 = rho[[2]]; surf1 = rsubdiv4[net1, mm, n, debug]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = ltriang[11,surf1,mm]]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = ltriang[11,surf2,mm]]; Print["Second patch done! "]; surf3 = rsubdiv4[theta1, mm, n, debug]; If[light === 1, surf3 = trianglight[surf3,mm], surf3 = ltriang[12,surf3,mm]]; Print["Third patch done! "]; surf4 = rsubdiv4[theta2, mm, n, debug]; If[light === 1, surf4 = trianglight[surf4,mm], surf4 = ltriang[12,surf4,mm]]; Print["Fourth patch done! "]; surf5 = rsubdiv4[rho1, mm, n, debug]; If[light === 1, surf5 = trianglight[surf5,mm], surf5 = ltriang[13,surf5,mm]]; Print["Fifth patch done! "]; surf6 = rsubdiv4[rho2, mm, n, debug]; If[light === 1, surf6 = trianglight[surf6,mm], surf6 = ltriang[13,surf6,mm]]; Print["Sixth patch done! "]; surf = Join[surf1, surf2]; surf = Join[surf, surf3]; surf = Join[surf, surf4]; surf = Join[surf, surf5]; surf = Join[surf, surf6]; Print["Ready to display! "]; If[fcnet === 1, image = {surf, controlnet3D[{net},mm]}, image = surf ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* Computes the new nets net1, net2, net3, in order to compute *) (* theta1, theta2, rho1, rho2 *) net123[{oldnet__}, mm_, debug_] := Block[ {net = {oldnet}, net0, net1, net2, net3, res, reftrig1, reftrig2, reftrig3, stop}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; stop = testm[net, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[net]; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; net3 = newcnet3[net0,mm,reftrig3]; Print["New Net 3: ", net3]; res = Join[{net1}, {net2}]; res = Join[res, {net3}]; res ]; (* Computes the new nets phi1, phi2, phi3 *) phinets[{inlnet__}, mm_, debug_] := Block[ {net = {inlnet}, anet, phi1, phi2, phi3, pt, i, j, res}, anet = convtomat[net,mm]; phi1 = {}; Do[ Do[ pt = anet[[i + 1, j + 1]]; If[OddQ[i], pt = -pt]; phi1 = Append[phi1, pt], {j, 0, mm - i} ], {i, 0, mm} ]; anet = convtomat[net,mm]; phi2 = {}; Do[ Do[ pt = anet[[i + 1, j + 1]]; If[OddQ[j], pt = -pt]; phi2 = Append[phi2, pt], {j, 0, mm - i} ], {i, 0, mm} ]; anet = convtomat[net,mm]; phi3 = {}; Do[ Do[ pt = anet[[i + 1, j + 1]]; If[OddQ[mm - i - j], pt = -pt]; phi3 = Append[phi3, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net phi1: ", phi1]; Print["New Net phi2: ", phi2]; Print["New Net phi3: ", phi3]; res = Join[{phi1}, {phi2}]; res = Join[res, {phi3}]; res ]; (* Computes the new nets theta1 and theta2 *) thetanet[{inlnet__}, mm_, debug_] := Block[ {net1, net2, net3, lnet = {inlnet}, anet, theta1, theta2, pt, i, j, k, res}, net1 = lnet[[1]]; net2 = lnet[[2]]; net3 = lnet[[3]]; anet = convtomat[net2,mm]; theta1 = {}; Do[ Do[ pt = anet[[j + 1, mm - i - j + 1]]; If[OddQ[i + j], pt = -pt]; theta1 = Append[theta1, pt], {j, 0, mm - i} ], {i, 0, mm} ]; anet = convtomat[net3,mm]; theta2 = {}; Do[ Do[ pt = anet[[i + 1, j + 1]]; If[OddQ[mm - i - j], pt = -pt]; theta2 = Append[theta2, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net theta1: ", theta1]; Print["New Net theta2: ", theta2]; res = Join[{theta1}, {theta2}]; res ]; (* Computes the new nets rho1 and rho2 *) rhonet[{inlnet__}, mm_, debug_] := Block[ {net1, net2, net3, lnet = {inlnet}, anet, rho1, rho2, pt, i, j, k, res}, net1 = lnet[[1]]; net2 = lnet[[2]]; net3 = lnet[[3]]; anet = convtomat[net3,mm]; rho1 = {}; Do[ Do[ pt = anet[[j + 1, mm - i - j + 1]]; If[OddQ[j], pt = -pt]; rho1 = Append[rho1, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net rho1: ", rho1]; anet = convtomat[net1,mm]; rho2 = {}; Do[ Do[ pt = anet[[mm - i - j + 1, i + 1]]; If[OddQ[mm - j], pt = -pt]; rho2 = Append[rho2, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net rho2: ", rho2]; res = Join[{rho1}, {rho2}]; res ]; (* To test the computation of nets theta1, theta2, rho1, rho2 *) testnet[{oldnet__}, mm_, debug_] := Block[ {net = {oldnet}, lnet, theta, rho, theta1, theta2, theta2b, rho1, rho2, rho2b, res, reftrig1, reftrig2, reftrig3, reftrig4}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; reftrig4 = {{-1, 1, 1}, {0, 0, 1}, {0, 1, 0}}; lnet = net123[net, mm, debug]; theta = thetanet[lnet, mm, debug]; theta1 = theta[[1]]; theta2 = theta[[2]]; theta2b = newcnet3[theta1,mm,reftrig4]; Print["theta2: ", theta2]; Print["theta2b: ", theta2b]; rho = rhonet[lnet, mm, debug]; rho1 = rho[[1]]; rho2 = rho[[2]]; rho2b = newcnet3[rho1,mm,reftrig4]; Print["rho2: ", rho2]; Print["rho2b: ", rho2b]; res = Join[{theta1}, {theta2}]; res = Join[res, {rho1}]; res = Join[res, {rho2}]; res ]; (* Computes the four new nets involved in showing all the patches *) (* and displays the closed surface *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) render3to6[{net__}, mm_, n_, light_, fcnet_, lx_, mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, oldnet, surf1, surf2, surf3, surf4, surf5, surf6, surf, image, stop, theta, rho, lnet, reftrig1, reftrig2, reftrig3, pt, i, j, k}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; lnet = net123[{net}, mm, debug]; theta = thetanet[lnet, mm, debug]; theta1 = theta[[1]]; theta2 = theta[[2]]; rho = rhonet[lnet, mm, debug]; rho1 = rho[[1]]; rho2 = rho[[2]]; surf3 = rsubdiv4[theta1, mm, n, debug]; If[debug === 20, Print["Third patch: ", surf3] ]; If[light === 1, surf3 = trianglight[surf3,mm], surf3 = ltriang[light,surf3,mm]]; Print["Third patch done! "]; surf4 = rsubdiv4[theta2, mm, n, debug]; If[debug === 20, Print["Fourh patch: ", surf4] ]; If[light === 20, surf4 = trianglight[surf4,mm], surf4 = ltriang[light,surf4,mm]]; Print["Fourth patch done! "]; surf5 = rsubdiv4[rho1, mm, n, debug]; If[debug === 20, Print["Fifth patch: ", surf5] ]; If[light === 20, surf5 = trianglight[surf5,mm], surf5 = ltriang[light,surf5,mm]]; Print["Fifth patch done! "]; surf6 = rsubdiv4[rho2, mm, n, debug]; If[debug === 20, Print["Sixth patch: ", surf6] ]; If[light === 20, surf6 = trianglight[surf6,mm], surf6 = ltriang[light,surf6,mm]]; Print["Sixth patch done! "]; surf = Join[surf3, surf4]; surf = Join[surf, surf5]; surf = Join[surf, surf6]; Print["Ready to display! "]; If[fcnet === 20, image = {surf, controlnet3D[{net},mm]}, image = surf ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* Computes the new nets 3 and 4 involved in showing all the patches *) (* and displays the closed surface *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) render34[{net__}, mm_, n_, light_, fcnet_, lx_, mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, lnet, surf1, surf2, surf3, surf4, surf5, surf6, surf, image, stop, theta, i, j, k}, stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; lnet = net123[{net}, mm, debug]; theta = thetanet[lnet, mm, debug]; theta1 = theta[[1]]; theta2 = theta[[2]]; surf3 = rsubdiv4[theta1, mm, n, debug]; If[debug === 20, Print["Third patch: ", surf3] ]; If[light === 1, surf3 = trianglight[surf3,mm], surf3 = ltriang[light,surf3,mm]]; Print["Third patch done! "]; surf4 = rsubdiv4[theta2, mm, n, debug]; If[debug === 20, Print["Fourth patch: ", surf4] ]; If[light === 1, surf4 = trianglight[surf4,mm], surf4 = ltriang[light,surf4,mm]]; Print["Fourth patch done! "]; surf = Join[surf3, surf4]; Print["Ready to display! "]; If[fcnet === 1, image = {surf, controlnet3D[{net},mm]}, image = surf ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* Computes the new nets 5 and 6 involved in showing all the patches *) (* and displays the closed surface *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) render56[{net__}, mm_, n_, light_, fcnet_, lx_, mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, oldnet, surf1, surf2, surf3, surf4, surf5, surf6, surf, image, stop, rho, lnet, reftrig1, reftrig2, reftrig3, pt, i, j, k}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; lnet = net123[{net}, mm, debug]; rho = rhonet[lnet, mm, debug]; rho1 = rho[[1]]; rho2 = rho[[2]]; surf5 = rsubdiv4[rho1, mm, n, debug]; If[debug === 20, Print["Fifth patch: ", surf5] ]; If[light === 1, surf5 = trianglight[surf5,mm], surf5 = ltriang[light,surf5,mm]]; Print["Fifth patch done! "]; surf6 = rsubdiv4[rho2, mm, n, debug]; If[debug === 20, Print["Sixth patch: ", surf6] ]; If[light === 1, surf6 = trianglight[surf6,mm], surf6 = ltriang[light,surf6,mm]]; Print["Sixth patch done! "]; surf = Join[surf5, surf6]; Print["Ready to display! "]; If[fcnet === 1, image = {surf, controlnet3D[{net},mm]}, image = surf ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; testnet[{oldnet__}, mm_, debug_] := Block[ {net = {oldnet}, lnet, theta, rho, theta1, theta2, theta2b, rho1, rho2, rho2b, res, reftrig1, reftrig2, reftrig3, reftrig4}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; reftrig4 = {{-1, 1, 1}, {0, 0, 1}, {0, 1, 0}}; lnet = net123[net, mm, debug]; theta = thetanet[lnet, mm, debug]; theta1 = theta[[1]]; theta2 = theta[[2]]; theta2b = newcnet3[theta1,mm,reftrig4]; Print["theta2: ", theta2]; Print["theta2b: ", theta2b]; rho = rhonet[lnet, mm, debug]; rho1 = rho[[1]]; rho2 = rho[[2]]; rho2b = newcnet3[rho1,mm,reftrig4]; Print["rho2: ", rho2]; Print["rho2b: ", rho2b]; res = Join[{theta1}, {theta2}]; res = Join[res, {rho1}]; res = Join[res, {rho2}]; res ]; (* Computes the new net 3 involved in showing all the patches *) (* and displays the closed surface *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) render3[{net__}, mm_, n_, light_, fcnet_, lx_, mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, net3, lnet, anet, theta1, theta2, rho1, rho2, oldnet, surf1, surf2, surf3, surf4, surf5, surf6, surf, image, stop, theta, reftrig1, reftrig2, reftrig3, pt, i, j, k}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; lnet = net123[{net}, mm, debug]; theta = thetanet[lnet, mm, debug]; theta1 = theta[[1]]; theta2 = theta[[2]]; net0 = maptohat[{net}]; oldnet = net0; surf3 = rsubdiv4[theta1, mm, n, debug]; If[debug === 20, Print["Third patch: ", surf3] ]; If[light === 1, surf3 = trianglight[surf3,mm], surf3 = ltriang[light,surf3,mm]]; Print["Third patch done! "]; Print["Ready to display! "]; If[fcnet === 1, image = {surf3, controlnet3D[{net},mm]}, image = surf3 ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* Computes the new net 4 involved in showing all the patches *) (* and displays the closed surface *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) render4[{net__}, mm_, n_, light_, fcnet_, lx_, mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, oldnet, surf1, surf2, surf3, surf4, surf5, surf6, surf, image, stop, theta, lnet, reftrig1, reftrig2, reftrig3, pt, i, j, k}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; lnet = net123[{net}, mm, debug]; theta = thetanet[lnet, mm, debug]; theta1 = theta[[1]]; theta2 = theta[[2]]; surf4 = rsubdiv4[theta2, mm, n, debug]; If[debug === 20, Print["Fourth patch: ", surf4] ]; If[light === 1, surf4 = trianglight[surf4,mm], surf4 = ltriang[light,surf4,mm]]; Print["Fourth patch done! "]; Print["Ready to display! "]; If[fcnet === 1, image = {surf4, controlnet3D[{net},mm]}, image = surf4 ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* Computes the new net 5 involved in showing all the patches *) (* and displays the closed surface *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) render5[{net__}, mm_, n_, light_, fcnet_, lx_, mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, oldnet, surf1, surf2, surf3, surf4, surf5, surf6, surf, image, stop, rho, lnet, reftrig1, reftrig2, reftrig3, pt, i, j, k}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; lnet = net123[{net}, mm, debug]; rho = rhonet[lnet, mm, debug]; rho1 = rho[[1]]; rho2 = rho[[2]]; surf5 = rsubdiv4[rho1, mm, n, debug]; If[debug === 20, Print["Fifth patch: ", surf5] ]; If[light === 1, surf5 = trianglight[surf5,mm], surf5 = ltriang[light,surf5,mm]]; Print["Fifth patch done! "]; Print["Ready to display! "]; If[fcnet === 1, image = {surf5, controlnet3D[{net},mm]}, image = surf5 ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* Computes the new net 6 involved in showing all the patches *) (* and displays the closed surface *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) render6[{net__}, mm_, n_, light_, fcnet_, lx_, mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, oldnet, surf1, surf2, surf3, surf4, surf5, surf6, surf, image, stop, rho, lnet, reftrig1, reftrig2, reftrig3, pt, i, j, k}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; lnet = net123[{net}, mm, debug]; rho = rhonet[lnet, mm, debug]; rho1 = rho[[1]]; rho2 = rho[[2]]; surf6 = rsubdiv4[rho2, mm, n, debug]; If[debug === 20, Print["Sixth patch: ", surf6] ]; If[light === 1, surf6 = trianglight[surf6,mm], surf6 = ltriang[light,surf6,mm]]; Print["Sixth patch done! "]; Print["Ready to display! "]; If[fcnet === 1, image = {surf6, controlnet3D[{net},mm]}, image = surf6 ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* Computes the six new nets involved in showing all the patches *) (* and displays the closed surface *) (* fractal version *) (* if light = 1, uses light triangulation trianglight, *) (* else fractriang *) (* if fcnet = 1, displays control net, otherwise skip it *) rfrac6[{net__}, mm_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, oldnet, surf1, surf2, surf3, surf4, surf5, surf6, surf, image, stop, reftrig1, reftrig2, reftrig3, pt, i, j, k}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; lnet = net123[{net}, mm, debug]; net1 = lnet[[1]]; net2 = lnet[[2]]; theta = thetanet[lnet, mm, debug]; theta1 = theta[[1]]; theta2 = theta[[2]]; rho = rhonet[lnet, mm, debug]; rho1 = rho[[1]]; rho2 = rho[[2]]; surf1 = rfracdiv4[net1, mm, n]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = fractriang[light,surf1,mm]]; Print["First patch done! "]; surf2 = rfracdiv4[net2, mm, n]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = fractriang[light,surf2,mm]]; Print["Second patch done! "]; surf3 = rfracdiv4[theta1, mm, n]; If[light === 1, surf3 = trianglight[surf3,mm], surf3 = fractriang[light,surf3,mm]]; Print["Third patch done! "]; surf4 = rfracdiv4[theta2, mm, n]; If[light === 1, surf4 = trianglight[surf4,mm], surf4 = fractriang[light,surf4,mm]]; Print["Fourth patch done! "]; surf5 = rfracdiv4[rho1, mm, n]; If[light === 1, surf5 = trianglight[surf5,mm], surf5 = fractriang[light,surf5,mm]]; Print["Fifth patch done! "]; surf6 = rfracdiv4[rho2, mm, n]; If[light === 1, surf6 = trianglight[surf6,mm], surf6 = fractriang[light,surf6,mm]]; Print["Sixth patch done! "]; surf = Join[surf1, surf2]; surf = Join[surf, surf3]; surf = Join[surf, surf4]; surf = Join[surf, surf5]; surf = Join[surf, surf6]; Print["Ready to display! "]; If[fcnet === 1, image = {surf, controlnet3D[{net},mm]}, image = surf ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* Computes the six new nets involved in showing all the patches *) (* and displays the closed surface *) (* fractal version *) (* if light = 1, uses light triangulation trianglight, *) (* This version colors patches 1,2 in blue, 3,4 in red, 5,6 in green *) (* if fcnet = 1, displays control net, otherwise skip it *) rfrac3col[{net__}, mm_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_] := Block[ {net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, oldnet, surf1, surf2, surf3, surf4, surf5, surf6, surf, image, stop, reftrig1, reftrig2, reftrig3, pt, i, j, k}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; lnet = net123[{net}, mm, debug]; net1 = lnet[[1]]; net2 = lnet[[2]]; theta = thetanet[lnet, mm, debug]; theta1 = theta[[1]]; theta2 = theta[[2]]; rho = rhonet[lnet, mm, debug]; rho1 = rho[[1]]; rho2 = rho[[2]]; surf1 = rfracdiv4[net1, mm, n]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = fractriang[11,surf1,mm]]; Print["First patch done! "]; surf2 = rfracdiv4[net2, mm, n]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = fractriang[11,surf2,mm]]; Print["Second patch done! "]; surf3 = rfracdiv4[theta1, mm, n]; If[light === 1, surf3 = trianglight[surf3,mm], surf3 = fractriang[12,surf3,mm]]; Print["Third patch done! "]; surf4 = rfracdiv4[theta2, mm, n]; If[light === 1, surf4 = trianglight[surf4,mm], surf4 = fractriang[12,surf4,mm]]; Print["Fourth patch done! "]; surf5 = rfracdiv4[rho1, mm, n]; If[light === 1, surf5 = trianglight[surf5,mm], surf5 = fractriang[13,surf5,mm]]; Print["Fifth patch done! "]; surf6 = rfracdiv4[rho2, mm, n]; If[light === 1, surf6 = trianglight[surf6,mm], surf6 = fractriang[13,surf6,mm]]; Print["Sixth patch done! "]; surf = Join[surf1, surf2]; surf = Join[surf, surf3]; surf = Join[surf, surf4]; surf = Join[surf, surf5]; surf = Join[surf, surf6]; Print["Ready to display! "]; If[fcnet === 1, image = {surf, controlnet3D[{net},mm]}, image = surf ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, running it twice *) (* using the reference triangles *) (* reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}, and *) (* reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}} *) subd2[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_, debug_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, oldnet}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; surf1 = rsubdiv4[net1, mm, n, debug]; surf1 = ltriang[debug,surf1,mm]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; surf2 = ltriang[debug,surf2,mm]; Print["Second patch done! "]; surf = Join[surf1, surf2]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, running it twice *) (* using the reference triangles *) (* reftrig1 = {{-a, b, 1 + a - b}, {-a, -b, 1 + a + b}, *) (* {a, b, 1 - a - b}}, and *) (* reftrig2 = {{a, -b, 1 -a + b}, {a, b, 1 - a - b}, *) (* {-a, -b, 1 + a + b}} *) new2net[{net__}, mm_, a_, b_, n_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, oldnet}, reftrig1 = {{-a, b, 1 + a - b}, {-a, -b, 1 + a + b}, {a, b, 1 - a - b}}; reftrig2 = {{a, -b, 1 - a + b}, {a, b, 1 - a - b}, {-a, -b, 1 + a + b}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; surf1 = rsubdiv4[net1, mm, n, debug]; surf1 = ltriang[debug,surf1,mm]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; surf2 = ltriang[debug,surf2,mm]; Print["Second patch done! "]; surf = Join[surf1, surf2]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, running it twice *) (* as a square patch, using the reference triangles *) (* reftrig1 = {{c, b, 1 - b - c}, {a, d, 1 - a - d}, {a, b, 1 - a - b}}; *) (* reftrig2 = {{c, b, 1 - b - c}, {a, d, 1 - a - d}, {c, d, 1 - c - d}}; *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* This is wrt the square [a, b] x [c, d] *) (* This versions lists c, a, d, b, in a more intuitive order! *) (* Also, when light = 20, patch 1 is blue and patch 2 is red *) rendersq[{net__}, mm_, c_, a_, d_, b_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, oldnet}, (* reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; *) (* reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, {c, d, 1 - c - d}}; *) reftrig1 = {{c, b, 1 - b - c}, {a, d, 1 - a - d}, {a, b, 1 - a - b}}; reftrig2 = {{c, b, 1 - b - c}, {a, d, 1 - a - d}, {c, d, 1 - c - d}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; surf1 = rsubdiv4[net1, mm, n, debug]; If[light === 1, surf1 = trianglight[surf1,mm], If[light === 20, surf1 = ltriang[11,surf1,mm], surf1 = ltriang[light,surf1,mm] ] ]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; If[light === 1, surf2 = trianglight[surf2,mm], If[light === 20, surf2 = ltriang[12,surf2,mm], surf2 = ltriang[light,surf2,mm] ] ]; Print["Second patch done! "]; surf = Join[surf1, surf2]; Print["Ready to display! "]; If[fcnet === 1, Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction], Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction] ]; ]; (* To display the result of n subdivision steps using subdivrec4 *) (* running it twice as a square patch, using the reference *) (* triangles *) (* reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, *) (* {a, b, 1 - a - b}}, and *) (* reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, *) (* {c, d, 1 - c - d}} *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* This is wrt the square [c, a] x [d, b] *) (* This versions lists c, a, d, b, in a more intuitive order! *) rendersq2[{net__}, mm_, c_, a_, d_, b_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, oldnet}, reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, {c, d, 1 - c - d}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; surf1 = subdivrec4[net1, mm, n]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = ltriang[light,surf1,mm]]; Print["First patch done! "]; surf2 = subdivrec4[net2, mm, n]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = ltriang[light,surf2,mm]]; Print["Second patch done! "]; surf = Join[surf1, surf2]; Print["Ready to display! "]; If[fcnet === 1, Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction], Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction] ]; ]; (* To display the result of n subdivision steps using subdivrec4 *) (* running it twice as a square patch, using the reference *) (* triangles *) (* reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, *) (* {a, b, 1 - a - b}}, and *) (* reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, *) (* {c, d, 1 - c - d}} *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* This is wrt the square [c, a] x [d, b] *) (* This versions does not print the bounding box *) rendersq2b[{net__}, mm_, c_, a_, d_, b_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, oldnet}, reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, {c, d, 1 - c - d}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; surf1 = subdivrec4[net1, mm, n]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = ltriang[light,surf1,mm]]; Print["First patch done! "]; surf2 = subdivrec4[net2, mm, n]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = ltriang[light,surf2,mm]]; Print["Second patch done! "]; surf = Join[surf1, surf2]; Print["Ready to display! "]; If[fcnet === 1, Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> False, Boxed -> False, lightkind[light], DisplayFunction -> $DisplayFunction], Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> False, Boxed -> False, lightkind[light], DisplayFunction -> $DisplayFunction] ]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, running it twice *) (* as a square patch, using the reference triangles *) (* reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, *) (* {a, b, 1 - a - b}}, and *) (* reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, *) (* {c, d, 1 - c - d}} *) (* This version uses the input net, and also computes *) (* the symmetric net wrt the plane x = 0, and renders the whole *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* This is wrt the square [c, a] x [d, b] *) (* This versions lists c, a, d, b, in a more intuitive order! *) rendersymx[{net__}, mm_, c_, a_, d_, b_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, oldnet, sxnet1, sxnet2, xsurf1, xsurf2}, reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, {c, d, 1 - c - d}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; sxnet1 = symmx[net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; sxnet2 = symmx[net2]; surf1 = rsubdiv4[net1, mm, n, debug]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = ltriang[light,surf1,mm]]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = ltriang[light,surf2,mm]]; Print["Second patch done! "]; surf = Join[surf1, surf2]; xsurf1 = rsubdiv4[sxnet1, mm, n, debug]; If[light === 1, xsurf1 = trianglight[xsurf1,mm], xsurf1 = ltriang[light,xsurf1,mm]]; surf = Join[surf, xsurf1]; Print["Symmetric (x -> -x) of first patch done! "]; xsurf2 = rsubdiv4[sxnet2, mm, n, debug]; If[light === 1, xsurf2 = trianglight[xsurf2,mm], xsurf2 = ltriang[light,xsurf2,mm]]; surf = Join[surf, xsurf2]; Print["Symmetric (x -> -x) of second patch done! "]; Print["Ready to display! "]; If[fcnet === 1, Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction], Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction] ]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, running it twice *) (* as a square patch, using the reference triangles *) (* reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, *) (* {a, b, 1 - a - b}}, and *) (* reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, *) (* {c, d, 1 - c - d}} *) (* This version uses the input net, and also computes *) (* the symmetric net wrt the plane y = 0, and renders the whole *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* This is wrt the square [c, a] x [d, b] *) (* This versions lists c, a, d, b, in a more intuitive order! *) rendersymy[{net__}, mm_, c_, a_, d_, b_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, oldnet, sxnet1, sxnet2, xsurf1, xsurf2}, reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, {c, d, 1 - c - d}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; sxnet1 = symmy[net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; sxnet2 = symmy[net2]; surf1 = rsubdiv4[net1, mm, n, debug]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = ltriang[light,surf1,mm]]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = ltriang[light,surf2,mm]]; Print["Second patch done! "]; surf = Join[surf1, surf2]; xsurf1 = rsubdiv4[sxnet1, mm, n, debug]; If[light === 1, xsurf1 = trianglight[xsurf1,mm], xsurf1 = ltriang[light,xsurf1,mm]]; surf = Join[surf, xsurf1]; Print["Symmetric (y -> -y) of first patch done! "]; xsurf2 = rsubdiv4[sxnet2, mm, n, debug]; If[light === 1, xsurf2 = trianglight[xsurf2,mm], xsurf2 = ltriang[light,xsurf2,mm]]; surf = Join[surf, xsurf2]; Print["Symmetric (y -> -y) of second patch done! "]; Print["Ready to display! "]; If[fcnet === 1, Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction], Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction] ]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, running it twice *) (* as a square patch, using the reference triangles *) (* reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, *) (* {a, b, 1 - a - b}}, and *) (* reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, *) (* {c, d, 1 - c - d}} *) (* This version uses the input net, and also computes *) (* the symmetric net wrt the plane z = 0, and renders the whole *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* This is wrt the square [c, a] x [d, b] *) (* This versions lists c, a, d, b, in a more intuitive order! *) rendersymz[{net__}, mm_, c_, a_, d_, b_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, oldnet, sxnet1, sxnet2, xsurf1, xsurf2}, reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, {c, d, 1 - c - d}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; sxnet1 = symmz[net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; sxnet2 = symmz[net2]; surf1 = rsubdiv4[net1, mm, n, debug]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = ltriang[light,surf1,mm]]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = ltriang[light,surf2,mm]]; Print["Second patch done! "]; surf = Join[surf1, surf2]; xsurf1 = rsubdiv4[sxnet1, mm, n, debug]; If[light === 1, xsurf1 = trianglight[xsurf1,mm], xsurf1 = ltriang[light,xsurf1,mm]]; surf = Join[surf, xsurf1]; Print["Symmetric (z -> -z) of first patch done! "]; xsurf2 = rsubdiv4[sxnet2, mm, n, debug]; If[light === 1, xsurf2 = trianglight[xsurf2,mm], xsurf2 = ltriang[light,xsurf2,mm]]; surf = Join[surf, xsurf2]; Print["Symmetric (z -> -z) of second patch done! "]; Print["Ready to display! "]; If[fcnet === 1, Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction], Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction] ]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, running it twice *) (* as a square patch, using the reference triangles *) (* reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, *) (* {a, b, 1 - a - b}}, and *) (* reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, *) (* {c, d, 1 - c - d}} *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* This is wrt the square [c, a] x [d, b] *) sqrender[{net__}, mm_, a_, b_, c_, d_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, oldnet}, reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, {c, d, 1 - c - d}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; surf1 = rsubdiv4[net1, mm, n, debug]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = ltriang[light,surf1,mm]]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = ltriang[light,surf2,mm]]; Print["Second patch done! "]; surf = Join[surf1, surf2]; Print["Ready to display! "]; If[fcnet === 1, Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction], Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction] ]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, running it twice *) (* as a square patch, using the reference triangles *) (* reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, *) (* {a, b, 1 - a - b}}, and *) (* reftrig2 = {{a, d, 1 - a - d}, {c, d, 1 - c - d}, *) (* {a, b, 1 - a - b} (* Warning: this is NOT the same as reftrig2 in sqrender the last two points are swapped *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* This is wrt the square [c, a] x [d, b] *) revsqrender[{net__}, mm_, a_, b_, c_, d_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, oldnet}, reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; reftrig2 = {{a, d, 1 - a - d}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; surf1 = rsubdiv4[net1, mm, n, debug]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = ltriang[light,surf1,mm]]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = ltriang[light,surf2,mm]]; Print["Second patch done! "]; surf = Join[surf1, surf2]; Print["Ready to display! "]; If[fcnet === 1, Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction], Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction] ]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, fractal version, *) (* running it twice using the reference triangles *) (* reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}, and *) (* reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}} *) fractal2[{net__}, mm_, n_, lx_,mx_,ly_,my_,lz_,mz_,light_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, oldnet}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; surf1 = rfracdiv4[net1, mm, n]; surf1 = fractriang[light,surf1,mm]; Print["First patch done! "]; surf2 = rfracdiv4[net2, mm, n]; surf2 = fractriang[light,surf2,mm]; Print["Second patch done! "]; surf = Join[surf1, surf2]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, running it twice *) (* as a square patch, using the reference triangles *) (* reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, *) (* {a, b, 1 - a - b}}, and *) (* reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, *) (* {c, d, 1 - c - d}} *) (* if light = 1, uses light triangulation trianglight, *) (* else fractriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* This is wrt the square [c, a] x [d, b] *) (* Warning: the frames are NOT those used in rendersq *) sqfractal[{net__}, mm_, a_, b_, c_, d_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, image, stop, oldnet}, reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, {c, d, 1 - c - d}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; surf1 = rfracdiv4[net1, mm, n]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = fractriang[light,surf1,mm]]; Print["First patch done! "]; surf2 = rfracdiv4[net2, mm, n]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = fractriang[light,surf2,mm]]; Print["Second patch done! "]; surf = Join[surf1, surf2]; Print["Ready to display! "]; If[fcnet === 1, Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, lightkind[light], AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction], Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, lightkind[light], AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction] ]; ]; (* To display a single call to the de Casteljau algorithm *) (* returns 3 control nets *) (* if light = 1, uses light triangulation trianglight, *) (* else fractriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) decasone[{net__}, mm_, {a__}, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_, debug_] := Block[ {net2, surf, image, stop, oldnet, pt = {a}}, stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net2 = maptohat[{net}]; oldnet = net2; Print["Input Net: ", net2]; surf = sdecas3[net2, mm, pt]; If[light === 1, surf = trianglight[surf,mm], surf = fractriang[light,surf,mm]]; Print["Ready to display! "]; If[fcnet === 1, image = {surf, controlnet3D[{net},mm]}, image = surf ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* To display a single call to the de Casteljau algorithm *) (* returns the layers of the pyramid *) (* if light = 1, uses light triangulation trianglight, *) (* else fractriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) decaspyr[{net__}, mm_, {a__}, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_, debug_] := Block[ {net2, surf1, pnet, lnet, fpoint, trig, image, stop, oldnet, pt = {a}}, stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net2 = maptohat[{net}]; oldnet = net2; Print["Input Net: ", net2]; surf1 = decas3a[net2, mm, pt]; fpoint = Last[surf1]; surf1 = Drop[surf1,-1]; pnet = Last[surf1]; lnet = Join[pnet,fpoint]; (* Print["lnet = ",lnet]; *) lnet = projnet[lnet]; trig = {Line[{lnet[[1]],lnet[[4]]}], Line[{lnet[[2]],lnet[[4]]}], Line[{lnet[[3]],lnet[[4]]}]}; (* Print["trig = ",trig]; *) If[light === 1, surf1 = trianglight[surf1,mm], surf1 = ltriang[light,surf1,mm]]; surf1 = Append[surf1, {RGBColor[1, 0, 0], trig}]; (* Print["surf1 3 =", surf1]; *) Print["Ready to display! "]; If[fcnet === 1, image = {surf1, controlnet3D[{net},mm]}, image = surf1 ]; Show[ Graphics3D[{Thickness[0.008], PointSize[0.1], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, first using a new *) (* reference triangle *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) trender[{net__}, mm_, {newtrig__}, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_, debug_] := Block[ {net2, surf, image, stop, oldnet}, stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net2 = maptohat[{net}]; oldnet = net2; net2 = newcnet3[net2,mm,{newtrig}]; Print["New Net: ", net2]; surf = rsubdiv4[net2, mm, n, debug]; If[light === 1, surf = trianglight[surf,mm], surf = ltriang[light,surf,mm]]; Print["Ready to display! "]; If[fcnet === 1, image = {surf, controlnet3D[{net},mm]}, image = surf ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using rfracdiv4 *) (* i.e. splitting into 4 patches recursively, first using a new *) (* reference triangle *) (* if light = 1, uses light triangulation trianglight, *) (* else fractriang *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* fractal version *) tfractal[{net__}, mm_, {newtrig__}, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_] := Block[ {net2, surf, image, stop, oldnet}, stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net2 = maptohat[{net}]; oldnet = net2; net2 = newcnet3[net2,mm,{newtrig}]; Print["New Net: ", net2]; surf = rfracdiv4[net2, mm, n]; If[light === 1, surf = trianglight[surf,mm], surf = fractriang[light,surf,mm]]; Print["Ready to display! "]; If[fcnet === 1, image = {surf, controlnet3D[{net},mm]}, image = surf ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* checks for collinearity of three points in the plane *) collin[a__, b__, c__] := Block[ {a1, a2, b1, b2, c1, c2, d1, d2, res}, a1 = a[[1]]; a2 = a[[2]]; b1 = b[[1]]; b2 = b[[2]]; c1 = c[[1]]; c2 = c[[2]]; d1 = (c1 - a1)*(b2 - a2); d2 = (c2 - a2)*(b1 - a1); If[d1 === d2, res = 1, res = 0]; res ]; (* checks to see whether a vector is almost a zero vector *) nearzero[a_] := Block[ {stop, i, d, l}, l = Length[a]; i = 1; stop = 1; While[stop === 1 && i <= l, d = a[[i]]; If[Abs[d] < 10^(-16), i = i + 1, stop = 0]; ]; If[stop === 1, Print["*** Warning, Near Zero Vector: ", a, " ***"]]; stop ]; (* computes a control net for a curve on the surface, *) (* given two points a = (a1, a2, a3) and b = (b1, b2, b3) *) (* in the parameter plane, wrt the reference triangle (r, s, t) *) (* We use newcnet3 to get the control points of the *) (* curve segment defined over [0, 1], with end points a, b. *) curvcpoly[{net__}, m_, a_, b_] := Block[ {cc = {net}, res, net1, trigr, trigs, trigt, i, r, s, t, aa, bb}, (r = {1, 0, 0}; s = {0, 1, 0}; t = {0, 0, 1}; aa = a; bb = b; trigr = {r, bb, aa}; trigs = {s, bb, aa}; trigt = {t, bb, aa}; If[collin[aa,bb,r] === 0, net1 = newcnet3[cc, m, trigr], If[collin[aa,t,r] === 1, net1 = newcnet3[cc, m, trigs], net1 = newcnet3[cc, m, trigt] ] ]; res = {}; Do[ res = Append[res, net1[[i]]], {i, 1, m + 1} ]; res ) ]; (* To render a closed curve on a surface, specified by a control *) (* polygon computed by curvcpoly from two points on the surface *) (* specified by points a, b, in the parameter plane. *) clocurv3D[{net__},mm_,n_,a_,b_,lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {cpoly = {net}, l, m, pol, curv1, curv2, pol1, pol2, r, s, stop, image}, stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; r = -1; s = 1; l = -1; m = 1; pol = maptohat[cpoly]; pol1 = curvcpoly[pol, mm, a, b]; Print["New Cont. Poly.: ", pol1]; pol2 = negpoly[pol1]; Print["Dual Cont. Poly.: ", pol2]; curv1 = subdiv[pol1, l, m, n]; Print["First curve segment done! "]; If[debug === 20, Print["First curve: ", curv1] ]; curv1 = tolineseg[curv1]; curv2 = subdiv[pol2, l, m, n]; If[debug === 20, Print["Second curve: ", curv2] ]; curv2 = tolineseg[curv2]; Print["Second curve segment done! "]; Print["Ready to display! "]; image = {{RGBColor[0,1,0], curv1}, {RGBColor[1,0,0], curv2}}; Show[Graphics3D[{Thickness[0.001], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction] ]; (* Functions needed to render a closed surface by closed curves *) (* corresponding to the u-curves and the v-curves *) uvcurves[{net__},mm_,n_,ni_,debug_] := Block[ {pol = {net}, l, m, curv1, curv2, pol1, pol2, r, s, a, b, x, y, i, surf1, surf2, num, res}, r = -1; s = 1; l = -1; m = 1; surf1 = {}; surf2 = {}; If[ni === 0, num = 1, num = ni]; Do[ x = -1 + i/num; a = {x, -1, 2 - x}; b = {x, 1, -x}; pol1 = curvcpoly[pol, mm, a, b]; pol2 = negpoly[pol1]; curv1 = subdiv[pol1, l, m, n]; If[debug === 20, Print["First curve: ", curv1] ]; curv1 = tolineseg[curv1]; surf1 = Join[surf1,curv1]; curv2 = subdiv[pol2, l, m, n]; If[debug === 20, Print["Second curve: ", curv2] ]; curv2 = tolineseg[curv2]; surf1 = Join[surf1,curv2]; Print["u-Curve ",i," done! "], {i, 0, 2*ni} ]; Print["u-Curves done! "]; Do[ y = -1 + i/num; a = {-1, y, 2 - y}; b = {1, y, -y}; pol1 = curvcpoly[pol, mm, a, b]; pol2 = negpoly[pol1]; curv1 = subdiv[pol1, l, m, n]; If[debug === 20, Print["First curve: ", curv1] ]; curv1 = tolineseg[curv1]; surf2 = Join[surf2,curv1]; curv2 = subdiv[pol2, l, m, n]; If[debug === 20, Print["Second curve: ", curv2] ]; curv2 = tolineseg[curv2]; surf2 = Join[surf2,curv2]; Print["v-Curve ", i, " done! "], {i, 0, 2*ni} ]; Print["v-Curves done! "]; res = Join[{surf1}, {surf2}]; res ]; (* To render a closed surface by closed curves *) (* corresponding to the u-curves and the v-curves *) (* wrt to square [-1, 1] x [-1, 1] *) lrender[{net__},mm_,n_,ni_,lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {cpoly = {net}, pol, curv1, curv2, pol1, pol2, stop, image, a, b, x, y, i, surf1, surf2, surf}, stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; pol = maptohat[cpoly]; surf1 = {}; surf2 = {}; surf = uvcurves[pol,mm,n,ni,debug]; surf1 = surf[[1]]; surf2 = surf[[2]]; Print["Ready to display!"]; image = {{RGBColor[0,0,1], surf1}, {RGBColor[1,0,0], surf2}}; Show[Graphics3D[{Thickness[0.001], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction] ]; (* To render a closed surface by closed curves *) (* corresponding to the u-curves and the v-curves *) (* This version computes closed curves for the nets net0, theta1, rho1 *) lfrender[{net__},mm_,n_,ni_,lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {cpoly = {net}, stop, net0, net1, net2, net3, anet, theta1, theta2, rho1, rho2, oldnet, surf1, surf2, surf3, surf4, surf5, surf6, surfu, surfv, surf, image, reftrig1, reftrig2, reftrig3, pt, i, j, k}, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}; reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}; reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}; stop = testm[cpoly, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[cpoly]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; anet = convtomat[net2,mm]; theta1 = {}; Do[ Do[ pt = anet[[j + 1, mm - i - j + 1]]; If[OddQ[i + j], pt = -pt]; theta1 = Append[theta1, pt], {j, 0, mm - i} ], {i, 0, mm} ]; net3 = newcnet3[net0,mm,reftrig3]; Print["New Net 3: ", net3]; Print["New Net theta1: ", theta1]; anet = convtomat[net3,mm]; rho1 = {}; Do[ Do[ pt = anet[[j + 1, mm - i - j + 1]]; If[OddQ[j], pt = -pt]; rho1 = Append[rho1, pt], {j, 0, mm - i} ], {i, 0, mm} ]; Print["New Net rho1: ", rho1]; surf = uvcurves[net0,mm,n,ni,debug]; surf1 = surf[[1]]; surf2 = surf[[2]]; Print["Curves for patch 1,2 done!"]; surf = uvcurves[theta1,mm,n,ni,debug]; surf3 = surf[[1]]; surf4 = surf[[2]]; Print["Curves for patch 3,4 done!"]; surf = uvcurves[rho1,mm,n,ni,debug]; surf5 = surf[[1]]; surf6 = surf[[2]]; Print["Curves for patch 5,6 done!"]; surfu = Join[surf1, surf3]; surfu = Join[surfu, surf5]; surfv = Join[surf2, surf4]; surfv = Join[surfv, surf6]; Print["Ready to display!"]; image = {{RGBColor[0,0,1], surfu}, {RGBColor[1,0,0], surfv}}; Show[Graphics3D[{Thickness[0.001], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction] ]; (* To render a curve on a surface, specified by a control *) (* polygon computed by curvcpoly from two points on the surface *) (* specified by points a, b, in the parameter plane. *) scurv3D[{net__},mm_,n_,a_,b_,lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {cpoly = {net}, pol, curv1, pol1, r, s, stop, image}, stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; r = 0; s = 1; pol = maptohat[cpoly]; pol1 = curvcpoly[pol, mm, a, b]; Print["New Cont. Poly.: ", pol1]; curv1 = subdiv[pol1, r, s, n]; Print["Curve segment done! "]; If[debug === 20, Print["Curve: ", curv1] ]; curv1 = tolineseg[curv1]; Print["Ready to display! "]; image = {RGBColor[1,0,0], curv1}; Show[Graphics3D[{Thickness[0.001], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction] ]; (* To render a closed curve on a surface, specified by a control *) (* polygon computed by curvcpoly, together with the surface, *) (* from two points on the surface *) (* specified by points pta, ptb, in the parameter plane. *) (* This is wrt the square [c, a] x [d, b] *) sqcurv[{net__}, mm_, a_, b_, c_, d_, n_, nn_, pta_, ptb_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_, debug_] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, l, m, r, s, pol1, pol2, curv1, curv2, np1, image1, image2, oldnet}, reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, {c, d, 1 - c - d}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; surf1 = rsubdiv4[net1, mm, n, debug]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = ltriang[light,surf1,mm]]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = ltriang[light,surf2,mm]]; Print["Second patch done! "]; surf = Join[surf1, surf2]; r = -1; s = 1; l = -1; m = 1; pol1 = curvcpoly[net0, mm, pta, ptb]; np1 = pol1; Print["New Cont. Poly.: ", np1]; pol2 = negpoly[np1]; Print["Dual Cont. Poly.: ", pol2]; curv1 = subdiv[pol1, l, m, nn]; If[debug === 20, Print["First curve: ", curv1] ]; curv1 = tolineseg[curv1]; Print["First curve segment done! "]; curv2 = subdiv[pol2, l, m, nn]; If[debug === 20, Print["Second curve: ", curv2] ]; curv2 = tolineseg[curv2]; Print["Second curve segment done! "]; Print["Ready to display! "]; If[fcnet === 1, image1 = {Thickness[0.0006], surf}, image1 = {Thickness[0.0006], controlnet3D[{net},mm], surf} ]; image2 = {{RGBColor[0,0,1], curv1}, {RGBColor[1,0,0], curv2}}; Show[Graphics3D[{image1, image2}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* To build a full rectangular grid of edges from a rectangular net of degree (p, q) *) recgrid[{cnet__},p_,q_] := Block[ {cc = {cnet}, i, j, lseg = {}, edge}, ( lseg = {}; (* first, links horizontal line segments *) Do[ Do[ edge= {cc[[i + 1, j + 1]], cc[[i + 1, j + 2]]}; lseg = Append[lseg, edge], {j, 0, q - 1} ], {i, 0, p} ]; (* Nest, links vertical line segments *) Do[ Do[ edge= {cc[[i + 1, j + 1]], cc[[i + 2, j + 1]]}; lseg = Append[lseg, edge], {i, 0, p - 1} ], {j, 0, q} ]; lseg ) ]; (* To build a tiling by solid squares from a rectangular net of degree (p, q) *) solrecgrid[{cnet__},p_,q_] := Block[ {cc = {cnet}, i, j, lseg = {}, square}, ( lseg = {}; Do[ Do[ square = {cc[[i, j]], cc[[i, j + 1]], cc[[i + 1, j + 1]], cc[[i + 1, j]]}; lseg = Append[lseg, square], {j, 1, q} ], {i, 1, p} ]; lseg ) ]; (* To build a rectangular grid only containing the boundary points *) (* from a rectangular net of degree (p, q) *) recgridlight[{cnet__},p_,q_] := Block[ {cc = {cnet}, i, j, lseg = {}, edge}, ( lseg = {}; (* first, links horizontal line segments, bottom line *) Do[ edge= {cc[[1, j + 1]], cc[[1, j + 2]]}; lseg = Append[lseg, edge], {j, 0, q - 1} ]; (* next, links horizontal line segments, top line *) Do[ edge= {cc[[p + 1, j + 1]], cc[[p + 1, j + 2]]}; lseg = Append[lseg, edge], {j, 0, q - 1} ]; (* next, links vertical line segments, left line *) Do[ edge= {cc[[i + 1, 1]], cc[[i + 2, 1]]}; lseg = Append[lseg, edge], {i, 0, p - 1} ]; (* next, links vertical line segments, right line *) Do[ edge= {cc[[i + 1, q + 1]], cc[[i + 2, q + 1]]}; lseg = Append[lseg, edge], {i, 0, p - 1} ]; lseg ) ]; (* Calls segltriang if switch === -1 (full triangulation by line segments) and segrecmesh, otherwise if switch === 4, triangulation by solid colored rectangles with 4 colors, else recmeshsolid, default lighting or 2 colors Only called when switch =!= 1 *) (* When switch === 3, do not create vertices *) (* When switch > 10, unique color, blue = 11, red = 12, green = 13, yellow = 14 *) recmesh[switch_,{cnet__},p_,q_] := Block[ {cc = {cnet}, res}, ( If[switch === -1, res = segrecmesh[cc,p,q], If[switch === 4, res = solrecmesh[cc,p,q], If[switch > 10, res = recmeshcolor[cc,p,q,switch], res = recmeshsolid[switch,cc,p,q] ] ] ]; res ) ]; (* To build a list of Opaque rectangles from a list of quadruples of nodes *) solbuildrecmesh[{cnet__}] := Block[ {cc = {cnet}, ll, square, i, lseg = {}, res}, (ll = Length[cc]; Do[ square = cc[[i]]; If[Length[square] === 4 , lseg = Append[lseg, Polygon[square]], lseg = Append[lseg, {RGBColor[1,0,0], Point[square]}] ], {i, 1, ll} ]; lseg ) ]; (* Builds a list of line segments for the union of lists of *) (* rectangular nets. This version calls recgrid *) (* Also adds the four corners of each rectangle *) (* uses different colors on four subpatches *) segrecmesh[{cnet__},p_,q_] := Block[ {cc = {cnet}, ll, l, i, j, lineseg = {}, linelis, newlist, edge, pt, res}, ( ll = Length[cc]; l = ll/4; Do[ Do[ edgelist= projnet[cc[[j]]]; newlist = makerecnet[edgelist, p, q]; linelis = recgrid[newlist, p, q]; edge = linelis[[1]]; pt = First[edge]; (* adds lower left corner *) linelis = Append[linelis,pt]; edge = linelis[[q]]; pt = Last[edge]; (* adds lower right corner *) linelis = Append[linelis,pt]; edge = linelis[[q*p + 1]]; pt = First[edge]; (* adds upper left corner *) linelis = Append[linelis,pt]; edge = linelis[[q*(p + 1)]]; pt = Last[edge]; (* adds upper right corner *) linelis = Append[linelis,pt]; linelis = buildmesh[linelis]; If[j === 1, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[j === 2, linelis = Prepend[linelis, RGBColor[1, 0.8, 0]]]; (* orange *) If[j === 3, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) If[j === 4, linelis = Prepend[linelis, RGBColor[0, 0.7, 0]]]; (* green *) lineseg = Join[lineseg, linelis], {j, 1, 4} ]; cc = Drop[cc, 4], {i, 1, l} ]; lineseg ) ]; (* Builds a list of solid squares for the union of lists of *) (* rectangular nets. This version calls solrecgrid *) (* Also adds the four corners of each rectangle *) recmeshsolid[switch_,{cnet__},p_,q_] := Block[ {cc = {cnet}, ll, i, lineseg = {}, linelis, newlist, edge, pt, res}, ( ll = Length[cc]; Do[ edgelist= projnet[cc[[i]]]; newlist = makerecnet[edgelist, p, q]; linelis = solrecgrid[newlist, p, q]; If[switch =!= 3, edge = linelis[[1]]; pt = First[edge]; (* adds lower left corner *) linelis = Append[linelis,pt] ]; If[switch =!= 3, edge = linelis[[q]]; pt = edge[[2]]; (* adds lower right corner *) linelis = Append[linelis,pt] ]; If[switch =!= 3, edge = linelis[[q*(p-1) + 1]]; pt = Last[edge]; (* adds upper left corner *) linelis = Append[linelis,pt] ]; If[switch =!= 3, edge = linelis[[q*p]]; pt = edge[[3]]; (* adds upper right corner *) linelis = Append[linelis,pt] ]; linelis = solbuildrecmesh[linelis]; lineseg = Join[lineseg, linelis], {i, 1, ll} ]; lineseg ) ]; (* Builds a list of solid squares for the union of lists of *) (* rectangular nets. This version calls solrecgrid *) (* Also adds the four corners of each rectangle *) (* uses different colors on four subpatches *) solrecmesh[{cnet__},p_,q_] := Block[ {cc = {cnet}, ll, l, i, j, lineseg = {}, linelis, newlist, edge, pt, res}, ( ll = Length[cc]; l = ll/4; Do[ Do[ edgelist= projnet[cc[[j]]]; newlist = makerecnet[edgelist, p, q]; linelis = solrecgrid[newlist, p, q]; edge = linelis[[1]]; pt = First[edge]; (* adds lower left corner *) linelis = Append[linelis,pt]; edge = linelis[[q]]; pt = edge[[2]]; (* adds lower right corner *) linelis = Append[linelis,pt]; edge = linelis[[q*(p-1) + 1]]; pt = Last[edge]; (* adds upper left corner *) linelis = Append[linelis,pt]; edge = linelis[[q*p]]; pt = edge[[3]]; (* adds upper right corner *) linelis = Append[linelis,pt]; linelis = solbuildrecmesh[linelis]; If[j === 1, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[j === 2, linelis = Prepend[linelis, RGBColor[0, 1, 0]]]; (* green *) If[j === 3, linelis = Prepend[linelis, RGBColor[1, 1, 0]]]; (* yellow *) If[j === 4, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) lineseg = Join[lineseg, linelis], {j, 1, 4} ]; cc = Drop[cc, 4], {i, 1, l} ]; lineseg ) ]; (* Builds a list of solid squares for the union of lists of *) (* rectangular nets. This version calls solrecgrid *) (* assigns a unique color to all subpatches *) (* Either blue = 11, red = 12, green = 13, or yellow = 14 *) recmeshcolor[{cnet__},p_,q_,col_] := Block[ {cc = {cnet}, ll, i, lineseg = {}, linelis, newlist, edge, pt, res, color}, ( color = col; ll = Length[cc]; Do[ edgelist= projnet[cc[[i]]]; newlist = makerecnet[edgelist, p, q]; linelis = solrecgrid[newlist, p, q]; linelis = solbuildrecmesh[linelis]; If[color === 11, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[color === 12, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) If[color === 13, linelis = Prepend[linelis, RGBColor[0, 1, 0]]]; (* green *) If[color === 14, linelis = Prepend[linelis, RGBColor[1, 1, 0]]]; (* yellow *) lineseg = Join[lineseg, linelis], {i, 1, ll} ]; lineseg ) ]; (* Builds a list of line segments for the union of lists of *) (* rectangular nets. This version calls recgridlight *) (* Also adds the four corners of each rectangle *) (* uses different colors on four subpatches *) recmeshlight[{cnet__},p_,q_] := Block[ {cc = {cnet}, ll, l, i, j, lineseg = {}, linelis, newlist, edge, pt}, ( ll = Length[cc]; l = ll/4; Do[ Do[ edgelist= projnet[cc[[j]]]; newlist = makerecnet[edgelist, p, q]; linelis = recgridlight[newlist, p, q]; edge = linelis[[1]]; pt = First[edge]; (* adds lower left corner *) linelis = Append[linelis,pt]; edge = linelis[[q]]; pt = Last[edge]; (* adds lower right corner *) linelis = Append[linelis,pt]; edge = linelis[[q + 1]]; pt = First[edge]; (* adds upper left corner *) linelis = Append[linelis,pt]; edge = linelis[[2*q]]; pt = Last[edge]; (* adds upper right corner *) linelis = Append[linelis,pt]; linelis = buildmesh[linelis]; If[j === 1, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[j === 2, linelis = Prepend[linelis, RGBColor[1, 0.8, 0]]]; (* orange *) If[j === 3, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) If[j === 4, linelis = Prepend[linelis, RGBColor[0, 0.7, 0]]]; (* green *) lineseg = Join[lineseg, linelis], {j, 1, 4} ]; cc = Drop[cc, 4], {i, 1, l} ]; lineseg ) ]; (* To make a rectangular net of degree (p, q) from a list of length (p + 1)x(q + 1) *) makerecnet[{net__}, p_, q_] := Block[ {oldnet = {net}, newnet = {}, row, pt, i, j, n}, n = Length[oldnet]; Do[row = {}; Do[ pt = oldnet[[i*(q + 1) + j + 1]]; row = Append[row,pt], {j, 0, q} ]; newnet = Append[newnet, row], {i, 0, p} ]; newnet ]; (* De Casteljau algorithm with subdivision into four rectangular nets *) (* for a rectangular net of degree (p, q) *) recdecas[{oldnet__}, p_, q_, r1_, s1_, r2_, s2_, u_, v_, debug_] := Block[ {net = {oldnet}, temp, newneta, newnetb, newnet1, newnet2, newnet}, newneta = {}; newnetb = {}; Print[" p = ", p, ", q = ", q, ", in recdecas"]; temp = vrecdecas[net, p, q, r2, s2, v, debug]; newneta = temp[[1]]; newnetb = temp[[2]]; newnet1 = urecdecas[newneta, p, q, r1, s1, u, debug]; newnet2 = urecdecas[newnetb, p, q, r1, s1, u, debug]; newnet = Join[newnet1, newnet2]; If[debug === 20, Print[" Subdivided nets: ", newnet]]; newnet ]; (* De Casteljau algorithm with subdivision into two rectangular nets *) (* for a rectangular net of degree (p, q). Subdivision along the u-curves *) urecdecas[{oldnet__}, p_, q_, r1_, s1_, u_, debug_] := Block[ {net = {oldnet}, i, j, temp, bistar, row1, row2, pt, newnet1, newnet2, newnet}, bistar = {}; row1 = {}; row2 = {}; newnet1 = {}; newnet2 = {}; Do[ bistar = {}; row1 = {}; row2 = {}; Do[ pt = net[[(q + 1)*i + j + 1]]; bistar = Append[bistar, pt], {i, 0, p} ]; If[debug === 20, Print[" bistar: ", bistar]]; temp = subdecas[bistar, r1, s1, u]; row1 = temp[[1]]; row2 = temp[[2]]; newnet1 = Join[newnet1, {row1}]; newnet2 = Join[newnet2, {row2}], {j, 0, q} ]; newnet1 = rectrans[newnet1]; newnet2 = rectrans[newnet2]; newnet = Join[{newnet1}, {newnet2}]; If[debug === 20, Print[" Subdivided nets: ", newnet]]; newnet ]; (* Converts a matrix given as a list of columns into a linear list of rows *) rectrans[{oldnet__}] := Block[ {net = {oldnet}, i, j, pp, qq, row, pt, newnet}, row = {}; newnet = {}; qq = Length[net]; pp = Length[net[[1]]]; Do[ Do[ row = net[[j]]; pt = row[[i]]; newnet = Append[newnet, pt], {j, 1, qq} ], {i, 1, pp} ]; newnet ]; (* De Casteljau algorithm with subdivision into two rectangular nets *) (* for a rectangular net of degree (p, q). Subdivision along the v-curves *) vrecdecas[{oldnet__}, p_, q_, r2_, s2_, v_, debug_] := Block[ {net = {oldnet}, i, j, temp, bstarj, row1, row2, pt, newneta, newnetb, newnet}, bistar = {}; row1 = {}; row2 = {}; newneta = {}; newnetb = {}; Do[ bstarj = {}; row1 = {}; row2 = {}; Do[ pt = net[[(q + 1)*i + j + 1]]; bstarj = Append[bstarj, pt], {j, 0, q} ]; If[debug === 20, Print[" bstarj: ", bstarj]]; temp = subdecas[bstarj, r2, s2, v]; row1 = temp[[1]]; row2 = temp[[2]]; newneta = Join[newneta, row1]; newnetb = Join[newnetb, row2], {i, 0, p} ]; newnet = Join[{newneta}, {newnetb}]; If[debug === 20, Print[" Subdivided nets: ", newnet]]; newnet ]; (* Computes a new rectangular net of degree (p, q) *) (* wrt to frames (a, b) and (c, d) on the affine line *) (* In the hat space *) recnewnet[{oldnet__}, p_, q_, a_, b_, c_, d_, debug_] := Block[ {net = {oldnet}, temp, newnet1, newnet2, newnet3, newnet4, rnet, r1, s1, r2, s2, ll, i}, r1 = 0; s1 = 1; r2 = 0; s2 = 1; newnet1 = {}; newnet2 = {}; newnet3 = {}; newnet4 = {}; If[d =!= r2, temp = vrecdecas[net, p, q, r2, s2, d, debug]; newnet1 = temp[[1]], newnet1 = net ]; If[debug === 20, Print[" newnet1: ", newnet1]]; If[b =!= r1, temp = urecdecas[newnet1, p, q, r1, s1, b, debug]; newnet2 = temp[[1]], newnet2 = newnet1 ]; If[debug === 20, Print[" newnet2: ", newnet2]]; If[d =!= r2, temp = vrecdecas[newnet2, p, q, r2, d, c, debug]; newnet3 = temp[[2]], temp = vrecdecas[newnet2, p, q, r2, s2, c, debug]; newnet3 = temp[[1]] ]; If[debug === 20, Print[" newnet3: ", newnet3]]; If[b =!= r1, temp = urecdecas[newnet3, p, q, r1, b, a, debug]; newnet4 = temp[[2]], temp = urecdecas[newnet3, p, q, r1, s1, a, debug]; (* needs to reverse the net in this case *) rnet = temp[[1]]; newnet4 = {}; ll = Length[rnet]; Do[ newnet4 = Prepend[newnet4, rnet[[i]]], {i, 1, ll} ] ]; If[debug === 20, Print[" newnet4: ", newnet4]]; newnet4 ]; (* Computes a new rectangular net of degree (p, q) *) (* wrt to frames (a, b) and (c, d) on the affine line *) (* down in the affine space *) newrecnet[{oldnet__}, p_, q_, a_, b_, c_, d_, debug_] := Block[ {net = {oldnet}, newnet1, newnet2, newnet}, newnet1 = maptohat[net]; newnet2 = recnewnet[newnet1, p, q, a, b, c, d, debug]; newnet = piomeg[newnet2]; newnet ]; (* performs a subdivision step on each rectangular net in a list *) (* using recdecas4, i.e., splitting into four subpatches *) recitersub4[{net__}, p_, q_, r1_, s1_, r2_, s2_, debug_] := Block[ {cnet = {net}, lnet = {}, l, i, u, v}, (l = Length[cnet]; u = 1/2; v = 1/2; Do[ lnet = Join[lnet, recdecas[cnet[[i]], p, q, r1, s1, r2, s2, u, v, debug]] , {i, 1, l} ]; lnet ) ]; (* performs n subdivision steps using recitersub4, i.e., splits *) (* a rectangular patch into 4 subpatches *) recsubdiv4[{net__}, p_, q_, n_, debug_] := Block[ {newnet = {net}, i, r1, s1, r2, s2}, (r1 = 0; s1 = 1; r2 = 0; s2 = 1; newnet = {newnet}; Do[ newnet = recitersub4[newnet, p, q, r1, s1, r2, s2, debug], {i, 1, n} ]; newnet ) ]; (* prepares rectangular control net *) recontrolnet3D[{net__}, p_, q_] := Block[ {res, net2, newlist, linelis}, net2 = strip[{net}]; newlist = makerecnet[net2, p, q]; linelis = recgrid[newlist, p, q]; res = buildmesh[linelis]; res = Prepend[res, RGBColor[0,1,0]]; res ]; (* To display the result of n subdivision steps using recsubdiv4 *) (* i.e. splitting the reference rectangle into 4 subrectangles *) (* if fcnet = 1, displays control net, otherwise skip it *) (* if light = 1, uses light rectangular grid recmeshlight, *) (* else recmesh *) recrender[{net__}, p_, q_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0 = {net}, net2, surf0, surf, image, stop}, stop = testpq[net0, p, q]; If[stop === 1, Return["*** Unable to run ***"]]; Print[" Input rectangular net: ", net0]; net2 = maptohat[net0]; surf0 = recsubdiv4[net2, p, q, n, debug]; If[debug === 20, Print["surf = ", surf0]]; If[light === 1, surf = recmeshlight[surf0,p,q], surf = recmesh[light,surf0,p,q]]; If[fcnet === 1, image = {surf, recontrolnet3D[net0, p, q]}, image = surf ]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n subdivision steps using recsubdiv4 *) (* i.e. splitting the reference rectangle into 4 subrectangles *) (* if fcnet = 1, displays control net, otherwise skip it *) (* if light = 1, uses light rectangular grid recmeshlight, *) (* else recmesh *) (* This version first recomputes a net wrt [a, b] X [c, d] *) (* when debug = 99, skip maptohat *) recsqrender[{net__}, p_, q_, n_, a_, b_, c_, d_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0 = {net}, net1, net2, net3, surf0, surf, image, stop}, stop = testpq[net0, p, q]; If[stop === 1, Return["*** Unable to run ***"]]; Print[" Input rectangular net: ", net0]; If[debug === 99, net1 = net0, net1 = maptohat[net0]]; net2 = recnewnet[net1, p, q, a, b, c, d, debug]; Print["New net = ", net2]; surf0 = recsubdiv4[net2, p, q, n, debug]; If[debug === 20, Print["surf = ", surf0]]; If[light === 1, surf = recmeshlight[surf0,p,q], surf = recmesh[light,surf0,p,q]]; If[fcnet === 1, image = {surf, recontrolnet3D[net0, p, q]}, image = surf ]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* computes the symmetric of a net under x -> -x, y -> -y *) symmxy[{cnet__}] := Block[ {net = {cnet}, l, j, newnet, pt, newpt}, ( newnet = {}; l = Length[net]; Do[ pt = net[[j]]; newpt = {-pt[[1]], -pt[[2]], pt[[3]], pt[[4]]}; newnet = Join[newnet, {newpt}], {j, 1, l} ]; newnet ) ]; (* computes the symmetric of a net under x -> -x (wrt to plane x = 0 *) symmx[{cnet__}] := Block[ {net = {cnet}, l, j, newnet, pt, newpt}, ( newnet = {}; l = Length[net]; Do[ pt = net[[j]]; newpt = {-pt[[1]], pt[[2]], pt[[3]], pt[[4]]}; newnet = Join[newnet, {newpt}], {j, 1, l} ]; newnet ) ]; (* computes the symmetric of a net under y -> -y (wrt to plane y = 0 *) symmy[{cnet__}] := Block[ {net = {cnet}, l, j, newnet, pt, newpt}, ( newnet = {}; l = Length[net]; Do[ pt = net[[j]]; newpt = {pt[[1]], -pt[[2]], pt[[3]], pt[[4]]}; newnet = Join[newnet, {newpt}], {j, 1, l} ]; newnet ) ]; (* computes the symmetric of a net under z -> -z (wrt to plane z = 0 *) symmz[{cnet__}] := Block[ {net = {cnet}, l, j, newnet, pt, newpt}, ( newnet = {}; l = Length[net]; Do[ pt = net[[j]]; newpt = {pt[[1]], pt[[2]], -pt[[3]], pt[[4]]}; newnet = Join[newnet, {newpt}], {j, 1, l} ]; newnet ) ]; (* 4D version !!! Use ltriang4ss *) (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, running it twice *) (* as a square patch, using the reference triangles *) (* reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, *) (* {a, b, 1 - a - b}}, and *) (* reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, *) (* {c, d, 1 - c - d}} *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* This is wrt the square [c, a] x [d, b] *) (* This versions lists c, a, d, b, in a more intuitive order! *) render4sq[{net__}, mm_, c_, a_, d_, b_, n_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_,{aaa__},{bbb__}] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, oldnet, newbase, pa = {aaa}, pb = {bbb}, H}, reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, {c, d, 1 - c - d}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; H = pb - pa; Print[" H = ", H]; newbase = basis4D[H]; Print[" newbase = ", newbase]; surf1 = rsubdiv4[net1, mm, n, debug]; surf1 = ltriang4ss[newbase,pa,pb,surf1,mm]; Print["First patch done! "]; surf2 = rsubdiv4[net2, mm, n, debug]; surf2 = ltriang4ss[newbase,pa,pb,surf2,mm]; Print["Second patch done! "]; surf = Join[surf1, surf2]; Print["Ready to display! "]; If[fcnet === 1, Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction], Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction] ]; ]; (* for 4D version *) ltriang4ss[{newb__},{aaa__},{bbb__},{cnet__},oldm_] := Block[ {cc = {cnet}, ll, l, i, j, lineseg = {}, edglist, linelis, m, pa = {aaa}, pb = {bbb}, edge, pt, nnb = {newb}}, ( ll = Length[cc]; l = ll/4; Do[ Do[ m = netsize[cc[[j]]]; edgelist= projnet4D[pa,pb,cc[[j]],nnb]; linelis = soltriangul[edgelist,m]; edge = Last[linelis]; pt = edge[[2]]; (* top corner *) linelis = Append[linelis,pt]; edge = First[linelis]; pt = First[edge]; (* lower left corner *) linelis = Append[linelis,pt]; edge = linelis[[m]]; pt = Last[edge]; (* lower right corner *) linelis = Append[linelis,pt]; linelis = solbuildmesh[linelis]; If[j === 1, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[j === 2, linelis = Prepend[linelis, RGBColor[1, 1, 0]]]; (* yellow *) If[j === 3, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) If[j === 4, linelis = Prepend[linelis, RGBColor[0, 1, 0]]]; (* green *) lineseg = Join[lineseg, linelis], {j, 1, 4} ]; cc = Drop[cc, 4], {i, 1, l} ]; lineseg ) ]; (* 4D version !!! Use ltriang5D *) (* performs the change of basis before subdividing *) (* Using quaternion multiplication ! *) (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, running it twice *) (* as a square patch, using the reference triangles *) (* reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, *) (* {a, b, 1 - a - b}}, and *) (* reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, *) (* {c, d, 1 - c - d}} *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* This is wrt the square [c, a] x [d, b] *) (* This versions lists c, a, d, b, in a more intuitive order! *) (* For central projection of center a, aflag = 10, *) (* Otherwise parallel projection using vector bbb - aaa *) (* For two patches, one blue one red, use light = 2 (<> 4 and <= 10) *) (* For 4 colors, use light = 4 *) (* For one solid color, use light = 11, 12, 13 (blue, red, green) *) render5qsq[{net__}, mm_, c_, a_, d_, b_, n_, light_, fcnet_, aflag_, lx_,mx_,ly_,my_,lz_,mz_,debug_,{aaa__},{bbb__}] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, oldnet, newbase, pa = {aaa}, pb = {bbb}, h, nh, H, HC, movnet1, movnet2, flag, qa, qb, qc, qd, N1, norm1}, reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, {c, d, 1 - c - d}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; flag = aflag; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; h = pb - pa; qa = h[[1]]; qb = h[[2]]; qc = h[[3]]; qd = h[[4]]; nh = {qd, qc, -qb, qa}; N1 = qa^2 + qb^2 + qc^2 + qd^2; norm1 = N[Sqrt[N1]]; nh = nh/norm1; movnet1 = quatmotion[pa, nh, net1]; Print["Moved Net 1: ", movnet1]; H = mulquat[nh, h]; Print["H = ", H]; surf1 = rsubdiv4[movnet1, mm, n, debug]; If[light === 4 || light > 10, surf1 = ltriang5D[H,surf1,mm,flag,light], surf1 = ltriang5D[H,surf1,mm,flag,11] ]; Print["First patch done! "]; movnet2 = quatmotion[pa, nh, net2]; Print["Moved Net 2: ", movnet2]; surf2 = rsubdiv4[movnet2, mm, n, debug]; If[light === 4 || light > 10, surf2 = ltriang5D[H,surf2,mm,flag,light], surf2 = ltriang5D[H,surf2,mm,flag,12] ]; Print["Second patch done! "]; surf = Join[surf1, surf2]; Print["Ready to display! "]; If[fcnet === 1, Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction], Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction] ]; ]; (* 4D version !!! Use ltriang5D *) (* performs the change of basis before subdividing *) (* To display the result of n subdivision steps using rsubdiv4 *) (* i.e. splitting into 4 patches recursively, running it twice *) (* as a square patch, using the reference triangles *) (* reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, *) (* {a, b, 1 - a - b}}, and *) (* reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, *) (* {c, d, 1 - c - d}} *) (* if fcnet = 1, displays control net, otherwise skip it *) (* *) (* This is wrt the square [c, a] x [d, b] *) (* This versions lists c, a, d, b, in a more intuitive order! *) (* For central projection of center a, aflag = 10, *) (* Otherwise parallel projection using vector bbb - aaa *) (* For two patches, one blue one red, use light = 2 (<> 4 and <= 10) *) (* For 4 colors, use light = 4 *) (* For one solid color, use light = 11, 12, 13 (blue, red, green) *) render5sq[{net__}, mm_, c_, a_, d_, b_, n_, light_, fcnet_, aflag_, lx_,mx_,ly_,my_,lz_,mz_,debug_,{aaa__},{bbb__}] := Block[ {net0, net1, net2, surf1, surf2, surf, reftrig1, reftrig2, stop, oldnet, newbase, pa = {aaa}, pb = {bbb}, H, h, movnet1, movnet2, TP, flag}, reftrig1 = {{c, b, 1 - b - c}, {c, d, 1 - c - d}, {a, b, 1 - a - b}}; reftrig2 = {{a, d, 1 - a - d}, {a, b, 1 - a - b}, {c, d, 1 - c - d}}; stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; flag = aflag; net1 = newcnet3[net0,mm,reftrig1]; Print["New Net 1: ", net1]; net2 = newcnet3[net0,mm,reftrig2]; Print["New Net 2: ", net2]; TP = makemotion[pa,pb]; movnet1 = rigmotion[pa, net1, TP]; Print["Moved Net 1: ", movnet1]; h = pb - pa; H = makeH[h, TP]; Print["H = ", H]; surf1 = rsubdiv4[movnet1, mm, n, debug]; If[light === 4 || light > 10, surf1 = ltriang5D[H,surf1,mm,flag,light], surf1 = ltriang5D[H,surf1,mm,flag,11] ]; Print["First patch done! "]; movnet2 = rigmotion[pa, net2, TP]; Print["Moved Net 2: ", movnet2]; surf2 = rsubdiv4[movnet2, mm, n, debug]; If[light === 4 || light > 10, surf2 = ltriang5D[H,surf2,mm,flag,light], surf2 = ltriang5D[H,surf2,mm,flag,12] ]; Print["Second patch done! "]; surf = Join[surf1, surf2]; Print["Ready to display! "]; If[fcnet === 1, Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf, controlnet3D[{net},mm]}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction], Show[ Graphics3D[{Thickness[0.0006], facecolors[light], surf}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction] ]; ]; (* for 4D version with projnet5D *) ltriang5ss[{bbb__},{cnet__},oldm_,aflag_] := Block[ {cc = {cnet}, ll, l, i, j, lineseg = {}, edglist, linelis, m, pb = {bbb}, edge, pt, flag}, ( ll = Length[cc]; l = ll/4; flag = aflag; Do[ Do[ m = netsize[cc[[j]]]; edgelist= projnet5D[pb,cc[[j]],flag]; linelis = soltriangul[edgelist,m]; edge = Last[linelis]; pt = edge[[2]]; (* top corner *) linelis = Append[linelis,pt]; edge = First[linelis]; pt = First[edge]; (* lower left corner *) linelis = Append[linelis,pt]; edge = linelis[[m]]; pt = Last[edge]; (* lower right corner *) linelis = Append[linelis,pt]; linelis = solbuildmesh[linelis]; If[j === 1, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[j === 2, linelis = Prepend[linelis, RGBColor[1, 1, 0]]]; (* yellow *) If[j === 3, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) If[j === 4, linelis = Prepend[linelis, RGBColor[0, 1, 0]]]; (* green *) lineseg = Join[lineseg, linelis], {j, 1, 4} ]; cc = Drop[cc, 4], {i, 1, l} ]; lineseg ) ]; (* Builds a list of solid triangles for the union of lists of nets *) (* This version calls soltriangul *) (* Also adds the three corners of *) (* each pyramidal net as list elements *) (* assigns a unique color to all subpatches *) (* Either blue = 11, red = 12, or green = 13 *) (* 4D version works with projnet5D *) ltriangcolor5[{bbb__},{cnet__},oldm_,aflag_,col_] := Block[ {cc = {cnet}, ll, b = {bbb}, i, lineseg = {}, edglist, linelis, m, color}, ( color = col; ll = Length[cc]; Do[ m = netsize[cc[[i]]]; edgelist= projnet5D[b,cc[[i]],aflag]; linelis = soltriangul[edgelist,m]; linelis = solbuildmesh[linelis]; If[color === 11, linelis = Prepend[linelis, RGBColor[0, 0, 1]]]; (* blue *) If[color === 12, linelis = Prepend[linelis, RGBColor[1, 0, 0]]]; (* red *) If[color === 13, linelis = Prepend[linelis, RGBColor[0, 1, 0]]]; (* green *) lineseg = Join[lineseg, linelis], {i, 1, ll} ]; lineseg ) ]; (* Calls ltriang5ss if col === 4 otherwise ltriangcolor5 *) (* When col > 10, unique color, blue = 11, red = 12, green = 13 *) ltriang5D[{bbb__},{cnet__},oldm_,aflag_,col_] := Block[ {cc = {cnet}, res, b = {bbb}}, ( If[col === 4, res = ltriang5ss[b,cc,oldm,aflag], res = ltriangcolor5[b,cc,oldm,aflag,col] ]; res ) ]; (* To project a net in the hat space back onto the affine space With projection on a hyperplane passing through b and orthogonal to b - a, 4D version, *) projnet4D[{a__},{b__},{net__},{newbase__}] := Block[ {snet = {net}, newnet = {}, A = {a}, B = {b}, pt, h, j, l2, X, Y, TP = {newbase}, xx, yy}, l2 = Length[snet]; Do[ pt = snet[[j]]; h = projpt[pt]; X = centproj[A,B,h]; xx = X - A; yy = TP . xx; (* Print["yy = ", yy]; *) Y = Drop[yy,-1]; newnet = Append[newnet,Y], {j, 1, l2} ]; newnet ]; (* To project a net in the hat space back onto the affine space With projection on a hyperplane passing through b and orthogonal to b - 0, 4D version, *) (* central projection of center a if aflag = 10, *) (* otherwise, parallel projection using vector h *) (* This version is designed to work with basis5D *) projnet5D[{h__}, {net__}, aflag_] := Block[ {snet = {net}, newnet = {}, H = {h}, w, pt, npt, j, l2, Y, YY, flag}, l2 = Length[snet]; flag = aflag; Do[ pt = snet[[j]]; npt = projpt[pt]; If[flag === 10, w = npt[[4]]; (* Print["w = ", w]; *) If[w === 0, Print["w = 0 in projnet5; npt = ", npt]; w = 10^(-6)]; YY = (H[[4]]*npt)/w, YY = npt ]; Y = Drop[YY,-1]; newnet = Append[newnet,Y], {j, 1, l2} ]; (* Print[" newnet in projnet5D = ",newnet]; *) newnet ]; (* To apply the rigid motion defined by newbase to a vector *) (* vv a 5D vector *) (* This version is designed to work with basis5D *) rmotionv[{vv__},{newbase__}] := Block[ {v = {vv}, TP = {newbase}, yy}, yy = TP . v; (* Print["yy = ", yy]; *) yy ]; (* To apply the rigid motion defined by newbase and the new origin a, to a control net *) (* This version is designed to work with basis5D *) rigmotion[{a__},{net__},{newbase__}] := Block[ {snet = {net}, newnet = {}, A = {a}, pt, pta, j, l2, Y, AA, w, matrix = {newbase}}, l2 = Length[snet]; Do[ pt = snet[[j]]; AA = Append[A,0]; w = Last[pt]; If[w =!= 0, pta = pt - AA*w, pta = pt - AA]; Y = rmotionv[pta,matrix]; (* Print["Y = ", Y]; *) newnet = Append[newnet,Y], {j, 1, l2} ]; newnet ]; (* To apply the rigid motion defined by newbase to a vector *) (* vv a 5D vector *) (* This version is designed to work with basis5D *) (* and uses quaternion multiplication *) qmotion[vv_, h_] := Block[ {v = vv, H = h, w, pt, newpt, yy}, w = Last[v]; pt = Drop[v, -1]; newpt = mulquat[H, pt]; yy = Append[newpt, w]; (* Print["yy = ", yy]; *) yy ]; (* To apply the rigid motion defined by newbase and the new origin a, to a control net *) (* This version is designed to work with basis5D *) (* and uses quaternion multiplication *) quatmotion[a_, h_, {net__}] := Block[ {snet = {net}, newnet = {}, H = h, pt, pta, j, l2, Y, A = a, AA, w, matrix = {newbase}}, l2 = Length[snet]; Do[ pt = snet[[j]]; AA = Append[A,0]; w = Last[pt]; If[w =!= 0, pta = pt - AA*w, pta = pt - AA]; Y = qmotion[pta,H]; (* Print["Y = ", Y]; *) newnet = Append[newnet,Y], {j, 1, l2} ]; newnet ]; (* To make the rigid motion defined by the new origin a and the hyperplane passing trough b and orthogonal to b - a *) makemotion[{a__},{b__}] := Block[ {A = {a}, B = {b}, H, TP}, H = B - A; TP = basis5D[H]; Print[" newbase = ", TP]; TP ]; (* To make the new vector H after applying the rigid motion specified by a and matrix *) makeH[{h__},{matrix__}] := Block[ {H = {h}, Ha, TP = {matrix}, nBB}, Ha = Append[H, 0]; nBB = rmotionv[Ha,TP]; nBB = Drop[nBB, -1]; nBB ]; (* Given a vector h = (a, b, c, d), finds an orthonormal basis *) (* using the vectors involved in the matrix representation *) (* of the quaternion (a, b, c, d). The transpose of *) (* the change of basis matrix has these vectors as its rows. *) basis4D[{h__}] := Block[ {H = {h}, a, b, c, d, res = {}, N1, norm1, e1, e2, e3, e4, e5}, a = H[[1]]; b = H[[2]]; c = H[[3]]; d = H[[4]]; N1 = a^2 + b^2 + c^2 + d^2; norm1 = N[Sqrt[N1]]; e1 = {-b, a, d, -c}; e1 = e1/norm1; (* Print["e1 = ", e1]; *) e2 = {-c, -d, a, b}; e2 = e2/norm1; (* Print["e2 = ", e2]; *) e3 = {-d, c, -b, a}; e3 = e3/norm1; (* Print["e3 = ", e3]; *) e4 = {a, b, c, d}; e4 = e4/norm1; (* Print["e4 = ", e4]; *) res = {e1, e2, e3, e4}; res ]; (* Given a vector h = (a, b, c, d), finds an orthonormal basis *) (* using the vectors involved in the matrix representation *) (* of the quaternion (a, b, c, d). The transpose of *) (* the change of basis matrix has these vectors as its rows. *) (* This version is designed for the hat space *) basis5D[{h__}] := Block[ {H = {h}, a, b, c, d, res = {}, N1, norm1, e1, e2, e3, e4, e5}, a = H[[1]]; b = H[[2]]; c = H[[3]]; d = H[[4]]; N1 = a^2 + b^2 + c^2 + d^2; norm1 = N[Sqrt[N1]]; e1 = {-b, a, d, -c, 0}; e1 = e1/norm1; (* Print["e1 = ", e1]; *) e2 = {-c, -d, a, b, 0}; e2 = e2/norm1; (* Print["e2 = ", e2]; *) e3 = {-d, c, -b, a, 0}; e3 = e3/norm1; (* Print["e3 = ", e3]; *) e4 = {a, b, c, d, 0}; e4 = e4/norm1; (* Print["e4 = ", e4]; *) e5 = {0, 0, 0, 0, 1}; res = {e1, e2, e3, e4, e5}; res ]; (* computes the central projection of center a onto the hyperplane through b and orthogonal to ab *) centproj[{a__},{b__},{x__}] := Block[ {A = {a}, B = {b}, l1, l2, i, nab, dotab, result, pt, xx = {x}}, l1 = Length[A]; l2 = Length[B]; (* Print["A = ",A]; Print["B = ",B]; *) nab = 0; Do[ nab = nab + (B[[i]] - A[[i]])^2, {i, 1, l1} ]; (* Print["nab = ",nab]; *) If[x === a, Print[" ** Error, x = a"]; Return[{1, 0, 0, 0}]]; dotab = 0; Do[ dotab = dotab + (B[[i]] - A[[i]])*(xx[[i]] - A[[i]]), {i, 1, l1} ]; (* Print["dotab = ",dotab]; *) result = {}; Do[ result = Append[result, A[[i]] + (nab*(xx[[i]] - A[[i]]))/dotab], {i, 1, l1} ]; result ]; (* Computes the three new nets involved in showing all the patches *) (* method mapping a triangle using the mappings phi1, phi2, phi3 *) (* and displays the closed surface *) (* if debug = -55, call newcnet3 with reftrig3 *) (* if light = 1, uses light triangulation trianglight, *) (* else ltriang *) (* if fcnet = 1, displays control net, otherwise skip it *) renderfour[{net__}, mm_, n_, light_, fcnet_, lx_, mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0, net1, anet, phi1, phi2, phi3, oldnet, surf1, surf2, surf3, surf4, surf, image, stop, phi, reftrig1 = {{-1, 1, 1}, {-1, -1, 3}, {1, 1, -1}}, reftrig2 = {{1, -1, 1}, {1, 1, -1}, {-1, -1, 3}}, reftrig3 = {{-1, 1, 1}, {1, 1, -1}, {1, -1, 1}}, pt, i, j, k}, stop = testm[{net}, mm]; If[stop === 1, Return["*** Unable to Run ***"]]; net0 = maptohat[{net}]; oldnet = net0; If[debug === -55, net1 = newcnet3[net0,mm,reftrig3], net1 = net0 ]; Print["New Net 1: ", net1]; phi = phinets[net1, mm, debug]; phi1 = phi[[1]]; phi2 = phi[[2]]; phi3 = phi[[3]]; surf1 = rsubdiv4[net1, mm, n, debug]; If[light === 1, surf1 = trianglight[surf1,mm], surf1 = ltriang[11,surf1,mm]]; Print["First patch done! "]; surf2 = rsubdiv4[phi1, mm, n, debug]; If[light === 1, surf2 = trianglight[surf2,mm], surf2 = ltriang[12,surf2,mm]]; Print["Second patch done! "]; surf3 = rsubdiv4[phi2, mm, n, debug]; If[light === 1, surf3 = trianglight[surf3,mm], surf3 = ltriang[13,surf3,mm]]; Print["Third patch done! "]; surf4 = rsubdiv4[phi3, mm, n, debug]; If[light === 1, surf4 = trianglight[surf4,mm], surf4 = ltriang[14,surf4,mm]]; Print["Fourth patch done! "]; surf = Join[surf1, surf2]; surf = Join[surf, surf3]; surf = Join[surf, surf4]; Print["Ready to display! "]; If[fcnet === 1, image = {surf, controlnet3D[{net},mm]}, image = surf ]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ]; (* Compute three new rectangular nets *) (* corresponding to the projectivities phi1, phi2, phi3 *) recphinet[{oldnet__}, p_, q_, debug_] := Block[ {net = {oldnet}, i, j, pt, phi1, phi2, phi3, phi}, phi1 = {}; phi2 = {}; phi3 = {}; Do[ Do[ pt = net[[(q + 1)*i + j + 1]]; If[OddQ[p - i], pt = -pt]; phi1 = Append[phi1, pt], {i, 0, p} ], {j, 0, q} ]; Do[ Do[ pt = net[[(q + 1)*i + j + 1]]; If[OddQ[q - j], pt = -pt]; phi2 = Append[phi2, pt], {i, 0, p} ], {j, 0, q} ]; Do[ Do[ pt = net[[(q + 1)*i + j + 1]]; If[OddQ[p + q - i - j], pt = -pt]; phi3 = Append[phi3, pt], {i, 0, p} ], {j, 0, q} ]; phi = Join[{phi1}, {phi2}]; phi = Join[phi, {phi3}]; Print[" rectangular nets phi1, ph2, ph3: ", phi]; phi ]; (* To display a closed rational surface defined by a rectangular net *) (* uses four patches: blue, red, green, yellow *) (* if fcnet = 1, displays control net, otherwise skip it *) (* if light = 1, uses light rectangular grid recmeshlight, *) (* else recmesh *) (* This version first recomputes a net wrt [a, b] X [c, d] *) (* when debug = 99, skip maptohat *) recfrender[{net__}, p_, q_, n_, a_, b_, c_, d_, light_, fcnet_, lx_,mx_,ly_,my_,lz_,mz_,debug_] := Block[ {net0 = {net}, net1, net2, surf0, surf1, surf2, surf3, phi, phi1, phi2, phi3, surf, image, stop}, stop = testpq[net0, p, q]; If[stop === 1, Return["*** Unable to run ***"]]; Print[" Input rectangular net: ", net0]; If[debug === 99, net1 = net0, net1 = maptohat[net0]]; net2 = recnewnet[net1, p, q, a, b, c, d, debug]; Print["New net = ", net2]; surf0 = recsubdiv4[net2, p, q, n, debug]; Print["patch 1 done"]; If[debug === 20, Print["surf = ", surf0]]; If[light === 1, surf0 = recmeshlight[surf0,p,q], surf0 = recmesh[11,surf0,p,q]]; phi = recphinet[net2, p, q, debug]; phi1 = phi[[1]]; phi2 = phi[[2]]; phi3 = phi[[3]]; surf1 = recsubdiv4[phi1, p, q, n, debug]; Print["patch 2 done"]; If[debug === 20, Print["surf = ", surf1]]; If[light === 1, surf1 = recmeshlight[surf1,p,q], surf1 = recmesh[12,surf1,p,q]]; surf2 = recsubdiv4[phi2, p, q, n, debug]; Print["patch 3 done"]; If[debug === 20, Print["surf = ", surf2]]; If[light === 1, surf2 = recmeshlight[surf2,p,q], surf2 = recmesh[13,surf2,p,q]]; surf3 = recsubdiv4[phi3, p, q, n, debug]; Print["patch 4 done"]; If[debug === 20, Print["surf = ", surf3]]; If[light === 1, surf3 = recmeshlight[surf3,p,q], surf3 = recmesh[14,surf3,p,q]]; surf = Join[{surf0}, {surf1}]; surf = Join[surf, {surf2}]; surf = Join[surf, {surf3}]; If[fcnet === 1, image = {surf, recontrolnet3D[net0, p, q]}, image = surf ]; Print["Ready to display! "]; Show[ Graphics3D[{Thickness[0.0006], facecolors[light], image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, lightkind[light], DisplayFunction -> $DisplayFunction]; ];