The short story:
I would like to ask for the logic behind the planarize component in the Ngon plug in.

The longer story:
I started recreating these days an alogrithm for planarzing hexagonal meshes I found here: https://webspace.science.uu.nl/~vaxma001/
in this article: Dupin Meshing: A Parameterization Approach to Planar Hex-Dominant Meshing

My knowledge of meshes is really small and I just got overhelmed when I saw these articles.So glad I found them.
The process is shown in this picture:

In this algorithm the final hexagons are reshaped by the base surface curvature which is obtained in step(a).

Recreating the hexagon mesh(my first nGon andventure in c#) I saw there is a planarize component in the fantastic Ngon Plug-in which planarizes hexagons without underlaying surface.

I would like to ask what the logic here is, since I have no clue.

2019 Download
Petras Vestartas, Nicolas Rogeau, Julien Gamerro, Yves Weinand
Modelling Workflow for Segmented Timber Shells using Wood-wood Connections
Design Modelling Symposium (DMS) Berlin 2019.

The planarization is probably the most basic one, when polygon faces are projected to average plane and then their vertices are combined together. This is run several times.
I used this method to planarize extruded mesh edges:

Other method involves tangent plane, when neighbour planes are intersected to computed planar polygon. Which is also constrained to particular geometries (mostly spheres, and negative curvature cases, so flat parts of geometry is always an issue):

I would suggest also to look at Kangaroo2 methods.
I am aware of the paper that you reference, but I doubt there is a close implementation in Rhino.
I believe it is also connected to remeshing and particular geometries. Just look at the naked edges of the geometry shown your figure, they are cut at the border. Meaning there is a surface parametrization to remesh underlying quad patches.

I think @DanielPiker mentioned that hexagons can be represented by two quads so two circles could be inscribed in this hexagon must be planar. There was a lot of discussion on the forum:

There is a general misconception about planarization, that you can put whatever tessellation and solver will give you a planar mesh. It is really connected to underlying curvature and subdivision and often 4 sided nurbs subdivision would not planarize well because the mesh cells are not oriented according to curvature as is described in the paper.

Thank you very much for taking the time to answer so extensively!Really apreciated!

Your paper is just perfect for me, since the algorithm is easier to understand for me and I think you explained it very well.Congratulations on that. Your work has been an inspiration for me since several years.

I saw the tutorials on K2, also for evolution tools, but my incentive is to learn more about meshes and planarisation, tthat´s why I wanted to code it.

The discussion of the Landesgartenschau originally piqued my interest. There is where I saw first the vector field for the surface curvature, which I was able to recreate thanks to a discussion from you and Daniel Piker.
What I still not understand completly in the Landesgartenschau algorithm is the how the agents are influenced by the forcefield, I assume they are guided by some amount of the vectors they pass by.

I searched some detailed explanations to recreate the process of this pavilion, that´s why I ended up in the paper I posted.

Also your physical models are really nice, my aim is to make a physical model after recreating the planarisation algorithm.

May I ask one more thing: In your code , how do I obtain the weigths array/ the weigth for each point(is it the amaount of faces it belongs?)

It seems in the code I just use weight 1 for all vertices. If you want you can explore this further:
For the paper you are talking about you can ask Long or try to read his paper.

public void Solve2(ref Mesh mesh, int[] allNGonVArray, ref double deviation, int[][] faceMap, bool[] nakedStatus, ref bool[] planar, double tangentialT)//ref Plane[] NGonPlanes,
{
if (tangentialT > 0)
{
Point3d[] averageP = new Point3d[allNGonVArray.Length];
Vector3d[] averagePVecs = new Vector3d[allNGonVArray.Length];
Vector3d[][] averagePVecs_ = new Vector3d[allNGonVArray.Length][];
int[] averagePCount = new int[allNGonVArray.Length];
for (int i = 0; i < allNGonVArray.Length; i++)
{
//Position center point
//Vec array for surrounding vertices
int n = conV[i].Count + 1;
averagePVecs_[i] = new Vector3d[n];
Point3d Avg = new Point3d();
var Vecs = new Vector3d[n - 1];
for (int j = 1; j < n; j++)
{
Avg += mesh.TopologyVertices[conV[i][j - 1]]; //sum up all vertices
Vecs[j - 1] = mesh.TopologyVertices[conV[i][j - 1]] - mesh.TopologyVertices[allNGonVArray[i]]; //get vector to center point
}
double Inv = 1.0 / (n - 1);
Avg *= Inv;
//Vector3d Smooth = 0.5 * (Avg - p[PIndex[0]].Position);
Vector3d Smooth = 0.5 * (Avg - mesh.TopologyVertices[allNGonVArray[i]]);
Vector3d Normal = new Vector3d();
for (int j = 0; j < Vecs.Length; j++)
{
Normal += Vector3d.CrossProduct(Vecs[j], Vecs[(j + 1) % Vecs.Length]);
}
Normal.Unitize();
Smooth -= Normal * (Normal * Smooth);
//Assign values
averagePVecs_[i][0] = Smooth;
averagePVecs[i] += Smooth;
averagePCount[i] += 1;
//Rhino.RhinoApp.WriteLine(Smooth.ToString());
Smooth *= -Inv;
for (int j = 1; j < n; j++)
{
int localID = Array.IndexOf(allNGonVArray, conV[i][j - 1]);
averagePVecs_[i][j] = Smooth;
averagePVecs[localID] = Smooth;
averagePCount[localID] += 1;
//Rhino.RhinoApp.WriteLine(Smooth.ToString());
}
}
//Assign vector one by one for currect point and neighbours
for (int i = 0; i < allNGonVArray.Length; i++)
{
if (this.nakedStatusSimple[allNGonVArray[i]])
continue;
mesh.TopologyVertices[allNGonVArray[i]] += (Vector3f)(averagePVecs_[i][0] / averagePCount[i] * tangentialT);
for (int j = 1; j < averagePVecs_[i].Length; j++)
{
int localID = Array.IndexOf(allNGonVArray, conV[i][j - 1]);
//Rhino.RhinoApp.WriteLine((averagePVecs_[i][j] / averagePCount[i]).ToString());
mesh.TopologyVertices[allNGonVArray[localID]] += (Vector3f)(averagePVecs_[i][j] / averagePCount[i] * tangentialT);
}
}
}
double currentDeviation = 0;
Point3d[] allPts = new Point3d[allNGonVArray.Length];
double[] allW = new double[allNGonVArray.Length];
//Get vector per face
for (int i = 0; i < faceV.Length; i++) {
Point3d[] facePts = new Point3d[faceV[i].Length];
for (int j = 0; j < faceV[i].Length; j++)
facePts[j] = mesh.TopologyVertices[faceV[i][j]];
Plane.FitPlaneToPoints(facePts, out Plane plane, out double dev);
currentDeviation += dev;
for (int j = 0; j < faceV[i].Length; j++) {
int localID = allNGonVArrayFlip[faceV[i][j]];
Point3d cp = plane.ClosestPoint(mesh.TopologyVertices[faceV[i][j]]);
allW[localID] += 1;
allPts[localID] += cp;
}
}
for (int i = 0; i < allNGonVArray.Length; i++)
{
mesh.TopologyVertices[allNGonVArray[i]] = (Point3f)(allPts[i] / allW[i]);
deviation = currentDeviation;/// faceV.Length;
}
}