BeginPackage["knifepackage`"]; (* * Copyright (c) 1991, The Geometry Center * University of Minnesota * 1300 South Second Street * Minneapolis, MN 55454 * * email address: software@geom.umn.edu * * This software is copyrighted as noted above. It may be freely copied, * modified, and redistributed, provided that the copyright notice is * preserved on all copies. * * There is no warranty or other guarantee of fitness for this software, * it is provided solely "as is". Bug reports or fixes may be sent * to the authors, who may or may not act on them as they desire. * * You may not include this software in a program or other software product * without supplying the source, or without informing the end-user that the * source is available for no extra charge. * * If you modify this software, you should include a notice giving the * name of the person performing the modification, the date of modification, * and the reason for such modification. * * The National Science and Technology Research Center for * Computation and Visualization of Geometric Structures *) dimensionBlade::usage = "dimensionBlade[inequal, iconst] will splice the object defined by the given inequalities by the splicing plane, and show a graph of the 3 dimensional intersection as the plane travels along its perpendicular. Inequalities should be given as lists of coordinates, the first being the plane's normal and the second the constant the normal dotted with {x,y,z(,w)} equals: {{{{ineq11},cons11},{{ineq12},cons12}},{{{ineq21}, cons21},{{ineq22}, cons22}},etc} Note that a union of inequalities can be used to describe the object. Equalities can be used as well, but entered as an option. Type ?Equalities for more information."; findExtremal::usage = "findExtremal[equal, econst, inequal, iconst] finds the extremum points of the intersection of the equalities with the inequalities."; Equalities::usage = "Equalities is an option for dimensionBlade that allows the user to use equalities in addition to inequalities to describe an object. Equalities are entered in the exact format as inequalities as described under ?dimensionBlade. Note that if one object in a union needs an equality, each object has a position in the list which must be filled in as {} if it uses no equalities."; SplicingPlane::usage = "SplicingPlane is an option for dimensionBlade which allows the user to specify which hyperplane (or plane in the three dimensional case) to 'splice' the object with. It is the intersection with this plane which will be shown on screen. Default setting is {0,0,0,1}."; CompoundReflection::usage = "CompoundReflection is an option for dimensionBlade which allows the user to specify whether he wants to reflect an object over different hyperplanes of symmetry independently of each other or allow compound reflections (filling more space faster). Default setting is True."; IntersectionRange::usage = "IntersectionRange is an option for dimensionBlade that allows the user to specify what constant values his splicing hyperplane should fluctuate between. Default setting is {-1, 1}."; StepSize::usage = "StepSize is an option for dimensionBlade which allows the user to specify what the increment size should be for the constant of the splicing hyperplane. Default setting is .2."; ViewingPoint::usage = "ViewingPoint is an option for dimensionBlade which allows the user to specify the ViewPoint from which cross sections should be viewed. Default setting is {5,5,5}."; RoundOff::usage = "RoundOff is an option for dimensionBlade which allows the user to specify the precision that should be used when rounding off solutions to the linear equations solved to find vertices. Default value is 100, meaning solutions are solved to the nearest hundredth."; ShadedFaces::usage = "ShadedFaces is an option for dimensionBlade that allows the user to specify whether the faces of the cross section should be shaded in, or whether only the lattice should be shown. Default value is True (False implies only lattice is shown)."; XRange::usage = "XRange is an option for dimensionBlade that allows the user to specify what range of x-coordinates should be shown in the graph of the cross section. Default setting is {-1,1}."; YRange::usage = "YRange is an option for dimensionBlade that allows the user to specify what range of y-coordinates should be shown in the graph of the cross section. Default setting is {-1,1}."; ZRange::usage = "ZRange is an option for dimensionBlade that allows the user to specify what range of z-coordinates should be shown in the graph of the cross section. Default setting is {-1,1}."; ShowObject::usage = "ShowObject is an option for dimensionBlade that allows the user to specify whether or not the lattice for the entire object should be shown as well as the cross section. Default setting is True."; Reflection::usage = "Reflection is an option for dimensionBlade that allows the user to specify which hyperplanes to reflect over. A list of one or more hyperplanes should be entered in the form: {{{normal1},constant1},{{normal2},constant2},etc.} Default setting is {}, which signifies no reflections." Lattices::usage = "Lattices is an option for dimensionBlade that allows the user to specify whether the edges, or lattice, of the intersection should be graphed. Default setting is true."; ShowBox::usage = "ShowBox is an option for dimensionBlade that allows the user to specify whether to place a box around the area being shown in the graph. Default setting is false."; PrintPolyData::usage = "PrintPolyData is an option for dimensionBlade which allows the user to print out the data being used for the faces of the images being produced. Default setting is false."; PrintTotalData::usage = "PrintTotalData is an option for dimensionBlade which allows the user to print out the raw data being used to produce all output (in the image). Default setting is false."; ShowOnlyProjection::usage = "ShowOnlyProjection is an option for dimensionBlade which allows the user to show only the entire projection, no cross sections, on screen. Default setting is false."; OverLayEdges::usage = "Overlay Edges is an option for dimensionBlade which allows the user to overlay all edges from all slices onto one graph. Default setting is false."; OverLayVertices::usage = "OverLayVertices is an option for dimensionBlade which allows the uswer to overlay all vertices from all slices onto a single graph, thus mapping out the edges of the complete object. Default setting is false."; InstantPicture::usage = "InstantPicture is an option for dimensionBlade which allows the user to show graphs as they are calculated, rather than only seeing them after whole series is calculated. Default setting is false."; AllPictures::usage = "AllPictures is an option for dimensionBlade which allows the user to show the graphs after all calculations for all pictures are finished. This is handy when doing OverLayVertices as all vertices will be shown in all pictures. Default setting is false."; GraphicsList::usage = "GraphicsList is an option for dimensionBlade which causes the function to return a list of all graphics produced, saving all output. Note that to use this option the entered object must contain no equalities.Default setting is false."; Vertices::usage = "Vertices is an option for dimensionBlade that allows the user to enter objects in lists of vertex points. Default value is False." Dual::usage = "Dual is an option for dimensionBlade that allows the user to see the dual of the object described by the entered inequalites. Note that to use this option the entered object must contain no equalities. Default value is False." ObjectData::usage = "ObjectData is an option for dimensionBlade which allows the user to enter equations for the object which all slices should be imbedded in. User defined objects can only be composed of inequalities. Default value is the object initially input."; PrintDualData::usage = "PrintDualData is an option for dimensionBlade which allows the user to print out the equations of the boundaries for the dual being computed. Default value is False."; Begin["Private`"]; pi = N[Pi]; Options[dimensionBlade] := { Equalities -> {}, Vertices -> False, Dual -> False, SplicingPlane -> {0,0,0,1}, Reflection -> {}, CompoundReflection -> True, ShowObject -> True, ObjectData -> {}, ShadedFaces -> True, Lattices -> True, IntersectionRange -> {-1, 1}, StepSize -> .2, XRange -> {-1,1}, YRange -> {-1,1}, ZRange -> {-1,1}, ViewingPoint -> {5,5,5}, RoundOff -> 100, ShowBox -> False, PrintPolyData -> False, PrintTotalData -> False, PrintDualData -> False, OverLayVertices -> False, OverLayEdges -> False, ShowOnlyProjection -> False, InstantPicture -> True, AllPictures -> False, GraphicsList -> False }; dimensionBlade[inequala_List, options___] := Block[{index, data, predata, polyData, objectData, apart = False, inequal = inequala, splicingPlane, pos = 0, obeq, dimension4, compoundReflection, min, max, stepSize,test, equal,elist, viewingPoint, roundOff, shadedFaces, xRange, yRange, zRange, dual, reflection, totalData={}, lattices, showBox, showObject, printPolyData, printTotalData, totalPolyData = {}, overLayEdges, overLayVertices, totalEdgeData = {}, totalVertexData = {}, showOnlyProjection, vertices, instantPicture, allPictures, roundOff, glist = {}, graphicsList, g}, {splicingPlane, compoundReflection, {min, max}, stepSize, equal, viewingPoint, roundOff, shadedFaces, xRange, yRange, zRange, reflection, lattices, showBox, showObject, printPolyData, printTotalData, showOnlyProjection, overLayEdges, overLayVertices, instantPicture, allPictures, roundOff, graphicsList, dual, vertices, objectData, printDualData} = {SplicingPlane, CompoundReflection, IntersectionRange, StepSize, Equalities, ViewingPoint, RoundOff, ShadedFaces, XRange, YRange, ZRange, Reflection, Lattices, ShowBox, ShowObject, PrintPolyData, PrintTotalData, ShowOnlyProjection, OverLayEdges, OverLayVertices, InstantPicture, AllPictures, RoundOff, GraphicsList, Dual, Vertices, ObjectData, PrintDualData} /. {options} /. Options[dimensionBlade]; {dual, vertices} = {Dual,Vertices} /. {options} /. Options[dimensionBlade2]; Share[]; dimension4 = (First[Flatten[Map[Length, inequal, {-2}]]]==4); If[!dimension4 && (Length[splicingPlane] == 4), splicingPlane = {0,0,1}]; If[vertices, inequal = findDual[Transpose[{#, Map[(#.#^.5)&, #]}], roundOff]& /@inequal]; If[dual, inequal = findDual[#, roundOff]& /@ inequal; If[printDualData, Print["Dual's boundaries:"]; Print@inequal;Print[]]]; If[reflection != {}, If[compoundReflection, inequal = reflectTogether[reflection, inequal]; If[equal != {}, equal = reflectTogether[reflection, equal]], inequal = reflectApart[reflection, inequal]; If[equal != {}, equal = reflectTogether[reflection, equal]]]]; If[showObject || showOnlyProjection, If[objectData == {}, objectData = (totalData = inequal); obeq = equal, obeq = {}]; objectData = Flatten[Map[( pos++; If[(obeq == {}) || (obeq[[pos]] == {}), findExtremal[{{},{}},Transpose[#], roundOff, 1], findExtremal[elist=Transpose[obeq[[pos]]], Transpose[#], roundOff, Length@NullSpace@First@elist]])&, objectData], -4]; If[dimension4 && (objectData != {}), objectData = makeList[rotateFourth[splicingPlane, objectData]]]]; If[showOnlyProjection, If[totalData != inequal, pos = 0; totalData = Flatten[Map[( pos++; If[(equal == {}) || (equal[[pos]] == {}), findExtremal[{{},{}},Transpose[#], roundOff, 1], findExtremal[elist=Transpose[equal[[pos]]], Transpose[#], roundOff, Length@NullSpace@First@elist]])&, inequal], -4]; If[dimension4 && (totalData != {}), totalData = makeList[rotateFourth[splicingPlane, totalData]]], totalData = objectData]; glist = {makeGraph[totalEdgeData, totalData, shadedFaces, lattices, showObject, If[objectData != totalData, objectData, {}], overLayEdges, overLayVertices, showBox, viewingPoint, xRange, yRange, zRange, totalVertexData]}, Do[ pos = 0; totalData = Select[Flatten[Map[( If[equal == {}, predata = findExtremal[elist={{splicingPlane}, {index}}, Transpose[#], roundOff, Length@NullSpace@First@elist], predata = findExtremal[elist=Transpose[Append[equal[[++pos]], {splicingPlane, index}]], Transpose[#], roundOff, Length@NullSpace@First@elist]]; While[Depth[predata]<6, predata = {predata}]; If[((Flatten@Union@predata != {}) && dimension4), data = makeList[rotateFourth[splicingPlane, predata]], data = predata])&, inequal],-4], #!={}&]; If[Flatten[Union[totalData]] != {}, AppendTo[totalEdgeData, totalData]]; totalVertexData = Union[Flatten[totalEdgeData, -2]]; If[instantPicture, g=makeGraph[totalEdgeData, totalData, shadedFaces, lattices, showObject, objectData, overLayEdges, overLayVertices, showBox, viewingPoint, xRange, yRange, zRange, totalVertexData]]; If[graphicsList, AppendTo[glist, g]], {index, min, max, stepSize}]; If[allPictures, glist = Map[(makeGraph[totalEdgeData, totalData, shadedFaces, lattices, showObject, objectData, overLayEdges, overLayVertices, showBox, viewingPoint, xRange, yRange, zRange, totalVertexData])&, totalEdgeData]]]; If[graphicsList, Return[glist]]]; makeGraph[totalEdgeData_, totalData_, shadedFaces_, lattices_, showObject_, objectData_, overLayEdges_, overLayVertices_, showBox_, viewingPoint_, xRange_, yRange_, zRange_, totalVertexData_] := Block[{polyData,g}, polyData = Map[goodSols, Map[polygonalize, totalData, {-4}], {-3}]; If[printPolyData, Print@"Polygon Data:"; Print[polyData]; Print[]]; If[printTotalData, Print@"Total Data:";Print@totalData;Print[]]; Show[g = Graphics3D[ {EdgeForm[], If[shadedFaces, Map[(If[Length[Level[#, {-2}]] > 2, Polygon[#], EdgeForm[]])&, polyData, {-3}], EdgeForm[]], If[lattices, Map[Line, totalData, {-3}],EdgeForm[]], If[showObject, Map[Line, objectData, {-3}], EdgeForm[]], If[overLayEdges, Map[(If[#!={},Line@#, EdgeForm[]])&, totalEdgeData, {-3}], EdgeForm[]], If[overLayVertices, Map[(If[#!={}, Point@#, EdgeForm[]])&, totalVertexData, {-2}], EdgeForm[]]}, Boxed -> showBox, ViewPoint -> viewingPoint, PlotRange -> {xRange,yRange,zRange}]]]; findDual[object_List, roundOff_] := Block[{inequal={}, dual, vert}, Map[AppendTo[inequal, #]&, findExtremal[{{},{}}, Transpose@object, roundOff, 1, True], {-2}]; inequal = Union[inequal]; Return[Transpose[{inequal, Map[(#.#)^.5&, inequal]}]]]; findExtremal[{equal_List, econst_List}, {inequal_List, iconst_List}, roundOff_, dim_, vertices_:False] := Block[{pos=0}, If[dim==0, Return@solvePoint[equal, econst, inequal, iconst, roundOff], If[vertices && (equal != {}), Return@Union[Select[Map[getAnswer[equal, econst, inequal, iconst, #, roundOff, vertices, dim, pos]&, Drop[inequal, ++pos]], # != {} &]], Return@Union[Select[Map[getAnswer[equal, econst, inequal, iconst, #, roundOff, vertices, dim, ++pos]&, inequal], # != {} &]]]]]; getAnswer[equal_List, econst_List, inequal_List, iconst_List, elem1_, roundOff_, vertices_, old_, pos_] := Block[{elem2, index, goon, new}, {goon, new} = lessDimension[elem2 = Prepend[equal, elem1], old]; If[goon, Return[findExtremal[{elem2, Prepend[econst, iconst[[pos]]]}, If[!vertices, {Drop[inequal, {pos,pos}], Drop[iconst, {pos,pos}]}, {inequal, iconst}], roundOff, new, vertices]], Return[{}]]]; solvePoint[equal_List, econst_List, inequal_List, iconst_List, roundOff_] := Block[{segment}, segment = Release[Chop[LinearSolve[equal, econst]]]; If[ And[inequal != {}, segment != {}, And @@(LessEqual @@#& /@ Transpose[{inequal.Flatten@segment, iconst+1/roundOff}])], Return[N[Round[roundOff segment] / roundOff]], Return[If[inequal == {}, segment, {}]]]]; lessDimension[equal_List, old_] := Block[{dim}, dim = Length@First@equal-Length@Select[RowReduce@equal, #.#!=0&]; {(dim < old)||(Length[equal]==1), dim}]; reflect[{vector_, const_}, mat_List] := Block[{temp = N[(const / (vector.vector)^.5)], n, n2, c, tran, pos = 0, mat2 = mat}, Map[( n = First[#]; n2 = n-2(n.vector)/(vector.vector)vector; c = Last[#];pos++; tran = const(n.vector)/(vector.vector); mat2[[pos]] = {n2, c - 2tran})&, mat2, {-3}]; Return[mat2]]; reflectApart[reflectionLines_List, inequal_List] := Block[{newInequal = inequal, pos, mat = {}, current}, Map[( pos = 0; current = #; Map[( pos++; mat = #; AppendTo[newInequal, mat = reflect[current, mat]])&, inequal]; )&, reflectionLines]; Return[newInequal]]; reflectTogether[reflectionLines_, inequal_List] := Block[{num = Length[reflectionLines], newInequal = inequal, n, mat = {}}, n = num; Nest[( Map[( mat = #; AppendTo[newInequal, mat = reflect[reflectionLines[[n]], mat]])&, newInequal]; n-- )&, newInequal, num]; Return[newInequal]]; goodSols[faces_List] := Select[faces, Length[#] >= 3 &]; polygonalize[vertices_List] := Block[{elem1, elem2, dummy = Select[vertices, Length[#] == 2 &], pos, newList = {}}, If[dummy != {}, newList = {First[First[dummy]]}; While[dummy != {}, elem1 = Last[newList]; elem2 = First[Select[dummy, MemberQ[#, elem1]&]]; If[!MemberQ[newList, Flatten[Select[elem2, # != elem1 &],1]], AppendTo[newList, Flatten[Select[elem2, # != elem1 &],1]]]; dummy = Drop[dummy, {(pos = First[Flatten[Position[dummy, elem2]]]), pos}]]]; Return[newList]]; makeList[data_List] := Map[Take[#, 3] &, data, {-2}]; rotateFourth[{a_,b_,c_,d_}, vertices_List] := Block[{omega, theta, delta, rotation}, If[b != 0, delta = N[ArcTan[a/b]], If[a != 0, delta = pi/2, delta = 0]]; If[c != 0, theta = N[ArcTan[(a Sin[delta] + b Cos[delta]) / c]], If[({a,b} != {0,0}), theta = pi/2, theta = 0]]; If[d != 0, omega = N[ArcTan[((a Sin[delta] + b Cos[delta]) Sin[theta] + c Cos[theta]) / d]], If[{a,b,c} != {0,0,0}, omega = pi/2, omega = 0]]; rotation = Release[ {{1,0,0,0}, {0,1,0,0}, {0,0,Cos[omega],-Sin[omega]}, {0,0,Sin[omega],Cos[omega]}}. {{1,0,0,0}, {0,Cos[theta],-Sin[theta],0}, {0,Sin[theta],Cos[theta],0}, {0,0,0,1}}. {{Cos[delta],-Sin[delta],0,0}, {Sin[delta],Cos[delta],0,0}, {0,0,1,0}, {0,0,0,1}}]; Chop[Map[rotation.# &, vertices, {-2}]]]; End[]; EndPackage[];