At the weekend I was looking at this paper The curvature line parametrization from circular
nets on a surface, by Bobenko and Tsarev
…and they describe what I thought was an interesting purely geometric construction for finding principal curvature directions on a mesh.
I tried it out in Grasshopper, and it seems to work fairly well. I’m just sharing it here in case it is of use or interest to anyone else. It may not be the fastest or most robust method for doing this, but I thought it was neat how simple it is, and that this is possible with only native components and no scripting:

Having an idea of the principal curvatures can be useful when you are trying to build a planar quad mesh or developable strip discretization of the surface.

The basic idea is:
Take a small circle through 3 points on the mesh. It will intersect the mesh in at least one more point. Take the diagonals connecting these 4 points, and the bisectors of these are the principal directions.
The ideal size of circle to use depends on the resolution of the mesh. I also noticed it generally works better if you subdivide the mesh first with Weaverbird to get something denser and smoother.

If I recall, there’s a triangulation step before the intersection. The asymmetry might be coming from this, and if so, making sure the triangulation directions are symmetrical should get rid of it.

As for getting continuous curves, certainly one could script something to trace paths using standard vector field integration techniques, with this circle intersection method being used for the sampling.

I think though there could also be a better way by creating a custom 6dof goal in kangaroo that would make a polyline slide around on the mesh to align itself with the principal curvature, based on alignment of one base vector with the mesh normal, and one staying coplanar with those of the adjacent frames, to remove torsion.
I’d have to test this out, but I think it could allow some interesting options for controlling the spacing.

The simplest way of integrating a vector field on a mesh like this is Euler’s method, which is just:

Choose a starting point on the mesh

Evaluate the principal curvature vectors at this point.

Pick one of these 4 directions (+/- for each curvature)

Move the point by this vector multiplied by some small step size h

Project it back onto the mesh

Return to step 2

In step 3, after the first iteration, you’ll need to choose between the options according to which is closest to the previous direction of motion.
Euler’s method is prone to drift, but you can always choose a smaller h to reduce this.
If you need something that lets you take bigger steps while losing less accuracy, you can look at other methods that make multiple evaluations per step, like Verlet or Runge-Kutta.

@DanielPiker I tried a simple case of hypar, when the subdivision is rough, the approximation fails, when it is large, it is good, except naked parts. It is also closely connected to triangulation.

In computer graphics it is usually refereed to Laplace equation to somehow approximate the smooth version of the mesh.

Do you know how to connect Laplace equation to make it work on low res version?

Do you mean like Laplacian smoothing? Yes, if you sample the curvature at some distribution of points, and connect them (could be just by proximity), you can adjust the frame at each point towards the average of it’s neighbours to smooth the field.
Again the complication here arises in choosing how to match the 4 directions between frames in this averaging. Especially if you want it to be rigorous around any singularities in the field.
I have some old scripts for this. I’ll try and dig them out later.

Yes laplace smoothing, I constantly bump the equation below.
I am wondering if the difference between triangulation would not influence the result at specific points, if having this smoothing equation. But I have no idea how to connect the two.

As a more general view, think of Laplacian as a technique that operates with the differential around a point, that is, how the neighboring domain (infinitesimally close) changes with respect to a position. That temperature tends to equilibrium mathematically means that its laplacian tends to zero (minimizing the difference of each point with its neighboring space). In meshes or more generally in graphs, the differential is approximated with the neighboring nodes. In the case of smoothing, the position at the next interaction is a point between the mean of the neighbors and the current position.

Different topology (triangulation) will give a different mesh when smoothed (slightly different).

Yes, using the cotan weights is a way of making the smoothing less dependent on the quality of the triangulation, compared to just taking the simple unweighted average.
Remember though that here we can choose the points at which we take our initial curvature samples. The closer the distribution of these points is to an equilateral grid, the less difference there is between cotan weighted vs uniform weighted.

looks like it works differently in saddle shaped areas and convex areas,
furthermore sometimes it flips the direction of the axis arbitrarily what could it be?

I think I got it, the problem is due to the fact that sometimes the circles do not touch the mesh in 3 points and that messes up the algorithm, I’v solved it by rebuilding the circles in a different way