1. Problem statement

An old answer of mine: Computing Gaussian curvature shows how to color parametric curves and surfaces according to their curvature. However, the methods used there will fail for non-differentiable and very complicated functions. A workaround, not perfect, but the only one I know, is to color such surfaces according to their polygon sizes. Polygon areas decrease in zones of high curvature. I know how to use this "area method" in my Blender/Python application, but don't know how to transfer it to Mathematica.

2. Curvature of a Torus (formula method)

GaussianCurvature[f_] :=
 With[{dfu = D[f, u], dfv = D[f, v]}, 
  Simplify[(Det[{D[dfu, u], dfu, dfv}] * Det[{D[dfv, v], dfu, dfv}] - 
      Det[{D[f, u, v], dfu, dfv}]^2) / (dfu.dfu * dfv.dfv - (dfu.dfv)^2)^2]]

torus = {(2 + Cos[v]) Cos[u], (2 + Cos[v]) Sin[u], Sin[v]};

curvature = GaussianCurvature[torus]

Cos[v] / (2 + Cos[v])

range = Last[PlotRange /. AbsoluteOptions[Plot3D[curvature, {u, 0, 2 Pi}, {v, 0, 2 Pi}], PlotRange]]

{-1., 0.333333}

Legended[
 ParametricPlot3D[torus, {u, 0, 2 Pi}, {v, 0, 2 Pi},
  ColorFunction ->
   Function[{x, y, z, u, v},ColorData["TemperatureMap"][Rescale[curvature, range]]],
  ColorFunctionScaling -> False,
  ImageSize -> 400,
  Mesh -> False,
  PlotPoints -> 70],
 BarLegend[{"TemperatureMap", range}, Automatic]]

3. Curvature of a Torus with Blender/Python (area method)

The Blender 3.6.5 Python API has a method, calc_area() which returns a list of polygon sizes (basically with one line of code):

p.list = [f.calc_area() for f in bm.faces]

("Faces" are rectangular polygons defined by 4 points). A rescaled list of these sizes is then used to index the color table and color the polygons. We can see that the Blender result is very similar to the above Mathematica plot:

4. The Klein bottle

KleinBottle[u_, v_] :=
 Module[{a, b, c, x, y, z},
  a = 6 Cos[u] (1 + Sin[u]);
  b = 16 Sin[u];
  c = 5 (1 - Cos[u]/2);
  x = If[Pi < u <= 2 Pi, a + c Cos[v + Pi], a + c Cos[u] Cos[v]];
  y = If[Pi < u <= 2 Pi, b, b + c Sin[u] Cos[v]];
  z = c Sin[v];
  {x, y, z}]

klein =
 ParametricPlot3D[
  KleinBottle[u, v], {u, 0, 2 Pi}, {v, 0, 2 Pi},
  PlotPoints -> 60]

Because it is not possible to apply the curvature formula to this complicated surface, I show how it would look like in Blender using the "polygon area method" (with a logarithmic scale):

5. Question

I know how to extract the polygons, f.e. the first polygon of the Klein bottle:

Cases[Normal[klein], _Polygon, Infinity][[1]] // PolygonCoordinates

{{8.80683, 0.984612, 0.133244}, {9.10008, 1.96794, 2.67745*10^-7}, {9.08591, 1.96642, 0.267239}}

But I don't know how to proceed from there to calculate the polygon sizes and color the bottle similar to the last screenshot.