(*******************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which
typically will have the same name as this file except ending in
".nb" instead of ".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent
to using the Evaluate Initialization Cells menu command in the front
end.

DO NOT EDIT THIS FILE.  This entire file is regenerated
automatically each time the parent Notebook file is saved in the
Mathematica front end.  Any changes you make to this file will be
overwritten.
***********************************************************************)



SetDirectory["./"];







ufile = "u.txt";



chopTolerance = 10^(-8);
xyMaxTolerance = 10^(-3);
changeSign =False;

zeroColor = RGBColor[1,0,0];
zeroThickness = Thickness[.004];

posColor =RGBColor[1,1,1];
posThickness = Thickness[.004];

negColor = RGBColor[0,0,0];
negThickness = Thickness[.004];

zeroRegionColor = RGBColor[.5,.5,.5];
posRegionColor = RGBColor[0,0,0];
negRegionColor = RGBColor[1,1,1];
radiusofDisk = 0.005;
bdyLevel = 5 ;
bdyThickness = Thickness[0.004];
imgSize = 400;

nContoursMax5 = 5;
nContoursMax4 = 3;



zListUnchopped = Import[ufile, "List"];
numPts = Length[zListUnchopped];
Which[numPts \[Equal] 1261, 
    Lev = 4; nContoursMax = nContoursMax4, 
    numPts \[Equal] 11605, 
    Lev = 5; nContoursMax = nContorusMax5,
    True, Print["Only level 4 and 5 are supported."]];






\!\(\(bdyfile\  = \ "\<jboundary\>" <> 
        ToString[bdyLevel] <> "\<.txt\>";\)\[IndentingNewLine]
  \(bdyList\  = Import[bdyfile, \ "\<Table\>"];\)\[IndentingNewLine]
  \(bdyList\  = \ bdyList\ /\ \@3. ;\)\[IndentingNewLine]
  \(Print["\<the bdyList has length \>", \ 
      Length[bdyList]\ ];\)\[IndentingNewLine]
  \(bdyPlot\  = \ Graphics[{bdyThickness, Line[bdyList]}];\)\)



\!\(\(lineSegmentQ[z1_, z2_, \ z3_]\  := \ 
      Module[{}, \[IndentingNewLine]If[z1 < 0\  && z2 < 0 && z3 < 0, \ 
          Return[False]\ ]; \[IndentingNewLine]If[
          z1 > 0\  && z2 > 0 && z3 > 0, \ 
          Return[False]\ ]; \[IndentingNewLine]If[
          z1 \[Equal] 0\  && z2*z3 > 0, \ 
          Return[False]]; \[IndentingNewLine]If[
          z2 \[Equal] 0\  && z3*z1 > 0, \ 
          Return[False]]; \[IndentingNewLine]If[
          z3 \[Equal] 0\  && z1*z2 > 0, \ 
          Return[False]]; \[IndentingNewLine]If[
          z1 \[Equal] 0\  && z2 \[Equal] 0 && z3 == 0, \ 
          Return[False]]; \[IndentingNewLine]True];\)\[IndentingNewLine]
  \(contourSegment[v1_, \ v2_, \ v3_, \ z1_, \ z2_, \ 
        z3_]\  := \ \[IndentingNewLine]Module[{endPoints}, \
\[IndentingNewLine]If[z1 < 0\  && z2 < 0 && z3 < 0, \ 
          Return[Null]\ ]; \[IndentingNewLine]If[
          z1 > 0\  && z2 > 0 && z3 > 0, \ 
          Return[Null]\ ]; \[IndentingNewLine]If[
          z1 \[Equal] 0\  && z2*z3 > 0, \ 
          Return[Null]]; \[IndentingNewLine]If[
          z2 \[Equal] 0\  && z3*z1 > 0, \ 
          Return[Null]]; \[IndentingNewLine]If[
          z3 \[Equal] 0\  && z1*z2 > 0, \ 
          Return[Null]]; \[IndentingNewLine]If[
          z1 \[Equal] 0\  && z2 \[Equal] 0 && z3 == 0, \ 
          Return[Null]]; \[IndentingNewLine]\[IndentingNewLine]If[
          z1\  \[Equal] \ 0\  && \ z2\  \[Equal] \ 0\  && \ 
            z3 \[NotEqual] 0, \ 
          Return[Line[{v1, \ v2}]]]; \[IndentingNewLine]If[
          z2\  \[Equal] \ 0\  && \ z3\  \[Equal] \ 0\  && \ 
            z1 \[NotEqual] 0, \ 
          Return[Line[{v2, \ v3}]]]; \[IndentingNewLine]If[
          z3\  \[Equal] \ 0\  && \ z1\  \[Equal] \ 0\  && \ 
            z2 \[NotEqual] 0, \ 
          Return[Line[{v3, \ 
                v1}]]]; \[IndentingNewLine]\[IndentingNewLine]endPoints\  = \ \
{}; \[IndentingNewLine]If[z1\  \[Equal] \ 0\  && \ z2*z3 < 0, \ 
          endPoints\  = \ {v1}]; \[IndentingNewLine]If[
          z2\  \[Equal] \ 0\  && \ z3*z1 < 0, \ endPoints\  = \ {v2}]; 
        If[z3\  \[Equal] \ 0\  && \ z1*z2 < 0, \ 
          endPoints\  = \ {v3}]; \[IndentingNewLine]If[z1\ *\ z2\  < \ 0, \ 
          AppendTo[
            endPoints, \ \(z2*\ v1\  - \ z1*v2\)\/\(z2 - z1\)]\ ]; \
\[IndentingNewLine]If[z2\ *\ z3\  < \ 0, \ 
          AppendTo[endPoints, \ \(z3*\ v2\  - \ z2*v3\)\/\(z3 - z2\)]\ ]; 
        If[z3\ *\ z1\  < \ 0, \ 
          AppendTo[
            endPoints, \ \(z1*\ v3\  - \ z3*v1\)\/\(z1 - z3\)]\ ]; \
\[IndentingNewLine]Return[Line[endPoints]]\[IndentingNewLine]];\)\)

readxyFiles[dataFile_] := Module[{neighbors, xy},
      OpenRead[dataFile];
      neighborsList = {};
      xyList = {};
      ptLabel = Read[dataFile, Number];
      While[NumberQ[ptLabel],
        neighbors = {};
        Do[neighbor = Read[dataFile, Number];
          If[neighbor > 0, AppendTo[neighbors, neighbor]],{6}];
        AppendTo[neighborsList, neighbors];
        xy = Read[dataFile,{Number, Number}];
        AppendTo[xyList, xy];
        ptLabel = Read[dataFile, Number]
        ];
      numPts = Length[xyList];
      Close[dataFile];
      ];

\!\(\(hex[{x_, y_}, \ h_]\  := \ 
      Polygon[{{x + h\/2, \ y + h\/\(2 \@ 3\)}, \ {x\ , \ 
            y + \ h\/\@3}, \ {x - h\/2, \ y + h\/\(2 \@ 3\)}, \ {x - h\/2, \ 
            y - h\/\(2 \@ 3\)}, \ \ {x\ , \ y - \ h\/\@3}, \ {x + h\/2, \ 
            y - h\/\(2 \@ 3\)}, \ {x + h\/2, \ y + h\/\(2 \@ 3\)}}];\)\)



\!\(\(Clear[a, \ b, \ c, \ d, \ e, \ f, \ g, \ h];\)\[IndentingNewLine]
  \(data\  = {{0, 0, \ a}, \ {\(-1\), \ 0, \ b}, \ {\(-1\)/2, \ \@3/2, \ 
          c}, \ {\(-1\)/2, \ \(-\@3\)/2, \ d}, \ {1/2, \ \@3/2, \ 
          e}, \ {1/2, \ \(-\@3\)/2, \ f}, \ {1, \ 0, \ 
          g}};\)\[IndentingNewLine]
  \(fxy[a_, b_, \ c_, \ d_, \ e_, \ f_, \ g_]\  = \ 
      Chop[Fit[data, {x^2, \ y^2, \ x\ y, x, \ y, \ 1}, \ {x, 
            y}]];\)\[IndentingNewLine]
  \(soln\  = \ 
      Solve[{\[PartialD]\_x\ fxy[a, b, c, d, e, f, g] \[Equal] 
            0, \ \[PartialD]\_y\ fxy[a, b, c, d, e, f, g] \[Equal] 0}, \ {x, 
          y}];\)\[IndentingNewLine]
  \(critPt[a_, {b_, c_, d_, e_, f_, g_}, \ h_]\  = \ 
      h {x, y} /. 
        soln\_\(\(\[LeftDoubleBracket]\)\(1\)\(\[RightDoubleBracket]\)\);\)\)



h = 2/3^Lev;
Print["The plot will be made at level ", Lev];trianglefile = 
  "jtriangles"<>ToString[Lev]<>".txt";
triangleList = Import[trianglefile, "Table"];
Print["There are ", Length[triangleList], " triangles in the ", trianglefile,
    " file."];





xyFile= "jall.nbd.xy."<>ToString[Lev]<>".txt";
readxyFiles[xyFile]
Print["There are ", numPts, " grid points at level ", Lev,"."];



zList = Chop[zListUnchopped, chopTolerance];
If[changeSign, zList = - zList];



\!\(\(xyMaxList\  = \ {posColor};\)\[IndentingNewLine]
  \(xyMinList\  = \ {negColor};\)\[IndentingNewLine]
  Do[zi\  = \ 
      zList\_\(\(\[LeftDoubleBracket]\)\(i\)\(\[RightDoubleBracket]\)\); \
\[IndentingNewLine]znbrs\  = \ 
      zList\_\(\(\[LeftDoubleBracket]\)\(neighborsList\_\(\(\
\[LeftDoubleBracket]\)\(i\)\(\[RightDoubleBracket]\)\)\)\(\
\[RightDoubleBracket]\)\); \[IndentingNewLine]maxNbrs\  = \ 
      Max[znbrs]; \[IndentingNewLine]minNbrs\  = \ 
      Min[znbrs]; \[IndentingNewLine]If[
      zi\  > \ maxNbrs\  && \ maxNbrs\  > \ minNbrs\  && \ 
        Length[znbrs]\  \[Equal] \ 6, \[IndentingNewLine]cp\  = \ 
        xyList\_\(\(\[LeftDoubleBracket]\)\(i\)\(\[RightDoubleBracket]\)\)\ \ \
 + \ critPt[zi, \ znbrs, \ h]; \[IndentingNewLine]AppendTo[xyMaxList, \ 
        Disk[cp, \ radiusofDisk]]\[IndentingNewLine]]; \[IndentingNewLine]If[
      zi\  < \ minNbrs\  && \ maxNbrs\  > \ minNbrs\  && \ 
        Length[znbrs]\  \[Equal] \ 6, \[IndentingNewLine]cp\  = \ 
        xyList\_\(\(\[LeftDoubleBracket]\)\(i\)\(\[RightDoubleBracket]\)\)\ \ \
 + \ critPt[zi, \ znbrs, \ h]; \[IndentingNewLine]AppendTo[xyMinList, \ 
        Disk[cp, \ 
          radiusofDisk]]\[IndentingNewLine]];\[IndentingNewLine], {i, \ 
      numPts}]\[IndentingNewLine]
  \(numLocalMax\  = \ Length[xyMaxList] - 1;\)\[IndentingNewLine]
  \(numLocalMin\  = \ Length[xyMinList]\  - 1;\)\[IndentingNewLine]
  \(numLocalExtrema\  = \ numLocalMax\  + \ numLocalMin;\)\[IndentingNewLine]
  \(Print["\<I found \>", 
      numLocalMax, \ "\< local maximum points.\>"];\)\[IndentingNewLine]
  \(\(Print["\<I found \>", 
      numLocalMin, \ "\< local minimum points.\>"];\)\(\[IndentingNewLine]\)
  \)\)





zMax = Max[zList];
zMin = Min[zList];
supNorm = Max[zMax, -zMin];
If[numLocalExtrema \[Equal] 0, numLocalExtrema = 1];
nContours = Ceiling[nContoursMax/Sqrt[numLocalExtrema]];
Print["nContours = ", nContours];
contours = Table[supNorm*i/(nContours+1),{i, 1, nContours}];



\!\(\(posLineSegments\  = \ {posColor, \ posThickness};\)\[IndentingNewLine]
  \(zeroLineSegments\  = \ {zeroColor, zeroThickness};\)\[IndentingNewLine]
  \(negLineSegments\  = \ {negColor, \ negThickness};\)\[IndentingNewLine]
  \(Do[\[IndentingNewLine]vertexLabelList\  = \ 
        triangleList\_\(\(\[LeftDoubleBracket]\)\(tLabel\)\(\
\[RightDoubleBracket]\)\); \[IndentingNewLine]{v1, \ v2, \ v3}\  = \ 
        xyList\_\(\(\[LeftDoubleBracket]\)\(vertexLabelList\)\(\
\[RightDoubleBracket]\)\); \[IndentingNewLine]z123 = \({z1, \ z2, \ z3}\  = \ 
          zList\_\(\(\[LeftDoubleBracket]\)\(vertexLabelList\)\(\
\[RightDoubleBracket]\)\)\); \[IndentingNewLine]If[lineSegmentQ[z1, z2, z3], 
        AppendTo[zeroLineSegments, \ 
          contourSegment[v1, v2, v3, z1, z2, 
            z3]]]; \[IndentingNewLine]Do[\[IndentingNewLine]zValues\  = \ 
          z123\  - \ 
            contours\_\(\(\[LeftDoubleBracket]\)\(j\)\(\[RightDoubleBracket]\)\
\); \[IndentingNewLine]{z1, z2, z3}\  = \ zValues; \[IndentingNewLine]If[
          lineSegmentQ[z1, z2, z3], 
          AppendTo[posLineSegments, \ 
            contourSegment[v1, v2, v3, \ z1, z2, 
              z3]]]; \[IndentingNewLine]{z1, z2, z3}\  = \ 
          z123\  + \ 
            contours\_\(\(\[LeftDoubleBracket]\)\(j\)\(\[RightDoubleBracket]\)\
\); \[IndentingNewLine]If[lineSegmentQ[z1, z2, z3], 
          AppendTo[negLineSegments, \ 
            contourSegment[v1, v2, v3, z1, z2, z3]]], {j, \ 
          nContours}]\[IndentingNewLine], {tLabel, 1, 
        Length[triangleList]}];\)\[IndentingNewLine]
  \(ctPlot\  = \ {bdyPlot, 
        Graphics[{zeroLineSegments, posLineSegments, negLineSegments, \ 
            xyMinList, \ xyMaxList}]};\)\)



\!\(\(hexListPos\  = \ {posRegionColor};\)\[IndentingNewLine]
  \(hexListZero\  = \ {zeroRegionColor};\)\[IndentingNewLine]
  \(Do[zi\  = \ 
        zList\_\(\(\[LeftDoubleBracket]\)\(i\)\(\[RightDoubleBracket]\)\); \
\[IndentingNewLine]If[zi\  > \ 0, 
        AppendTo[hexListPos, 
          hex[xyList\_\(\(\[LeftDoubleBracket]\)\(i\)\(\[RightDoubleBracket]\)\
\), \ h]]]; \[IndentingNewLine]If[zi\  \[Equal] \ 0, \ 
        AppendTo[hexListZero, \ 
          hex[xyList\_\(\(\[LeftDoubleBracket]\)\(i\)\(\[RightDoubleBracket]\)\
\), \ h]]], \ {i, \ numPts}];\)\[IndentingNewLine]
  \(bgPlot\  = \ Graphics[{hexListPos, \ hexListZero}];\)\)



fig = Show[bgPlot, ctPlot,DisplayFunction\[Rule]Identity,
      ImageSize \[Rule]imgSize,AspectRatio\[Rule]Automatic];

Export["contourPlot.eps", fig];
Export["u.gif", fig];

