This SetDirectory command is needed in the PC version of Mathematica.  It is not an "initialization cell" so it does not get written to the .m file.

SetDirectory["./"];

This notebook writes a .m file which can be run non-interactively.  The notebook uses the "u.txt" file in the current directory and produces two figures, contourPlot.eps and contourPlot.gif.  It is assumed that several files are in the current directory.

The readxyFiles[dataFile_] function produces the global variables:
*neighborsList : List of lists of neighbors.  \(TraditionalForm\`neighborsList\_\(\([[\)\(5\)\(]]\)\)\) is a list of the neighbors of point 5.  The "0"s in the text file which indicate fewer than 6 neighbors are not included in this list.
*xyList
*numPts : length of any of the previous lists.

Other global varibles:
bdyList
bdyPlot
triangleList
lineSegmentQ[z1_,z2_, z3_]
contourSegment[v1_, v2_, v3_, z1_, z2_, z3_]

To Do:
If the point has less than 6 neighbors, plot a star of David rather than a hexagon.  This will avoid the white showing through near the boundary.
Print white and black triangles as required, in place of the zero contour which is now red.

Current Problems:  
There is a problem at the x points of the zero contour.  A triangle with one u value 0 (on a line of reflection symmetry) and two values of opposite sign on the same side of the line of symmetry automatically give a contour line segment which goes straight to the grid point on the line of reflection symmetry.

ufile  =  "u.txt";

The following cell has parameters that one might want to tweak.

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;

I might also want to make the boundary level (bdyLev) depend on level, like the nContoursMax values.

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

The first part produces the boundary by reading in the jboundaryL.txt files produced by make Region List.nb.  This text file has the points on the boundary listed in order.

The readxyzFiles[ ] function produces the global variables:
*neighborsList : List of lists of neighbors.  \(TraditionalForm\`neighborsList\_\(\([[\)\(5\)\(]]\)\)\) is a list of the neighbors of point 5.  The "0"s in the text file which indicate fewer than 6 neighbors are not included in this list.
*xyList
*zList
*numPts : length of any of the previous lists.

The factor of sqrt(3) below is there because it was produced when we were using a different size grid.

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

lineSegmentQ[z1_, z2_,  z3_]  :=  
        Module[{}, 
If[z1 < 0  && z2 < 0 && z3 < 0,  
            Return[False] ]; 
If[
            z1 > 0  && z2 > 0 && z3 > 0,  
            Return[False] ]; 
If[
            z1 == 0  && z2*z3 > 0,  
            Return[False]]; 
If[
            z2 == 0  && z3*z1 > 0,  
            Return[False]]; 
If[
            z3 == 0  && z1*z2 > 0,  
            Return[False]]; 
If[
            z1 == 0  && z2 == 0 && z3 == 0,  
            Return[False]]; 
True];
contourSegment[v1_,  v2_,  v3_,  z1_,  z2_,  
          z3_]  :=  
Module[{endPoints}, 
If[z1 < 0  && z2 < 0 && z3 < 0,  
            Return[Null] ]; 
If[
            z1 > 0  && z2 > 0 && z3 > 0,  
            Return[Null] ]; 
If[
            z1 == 0  && z2*z3 > 0,  
            Return[Null]]; 
If[
            z2 == 0  && z3*z1 > 0,  
            Return[Null]]; 
If[
            z3 == 0  && z1*z2 > 0,  
            Return[Null]]; 
If[
            z1 == 0  && z2 == 0 && z3 == 0,  
            Return[Null]]; 

If[
            z1  ==  0  &&  z2  ==  0  &&  
              z3 ≠ 0,  
            Return[Line[{v1,  v2}]]]; 
If[
            z2  ==  0  &&  z3  ==  0  &&  
              z1 ≠ 0,  
            Return[Line[{v2,  v3}]]]; 
If[
            z3  ==  0  &&  z1  ==  0  &&  
              z2 ≠ 0,  
            Return[Line[{v3,  
                  v1}]]]; 

endPoints  =  {}; 
If[z1  ==  0  &&  z2*z3 < 0,  
            endPoints  =  {v1}]; 
If[
            z2  ==  0  &&  z3*z1 < 0,  endPoints  =  {v2}]; 
          If[z3  ==  0  &&  z1*z2 < 0,  
            endPoints  =  {v3}]; 
If[z1 * z2  <  0,  
            AppendTo[
              endPoints,  ((z2* v1  -  z1*v2)/(z2 - z1))] ]; 
If[z2 * z3  <  0,  
            AppendTo[endPoints,  ((z3* v2  -  z2*v3)/(z3 - z2))] ]; 
          If[z3 * z1  <  0,  
            AppendTo[
              endPoints,  ((z1* v3  -  z3*v1)/(z1 - z3))] ]; 
Return[Line[endPoints]]
];

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 (Sqrt[ 3])))},  {x ,  
              y +  (h/(Sqrt[3]))},  {x - (h/2),  
              y + (h/(2 (Sqrt[ 3])))},  {x - (h/2),  
              y - (h/(2 (Sqrt[ 3])))},   {x ,  y -  (h/(Sqrt[3]))},  {x + (h/2),  
              y - (h/(2 (Sqrt[ 3])))},  {x + (h/2),  
              y + (h/(2 (Sqrt[ 3])))}}];

The following function makes a better position of the maximum points on the hexagonal grid.  I'll assume that the middle point is a maximum (with value a), and that the 6 neighboring points are all in the grid, with values b, c, d, e, f, and g. The order of the grid points is always the same, as listed in the data below.  I use Mathematica's ability to do a least square fit symbolically.  The best fit of the 6 coeffecients in the general quadratic function of 2 variables are found for a function with values at the center of a hexagon and the 6 vertices.

The result is a function crtitPt[uMax, uNbhList, h] that gives the position of the critical point of the quadratic fit. 

(With zero Neuman BCs, this would function would have to be modified, since the max/min points would often be on the boundary.)

Clear[a,  b,  c,  d,  e,  f,  g,  
        h];
data  = {{0, 0,  a},  {-1,  0,  b},  {-1/2,  (Sqrt[3])/2,  
            c},  {-1/2,  -(Sqrt[3])/2,  d},  {1/2,  (Sqrt[3])/2,  
            e},  {1/2,  -(Sqrt[3])/2,  f},  {1,  0,  
            g}};
fxy[a_, b_,  c_,  d_,  e_,  f_,  g_]  =  
        Chop[Fit[
            data, {x^2,  y^2,  x y, x,  y,  1},  {x, 
              y}]];
soln  =  
        Solve[{(∂_x) fxy[a, b, c, d, e, f, g] == 
              0,  (∂_y) fxy[a, b, c, d, e, f, g] == 
              0},  {x, y}];
critPt[a_, {b_, c_, d_, e_, f_, g_},  h_]  =  
        h {x, y} /. 
          (Subscript[soln, [[1]]]);

It takes about 13 seconds at level 5 and about 2 seconds at level 4 to read in the triangle files on Jim's PC.

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];

The next cell finds the local max and min points.  The test of Length[znbrs] == 6 is there because without it there are some false local extrema near the boundary.  This should work with Neumann solutions too: the local extrema on the boundary are just ignored.

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

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

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

In the future, there will be no hexListZero

hexListPos  =  {posRegionColor};
hexListZero  =  {zeroRegionColor};
Do[
        zi  =  
          (Subscript[zList, [[i]]]); 
If[zi  >  0, 
          AppendTo[hexListPos, 
            hex[(Subscript[xyList, [[i]]]),  h]]]; 
If[zi  ==  0,  
          AppendTo[hexListZero,  
            hex[(Subscript[xyList, [[i]]]),  h]]],  {i,  numPts}];
bgPlot  =  Graphics[{hexListPos,  hexListZero}];

If we want to put a text box in the figure we can use something like this:
titlePlot = 
    Graphics[Text[
        StyleForm[figTitle, FontSize -> txtSize], {.3,-.58},{0,-1}]];
        fig = Show[bgPlot, ctPlot, titlePlot, ImageSize -> 300];

fig  =  
        Show[bgPlot,  ctPlot, DisplayFunction -> Identity, 
          ImageSize  -> imgSize, 
          AspectRatio -> Automatic];

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