How can I build a 3D slab of hexagons or a thick hexagon mesh?

Does this do what you are looking for?

The following function takes a parameterization surface and a region region and tries to mesh it with hexagons of radius meshsize. Afterwards, it maps surface over it and creates a mesh of extruded hexagons of thickness thickness

ClearAll[hexhex]
hexhex[surface_, thickness_, meshsize_, region_] := 
 Module[{hex0, centers0, centers, m, n, shifts, planehex, 
   surfacenormal, normals, midlayer, toplayer, bottomlayer, topidx, 
   bottomidx, mantles, B, p0, regmemb},
  B = BoundingRegion[DiscretizeRegion[region]];
  p0 = RegionCentroid[B];
  regmemb = RegionMember[region];

  hex0 = Table[meshsize {Cos[Pi k/3.], Sin[Pi k/3.]}, {k, 0, 5}];
  shifts = MovingAverage[(2. meshsize) Table[ {Cos[Pi k/3.], Sin[Pi k/3.]}, {k, -1, 1}], 2];
  {m, n} = Max /@ Transpose[{
      Ceiling[Abs[LinearSolve[shifts\[Transpose], B[[2]] - p0]]],
      Ceiling[ Abs[LinearSolve[ shifts\[Transpose], {B[[1, 2]], B[[2, 1]]} - p0]]]
      }];
  centers0 = Plus[
    Flatten[Outer[List, Range[-m, m], Range[-n, n]], 1].shifts,
    ConstantArray[p0, (2 m + 1) (2 n + 1)]
    ];
  centers = Pick[centers0, regmemb /@ centers0];
  planehex = Outer[Plus, centers, hex0, 1];

  Quiet[Block[{X},
    surfacenormal = X \[Function] Evaluate[
       Normalize[
        Cross @@Transpose[D[surface[{X[[1]], X[[2]]}], {{X[[1]], X[[2]]}, 1}]]]
       ]]];

  midlayer = Map[surface, planehex, {2}];
  normals = Map[surfacenormal, planehex, {2}];
  toplayer = midlayer + 0.5 thickness normals;
  bottomlayer = midlayer - 0.5 thickness normals;

  topidx = Partition[Append[Range[6], 1], 2, 1];
  bottomidx = Reverse /@ topidx;
  mantles = Join[
    ArrayReshape[
     toplayer[[All, Flatten[topidx]]], {Length[toplayer], 6, 2, 3}],
    ArrayReshape[
     bottomlayer[[All, Flatten[bottomidx]]], {Length[bottomlayer], 6, 
      2, 3}],
    3
    ];
  {toplayer, Flatten[mantles, 1], bottomlayer}
  ]

Usage example:

surface = Quiet[Block[{X, z, f, g},
    f = z \[Function] 1;
    g = z \[Function] z;
    X \[Function] Evaluate[
      N@ComplexExpand[{
          Re[Integrate[f[z] (1 - g[z]^2)/2, z]],
          Re[Integrate[I f[z] (1 + g[z]^2)/2, z]],
          Re[Integrate[f[z] g[z], z]]
          } /. {z -> (X[[1]] + I X[[2]])}]
      ]
    ]];
meshsize = 0.025;
thickness = 0.1;
region = Disk[{0., 0.}, 1.6];
data = hexhex[surface, thickness, meshsize, region];
Graphics3D[{Specularity[White, 30], EdgeForm[{Thin, Black}],
  Darker@Darker@Red, Polygon[data[[1]]],
  Darker@Darker@Blue, Polygon[data[[2]]],
  Darker@Darker@Green, Polygon[data[[3]]]
  },
 Lighting -> "Neutral"
 ]

The planar region may be rather arbitrary, for example, we can use this sea star:

c = t \[Function] (2 + Cos[5 t])/3 {Cos[t], Sin[t]};
region = Module[{pts, edges, B},
   pts = Most@Table[c[t], {t, 0., 2. Pi, 2. Pi/2000}];
   edges = Append[Transpose[{Range[1, Length[pts] - 1], Range[2, Length[pts]]}], {Length[pts], 1}];
   BoundaryMeshRegion[pts, Line[edges]]
   ];
data = hexhex[surface, thickness, 0.5 meshsize, region];
Graphics3D[{Specularity[White, 30], EdgeForm[{Thin, Black}],
  Darker@Darker@Red, Polygon[data[[1]]],
  Darker@Darker@Blue, Polygon[data[[2]]],
  Darker@Darker@Green, Polygon[data[[3]]]
  },
 Lighting -> "Neutral"
 ]

Edit

Because DiscretizeRegion and RegionMember were introduced with version 10, I also provide the following function that takes a list of 6-tuples of points in the plane that represents the hexagon, maps them by the parameterization surface to $\mathbb{R}^3$, and extrudes them.

ClearAll[hexhex2]
hexhex2[surface_, thickness_, planehex_] := Module[{surfacenormal, normals, midlayer, toplayer, bottomlayer, topidx, bottomidx, mantles, B, p0, regmemb},
  Quiet[Block[{X}, 
    surfacenormal = X \[Function] Evaluate[ Normalize[ Cross @@ Transpose[D[surface[{X[[1]], X[[2]]}], {{X[[1]], X[[2]]}, 1}]]]]]];
  midlayer = Map[surface, planehex, {2}];
  normals = Map[surfacenormal, planehex, {2}];
  toplayer = midlayer + 0.5 thickness normals;
  bottomlayer = midlayer - 0.5 thickness normals;
  topidx = Partition[Append[Range[6], 1], 2, 1];
  bottomidx = Reverse /@ topidx;
  mantles = 
   Join[ArrayReshape[
     toplayer[[All, Flatten[topidx]]], {Length[toplayer], 6, 2, 3}], 
    ArrayReshape[
     bottomlayer[[All, Flatten[bottomidx]]], {Length[bottomlayer], 6, 
      2, 3}], 3];
  {toplayer, Flatten[mantles, 1], bottomlayer}];

Use it like this:

h2[x_, y_] := Table[N@{Cos[2 Pi k/6] + x, Sin[2 Pi k/6] + y}, {k, 6}]
planehex = Flatten[Table[h2[3 i + 3 ((-1)^j + 1)/4, Sqrt[3]/2 j], {i, 11}, {j, 12}], 1]; 
surface = X \[Function] {X[[1]], X[[2]], 0.};
thickness = 1;
data = hexhex2[surface, thickness, planehex];
Graphics3D[{Specularity[White, 30], EdgeForm[{Thin, Black}], 
  Darker@Darker@Red, Polygon[data[[1]]], Darker@Darker@Blue, 
  Polygon[data[[2]]], Darker@Darker@Green, Polygon[data[[3]]]}, 
 Lighting -> "Neutral"]

If you want a planar hexagon mesh, then this will do it:

face[{pt1_, pt2_}, h_: 1] := Polygon[{
   Append[pt1, 0],
   Append[pt2, 0],
   Append[pt2, h],
   Append[pt1, h]
   }]
top[pts_, h_: 1] := Polygon[Append[h] /@ pts]
bottom[pts_] := Polygon[Append[0] /@ pts]
hexagon[c_, h_: 1] := Module[{pts, hex},
  pts = Map[c + # &, CirclePoints[6]];
  hex = Prepend[#, Last[#]] &@pts;
  {top[hex, h], bottom[hex], face[#, h] & /@ Partition[hex, 2, 1]}
  ]

hexGrid3D[nx_, ny_, h_: 1] := Table[
   hexagon[{3 i + 3 ((-1)^j + 1)/4, Sqrt[3]/2 j}, 5],
   {i, nx}, {j, ny}
   ];

Graphics3D[
 hexGrid3D[10, 10, 3]
 ]

The argument h in hexGrid3D[nx, ny, h] is the height of the hexagons. nx and ny controls the number of elements in the x and y directions.