using FeMM.Grasshopper.Helpers;
using Grasshopper;
using Grasshopper.Kernel;
using Grasshopper.Kernel.Data;
using Grasshopper.Kernel.Parameters;
using Grasshopper.Kernel.Types;
using Rhino;
using Rhino.Geometry;
using Rhino.Geometry.Intersect;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Threading.Tasks;
namespace FeMM.Grasshopper.Components.Patterning
{
public class GeodesicOffsetsComponent : GH_Component
{
/// <summary>
/// Initializes a new instance of the GeodesicOffsetsComponent class.
/// </summary>
public GeodesicOffsetsComponent()
: base("Geodesic Offsets", "GO", "Generate the offsets of a list of geodesics line", CategoryNameConstants.CATEGORY_FEMM, CategoryNameConstants.SUBCATEGORY_PATTERNING)
{
}
/// <summary>
/// Registers all the input parameters for this component.
/// </summary>
protected override void RegisterInputParams(GH_InputParamManager pManager)
{
pManager.AddMeshParameter("Mesh", "M", "The mesh", GH_ParamAccess.item);
pManager.AddCurveParameter("Geodesics", "Gs", "The geodesics to copy", GH_ParamAccess.list);
int i = pManager.AddIntegerParameter("Side", "S", "The side of the geodesic line where to create the geodesic offset", GH_ParamAccess.item, 1);
Param_Integer par = (Param_Integer)pManager[i];
par.AddNamedValue("Left", -1);
par.AddNamedValue("Right", 1);
pManager.AddNumberParameter("Distance", "D", "The offset distance", GH_ParamAccess.item);
pManager.AddNumberParameter("Tolerance", "T", "The distance tolerance", GH_ParamAccess.item, 1);
pManager.AddNumberParameter("Resolution", "R", "The resolution for generating the geodesics lines", GH_ParamAccess.item, 0);
pManager.AddIntegerParameter("Iterations", "I", "The number of iteration for the geodesic calculation", GH_ParamAccess.item, 1000);
pManager.AddBooleanParameter("Vertical", "V", "Determine if the cutting plane is vertical or contained on starting geodesic", GH_ParamAccess.item, true);
}
/// <summary>
/// Registers all the output parameters for this component.
/// </summary>
protected override void RegisterOutputParams(GH_OutputParamManager pManager)
{
pManager.AddCurveParameter("Geodesics", "Gs", "The resulting geodesic curves", GH_ParamAccess.tree);
pManager.AddMeshParameter("Mesh", "Ms", "The mesh", GH_ParamAccess.list);
pManager.AddMeshParameter("Flatten Mesh", "Fl", "The resulting flatten mesh", GH_ParamAccess.list);
}
protected override void SolveInstance(IGH_DataAccess DA)
{
GH_Mesh mesh = null;
List<GH_Curve> geodesics = [];
int side = 0;
double dist = 0;
double tol = 0;
double res = 0;
int steps = 0;
bool isVertical = false;
if (!DA.GetData(0, ref mesh))
return;
if (!DA.GetDataList(1, geodesics))
return;
if (!DA.GetData(2, ref side))
return;
if (!DA.GetData(3, ref dist))
return;
if (!DA.GetData(4, ref tol))
return;
if (!DA.GetData(5, ref res))
return;
if (!DA.GetData(6, ref steps))
return;
if (!DA.GetData(7, ref isVertical))
return;
const double STEP_SIZE = 0.5;
const int MAX_ITERATIONS = 20;
try
{
List<Curve>[] results = new List<Curve>[geodesics.Count];
Mesh[] internalCuts = new Mesh[geodesics.Count];
Mesh[] flats = new Mesh[geodesics.Count];
// Move to origin
var boundingBox = mesh.Boundingbox;
var toOrigin = (Vector3d)(-boundingBox.Center);
var bufferMesh = mesh.Value.Duplicate() as Mesh;
bufferMesh.Translate(toOrigin);
// Fast flattening configuration
ARAPHelper.Configure(1000, 0, 0, 0, 0.1);
var (meshToCut, brepToCut) = PrepareCutMeshConcurrent(bufferMesh);
//for (int gIndex = 0; gIndex < geodesics.Count; gIndex++)
int parallelismDegree = Environment.ProcessorCount;
parallelismDegree -= Environment.ProcessorCount <= 8 ? 1 : 2;
var options = new ParallelOptions { MaxDegreeOfParallelism = parallelismDegree };
Parallel.For(0, geodesics.Count, options, (gIndex, state) =>
{
GH_Curve geodesic = geodesics[gIndex];
var gBufferMesh = mesh.Value.Duplicate() as Mesh;
//getting refPlane passing on geodesic
Curve gd = geodesic.Value.PullToMesh(gBufferMesh, RhinoDoc.ActiveDoc.ModelAbsoluteTolerance);
Point3d ini = gd.PointAtStart;
Point3d end = gd.PointAtEnd;
Plane refPlane;
if (isVertical)
{
Vector3d dir = end - ini;
dir.Unitize();
if (Math.Abs(dir * Vector3d.ZAxis) > 0.999)
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Error, $"{gIndex}° geodetic: Unable to use vertical with vertical geodesic");
return;
}
Vector3d normal = Vector3d.CrossProduct(dir, Vector3d.ZAxis);
refPlane = new Plane(ini, normal);
}
else
{
Point3d mid = gd.PointAtLength(0.5 * gd.GetLength());
MeshPoint mp = gBufferMesh.ClosestMeshPoint(mid, 0);
if (mp == null)
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Error, $"{gIndex}° geodetic: Unable to find the normal at mid point of geodesic");
return;
}
gBufferMesh.RebuildNormals();
Vector3d normal = gBufferMesh.FaceNormals[mp.FaceIndex];
refPlane = new Plane(ini, end, ini + 100 * normal);
}
bool reverse = false;
if (side == 1)
{
//Right, positive x
if (refPlane.Normal * Vector3d.XAxis < 0)
{
reverse = true;
}
}
else if (side == -1)
{
//Left, negative x
if (refPlane.Normal * Vector3d.XAxis > 0)
{
reverse = true;
}
}
else
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Error, $"{gIndex}° geodetic: NotSupportedException - Side must be Left or Right");
return;
}
if (reverse)
{
refPlane = new Plane(refPlane.Origin, -refPlane.Normal);
}
//finding max width
Point3d maxPoint = Point3d.Unset;
double maxDistance = double.MinValue;
Transform trans = Transform.ChangeBasis(Plane.WorldXY, refPlane);
for (int i = 0; i < gBufferMesh.Vertices.Count; i++)
{
Point3d p = (Point3d)gBufferMesh.Vertices[i];
Point3d localP = trans * p;
if (localP.Z > maxDistance)
{
maxDistance = localP.Z;
maxPoint = p;
}
}
if (maxDistance < 0.05 * dist)
{
//Too small residual length
return;
}
gd.Translate(toOrigin);
ini += toOrigin;
end += toOrigin;
maxPoint += toOrigin;
refPlane.Translate(toOrigin);
int iter = 1;
double planeShift = dist;
double lowerBound = -1;
double upperBound = -1;
do
{
//finding cutting curves
var cuttingCurves = new List<Curve> { gd };
bool exceedMaxDistance = planeShift >= maxDistance + RhinoDoc.ActiveDoc.ModelAbsoluteTolerance;
var result = new List<Curve>();
if (exceedMaxDistance)
{
result = null;
}
else
{
var shifted = new Plane(refPlane.Origin + refPlane.Normal * planeShift, refPlane.Normal);
Polyline[] pls = Intersection.MeshPlane(bufferMesh, shifted);
for (int i = 0; i < pls.Length; i++)
{
Polyline pl = pls[i];
if (res == 0)
res = pl.First.DistanceTo(pl.Last) / 20.0;
var mg = new Common.Geometry.MeshGeodesic(bufferMesh, pl.First, pl.Last, res, STEP_SIZE, steps);
result.Add(MeshGeodesicComponent.ProjectPolylineToMesh(mg.Polyline, bufferMesh));
}
cuttingCurves.AddRange(result);
}
//cutting
//foreach (Curve c in cuttingCurves)
//{
// RhinoDoc.ActiveDoc.Objects.AddCurve(c);
//}
Mesh[] cuts = CutMeshConcurrent(meshToCut, brepToCut, [.. cuttingCurves]);
//Finding the internal piece
Mesh internalCut = null;
for (int i = 0; i < cuts.Length; i++)
{
Mesh cut = cuts[i];
double maximumDistance = RhinoDoc.ActiveDoc.ModelAbsoluteTolerance * 10;
bool touchAll = true;
bool touched;
for (int j = 0; j < cuttingCurves.Count; j++)
{
Curve c = cuttingCurves[j];
touched = false;
for (int k = 0; k < cut.Vertices.Count; k++)
{
Point3d v = (Point3d)cut.Vertices[k];
if (c.ClosestPoint(v, out double t, maximumDistance))
{
//Point3d cp = c.PointAt(t);
//double d = cp.DistanceTo(v);
touched = true;
break;
}
}
if (!touched)
{
touchAll = false;
break;
}
}
if (touchAll)
{
if (exceedMaxDistance)
{
touched = false;
for (int p = 0; p < cut.Vertices.Count; p++)
{
Point3d v = (Point3d)cut.Vertices[p];
if (v.DistanceTo(maxPoint) < maximumDistance)
{
touched = true;
break;
}
}
//RhinoDoc.ActiveDoc.Objects.AddPoint(maxPoint);
if (!touched)
{
touchAll = false;
}
}
if (touchAll)
{
internalCut = cut;
break;
}
}
}
if (internalCut == null)
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Error, $"{gIndex}° geodetic: Unable to find a mesh internal to the cut");
return;
}
//internalCut.Reduce(Convert.ToInt32(internalCut.Faces.Count * 0.9), true, 5, false);
//RhinoDoc.ActiveDoc.Objects.AddMesh(internalCut);
//Flattening
Mesh flat = internalCut.DuplicateMesh();
if (!ARAPHelper.Flattening(internalCut, flat))
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Error, $"{gIndex}° geodetic: Unable to flat the mesh");
return;
}
//Finding the width of the mesh
//Finding the nodes on the first cut, which impose the direction
int firstNode = -1;
int lastNode = -1;
for (int i = 0; i < internalCut.Vertices.Count; i++)
{
var p3d = new Point3d(internalCut.Vertices[i]);
if (firstNode < 0 && p3d.DistanceTo(ini) < tol)
firstNode = i;
if (lastNode < 0 && p3d.DistanceTo(end) < tol)
lastNode = i;
if (firstNode > -1 && lastNode > -1)
break;
}
// Fallback if searching with tol fails.
if (firstNode < 0)
{
firstNode = 0;
var minDistFirstNode = new Point3d(internalCut.Vertices[0]).DistanceTo(end);
for (int i = 1; i < internalCut.Vertices.Count; i++)
{
var p3d = new Point3d(internalCut.Vertices[i]);
var distFirst = p3d.DistanceTo(end);
if (distFirst < minDistFirstNode)
{
firstNode = i;
minDistFirstNode = distFirst;
}
}
}
if (lastNode < 0)
{
lastNode = 0;
var minDistLastNode = new Point3d(internalCut.Vertices[0]).DistanceTo(end);
for (int i = 1; i < internalCut.Vertices.Count; i++)
{
var p3d = new Point3d(internalCut.Vertices[i]);
var distEnd = p3d.DistanceTo(end);
if (distEnd < minDistLastNode)
{
lastNode = i;
minDistLastNode = distEnd;
}
}
}
Vector3d flatDir = flat.Vertices[lastNode] - flat.Vertices[firstNode];
flatDir.Unitize();
Vector3d orto = Vector3d.CrossProduct(flatDir, Vector3d.ZAxis);
double posMin = double.MaxValue;
double posMax = double.MinValue;
for (int i = 0; i < flat.Vertices.Count; i++)
{
Point3d v = (Point3d)flat.Vertices[i];
double pos = new Vector3d(v) * orto;
posMin = Math.Min(posMin, pos);
posMax = Math.Max(posMax, pos);
}
double width = Math.Abs(posMax - posMin);
if (width < 1000 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance)
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Error, $"{gIndex}° geodetic: Dimensione del cut troppo piccola");
return;
}
//RhinoDoc.ActiveDoc.Objects.AddMesh(flat);
if (exceedMaxDistance && width < dist)
{
return;
}
else if (Math.Abs(width - dist) < tol)
{
//RhinoDoc.ActiveDoc.Objects.AddMesh(flat);
//RhinoDoc.ActiveDoc.Objects.AddLine(flat.Vertices[firstNode], flat.Vertices[lastNode]);
foreach (var curve in result)
curve.Translate(-toOrigin);
results[gIndex] = result;
internalCut.Translate(-toOrigin);
internalCuts[gIndex] = internalCut;
flat.Translate(-toOrigin);
flats[gIndex] = flat;
#if DEBUG
AddRuntimeMessage(GH_RuntimeMessageLevel.Remark, $"{gIndex}° geodetic: {iter} iterations");
#endif
break;
}
else if (iter > MAX_ITERATIONS)
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Error, $"{gIndex}° geodetic: Not converged after {MAX_ITERATIONS} iterations");
return;
}
else
{
if (width < dist)
{
lowerBound = planeShift;
}
else
{
upperBound = planeShift;
}
if (lowerBound > 0 && upperBound > 0)
{
planeShift = (lowerBound + upperBound) / 2;
}
else
{
planeShift = planeShift * dist / width;
}
iter++;
}
} while (true == true);
//}
});
DataTree<Curve> resultsTree = new();
for (int i = 0; i < geodesics.Count; i++)
if (results[i] != null)
resultsTree.AddRange(results[i], new GH_Path(i));
DA.SetDataTree(0, resultsTree);
DA.SetDataList(1, internalCuts.ToList());
DA.SetDataList(2, flats.ToList());
}
catch (Exception ex)
{
AddRuntimeMessage(GH_RuntimeMessageLevel.Error, ex.Message);
}
}
public static Mesh[] CutMesh(Mesh mesh, Curve[] geodesics)
{
Mesh bufferMesh = mesh.DuplicateMesh();
//Triangulation of faces
bufferMesh.Faces.ConvertQuadsToTriangles();
bufferMesh.Weld(Math.PI);
bufferMesh.UnifyNormals();
bufferMesh.Normals.ComputeNormals();
bufferMesh.FaceNormals.UnitizeFaceNormals();
bufferMesh.Normals.UnitizeNormals();
//Converting to brep
Brep brep = Brep.CreateFromMesh(bufferMesh, true);
var cutters = new List<Brep>();
//Creating cutters
for (int j = 0; j < geodesics.Length; j++)
{
Curve c = geodesics[j];
if (!c.TryGetPolyline(out Polyline pl))
{
throw new NotSupportedException("Section is not a polyline");
}
double height = 1d / 20d * c.GetLength();
var plus = new Polyline();
var minus = new Polyline();
bool add = true;
double minDistance = c.GetLength() / (pl.Count - 1) / 10;
Point3d prevAdded = Point3d.Unset;
for (int i = 0; i < pl.Count; i++)
{
MeshPoint mpt;
Vector3d normal;
//Avoiding first node and last node, singularities points
if (i == 0)
{
mpt = bufferMesh.ClosestMeshPoint(pl[1], 100 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance);
}
else if (i == pl.Count - 1)
{
mpt = bufferMesh.ClosestMeshPoint(pl[pl.Count - 2], 100 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance);
}
else
{
mpt = bufferMesh.ClosestMeshPoint(pl[i], 100 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance);
}
if (mpt == null)
{
add = false;
break;
}
if (prevAdded == Point3d.Unset || i == pl.Count - 1 || prevAdded.DistanceTo(pl[i]) > minDistance)
{
normal = bufferMesh.FaceNormals[mpt.FaceIndex];
plus.Add(pl[i] + normal * height);
minus.Add(pl[i] - normal * height);
prevAdded = pl[i];
//RhinoDoc.ActiveDoc.Objects.AddPoint(pl[i]);
}
}
//Extending the cutting lines
Curve plusCurve = plus.ToPolylineCurve();
Curve minusCurve = minus.ToPolylineCurve();
plusCurve = plusCurve.Extend(CurveEnd.Both, height, CurveExtensionStyle.Line);
minusCurve = minusCurve.Extend(CurveEnd.Both, height, CurveExtensionStyle.Line);
//RhinoDoc.ActiveDoc.Objects.AddCurve(plusCurve);
//RhinoDoc.ActiveDoc.Objects.AddCurve(minusCurve);
if (add)
{
// Loft sometimes does not work in Rhino 8.
//cutters.AddRange(Brep.CreateFromLoft(new List<Curve>() { plusCurve, minusCurve }, Point3d.Unset, Point3d.Unset, LoftType.Straight, false));
// Sweep seem to be a overkill for cut surfaces.
//var stw = new SweepTwoRail();
//cutters.AddRange(stw.PerformSweep(plusCurve, minusCurve, new LineCurve(plusCurve.PointAtStart, minusCurve.PointAtStart)));
if (plusCurve.PointAtStart.DistanceTo(minusCurve.PointAtStart) >
plusCurve.PointAtStart.DistanceTo(minusCurve.PointAtEnd))
{
minusCurve.Reverse();
}
plusCurve = plusCurve.DuplicateCurve();
minusCurve = minusCurve.DuplicateCurve();
plusCurve.Domain = new Interval(0, 1);
minusCurve.Domain = new Interval(0, 1);
Surface ruled = NurbsSurface.CreateRuledSurface(plusCurve, minusCurve);
if (ruled is not null)
{
var brepRuled = ruled.ToBrep();
if (brepRuled is not null)
cutters.Add(brepRuled);
}
}
}
////Moving points away from cutters to avoid problems
//for (int j = 0; j < cutters.Count; j++)
//{
// Brep b = cutters[j];
// for (int i = 0; i < bufferMesh.Vertices.Count; i++)
// {
// if (b.ClosestPoint(bufferMesh.Vertices[i], out Point3d closestPoint, out ComponentIndex ci, out double s, out double t, 100 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance,
// out Vector3d normal))
// {
// bufferMesh.Vertices.SetVertex(i, closestPoint + 100 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance * normal);
// }
// }
//}
// Rhino8:
// 1e-7 starts working
// 1e-8 does more
// 1e-10 sometimes does not converge
// 1e-12 goes into error
//
// Rhino7:
// 1e-3 work correctly
// 1e-8 sometimes does not converge
double geomTolerance = RhinoDoc.ActiveDoc.ModelAbsoluteTolerance;
if (RhinoApp.Version.Major >= 8)
geomTolerance *= 0.00001;
else
geomTolerance *= 0.1;
Brep[] split_breps = brep.Split(cutters, geomTolerance);
//RhinoDoc.ActiveDoc.Objects.AddBrep(brep);
//foreach (Brep cut in cutters)
//{
// RhinoDoc.ActiveDoc.Objects.AddBrep(cut);
//}
//foreach (Brep split in split_breps)
//{
// RhinoDoc.ActiveDoc.Objects.AddBrep(split);
//}
var result = new List<Mesh>();
for (int k = 0; k < split_breps.Length; k++)
{
Brep b = split_breps[k];
//
// https://developer.rhino3d.com/api/RhinoCommon/html/M_Rhino_Geometry_BrepFace_Split.htm
// https://developer.rhino3d.com/api/RhinoCommon/html/M_Rhino_Geometry_Mesh_CreateFromBrep_1.htm
//
// https://github.com/gradientspace/geometry3Sharp
//
//Testing if have an area
if (b.GetArea() < 1000 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance) continue;
var brep_mesh = new Mesh();
int i = 0;
foreach (BrepFace bf in b.Faces)
{
Curve ol = bf.OuterLoop.To3dCurve();
double t0 = ol.Domain.T0;
double t1 = ol.Domain.T1;
double t = ol.Domain.T0;
var ts = new List<double> { t };
while (ol.GetNextDiscontinuity(Continuity.C1_continuous, t0, t1, out t))
{
ts.Add(t);
t0 = t;
}
Curve[] segs = ol.Split(ts);
List<Point3d> points = [.. segs.Select(s => s.PointAtStart).Distinct()];
//Removing double points
for (int index = points.Count - 1; index >= 1; index--)
{
bool remove;
Point3d p;
remove = false;
p = points[index];
for (int oIndex = index - 1; oIndex >= 0; oIndex--)
{
Point3d op;
op = points[oIndex];
if (p.DistanceTo(op) < 10 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance)
{
remove = true;
break;
}
}
if (remove)
{
points.RemoveAt(index);
}
}
if (points.Count > 2)
{
// Check point colinearity
for (int j = points.Count - 1; j >= 0; j--)
{
Point3d p1, p2;
if (j == 0)
p1 = points[points.Count - 1];
else
p1 = points[j - 1];
if (j == points.Count - 1)
p2 = points[0];
else
p2 = points[j + 1];
Vector3d v1 = points[j] - p1;
Vector3d v2 = p2 - points[j];
double ang = Vector3d.VectorAngle(v1, v2);
if (ang < 0.01)
{
points.RemoveAt(j);
break;
}
}
if (points.Count > 2)
{
int n = brep_mesh.Vertices.Count;
brep_mesh.Vertices.AddVertices(points);
if (points.Count == 3)
{
brep_mesh.Faces.AddFace(n, n + 1, n + 2);
}
else if (points.Count == 4)
{
brep_mesh.Faces.AddFace(n, n + 1, n + 2, n + 3);
}
else
{
Polyline pl;
MeshFace[] mfs;
pl = new Polyline(points)
{
points[0]
};
mfs = pl.TriangulateClosedPolyline();
if (mfs != null)
{
for (int f = 0; f < mfs.Length; f++)
{
MeshFace fs = mfs[f];
if (fs.IsTriangle)
{
brep_mesh.Faces.AddFace(n + fs.A, n + fs.B, n + fs.C);
}
else
{
throw new NotSupportedException();
}
}
}
}
}
}
i++;
}
brep_mesh.Vertices.CombineIdentical(true, true);
brep_mesh.Vertices.Align(10 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance);
brep_mesh.Normals.ComputeNormals();
brep_mesh.Compact();
brep_mesh.Weld(Math.PI);
brep_mesh.Smooth(1, true, true, false, true, SmoothingCoordinateSystem.Object);
brep_mesh.Faces.ConvertQuadsToTriangles();
brep_mesh.CollapseFacesByArea(1000 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance, 0);
for (i = 0; i < 3; i++)
{
brep_mesh.CollapseFacesByEdgeLength(false, 100 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance);
}
brep_mesh.Weld(Math.PI);
brep_mesh.RebuildNormals();
brep_mesh.UnifyNormals();
brep_mesh.Normals.ComputeNormals();
if (brep_mesh.IsValid)
{
result.Add(brep_mesh);
}
}
Vector3d direction = geodesics[0].PointAtEnd - geodesics[0].PointAtStart;
direction.Unitize();
if (Math.Abs(direction * Vector3d.ZAxis) > 0.999)
{
direction = Vector3d.CrossProduct(direction, Vector3d.YAxis);
}
else
{
direction = Vector3d.CrossProduct(direction, Vector3d.ZAxis);
}
direction.Unitize();
if (direction * Vector3d.XAxis < 0)
{
direction *= -1;
}
result.Sort(new MeshComparer(direction));
return [.. result];
}
public static (Mesh, Brep) PrepareCutMeshConcurrent(Mesh mesh)
{
Mesh bufferMesh = mesh.DuplicateMesh();
//Triangulation of faces
bufferMesh.Faces.ConvertQuadsToTriangles();
bufferMesh.Weld(Math.PI);
bufferMesh.UnifyNormals();
bufferMesh.Normals.ComputeNormals();
bufferMesh.FaceNormals.UnitizeFaceNormals();
bufferMesh.Normals.UnitizeNormals();
//Converting to brep
Brep brep = Brep.CreateFromMesh(bufferMesh, true);
return (bufferMesh, brep);
}
public static Mesh[] CutMeshConcurrent(Mesh bufferMesh, Brep brep, Curve[] geodesics)
{
var cutters = new List<Brep>();
//Creating cutters
for (int j = 0; j < geodesics.Length; j++)
{
Curve c = geodesics[j];
if (!c.TryGetPolyline(out Polyline pl))
{
throw new NotSupportedException("Section is not a polyline");
}
double height = 1d / 20d * c.GetLength();
var plus = new Polyline();
var minus = new Polyline();
bool add = true;
double minDistance = c.GetLength() / (pl.Count - 1) / 10;
Point3d prevAdded = Point3d.Unset;
for (int i = 0; i < pl.Count; i++)
{
MeshPoint mpt;
Vector3d normal;
//Avoiding first node and last node, singularities points
if (i == 0)
{
mpt = bufferMesh.ClosestMeshPoint(pl[1], 100 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance);
}
else if (i == pl.Count - 1)
{
mpt = bufferMesh.ClosestMeshPoint(pl[pl.Count - 2], 100 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance);
}
else
{
mpt = bufferMesh.ClosestMeshPoint(pl[i], 100 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance);
}
if (mpt == null)
{
add = false;
break;
}
if (prevAdded == Point3d.Unset || i == pl.Count - 1 || prevAdded.DistanceTo(pl[i]) > minDistance)
{
normal = bufferMesh.FaceNormals[mpt.FaceIndex];
plus.Add(pl[i] + normal * height);
minus.Add(pl[i] - normal * height);
prevAdded = pl[i];
//RhinoDoc.ActiveDoc.Objects.AddPoint(pl[i]);
}
}
//Extending the cutting lines
Curve plusCurve = plus.ToPolylineCurve();
Curve minusCurve = minus.ToPolylineCurve();
plusCurve = plusCurve.Extend(CurveEnd.Both, height, CurveExtensionStyle.Line);
minusCurve = minusCurve.Extend(CurveEnd.Both, height, CurveExtensionStyle.Line);
//RhinoDoc.ActiveDoc.Objects.AddCurve(plusCurve);
//RhinoDoc.ActiveDoc.Objects.AddCurve(minusCurve);
if (add)
{
// Loft sometimes does not work in Rhino 8.
//cutters.AddRange(Brep.CreateFromLoft(new List<Curve>() { plusCurve, minusCurve }, Point3d.Unset, Point3d.Unset, LoftType.Straight, false));
// Sweep seem to be a overkill for cut surfaces.
//var stw = new SweepTwoRail();
//cutters.AddRange(stw.PerformSweep(plusCurve, minusCurve, new LineCurve(plusCurve.PointAtStart, minusCurve.PointAtStart)));
if (plusCurve.PointAtStart.DistanceTo(minusCurve.PointAtStart) >
plusCurve.PointAtStart.DistanceTo(minusCurve.PointAtEnd))
{
minusCurve.Reverse();
}
plusCurve = plusCurve.DuplicateCurve();
minusCurve = minusCurve.DuplicateCurve();
plusCurve.Domain = new Interval(0, 1);
minusCurve.Domain = new Interval(0, 1);
Surface ruled = NurbsSurface.CreateRuledSurface(plusCurve, minusCurve);
if (ruled is not null)
{
var brepRuled = ruled.ToBrep();
if (brepRuled is not null)
cutters.Add(brepRuled);
}
}
}
////Moving points away from cutters to avoid problems
//for (int j = 0; j < cutters.Count; j++)
//{
// Brep b = cutters[j];
// for (int i = 0; i < bufferMesh.Vertices.Count; i++)
// {
// if (b.ClosestPoint(bufferMesh.Vertices[i], out Point3d closestPoint, out ComponentIndex ci, out double s, out double t, 100 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance,
// out Vector3d normal))
// {
// bufferMesh.Vertices.SetVertex(i, closestPoint + 100 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance * normal);
// }
// }
//}
// Rhino8:
// 1e-7 starts working
// 1e-8 does more
// 1e-10 sometimes does not converge
// 1e-12 goes into error
//
// Rhino7:
// 1e-3 work correctly
// 1e-8 sometimes does not converge
double geomTolerance = RhinoDoc.ActiveDoc.ModelAbsoluteTolerance;
if (RhinoApp.Version.Major >= 8)
geomTolerance *= 0.00001;
else
geomTolerance *= 0.1;
Brep[] split_breps = brep.Split(cutters, geomTolerance);
//RhinoDoc.ActiveDoc.Objects.AddBrep(brep);
//foreach (Brep cut in cutters)
//{
// RhinoDoc.ActiveDoc.Objects.AddBrep(cut);
//}
//foreach (Brep split in split_breps)
//{
// RhinoDoc.ActiveDoc.Objects.AddBrep(split);
//}
var result = new List<Mesh>();
for (int k = 0; k < split_breps.Length; k++)
{
Brep b = split_breps[k];
//
// https://developer.rhino3d.com/api/RhinoCommon/html/M_Rhino_Geometry_BrepFace_Split.htm
// https://developer.rhino3d.com/api/RhinoCommon/html/M_Rhino_Geometry_Mesh_CreateFromBrep_1.htm
//
// https://github.com/gradientspace/geometry3Sharp
//
//Testing if have an area
if (b.GetArea() < 1000 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance) continue;
var brep_mesh = new Mesh();
int i = 0;
foreach (BrepFace bf in b.Faces)
{
Curve ol = bf.OuterLoop.To3dCurve();
double t0 = ol.Domain.T0;
double t1 = ol.Domain.T1;
double t = ol.Domain.T0;
var ts = new List<double> { t };
while (ol.GetNextDiscontinuity(Continuity.C1_continuous, t0, t1, out t))
{
ts.Add(t);
t0 = t;
}
Curve[] segs = ol.Split(ts);
List<Point3d> points = [.. segs.Select(s => s.PointAtStart).Distinct()];
//Removing double points
for (int index = points.Count - 1; index >= 1; index--)
{
bool remove;
Point3d p;
remove = false;
p = points[index];
for (int oIndex = index - 1; oIndex >= 0; oIndex--)
{
Point3d op;
op = points[oIndex];
if (p.DistanceTo(op) < 10 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance)
{
remove = true;
break;
}
}
if (remove)
{
points.RemoveAt(index);
}
}
if (points.Count > 2)
{
// Check point colinearity
for (int j = points.Count - 1; j >= 0; j--)
{
Point3d p1, p2;
if (j == 0)
p1 = points[points.Count - 1];
else
p1 = points[j - 1];
if (j == points.Count - 1)
p2 = points[0];
else
p2 = points[j + 1];
Vector3d v1 = points[j] - p1;
Vector3d v2 = p2 - points[j];
double ang = Vector3d.VectorAngle(v1, v2);
if (ang < 0.01)
{
points.RemoveAt(j);
break;
}
}
if (points.Count > 2)
{
int n = brep_mesh.Vertices.Count;
brep_mesh.Vertices.AddVertices(points);
if (points.Count == 3)
{
brep_mesh.Faces.AddFace(n, n + 1, n + 2);
}
else if (points.Count == 4)
{
brep_mesh.Faces.AddFace(n, n + 1, n + 2, n + 3);
}
else
{
Polyline pl;
MeshFace[] mfs;
pl = new Polyline(points)
{
points[0]
};
mfs = pl.TriangulateClosedPolyline();
if (mfs != null)
{
for (int f = 0; f < mfs.Length; f++)
{
MeshFace fs = mfs[f];
if (fs.IsTriangle)
{
brep_mesh.Faces.AddFace(n + fs.A, n + fs.B, n + fs.C);
}
else
{
throw new NotSupportedException();
}
}
}
}
}
}
i++;
}
brep_mesh.Vertices.CombineIdentical(true, true);
brep_mesh.Vertices.Align(10 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance);
brep_mesh.Normals.ComputeNormals();
brep_mesh.Compact();
brep_mesh.Weld(Math.PI);
brep_mesh.Smooth(1, true, true, false, true, SmoothingCoordinateSystem.Object);
brep_mesh.Faces.ConvertQuadsToTriangles();
brep_mesh.CollapseFacesByArea(1000 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance, 0);
for (i = 0; i < 3; i++)
{
brep_mesh.CollapseFacesByEdgeLength(false, 100 * RhinoDoc.ActiveDoc.ModelAbsoluteTolerance);
}
brep_mesh.Weld(Math.PI);
brep_mesh.RebuildNormals();
brep_mesh.UnifyNormals();
brep_mesh.Normals.ComputeNormals();
if (brep_mesh.IsValid)
{
result.Add(brep_mesh);
}
}
Vector3d direction = geodesics[0].PointAtEnd - geodesics[0].PointAtStart;
direction.Unitize();
if (Math.Abs(direction * Vector3d.ZAxis) > 0.999)
{
direction = Vector3d.CrossProduct(direction, Vector3d.YAxis);
}
else
{
direction = Vector3d.CrossProduct(direction, Vector3d.ZAxis);
}
direction.Unitize();
if (direction * Vector3d.XAxis < 0)
{
direction *= -1;
}
result.Sort(new MeshComparer(direction));
return [.. result];
}
private class MeshComparer : IComparer<Mesh>
{
private Vector3d _direction;
public MeshComparer(Vector3d direction)
{
_direction = direction;
}
public int Compare(Mesh x, Mesh y)
{
if (x == null && y == null)
{
return 0;
}
else if (x == null)
{
return 1;
}
else if (y == null)
{
return -1;
}
Point3d centerX = x.GetBoundingBox(false).Center;
Point3d centerY = y.GetBoundingBox(false).Center;
double posX = new Vector3d(centerX) * _direction;
double posY = new Vector3d(centerY) * _direction;
if (Math.Abs(posX - posY) < RhinoDoc.ActiveDoc.ModelAbsoluteTolerance)
{
return 0;
}
else if (posX < posY)
{
return -1;
}
else
{
return 1;
}
}
}
public override GH_Exposure Exposure => GH_Exposure.secondary;
/// <summary>
/// Provides an Icon for the component.
/// </summary>
protected override System.Drawing.Bitmap Icon => Properties.Resources.GeodesicOffsetsIcon;
/// <summary>
/// Gets the unique ID for this component. Do not change this ID after release.
/// </summary>
public override Guid ComponentGuid => new("c1e35ec9-caf9-4dcf-abc1-41b8dcd34a46");
}
}