diff --git a/.gitattributes b/.gitattributes index 2ab7d79..a6e5699 100644 --- a/.gitattributes +++ b/.gitattributes @@ -12,4 +12,8 @@ # Denote files that are truly binary and should not be modified. *.png binary *.jpg binary -*.bmp binary \ No newline at end of file +*.bmp binary + +# Grasshopper documents: .gh is binary (no diff/merge); .ghx is XML text (diffable). +*.gh binary +*.ghx text \ No newline at end of file diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..6b3094d --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,108 @@ +# CLAUDE.md + +This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository. + +## What this is + +Bertini_real is a numerical algorithm for decomposing real algebraic curves (dim 1) and +surfaces (dim 2) into cell complexes, using **Bertini 1** as the homotopy-continuation engine. +The C++ side produces the decomposition; the Python package (`python/`) parses the on-disk +output and renders/plots it. MATLAB code in `matlab_codes/` is legacy post-processing. + +Two C++ executables are built: `bertini_real` (computes the decomposition) and `sampler` +(refines/adaptively samples an existing decomposition for smoother output). + +## Building (C++) + +As of July 2025 the build is **CMake** (previously autotools), and it requires **Bertini 1 +version >= 1.7** (which itself moved to CMake — install location and headers changed). + +```sh +mkdir build && cd build +cmake ../ +make # use -j for parallelism +make install # may need sudo; installs bertini_real + sampler to bin/ +``` + +Dependencies (all must be findable by CMake): MPFR, GMP, `bertini-parallel` (>=1.7, found via +`find_package(bertini1 1.7 CONFIG)`), Boost (>=1.53, components `filesystem` + `timer`), MPI +(openmpi/mpich), plus tools CMake, Flex, Bison. **Bertini 1 must be compiled from source +against the same GMP/MPFR/MPI libraries.** Custom CMake find-modules live in `cmake/`. +`brconfig.h.in` is configured into `build/config.h`. Source-file lists are in `files.cmake` +(not globbed) — **add new `.cpp`/`.hpp` files there**, not just to disk. + +Clone recursively (`git clone --recursive`) to get the `matlab_codes/brakelab` submodule. + +## Python package + +```sh +cd python && pip install -e . # installs the `bertini_real` package +``` +Deps: matplotlib, trimesh, dill, algopy, sympy, scipy, networkx. `glumpy` (OpenGL rendering) +is optional and imported defensively in `__init__.py`. + +Typical interactive use, run **from inside a decomposition output folder**: +```python +import bertini_real +bertini_real.gather_and_plot() # data.gather() -> plot.plot() +``` +`data.gather()` reads the raw on-disk output into a `Curve` or `Surface` object (chosen by the +dimension in the directory name); `data.gather_and_save()` also dills it to a `BRdataN.pkl`. + +## Running / tests + +There is **no automated test suite or test runner**. `test/curve/*` and `test/surface/*` are +example systems, each a directory containing a Bertini `input` file (and sometimes Python +plot/assemble scripts). The manual workflow for any example: + +1. Run Bertini 1 on the `input` file with `tracktype: 1` to produce a numerical irreducible + decomposition and a `witness_data` file. +2. Run `bertini_real` in that directory; it consumes `input` + `witness_data`. If there are + multiple components it prompts for which to decompose. Optional flags: + `-sphere -pi `. +3. Optionally run `sampler` to refine the decomposition. +4. Use the Python package to plot the result. + +Output is plain-text files in a subfolder of cwd, written incrementally after each major stage +(so a crash still leaves the last good parsable state). Key files: `decomp`, `vertex_set`, +copies of `input` + `witness_data`; curves add `E.edge`; surfaces add `S.surf` plus curve +sub-decompositions in their own subfolders. The README and `manual/bertini_real_manual.pdf` +document these formats in full. + +## C++ architecture + +`src/` and `include/` mirror each other and are organized by subsystem (file lists in +`files.cmake`). `src/bertini_real.cpp` and `src/sampler/` hold the two executable `main`s; +everything else compiles into both (`common_src`). + +- **`bertini1/`** — `bertini_extensions`: the C++ bridge to Bertini 1's C structures/headers + (`bertini_headers.hpp`). +- **`nag/`** — numerical algebraic geometry core. `nid` (numerical irreducible decomposition), + `witness_set`, `system_randomizer`, and `nag/solvers/` (the homotopy solvers: `midpoint`, + `multilintolin`, `nullspace`, `sphere_intersection`, `postProcessing`, common `solver`). +- **`decompositions/`** — the top-level algorithms: `curve`, `surface`, base `decomposition`, + and `checkSelfConjugate` (real-vs-complex detection). +- **`cells/`** + **`containers/`** — the cell-complex data model: `vertex`/`edge`/`face`/`cell`, + held in `vertex_set` and `holders`. +- **`symbolics/`** — symbolic preprocessing: `derivative_systems`, `isosingular` (deflation), + `nullspace`, `slicing`, `sphere_intersection`. +- **`io/`** — `fileops`, terminal `color`, and the **Flex** parser `partitionParse.l` (compiled + to `partitionParse.yy.c` at build time; CMake `flex_target` with prefix `partitionParse`). +- top-level: `programConfiguration` (CLI flags / config), `parallelism` (MPI master/worker), + `double_odometer`, `limbo`. + +The program is **MPI-parallel** (head/worker model in `parallelism`). C++14. + +## Python architecture + +`python/bertini_real/` mirrors the C++ cell model in Python objects. `data/` does the parsing +(`gather*`), `curve`/`surface`/`edge`/`face`/`vertex`/`cell`/`decomposition` are the parsed +types, `parse/` reads the directory naming convention, `plot/`/`glumpyplotter`/`anaglypy` do +rendering, `dehomogenize/` handles projective coords, `paths/` and `util/` are helpers. + +## Notes + +- `documentation/` holds the Doxygen config (`bertini_real.doxy.config`) for the C++ docs at + doc.bertinireal.com/cpp; `python/docs/` is the Sphinx source for doc.bertinireal.com/python. +- GitHub Actions only mirrors pushes to an MPI GitLab (`.github/workflows/github-gitlab-sync.yml`); + there is no CI build/test. diff --git a/grasshopper/bertini_real/Capping.cs b/grasshopper/bertini_real/Capping.cs new file mode 100644 index 0000000..181f8c8 --- /dev/null +++ b/grasshopper/bertini_real/Capping.cs @@ -0,0 +1,106 @@ +using System; +using System.Collections.Generic; +using Rhino; +using Rhino.Geometry; + +namespace bertini_real +{ + /// + /// Shared helpers for the cap components (Sphere Caps, Flat Caps): finding the piece mesh's + /// own naked boundary loops that lie on the sphere, and adding non-degenerate triangles. + /// + internal static class Capping + { + // a tiny absolute distance below which two cap vertices are treated as the same point + public const double CoincidentTol = 1e-9; + + /// + /// Closed loops (ordered lists of TopologyVertex indices) of the mesh's naked boundary that + /// lie on the sphere. Sets if some on-sphere boundary did not + /// form clean degree-2 cycles. Works in topology space (coincident vertices merged); zero- + /// length edges collapse and are ignored. + /// + public static List> OnSphereBoundaryLoops(Mesh mesh, Point3d center, double radius, double tol, out bool unclean) + { + unclean = false; + var result = new List>(); + + var topo = mesh.TopologyVertices; + var edges = mesh.TopologyEdges; + + bool OnSphere(int tv) + { + Point3d p = topo[tv]; + return Math.Abs(p.DistanceTo(center) - radius) < tol; + } + + var adj = new Dictionary>(); + void Link(int u, int v) + { + if (!adj.TryGetValue(u, out var lu)) { lu = new List(); adj[u] = lu; } + if (!lu.Contains(v)) lu.Add(v); + } + + for (int e = 0; e < edges.Count; e++) + { + if (edges.GetConnectedFaces(e).Length != 1) continue; // not naked + IndexPair ip = edges.GetTopologyVertices(e); + int i = ip.I, j = ip.J; + if (i == j) continue; // degenerate (collapsed) edge + if (!OnSphere(i) || !OnSphere(j)) continue; + Link(i, j); + Link(j, i); + } + + if (adj.Count == 0) return result; + + var visited = new HashSet(); + foreach (int start in adj.Keys) + { + if (visited.Contains(start)) continue; + + var loop = new List(); + int prev = -1, cur = start; + bool clean = true; + + while (true) + { + visited.Add(cur); + loop.Add(cur); + + var nbrs = adj[cur]; + if (nbrs.Count != 2) { clean = false; break; } // junction or dead-end + + int next = nbrs[0] != prev ? nbrs[0] : nbrs[1]; + if (next == start) break; // closed the loop + if (visited.Contains(next)) { clean = false; break; } + prev = cur; + cur = next; + } + + if (clean && loop.Count >= 3) + result.Add(loop); + else + unclean = true; + } + + return result; + } + + /// + /// Add a triangle, skipping only TRULY degenerate ones (a collapsed edge -- two coincident + /// vertices); thin-but-valid triangles must be kept, or high-resolution caps develop holes. + /// + public static void AddTri(Mesh m, int i, int j, int k) + { + Point3d a = m.Vertices[i]; + Point3d b = m.Vertices[j]; + Point3d c = m.Vertices[k]; + if (a.DistanceTo(b) < CoincidentTol || + b.DistanceTo(c) < CoincidentTol || + a.DistanceTo(c) < CoincidentTol) + return; + m.Faces.AddFace(i, j, k); + } + } +} diff --git a/grasshopper/bertini_real/CurveReadGhJson.cs b/grasshopper/bertini_real/CurveReadGhJson.cs new file mode 100644 index 0000000..de7b771 --- /dev/null +++ b/grasshopper/bertini_real/CurveReadGhJson.cs @@ -0,0 +1,99 @@ +using System; +using Grasshopper; +using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; +using Rhino.Geometry; + +namespace bertini_real +{ + /// + /// Reads a self-contained standalone-curve export (br_gh_export.json) written by Python's + /// Curve.export_gh_json. Brings the vertices in as ONE unified set and exposes each curve + /// piece as a polyline referring to those vertices by index. + /// + public class CurveReadGhJson : GH_Component + { + public CurveReadGhJson() + : base("Curve Read GH JSON", "CurveReadJSON", + "Read a bertini_real curve export: one unified vertex set, curve pieces as polylines", + "bertini_real", "Curve") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddTextParameter("File Path", "F", "Path to br_gh_export.json (a curve export)", GH_ParamAccess.item); + pManager.AddNumberParameter("Scale", "Sc", "Uniform scale applied on import (about the world origin), so the curve comes in bigger without a Scale component", GH_ParamAccess.item, 1.0); + Params.Input[1].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddPointParameter("Vertices", "V", "The single unified vertex set; list index = global vertex id", GH_ParamAccess.list); + pManager.AddCurveParameter("Curves", "C", "One polyline per curve piece", GH_ParamAccess.tree); + pManager.AddTextParameter("Curve Types", "T", "Type tag per curve piece, parallel to Curves", GH_ParamAccess.tree); + pManager.AddIntegerParameter("Curve Indices", "CI", "Per curve piece: vertex indices into Vertices", GH_ParamAccess.tree); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + string path = ""; + if (!DA.GetData(0, ref path)) return; + + double scale = 1.0; + DA.GetData(1, ref scale); + + GhExport content; + try + { + content = GhJsonIO.Load(path); + } + catch (Exception e) + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "Failed to read/parse JSON: " + e.Message); + return; + } + + if (content == null || content.decomposition_type != "curve") + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "Not a curve export (decomposition_type != 'curve')."); + return; + } + + var verts = GhJsonIO.ToVertices(content); + if (scale != 1.0) + for (int i = 0; i < verts.Count; i++) + verts[i] = verts[i] * scale; + + var curves = new DataTree(); + var types = new DataTree(); + var curveIdx = new DataTree(); + + if (content.curve_pieces != null) + { + foreach (var cp in content.curve_pieces) + { + var branch = new GH_Path(cp.piece_index); + + PolylineCurve pl = GhJsonIO.ToPolyline(cp.vertex_indices, verts); + if (pl != null) curves.Add(pl, branch); + + types.Add(cp.type, branch); + if (cp.vertex_indices != null) curveIdx.AddRange(cp.vertex_indices, branch); + } + } + + DA.SetDataList(0, verts); + DA.SetDataTree(1, curves); + DA.SetDataTree(2, types); + DA.SetDataTree(3, curveIdx); + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("telephone.png"); + + public override Guid ComponentGuid + { + get { return new Guid("7C2A9F10-3B8E-4D55-A1C7-6E0B4F92D38B"); } + } + } +} diff --git a/grasshopper/bertini_real/GhJsonIO.cs b/grasshopper/bertini_real/GhJsonIO.cs new file mode 100644 index 0000000..4ced2b7 --- /dev/null +++ b/grasshopper/bertini_real/GhJsonIO.cs @@ -0,0 +1,98 @@ +using System; +using System.Collections.Generic; +using System.IO; +using System.Text.Json; +using Rhino.Geometry; + +namespace bertini_real +{ + /// + /// Shared, Rhino-light helpers for reading the self-contained br_gh_export.json and turning + /// its unified vertex set + index data into Rhino geometry. Kept separate from the GH + /// components so the mapping logic stays small and could be unit-tested later. + /// + internal static class GhJsonIO + { + public static GhExport Load(string path) + { + return JsonSerializer.Deserialize(File.ReadAllText(path)); + } + + /// The single unified vertex set, as Rhino points (index = global vertex id). + public static List ToVertices(GhExport content) + { + var verts = new List(); + if (content?.vertices == null) return verts; + foreach (var v in content.vertices) + { + double x = v != null && v.Length > 0 ? v[0] : 0.0; + double y = v != null && v.Length > 1 ? v[1] : 0.0; + double z = v != null && v.Length > 2 ? v[2] : 0.0; + verts.Add(new Point3d(x, y, z)); + } + return verts; + } + + /// + /// Build a mesh on the FULL shared vertex cloud, so triangle indices stay global and + /// coincident vertices across pieces are recognized as the same point (enabling joins, + /// closed-solid detection, exact curve/mesh intersections). Deliberately does NOT + /// cull/compact vertices -- that would reindex and destroy the shared identity. + /// Degenerate and out-of-range triangles are skipped without dropping vertices. + /// + public static Mesh BuildMesh(GhMesh g, List verts) + { + if (g?.triangles == null) return null; + + var mesh = new Mesh(); + mesh.Vertices.AddVertices(verts); + + int n = verts.Count; + int[] tri = g.triangles; + for (int t = 0; t + 2 < tri.Length; t += 3) + { + int a = tri[t], b = tri[t + 1], c = tri[t + 2]; + if (a < 0 || b < 0 || c < 0 || a >= n || b >= n || c >= n) continue; // out of range + if (a == b || b == c || a == c) continue; // degenerate + mesh.Faces.AddFace(a, b, c); + } + + mesh.Normals.ComputeNormals(); + return mesh; + } + + /// + /// Polyline through the shared vertices, by index. Returns null for <2 valid points + /// (e.g. a nodal singularity), which the caller skips. + /// + public static PolylineCurve ToPolyline(int[] indices, List verts) + { + if (indices == null) return null; + + var pts = new List(); + int n = verts.Count; + foreach (int i in indices) + if (i >= 0 && i < n) pts.Add(verts[i]); + + if (pts.Count < 2) return null; + return new PolylineCurve(pts); + } + + /// + /// The decomposition's bounding sphere as a closed Brep (ready for boolean / capping + /// operations), or null if the sphere data is missing or degenerate. + /// + public static Brep ToSphereBrep(GhSphere s) + { + if (s == null || s.center == null || s.radius <= 0.0) return null; + + double[] c = s.center; + var center = new Point3d( + c.Length > 0 ? c[0] : 0.0, + c.Length > 1 ? c[1] : 0.0, + c.Length > 2 ? c[2] : 0.0); + + return new Sphere(center, s.radius).ToBrep(); + } + } +} diff --git a/grasshopper/bertini_real/PlugNegative.cs b/grasshopper/bertini_real/PlugNegative.cs index fe4079e..6340045 100644 --- a/grasshopper/bertini_real/PlugNegative.cs +++ b/grasshopper/bertini_real/PlugNegative.cs @@ -60,7 +60,7 @@ protected override void SolveInstance(IGH_DataAccess DA) //retrieve inputs if (!DA.GetData(0, ref wireHoleDia)) return; if (!DA.GetData(1, ref plugFactor)) return; - if (!DA.GetData(22, ref socketLength)) return; + if (!DA.GetData(2, ref socketLength)) return; if (!DA.GetData(3, ref lengthOverage)) return; if (!DA.GetData(4, ref bodyOverlap))return; if(!DA.GetData(5, ref eps))return; diff --git a/grasshopper/bertini_real/PlugParts.cs b/grasshopper/bertini_real/PlugParts.cs index 477af41..23c60aa 100644 --- a/grasshopper/bertini_real/PlugParts.cs +++ b/grasshopper/bertini_real/PlugParts.cs @@ -4,52 +4,85 @@ /* - * A Sanity File to hold all the classes used in the code. - * - * Data used in transformConnectors to store data from the JSON - * PieceData hold data on individual piece objects which has been parsed in TransformConnectors - * The rest are used to create the parts of a plug in Connectors.cs or PositivePlugComponent.cs - * plugBody create the main section of the plug with the tapered top - * plugtabs create the cutout box or wedges for the plug - * - * + * A Sanity File to hold helper classes used in the code. + * + * GhExport / GhPiece / GhMesh / GhSphere / GhSingularities / GhEmbeddedCurve / GhCurvePiece + * are the DTOs for br_gh_export.json (see GhJsonIO and the *ReadGhJson components). + * PlugBody / PlugTabs build the parts of a plug in Connectors.cs. */ namespace bertini_real { /// - /// Organize and store data from br_surf_piece_data.json file - /// NOTE: The JSON no longer has piece_indicies, but now piece_names which are file name strings + /// DTOs for the self-contained br_gh_export.json written by Python's + /// Surface.export_gh_json / Curve.export_gh_json. Property names must match the JSON + /// keys exactly (System.Text.Json matches by name). The single `vertices` array is the + /// unified vertex set; meshes carry only triangle indices into it, and curves carry only + /// ordered vertex-index lists into it -- so meshes and embedded curves refer to the same + /// points in Rhino. /// - /// - public class Data + /// + /// + public class GhExport { - public string[] piece_names { get; set; } //this will need to change - public int[][] singularities_on_pieces { get; set; } - public double[][] sing_directions { get; set; } - public double[][] sing_locations { get; set; } - public int[][] parities { get; set; } + public int format_version { get; set; } + public string decomposition_type { get; set; } // "surface" | "curve" + public string source_directory { get; set; } + public int num_variables { get; set; } + public int vertex_count { get; set; } + public double[][] vertices { get; set; } // the unified set; each is [x,y,z] + public GhSphere sphere { get; set; } // bounding sphere of the decomposition + public GhSingularities singularities { get; set; } // nodal-singularity connector data (surface) + public bool is_sampled { get; set; } // surface only + public GhPiece[] pieces { get; set; } // surface only + public GhCurvePiece[] curve_pieces { get; set; } // curve only } - - /// - /// Store data parsed from the Data class by peice - /// - /// - public class PieceData + + public class GhSphere + { + public double[] center { get; set; } // [x,y,z] + public double radius { get; set; } + } + + public class GhSingularities + { + public string[] piece_names { get; set; } // per piece + public double[][] locations { get; set; } // per singularity [x,y,z] + public double[][] directions { get; set; } // per singularity [x,y,z] + public int[][] parities { get; set; } // per singularity: value per piece (-1/0/1) + public int[][] on_pieces { get; set; } // per piece: compact singularity indices + } + + public class GhPiece { - public string piece_name { get; set; } public int piece_index { get; set; } - // public int[] indices { get; set; } - public int[] singsOnPiece { get; set; } - public Vector3d[] directions { get; set; } - public Vector3d[] locations { get; set; } - public int[] parities { get; set; } + public int[] face_indices { get; set; } + public GhMesh mesh_smooth { get; set; } // null when not sampled + public GhMesh mesh_raw { get; set; } + public GhEmbeddedCurve[] curves { get; set; } + } - public PieceData() - { + public class GhMesh + { + public int[] triangles { get; set; } // flat ijk, index into GhExport.vertices + public int triangle_count { get; set; } + } - } + public class GhEmbeddedCurve + { + public string type { get; set; } // critical|sphere|singular|midslice|critslice|unknown + public string curve_name { get; set; } + public int[] vertex_indices { get; set; } // ordered, index into GhExport.vertices } + public class GhCurvePiece + { + public int piece_index { get; set; } + public string type { get; set; } // "standalone" + public string curve_name { get; set; } + public int[] vertex_indices { get; set; } + } + + /* Used to Create the different parts of a positive plug */ public class PlugBody { diff --git a/grasshopper/bertini_real/SurfaceBooleanPiece.cs b/grasshopper/bertini_real/SurfaceBooleanPiece.cs new file mode 100644 index 0000000..e73bdea --- /dev/null +++ b/grasshopper/bertini_real/SurfaceBooleanPiece.cs @@ -0,0 +1,177 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using Grasshopper; +using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; +using Grasshopper.Kernel.Types; +using Rhino.Geometry; + +namespace bertini_real +{ + /// + /// Applies an ordered sequence of mesh boolean operations to each piece solid -- a fold, + /// left to right. Per piece: start from the Solid mesh, then for each Tool in order, union + /// it (sign +1) or subtract it (sign <= 0) from the running result. Order matters: each step + /// acts on the result of the previous one (e.g. union body, then subtract the hole through it). + /// + /// Operations default to subtract, so feeding only negatives "just subtracts" them. Booleans + /// need closed solids: a non-closed Solid is a hard warning. Features may be meshes or Breps + /// (Breps are meshed first). ("Feature" in the solid-modeling sense: an ordered additive or + /// subtractive operation on a body -- a plug body is additive, a wire hole subtractive.) + /// + public class SurfaceBooleanPiece : GH_Component + { + public SurfaceBooleanPiece() + : base("Boolean Piece", "BoolPiece", + "Fold an ordered list of mesh boolean ops (sign +1 union / -1 subtract) onto each piece solid", + "bertini_real", "Surface") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddMeshParameter("Solid", "S", "Closed piece mesh per piece (from Close Piece)", GH_ParamAccess.tree); + pManager.AddGeometryParameter("Features", "F", "Geometry to boolean in, in order, per piece (meshes or Breps)", GH_ParamAccess.tree); + pManager.AddIntegerParameter("Operations", "O", "Sign per feature, parallel to Features: +1 union, -1 subtract. Defaults to subtract.", GH_ParamAccess.tree); + Params.Input[1].Optional = true; // no features -> pass the solid through unchanged + Params.Input[2].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddMeshParameter("Result", "R", "Boolean result per piece", GH_ParamAccess.tree); + pManager.AddBooleanParameter("Is Closed", "X", "Whether each piece's result is a closed solid", GH_ParamAccess.tree); + pManager.AddTextParameter("Report", "Rep", "Per piece: any step that failed", GH_ParamAccess.tree); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + GH_Structure solids; + GH_Structure features; + GH_Structure ops; + if (!DA.GetDataTree(0, out solids)) return; + DA.GetDataTree(1, out features); // optional; no features -> solid passes through + DA.GetDataTree(2, out ops); // optional; default sign is -1 (subtract) + + var resultTree = new DataTree(); + var closedTree = new DataTree(); + var reportTree = new DataTree(); + + for (int b = 0; b < solids.PathCount; b++) + { + GH_Path path = solids.get_Path(b); + + // starting solid(s) + var current = new List(); + foreach (var goo in solids.get_Branch(path)) + if (goo is GH_Mesh gm && gm.Value != null) + current.Add(gm.Value.DuplicateMesh()); + + if (current.Count == 0) continue; + if (current.Any(m => !m.IsClosed)) + AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, + $"Piece at path {path}: solid is not closed; boolean results are unreliable on open meshes."); + + // features + signs for this piece, in order + var featureMeshes = (features != null && features.PathExists(path)) ? ToMeshList(features.get_Branch(path)) : new List(); + var signs = SignsForPath(ops, path, featureMeshes.Count); + + for (int i = 0; i < featureMeshes.Count; i++) + { + Mesh feature = featureMeshes[i]; + bool union = signs[i] > 0; + + if (!feature.IsClosed) + AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, + $"Piece at path {path}: feature {i} is not a closed mesh; the boolean may do nothing. " + + "Check the feature is a closed solid and actually overlaps the piece."); + + Mesh[] next = union + ? Mesh.CreateBooleanUnion(current.Concat(new[] { feature })) + : Mesh.CreateBooleanDifference(current, new[] { feature }); + + if (next == null || next.Length == 0) + { + reportTree.Add($"step {i} ({(union ? "union" : "subtract")}) failed", path); + AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, + $"Piece at path {path}: boolean step {i} ({(union ? "union" : "subtract")}) failed; kept the previous result."); + continue; // keep the last good `current` + } + current = next.ToList(); + } + + foreach (var m in current) + resultTree.Add(m, path); + closedTree.Add(current.All(m => m != null && m.IsClosed), path); + } + + DA.SetDataTree(0, resultTree); + DA.SetDataTree(1, closedTree); + DA.SetDataTree(2, reportTree); + } + + /// Convert a branch of geometry goos to meshes (meshing any Breps). + private static List ToMeshList(System.Collections.IList branch) + { + var meshes = new List(); + foreach (var goo in branch) + { + if (goo is GH_Mesh gm && gm.Value != null) + { + meshes.Add(gm.Value); + continue; + } + GeometryBase geo = (goo as IGH_GeometricGoo)?.IsValid == true + ? GH_Convert.ToGeometryBase(goo) + : null; + if (geo is Mesh m) + { + meshes.Add(m); + } + else if (geo is Brep brep) + { + Mesh[] fromBrep = Mesh.CreateFromBrep(brep, MeshingParameters.Default); + if (fromBrep != null && fromBrep.Length > 0) + { + // CreateFromBrep returns one mesh per face; appending leaves the cap + // seams unwelded (an open mesh), which makes mesh booleans no-op. Weld + // coincident vertices, and if it's still open (an uncapped tube, or a + // non-conforming seam between lateral + caps), fill the holes so the + // cutter becomes a closed solid. + var combined = new Mesh(); + foreach (var mm in fromBrep) combined.Append(mm); + combined.Vertices.CombineIdentical(true, true); + if (!combined.IsClosed) combined.FillHoles(); + combined.RebuildNormals(); + combined.Compact(); + meshes.Add(combined); + } + } + } + return meshes; + } + + /// Signs for a piece's tools; missing/short entries default to -1 (subtract). + private static List SignsForPath(GH_Structure ops, GH_Path path, int count) + { + var signs = new List(); + System.Collections.IList branch = (ops != null && ops.PathExists(path)) ? ops.get_Branch(path) : null; + for (int i = 0; i < count; i++) + { + int s = -1; + if (branch != null && i < branch.Count && branch[i] is GH_Integer gi) + s = gi.Value; + signs.Add(s); + } + return signs; + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("lego.png"); + + public override Guid ComponentGuid + { + get { return new Guid("6F2C1A94-8E3D-4B57-9A0C-2D7E5B41F8C3"); } + } + } +} diff --git a/grasshopper/bertini_real/SurfaceClosePiece.cs b/grasshopper/bertini_real/SurfaceClosePiece.cs new file mode 100644 index 0000000..64eba4f --- /dev/null +++ b/grasshopper/bertini_real/SurfaceClosePiece.cs @@ -0,0 +1,119 @@ +using System; +using System.Collections.Generic; +using Grasshopper; +using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; +using Grasshopper.Kernel.Types; +using Rhino.Geometry; + +namespace bertini_real +{ + /// + /// Joins a surface piece mesh with its spherical cap mesh(es) into a single welded mesh, + /// merging coincident boundary vertices so the result is watertight where the cap meets the + /// piece. Reports whether each joined piece is a closed solid. (A piece bounded only by the + /// sphere closes fully; one that also abuts a singular curve stays open there until joined to + /// its neighbor -- so IsClosed tells you which pieces are complete solids.) + /// + /// The piece meshes and caps must be the same sampling (raw with raw, sampled with sampled); + /// since the caps are built from the piece's own boundary, feeding matching trees preserves that. + /// + public class SurfaceClosePiece : GH_Component + { + public SurfaceClosePiece() + : base("Close Piece", "ClosePiece", + "Join a surface piece with its sphere cap(s) into a welded, ideally closed, mesh", + "bertini_real", "Surface") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddMeshParameter("Meshes", "M", "Surface piece meshes (from Surface Read GH JSON)", GH_ParamAccess.tree); + pManager.AddMeshParameter("Caps", "C", "Sphere cap meshes per piece (from Sphere Caps)", GH_ParamAccess.tree); + pManager.AddNumberParameter("Tolerance", "T", "Vertex merge tolerance for welding cap to piece", GH_ParamAccess.item, 1e-6); + Params.Input[2].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddMeshParameter("Closed", "M", "Piece joined with its cap(s), welded", GH_ParamAccess.tree); + pManager.AddBooleanParameter("Is Closed", "X", "Whether the joined mesh is a closed solid", GH_ParamAccess.tree); + pManager.AddCurveParameter("Naked Edges", "N", "Remaining naked (unwelded/open) edges, for diagnosing why a piece isn't closed", GH_ParamAccess.tree); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + GH_Structure pieces; + GH_Structure caps; + if (!DA.GetDataTree(0, out pieces)) return; + DA.GetDataTree(1, out caps); // caps may legitimately be empty for some pieces + + double tol = 1e-6; + DA.GetData(2, ref tol); + + var outMesh = new DataTree(); + var outClosed = new DataTree(); + var outNaked = new DataTree(); + + for (int b = 0; b < pieces.PathCount; b++) + { + GH_Path path = pieces.get_Path(b); + + var combined = new Mesh(); + + foreach (var goo in pieces.get_Branch(path)) + if (goo is GH_Mesh gm && gm.Value != null) + combined.Append(gm.Value); + + if (caps.PathExists(path)) + foreach (var goo in caps.get_Branch(path)) + if (goo is GH_Mesh gc && gc.Value != null) + combined.Append(gc.Value); + + if (combined.Faces.Count == 0) + continue; + + // weld: merge coincident vertices (the shared cap/piece boundary), drop junk + combined.Vertices.CombineIdentical(true, true); + combined.Faces.CullDegenerateFaces(); + combined.Vertices.CullUnused(); + combined.Compact(); + combined.RebuildNormals(); + combined.UnifyNormals(); + + // UnifyNormals makes the faces mutually consistent but seeds from an arbitrary face, + // so a closed solid can come out uniformly inside-out. Mesh.Volume() is signed by + // normal orientation, so a negative volume means the normals point inward -- flip the + // whole mesh outward (matches the Python pipeline's trimesh.fix_normals()). + bool closed = combined.IsClosed; + if (closed && combined.Volume() < 0.0) + { + combined.Flip(true, true, true); + combined.RebuildNormals(); + } + Polyline[] naked = combined.GetNakedEdges() ?? Array.Empty(); + if (!closed) + AddRuntimeMessage(GH_RuntimeMessageLevel.Remark, + $"Piece at path {path} is not a closed solid: {naked.Length} naked edge loop(s) remain " + + "(see Naked Edges output; likely an unwelded cap seam, an uncapped singular boundary, or a missing cap)."); + + outMesh.Add(combined, path); + outClosed.Add(closed, path); + foreach (var pl in naked) + outNaked.Add(pl, path); + } + + DA.SetDataTree(0, outMesh); + DA.SetDataTree(1, outClosed); + DA.SetDataTree(2, outNaked); + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("lego.png"); + + public override Guid ComponentGuid + { + get { return new Guid("D9C4B6A8-2E51-4F70-A38C-6B1D90E5F273"); } + } + } +} diff --git a/grasshopper/bertini_real/SurfaceColorByFunction.cs b/grasshopper/bertini_real/SurfaceColorByFunction.cs new file mode 100644 index 0000000..edd4694 --- /dev/null +++ b/grasshopper/bertini_real/SurfaceColorByFunction.cs @@ -0,0 +1,163 @@ +using System; +using System.Collections.Generic; +using System.Drawing; +using Grasshopper; +using Grasshopper.GUI.Gradient; +using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; +using Grasshopper.Kernel.Expressions; +using Grasshopper.Kernel.Types; +using Rhino.Geometry; + +namespace bertini_real +{ + /// + /// Colors piece meshes by a scalar function of position, evaluated at every mesh vertex. The + /// Function is an expression in x, y, z (e.g. "x^2 + y^2 + z^2", "Sin(x)*z"); its values are + /// normalized over all the meshes (or to an explicit Domain) and mapped through a color + /// gradient onto the mesh vertex colors. + /// + /// Rhino meshes carry vertex colors, and Spread Pieces preserves them, so this can go before + /// or after spreading; preview the output meshes to see the coloring. + /// + public class SurfaceColorByFunction : GH_Component + { + public SurfaceColorByFunction() + : base("Color By Function", "ColorFn", + "Color piece meshes by a scalar function of x,y,z evaluated at each vertex", + "bertini_real", "Surface") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddMeshParameter("Meshes", "M", "Per-piece meshes to color (one branch per piece)", GH_ParamAccess.tree); + pManager.AddTextParameter("Function", "F", "Scalar expression in x, y, z (e.g. x^2+y^2+z^2)", GH_ParamAccess.item, "z"); + pManager.AddColourParameter("Colours", "Cs", "Gradient stops (>= 2); default is a blue->red spectrum", GH_ParamAccess.list); + pManager.AddIntervalParameter("Domain", "D", "Value range mapped onto the gradient; default = the data's min..max", GH_ParamAccess.item); + Params.Input[1].Optional = true; + Params.Input[2].Optional = true; + Params.Input[3].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddMeshParameter("Meshes", "M", "Colored meshes (vertex colors set), same tree structure", GH_ParamAccess.tree); + pManager.AddNumberParameter("Values", "V", "Per-vertex function values, parallel to each mesh's vertices", GH_ParamAccess.tree); + pManager.AddIntervalParameter("Domain", "D", "The value range used for the gradient", GH_ParamAccess.item); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + GH_Structure meshes; + if (!DA.GetDataTree(0, out meshes)) return; + + string function = "z"; + DA.GetData(1, ref function); + if (string.IsNullOrWhiteSpace(function)) + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "Function expression is empty."); + return; + } + + var colours = new List(); + DA.GetDataList(2, colours); + + Interval domain = Interval.Unset; + bool hasDomain = DA.GetData(3, ref domain); + + // pass 1: evaluate the function at every vertex, tracking the global min/max + var entries = new List<(GH_Path path, Mesh mesh, double[] vals)>(); + double gmin = double.MaxValue, gmax = double.MinValue; + var parser = new GH_ExpressionParser(); + + for (int b = 0; b < meshes.PathCount; b++) + { + GH_Path path = meshes.get_Path(b); + foreach (var goo in meshes.get_Branch(path)) + { + if (!(goo is GH_Mesh gm) || gm.Value == null) continue; + Mesh mesh = gm.Value; + int n = mesh.Vertices.Count; + var vals = new double[n]; + for (int i = 0; i < n; i++) + { + Point3d p = mesh.Vertices[i]; + parser.AddVariable("x", p.X); + parser.AddVariable("y", p.Y); + parser.AddVariable("z", p.Z); + double v; + try { v = parser.Evaluate(function)._Double; } + catch (Exception e) + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "Could not evaluate function: " + e.Message); + return; + } + vals[i] = v; + if (v < gmin) gmin = v; + if (v > gmax) gmax = v; + } + entries.Add((path, mesh, vals)); + } + } + + if (entries.Count == 0) return; + + double lo = hasDomain ? domain.Min : gmin; + double hi = hasDomain ? domain.Max : gmax; + double span = hi - lo; + if (Math.Abs(span) < 1e-15) span = 1.0; + + GH_Gradient gradient = (colours != null && colours.Count >= 2) ? BuildGradient(colours) : DefaultGradient(); + + // pass 2: color each mesh and emit + var outMesh = new DataTree(); + var outValues = new DataTree(); + + foreach (var (path, mesh, vals) in entries) + { + Mesh dup = mesh.DuplicateMesh(); + dup.VertexColors.Clear(); + for (int i = 0; i < vals.Length; i++) + { + double t = (vals[i] - lo) / span; + if (t < 0) t = 0; else if (t > 1) t = 1; + dup.VertexColors.Add(gradient.ColourAt(t)); + } + outMesh.Add(dup, path); + outValues.AddRange(vals, path); + } + + DA.SetDataTree(0, outMesh); + DA.SetDataTree(1, outValues); + DA.SetData(2, new Interval(lo, hi)); + } + + private static GH_Gradient BuildGradient(List colours) + { + var g = new GH_Gradient(); + int n = colours.Count; + for (int i = 0; i < n; i++) + g.AddGrip((double)i / (n - 1), colours[i]); + return g; + } + + private static GH_Gradient DefaultGradient() + { + var g = new GH_Gradient(); + g.AddGrip(0.00, Color.FromArgb(0, 0, 200)); + g.AddGrip(0.25, Color.FromArgb(0, 200, 200)); + g.AddGrip(0.50, Color.FromArgb(0, 200, 0)); + g.AddGrip(0.75, Color.FromArgb(220, 220, 0)); + g.AddGrip(1.00, Color.FromArgb(220, 0, 0)); + return g; + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("import.png"); + + public override Guid ComponentGuid + { + get { return new Guid("3C7F1E08-5B62-4A9D-91C4-7E2A60D5F3B1"); } + } + } +} diff --git a/grasshopper/bertini_real/SurfaceConnectorsToFeatures.cs b/grasshopper/bertini_real/SurfaceConnectorsToFeatures.cs new file mode 100644 index 0000000..622ac02 --- /dev/null +++ b/grasshopper/bertini_real/SurfaceConnectorsToFeatures.cs @@ -0,0 +1,120 @@ +using System; +using System.Collections.Generic; +using Grasshopper; +using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; +using Grasshopper.Kernel.Types; + +namespace bertini_real +{ + /// + /// Weaves the four connector trees from Surface Place Components into a single ordered + /// Features tree plus matching Operations signs, ready to drop into Boolean Piece -- so you + /// don't have to Merge/Weave/sign things by hand. + /// + /// Order, per piece: each connector as (positive +1, then negative -1), plugs then sockets -- + /// the pos, neg, pos, neg sequence. Wire only the negatives and you get an all-subtract + /// Features list (the short-term case). Output Features/Operations are parallel and + /// per-piece ({piece}). + /// + public class SurfaceConnectorsToFeatures : GH_Component + { + public SurfaceConnectorsToFeatures() + : base("Connectors To Features", "Conn2Feat", + "Weave plug/socket connectors into one ordered Features tree + Operations signs for Boolean Piece", + "bertini_real", "Surface") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddGeometryParameter("Plugs positive", "Plugs+", "Positive plugs per piece (union)", GH_ParamAccess.tree); + pManager.AddGeometryParameter("Plugs negative", "Plugs-", "Negative plugs per piece (subtract)", GH_ParamAccess.tree); + pManager.AddGeometryParameter("Sockets positive", "Sockets+", "Positive sockets per piece (union)", GH_ParamAccess.tree); + pManager.AddGeometryParameter("Sockets negative", "Sockets-", "Negative sockets per piece (subtract)", GH_ParamAccess.tree); + Params.Input[0].Optional = true; + Params.Input[1].Optional = true; + Params.Input[2].Optional = true; + Params.Input[3].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddGeometryParameter("Features", "F", "Ordered boolean features per piece (for Boolean Piece)", GH_ParamAccess.tree); + pManager.AddIntegerParameter("Operations", "O", "Sign per feature, parallel to Features: +1 union, -1 subtract", GH_ParamAccess.tree); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + GH_Structure plugPos, plugNeg, socketPos, socketNeg; + DA.GetDataTree(0, out plugPos); + DA.GetDataTree(1, out plugNeg); + DA.GetDataTree(2, out socketPos); + DA.GetDataTree(3, out socketNeg); + + // union of piece indices present across all four inputs + var pieceIndices = new SortedSet(); + foreach (var tree in new[] { plugPos, plugNeg, socketPos, socketNeg }) + CollectPieceIndices(tree, pieceIndices); + + var featuresTree = new DataTree(); + var opsTree = new DataTree(); + + foreach (int pieceIndex in pieceIndices) + { + GH_Path path = new GH_Path(pieceIndex); + + var pp = BranchItems(plugPos, path); + var pn = BranchItems(plugNeg, path); + var sp = BranchItems(socketPos, path); + var sn = BranchItems(socketNeg, path); + + // plugs: (positive +1, negative -1) per connector, then sockets the same way + Weave(pp, pn, path, featuresTree, opsTree); + Weave(sp, sn, path, featuresTree, opsTree); + } + + DA.SetDataTree(0, featuresTree); + DA.SetDataTree(1, opsTree); + } + + /// Emit (positive +1, negative -1) pairs in order; extras (unequal counts) follow with their own sign. + private static void Weave(List positives, List negatives, + GH_Path path, DataTree features, DataTree ops) + { + int n = Math.Max(positives.Count, negatives.Count); + for (int i = 0; i < n; i++) + { + if (i < positives.Count) { features.Add(positives[i], path); ops.Add(+1, path); } + if (i < negatives.Count) { features.Add(negatives[i], path); ops.Add(-1, path); } + } + } + + private static List BranchItems(GH_Structure tree, GH_Path path) + { + var items = new List(); + if (tree != null && tree.PathExists(path)) + foreach (var goo in tree.get_Branch(path)) + if (goo is IGH_GeometricGoo gg && gg.IsValid) + items.Add(gg); + return items; + } + + private static void CollectPieceIndices(GH_Structure tree, SortedSet into) + { + if (tree == null) return; + for (int b = 0; b < tree.PathCount; b++) + { + var idx = tree.get_Path(b).Indices; + into.Add(idx.Length > 0 ? idx[idx.Length - 1] : b); + } + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("transform.png"); + + public override Guid ComponentGuid + { + get { return new Guid("4B8E2D17-5C6A-49F3-A1B0-3E9D7C625A48"); } + } + } +} diff --git a/grasshopper/bertini_real/SurfaceFlatCaps.cs b/grasshopper/bertini_real/SurfaceFlatCaps.cs new file mode 100644 index 0000000..51d383e --- /dev/null +++ b/grasshopper/bertini_real/SurfaceFlatCaps.cs @@ -0,0 +1,164 @@ +using System; +using System.Collections.Generic; +using Grasshopper; +using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; +using Grasshopper.Kernel.Types; +using Rhino.Geometry; + +namespace bertini_real +{ + /// + /// Caps each surface piece's on-sphere boundary loop with a FLAT fan to the loop's centroid, + /// rather than a cap that hugs the sphere (Sphere Caps). The sphere itself still shows where + /// the piece was cut; this just fills the opening flat -- so a Bertini-computed cylinder gets + /// flat disk ends, and blocky decompositions get faceted flat caps. + /// + /// Same on-sphere boundary detection as Sphere Caps; the cap is built on the piece mesh's own + /// boundary vertices (so it welds watertight in Close Piece), with the apex at the loop + /// centroid. Resolution adds concentric (linearly interpolated) rings for a denser flat cap. + /// + public class SurfaceFlatCaps : GH_Component + { + public SurfaceFlatCaps() + : base("Flat Caps", "FlatCaps", + "Cap a surface piece's on-sphere boundary loops with a flat fan to the loop centroid", + "bertini_real", "Surface") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddMeshParameter("Meshes", "M", "Surface piece meshes (from Surface Read GH JSON)", GH_ParamAccess.tree); + pManager.AddBrepParameter("Sphere", "S", "Bounding sphere Brep (from Surface Read GH JSON)", GH_ParamAccess.item); + pManager.AddNumberParameter("Tolerance", "T", "Distance tolerance for testing whether a boundary vertex lies on the sphere", GH_ParamAccess.item, 1e-3); + pManager.AddIntegerParameter("Resolution", "R", "Concentric rings from boundary to centroid (1 = a single flat fan)", GH_ParamAccess.item, 1); + Params.Input[2].Optional = true; + Params.Input[3].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddMeshParameter("Caps", "C", "Flat cap mesh(es) per piece (one per on-sphere loop)", GH_ParamAccess.tree); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + GH_Structure meshes; + if (!DA.GetDataTree(0, out meshes)) return; + + Brep sphereBrep = null; + if (!DA.GetData(1, ref sphereBrep) || sphereBrep == null) + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "No sphere supplied."); + return; + } + + double tol = 1e-3; + DA.GetData(2, ref tol); + + int resolution = 1; + DA.GetData(3, ref resolution); + if (resolution < 1) resolution = 1; + + BoundingBox bb = sphereBrep.GetBoundingBox(true); + Point3d center = bb.Center; + double radius = bb.Diagonal.X / 2.0; + if (radius <= 0) + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "Sphere has non-positive radius."); + return; + } + + var caps = new DataTree(); + + for (int b = 0; b < meshes.PathCount; b++) + { + GH_Path path = meshes.get_Path(b); + + foreach (var goo in meshes.get_Branch(path)) + { + var gm = goo as GH_Mesh; + if (gm?.Value == null) continue; + Mesh mesh = gm.Value; + + List> loops = Capping.OnSphereBoundaryLoops(mesh, center, radius, tol, out bool unclean); + if (unclean) + AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, + $"Piece at path {path}: on-sphere boundary is not all clean closed loops; those parts were not capped."); + + foreach (var loop in loops) + { + Mesh cap = BuildFlatCap(mesh, loop, resolution); + if (cap != null && cap.Faces.Count > 0) + caps.Add(cap, path); + } + } + } + + DA.SetDataTree(0, caps); + } + + /// + /// Flat fan over a closed loop: apex at the loop centroid, with + /// concentric rings linearly interpolated from the (fixed, shared) boundary toward the + /// centroid. Ring 0 keeps the exact boundary vertices so a later join welds watertight. + /// + private static Mesh BuildFlatCap(Mesh mesh, List loop, int resolution) + { + var topo = mesh.TopologyVertices; + int n = loop.Count; + + var pts = new Point3d[n]; + double cx = 0, cy = 0, cz = 0; + for (int k = 0; k < n; k++) + { + pts[k] = topo[loop[k]]; + cx += pts[k].X; cy += pts[k].Y; cz += pts[k].Z; + } + var centroid = new Point3d(cx / n, cy / n, cz / n); + + int R = Math.Max(1, resolution); + var cap = new Mesh(); + + // ring 0: exact boundary; rings 1..R-1: linearly interpolated toward the centroid + for (int k = 0; k < n; k++) cap.Vertices.Add(pts[k]); + for (int r = 1; r < R; r++) + { + double t = (double)r / R; + for (int k = 0; k < n; k++) + cap.Vertices.Add(pts[k] + t * (centroid - pts[k])); + } + int apexIdx = cap.Vertices.Add(centroid); + + int Idx(int r, int k) => r * n + k; + + for (int r = 0; r < R - 1; r++) + { + for (int k = 0; k < n; k++) + { + int k2 = (k + 1) % n; + Capping.AddTri(cap, Idx(r, k), Idx(r, k2), Idx(r + 1, k2)); + Capping.AddTri(cap, Idx(r, k), Idx(r + 1, k2), Idx(r + 1, k)); + } + } + for (int k = 0; k < n; k++) + { + int k2 = (k + 1) % n; + Capping.AddTri(cap, Idx(R - 1, k), Idx(R - 1, k2), apexIdx); + } + + if (cap.Faces.Count == 0) return null; + cap.Normals.ComputeNormals(); + cap.Compact(); + return cap; + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("lego.png"); + + public override Guid ComponentGuid + { + get { return new Guid("1D4A7E62-9C58-4B03-8E71-2F60A9C4D5B8"); } + } + } +} diff --git a/grasshopper/bertini_real/SurfaceGroupByPiece.cs b/grasshopper/bertini_real/SurfaceGroupByPiece.cs new file mode 100644 index 0000000..05618a7 --- /dev/null +++ b/grasshopper/bertini_real/SurfaceGroupByPiece.cs @@ -0,0 +1,106 @@ +using System; +using System.Collections; +using System.Collections.Generic; +using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; +using Grasshopper.Kernel.Types; + +namespace bertini_real +{ + /// + /// Groups each piece's mesh with its connectors into one branch per piece. Everything + /// upstream is already keyed per piece by tree branch (Surface Read GH JSON meshes, + /// Surface Place Components connectors), so this just merges those trees branch-by-branch + /// -- no JSON file, no pieceID-string routing. Each output branch is: the mesh(es) for that + /// piece followed by its connectors. + /// + public class SurfaceGroupByPiece : GH_Component + { + public SurfaceGroupByPiece() + : base("Surface Group By Piece", "SurfGroupPiece", + "Group each piece's mesh with its connectors (one branch per piece)", + "bertini_real", "Surface") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddGeometryParameter("Meshes", "M", "Piece meshes, one branch per piece (from Surface Read GH JSON)", GH_ParamAccess.tree); + pManager.AddGeometryParameter("Plugs positive", "Plugs+", "Positive plugs per piece (from Surface Place Components)", GH_ParamAccess.tree); + pManager.AddGeometryParameter("Plugs negative", "Plugs-", "Negative plugs per piece", GH_ParamAccess.tree); + pManager.AddGeometryParameter("Sockets positive", "Sockets+", "Positive sockets per piece", GH_ParamAccess.tree); + pManager.AddGeometryParameter("Sockets negative", "Sockets-", "Negative sockets per piece", GH_ParamAccess.tree); + + Params.Input[1].Optional = true; + Params.Input[2].Optional = true; + Params.Input[3].Optional = true; + Params.Input[4].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddGeometryParameter("Pieces with connectors", "PCs", + "DataTree, one branch per piece: the mesh(es) followed by that piece's connectors.", + GH_ParamAccess.tree); + pManager.AddIntegerParameter("Piece Indices", "Is", "Piece index of each output branch, in branch order", GH_ParamAccess.list); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + GH_Structure meshes, plugsPos, plugsNeg, socketsPos, socketsNeg; + if (!DA.GetDataTree(0, out meshes)) return; + DA.GetDataTree(1, out plugsPos); + DA.GetDataTree(2, out plugsNeg); + DA.GetDataTree(3, out socketsPos); + DA.GetDataTree(4, out socketsNeg); + + // piece index -> geometry (mesh first because meshes are ingested first) + var grouped = new SortedDictionary>(); + + void ingest(GH_Structure tree) + { + if (tree == null) return; + for (int b = 0; b < tree.PathCount; b++) + { + GH_Path path = tree.get_Path(b); + int pieceIndex = path.Indices.Length > 0 ? path.Indices[path.Indices.Length - 1] : b; + + if (!grouped.TryGetValue(pieceIndex, out var items)) + { + items = new List(); + grouped[pieceIndex] = items; + } + foreach (var goo in tree.get_Branch(path)) + if (goo is IGH_GeometricGoo gg) + items.Add(gg); + } + } + + ingest(meshes); + ingest(plugsPos); + ingest(plugsNeg); + ingest(socketsPos); + ingest(socketsNeg); + + var outTree = new GH_Structure(); + var pieceIndices = new List(); + foreach (var kv in grouped) + { + GH_Path path = new GH_Path(kv.Key); + foreach (var goo in kv.Value) + outTree.Append(goo, path); + pieceIndices.Add(kv.Key); + } + + DA.SetDataTree(0, outTree); + DA.SetDataList(1, pieceIndices); + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("lego.png"); + + public override Guid ComponentGuid + { + get { return new Guid("9094345E-E59C-4E5D-B260-E24CFC11EFC3"); } + } + } +} diff --git a/grasshopper/bertini_real/SurfaceMeshMode.cs b/grasshopper/bertini_real/SurfaceMeshMode.cs new file mode 100644 index 0000000..0a1952b --- /dev/null +++ b/grasshopper/bertini_real/SurfaceMeshMode.cs @@ -0,0 +1,52 @@ +using System; +using Grasshopper.Kernel; + +namespace bertini_real +{ + /// + /// Maps an integer to a Surface Read GH JSON "Mesh Mode" string, so a numeric slider can pick + /// the mode: 0 = auto, 1 = smooth, 2 = raw. Feed the output into the reader's Mesh Mode input. + /// + public class SurfaceMeshMode : GH_Component + { + private static readonly string[] Modes = { "auto", "smooth", "raw" }; + + public SurfaceMeshMode() + : base("Mesh Mode", "MeshMode", + "Pick a Surface Read GH JSON Mesh Mode by index: 0 = auto, 1 = smooth, 2 = raw", + "bertini_real", "Surface") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddIntegerParameter("Index", "i", "0 = auto, 1 = smooth, 2 = raw", GH_ParamAccess.item, 0); + Params.Input[0].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddTextParameter("Mesh Mode", "MM", "Mode string for Surface Read GH JSON", GH_ParamAccess.item); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + int index = 0; + DA.GetData(0, ref index); + + int clamped = Math.Max(0, Math.Min(Modes.Length - 1, index)); + if (clamped != index) + AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, + $"Index {index} out of range [0..{Modes.Length - 1}]; using '{Modes[clamped]}'."); + + DA.SetData(0, Modes[clamped]); + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("import.png"); + + public override Guid ComponentGuid + { + get { return new Guid("2E5D8B30-7A19-4C62-9F41-0B6C3E8A75D2"); } + } + } +} diff --git a/grasshopper/bertini_real/SurfacePlaceComponents.cs b/grasshopper/bertini_real/SurfacePlaceComponents.cs index b9614e7..776345d 100644 --- a/grasshopper/bertini_real/SurfacePlaceComponents.cs +++ b/grasshopper/bertini_real/SurfacePlaceComponents.cs @@ -1,289 +1,204 @@ using System; using System.Collections.Generic; -using System.Configuration; -using System.IO; -using System.Linq; -using System.Runtime.CompilerServices; -using System.Text.Json; +using Grasshopper; using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; using Grasshopper.Kernel.Types; using Rhino.Geometry; -///Component which creates and places connector geometry at singularities on a surface -///Placement is determined by br_piece_data.json which can be created using bertini_real write_piece_data() -///This component accepts 4 different connector geometry and requires at least 1 to run + +///Component which creates and places connector geometry at singularities on a surface. +///Placement comes from the singularity data emitted by "Surface Read GH JSON" (Sing Locations, +///Sing Directions, Sing Parities, Sing On Pieces) -- one JSON, one reader, no second file. +///Accepts up to 4 connector geometries and requires at least 1 to run. namespace bertini_real { public class SurfacePlaceComponents : GH_Component { - /// - /// Initializes a new instance of the MyComponent1 class. - /// public SurfacePlaceComponents() : base("Surface Place Components", "SurfPlaceComps", - "Read the specs for a surface from json, and place components at singularities, etc", + "Place plug/socket connectors at singularities, driven by Surface Read GH JSON's singularity outputs", "bertini_real", "Surface") { } /// - /// Registers all the input parameters for this component. - /// They can be accessed in SolveInstance with DA.GetData() - /// These will appear on the side of the component in which they appear - /// can be accessed in the solve isntance in this order. - /// If you change their order here you MUST change the index used to access them in the SolveInstance - /// DO NOT change the order of these once published/finalized + /// Inputs. The four singularity inputs wire directly from "Surface Read GH JSON". + /// Do NOT reorder once published; the SolveInstance indexes by position. /// protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) { - ///File to get data from. should be br_complete.json which is generated using write_piece() in python - pManager.AddTextParameter("File Path", "F", "Path of json file with surface specs", GH_ParamAccess.item); + pManager.AddPointParameter("Sing Locations", "SL", "Singularity locations (from Surface Read GH JSON)", GH_ParamAccess.list); + pManager.AddVectorParameter("Sing Directions", "SD", "Singularity connector directions (from Surface Read GH JSON)", GH_ParamAccess.list); + pManager.AddIntegerParameter("Sing Parities", "SP", "Per singularity: parity (-1/0/1) on each piece (from Surface Read GH JSON)", GH_ParamAccess.tree); + pManager.AddIntegerParameter("Sing On Pieces", "SOP", "Per piece: indices of the singularities on it (from Surface Read GH JSON)", GH_ParamAccess.tree); - ///Size and location play to adjust connectors. No inputs required by user because it has a default value - pManager.AddNumberParameter("Size", "S", "Scale factor for components", GH_ParamAccess.item, 0.01); - - ///The connector Brep prefabs place. At least one is required for the component to run, but it does not matter which one so all should be optional + // connector geometry is placed at true size -- orient + translate only, no scaling. pManager.AddGeometryParameter("Plug Positive", "Plug+", "Plug positive geometry", GH_ParamAccess.item); pManager.AddGeometryParameter("Plug Negative", "Plug-", "Plug negative geometry", GH_ParamAccess.item); - pManager.AddGeometryParameter("Socket Positive", "Socket+", "Socket positive geometry", GH_ParamAccess.item); pManager.AddGeometryParameter("Socket Negative", "Socket-", "Socket negative geometry", GH_ParamAccess.item); - ///All the geometries should be optional. We check that there is at least 1 geo input in the SolveInstance - Params.Input[1].Optional = true; - Params.Input[2].Optional = true; - Params.Input[3].Optional = true; - Params.Input[4].Optional = true; + + Params.Input[4].Optional = true; // geometries are individually optional; we require >=1 Params.Input[5].Optional = true; + Params.Input[6].Optional = true; + Params.Input[7].Optional = true; } /// - /// Registers all the output parameters for this component. - /// can be set in the SolveInstance using DA.SetData() - /// Appear on the side of the component in the order which they are listed - /// Do NOT change their order once published/finalized - /// if the order is chang you MUST UPDATE their index in the SolveInstance + /// Outputs. One branch per piece for the connector trees and per-piece diagnostics. + /// Do NOT reorder once published. /// protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) { - ///Note: only the connector geometries with an inputted prefab should be sent to output - ///Output a list of the negative Brep connectors transformed to every singularity - pManager.AddGeometryParameter("Plugs positive", "Plugs+", "Pos geos transofrmed", GH_ParamAccess.list); - pManager.AddGeometryParameter("Plugs negative", "Plugs-", "negative transformed geos", GH_ParamAccess.list); - ///Output a list of the positive Brep connectors transformed to every singularity - pManager.AddGeometryParameter("Sockets positive", "Sockets+", "Pos geos transformed", GH_ParamAccess.list); - pManager.AddGeometryParameter("Sockets negative", "Sockets-", "negative transformed", GH_ParamAccess.list); - ///Output a list of the positive Brep connectors transformed to every singularity - pManager.AddTextParameter("Piece filenames", "Fs", "Piece filenames", GH_ParamAccess.list); + pManager.AddGeometryParameter("Plugs positive", "Plugs+", "Positive plugs transformed, per piece", GH_ParamAccess.tree); + pManager.AddGeometryParameter("Plugs negative", "Plugs-", "Negative plugs transformed, per piece", GH_ParamAccess.tree); + pManager.AddGeometryParameter("Sockets positive", "Sockets+", "Positive sockets transformed, per piece", GH_ParamAccess.tree); + pManager.AddGeometryParameter("Sockets negative", "Sockets-", "Negative sockets transformed, per piece", GH_ParamAccess.tree); + + pManager.AddIntegerParameter("Sing Indices per Piece", "SingInds/Piece", "Singularity indices for each piece", GH_ParamAccess.tree); + pManager.AddIntegerParameter("Sing Parities per Piece", "SingPars/Piece", "Singularity parities for each piece", GH_ParamAccess.tree); + pManager.AddVectorParameter("Sing Directions per Piece", "SingDirs/Piece", "Singularity directions for each piece", GH_ParamAccess.tree); + pManager.AddPointParameter("Sing Locations per Piece", "SingLocs/Piece", "Singularity locations for each piece", GH_ParamAccess.tree); } - /// - /// This is the method that actually does the work. - /// GOAL: Take br_piece_data.json as produced from write_piece_data() in bertini_real - /// Produce a pair of connectors at each singularity - /// - /// The DA object is used to retrieve from inputs and store in outputs. protected override void SolveInstance(IGH_DataAccess DA) { - /* Get all inputs - * Create some empty variable - * pass the data from the input parameters to the variables */ - - ///empty variables - Brep plugPos = new Brep(); - Brep plugNeg = new Brep(); - Brep socketPos= new Brep(); - Brep socketNeg = new Brep(); - Double size = 0.01; - Point3d locationPlay = new Point3d(); - List locVectors = new List(); - List dirVectors = new List(); - string jsonPath = ""; - - List transformedPosPlugs = new List(); - List transformedNegPlugs = new List(); - List transformedPosSockets = new List(); - List transformedNegSockets = new List(); - - List piece_filenames = new List(); - - ///The parameters are stored an array. - ///To set a variable to a parameter we need to reference the parameter by its index - ///Do NOT want to change the order of these once published/finalized - if (!DA.GetData(0, ref jsonPath)) return; - DA.GetData(1, ref size); - if (!DA.GetData(2, ref plugPos)) return; - if (!DA.GetData(3, ref plugNeg)) return; - if (!DA.GetData(5, ref socketNeg)) return; - if (!DA.GetData(4, ref socketPos)) return; - ///Error checking inputs. Including a RuntimeMessage in script will automaticall generate an 'o' output on the component + var locations = new List(); + var directions = new List(); + GH_Structure parities; + GH_Structure onPieces; + + DA.GetDataList(0, locations); + DA.GetDataList(1, directions); + if (!DA.GetDataTree(2, out parities)) return; + if (!DA.GetDataTree(3, out onPieces)) return; + + Brep plugPos = BrepFromInput(DA, 4); + Brep plugNeg = BrepFromInput(DA, 5); + Brep socketPos = BrepFromInput(DA, 6); + Brep socketNeg = BrepFromInput(DA, 7); + + // any subset of the four connector geometries is fine (none through all); we simply + // place whatever is supplied and leave the rest empty. + + var plugsPos = new GH_Structure(); + var plugsNeg = new GH_Structure(); + var socketsPos = new GH_Structure(); + var socketsNeg = new GH_Structure(); + var indicesPerPiece = new GH_Structure(); + var paritiesPerPiece = new GH_Structure(); + var dirsPerPiece = new GH_Structure(); + var locsPerPiece = new GH_Structure(); + + // one branch per piece, from the Sing On Pieces tree + for (int b = 0; b < onPieces.PathCount; b++) + { + GH_Path piecePath = onPieces.get_Path(b); + int pieceIndex = piecePath.Indices.Length > 0 ? piecePath.Indices[piecePath.Indices.Length - 1] : b; + + plugsPos.EnsurePath(piecePath); + plugsNeg.EnsurePath(piecePath); + socketsPos.EnsurePath(piecePath); + socketsNeg.EnsurePath(piecePath); + + foreach (var goo in onPieces.get_Branch(piecePath)) + { + if (!(goo is GH_Integer gi)) continue; + int singIndex = gi.Value; - ///The component should not run in there are no prefab geometries - if ((!plugNeg.IsValid && plugPos.IsValid) || (!socketNeg.IsValid && socketPos.IsValid)) { - this.AddRuntimeMessage(GH_RuntimeMessageLevel.Remark, "Only positive geos inputted, ensure matching negative connectors are placed before combining with piece!"); - } //Remind the user if they only have positive geometries inputted that they will need negative geos if they want to combine with piece - else if(!plugNeg.IsValid && !plugPos.IsValid && !socketNeg.IsValid && !socketPos.IsValid) { - this.AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "At least one of the plug/socket pos/neg geometries is invalid!"); - return; - } + if (singIndex < 0 || singIndex >= locations.Count || singIndex >= directions.Count) + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, $"Singularity index {singIndex} out of range."); + continue; + } - /* Read and Parse the JSON File into a Data Object (defined in PlugParts.cs) */ - string text = File.ReadAllText(jsonPath); - ///this parses the JSON by key. The Data class must have properties the same name as the keys in the JSON file - ///Should Eventually include some runtimeMessage error handeling - var content = JsonSerializer.Deserialize(text); - ///JSON file is structured: - // { - // "piece_names":['filename1.stl', 'filename2.stl', ...], - // "singularities_on_pieces"[[0,1],[2,3],[...], ...] - // "sing_directions": [[vector for sing 0], [vec for sing 1], ... [vec for sing -1]], - // "sing_locations": [[loc for sing 0], [loc for sing 1], ... [loc for sing -1]], - // "parities": [[parities for singularity 0 on its pieces], [for sing 1], ..., [for sing -1]] - // } - // - ///where each property is a list of N lists, where N is the number of pieces. - ///Each list in a property corresponds to the property of the piece - ///each piece is defined by the properties at the same index in each property list - ///ex. piece 2 has piece_indices[1], singularities_on_piece[1], sing_directions[1], etc - ///This is a silly way of using JSON because now we need to sort the JSON into each piece - ///How we make the JSON in bertini_real write_piece_data really should be rewritten to organize by Sing or Piece where each Sing (or piece) has properties - - /* parse JSON data Piece objects */ - ///list for all the pieces - List allPieces = new List(); - ///each piece is represented by a list of indices. the number of peices = length of piece_indices - for (int piece_index= 0; piece_index < content.piece_names.Length; piece_index++) { - - string pieceName = content.piece_names[piece_index]; + int parity = ParityFor(parities, singIndex, pieceIndex); + Vector3d direction = directions[singIndex]; + Point3d location = locations[singIndex]; - piece_filenames.Add(pieceName); - //⚠️I would like to try just pass content to PieceData and have it do the work for me! - PieceData newPiece = new PieceData(); - newPiece.piece_name = pieceName; - // newPiece.indices = content.piece_indices[pieceName]; - newPiece.singsOnPiece = content.singularities_on_pieces[piece_index]; - - //there are vectors for each sing on the piece, need to turn the vectors from vectors into lists - //also append the vectors to the direction and location vector lists - for (int j = 0; j this piece gets the plug; -1 -> it gets the socket + if (parity == 1) { - ///add the direction and location vectors for this plug to the list of all vectors and locations - ///these lists are now unused and can be deleted, but I am keeping them for debugging - dirVectors.Add(dirVect); - locVectors.Add(locVect); - - if (plugNeg.IsValid) { - ///Create a new plug at this location and add it to the plug list - transformedNegPlugs.Add(moveComponents(newPiece.piece_name,locationPlay,size,dirVect,locVect,plugNeg)); - } - if (plugPos.IsValid) { - ///Create a new plug at this location and add it to the plug list - transformedPosPlugs.Add(moveComponents(newPiece.piece_name, locationPlay, size, dirVect, locVect, plugPos)); - } + if (plugNeg != null) + plugsNeg.Append(new GH_Brep(moveComponents(pieceIndex.ToString(), direction, new Vector3d(location), plugNeg)), piecePath); + if (plugPos != null) + plugsPos.Append(new GH_Brep(moveComponents(pieceIndex.ToString(), direction, new Vector3d(location), plugPos)), piecePath); } - - //add sockets if negative parity - else if (content.parities[singIndex][piece_index] == -1) + else if (parity == -1) { - ///add the direction and location vectors for this socket to the list of all vectors and locations - ///these lists are now unused and can be deleted, but I am keeping them for debugging - dirVectors.Add(dirVect); - locVectors.Add(locVect); - if (socketNeg.IsValid) - { - ///Create a new plug at this location and add it to the plug list - transformedNegSockets.Add(moveComponents(newPiece.piece_name, locationPlay, size, dirVect, locVect, socketNeg)); - } - if (socketPos.IsValid) - { - ///Create a new plug at this location and add it to the plug list - transformedPosSockets.Add(moveComponents(newPiece.piece_name, locationPlay, size, dirVect, locVect, socketPos)); - } + if (socketNeg != null) + socketsNeg.Append(new GH_Brep(moveComponents(pieceIndex.ToString(), direction, new Vector3d(location), socketNeg)), piecePath); + if (socketPos != null) + socketsPos.Append(new GH_Brep(moveComponents(pieceIndex.ToString(), direction, new Vector3d(location), socketPos)), piecePath); } } } - - - /* Set output data - * Set to the list of Geos - * 0 - Out must be text - * 1 - negComponents List - * 2 - posCompoents List */ - DA.SetDataList(0, transformedPosPlugs); - DA.SetDataList(1, transformedNegPlugs); - DA.SetDataList(2, transformedPosSockets); - DA.SetDataList(3, transformedNegSockets); - DA.SetDataList(4, piece_filenames); + DA.SetDataTree(0, plugsPos); + DA.SetDataTree(1, plugsNeg); + DA.SetDataTree(2, socketsPos); + DA.SetDataTree(3, socketsNeg); + DA.SetDataTree(4, indicesPerPiece); + DA.SetDataTree(5, paritiesPerPiece); + DA.SetDataTree(6, dirsPerPiece); + DA.SetDataTree(7, locsPerPiece); + } + + /// Parity (-1/0/1) of a singularity on a given piece, from the per-singularity parity tree. + private static int ParityFor(GH_Structure parities, int singIndex, int pieceIndex) + { + GH_Path path = new GH_Path(singIndex); + if (!parities.PathExists(path)) return 0; + var branch = parities.get_Branch(path); + if (pieceIndex < 0 || pieceIndex >= branch.Count) return 0; + return (branch[pieceIndex] as GH_Integer)?.Value ?? 0; + } + + /// Pull a Brep from an item geometry input, or null if absent/unconvertible. + private static Brep BrepFromInput(IGH_DataAccess DA, int index) + { + IGH_GeometricGoo goo = null; + if (!DA.GetData(index, ref goo) || goo == null) return null; + var geo = GH_Convert.ToGeometryBase(goo); + return geo as Brep; } /// - /// Helper function which makes a new connector and places it at the singularity on the piece*/ + /// Make a new connector and place it at a singularity: orient to the direction + /// (phi about Y, theta about Z) and translate to the location. The connector keeps its + /// own (true) size -- no scaling here; size it upstream. /// - /// User input which changes the distance of the connector from the singularity - /// User input for scaling of the connector - /// Direction vector points from the center of the piece to the singularity - /// Location of the singularity where the connector belongs - /// Geometry of the connector prefab to be created - /// A new connector Brep at a singularity - private Brep moveComponents(string piece_name, Point3d locationPlay, double size, Vector3d direction, Vector3d location, Brep geo) { - ///create a new connector - Brep newConnector = geo.DuplicateBrep(); - - ///find our angles + private Brep moveComponents(string piece_name, Vector3d direction, Vector3d location, Brep geo) + { + Brep newConnector = geo.DuplicateBrep(); + double phi = Math.Acos(direction[2] / direction.Length); double theta = Math.Atan2(direction[1], direction[0]); - ///create some transformation matricies and then tranform the connector - var sf = Transform.Scale(Point3d.Origin + locationPlay, size); - var rf = Transform.Rotation(phi, Vector3d.YAxis, Point3d.Origin); - - newConnector.Transform(sf); + var rf = Transform.Rotation(phi, Vector3d.YAxis, Point3d.Origin); newConnector.Transform(rf); - - rf = Transform.Unset; //clear the rotation matrix to be reused - + rf = Transform.Rotation(theta, Vector3d.ZAxis, Point3d.Origin); newConnector.Transform(rf); var xf = Transform.Translation(location); newConnector.Transform(xf); - ///add user data so the piece the the connector is attached to can later be identified newConnector.SetUserString("pieceID", piece_name); - ///send back the connector return newConnector; - } - - /// - /// Provides an Icon for the component. - /// + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("transform.png"); - /// - /// Gets the unique ID for this component. Do not change this ID after release. - /// public override Guid ComponentGuid { get { return new Guid("E636CDFC-C219-49A7-999A-06E91DE10B94"); } } } } - -///How to Importing STL -///String stlPiecePath = stlPath + "\\br_piece_smooth_"+piece.indices[0]+ "-" +piece.indices[1]+"-" +piece.indices[2]+"_solid.stl"; -///DA.SetData(4, stlPiecePath); - -///ActiveDoc.Import(stlPiecePath); \ No newline at end of file diff --git a/grasshopper/bertini_real/SurfaceReadGhJson.cs b/grasshopper/bertini_real/SurfaceReadGhJson.cs new file mode 100644 index 0000000..f96d8c0 --- /dev/null +++ b/grasshopper/bertini_real/SurfaceReadGhJson.cs @@ -0,0 +1,196 @@ +using System; +using System.Collections.Generic; +using Grasshopper; +using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; +using Rhino.Geometry; + +namespace bertini_real +{ + /// + /// Reads a self-contained surface export (br_gh_export.json) written by Python's + /// Surface.export_gh_json. Brings the vertices in as ONE unified set, then exposes each + /// nonsingular piece as a mesh and the curve pieces embedded on it as polylines -- all + /// referring to the same vertices by index. + /// + public class SurfaceReadGhJson : GH_Component + { + public SurfaceReadGhJson() + : base("Surface Read GH JSON", "SurfReadJSON", + "Read a bertini_real surface export: one unified vertex set, pieces as meshes, embedded curves as polylines", + "bertini_real", "Surface") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddTextParameter("File Path", "F", "Path to br_gh_export.json (a surface export)", GH_ParamAccess.item); + pManager.AddTextParameter("Mesh Mode", "MM", "auto | smooth | raw (auto = smooth when sampled, else raw)", GH_ParamAccess.item, "auto"); + pManager.AddNumberParameter("Scale", "Sc", "Uniform scale applied on import (about the world origin), so the whole model -- meshes, curves, sphere, singularities -- comes in bigger without a Scale component", GH_ParamAccess.item, 1.0); + Params.Input[1].Optional = true; + Params.Input[2].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddPointParameter("Vertices", "V", "The single unified vertex set; list index = global vertex id", GH_ParamAccess.list); + pManager.AddMeshParameter("Meshes", "M", "One mesh per nonsingular piece, built on the shared vertices", GH_ParamAccess.tree); + pManager.AddIntegerParameter("Mesh Faces", "MF", "Per piece: flat triangle indices into Vertices", GH_ParamAccess.tree); + pManager.AddCurveParameter("Curves", "C", "Embedded curve pieces per surface piece", GH_ParamAccess.tree); + pManager.AddTextParameter("Curve Types", "T", "Type tag per curve (critical/sphere/singular/midslice/critslice), parallel to Curves", GH_ParamAccess.tree); + pManager.AddIntegerParameter("Curve Indices", "CI", "Per curve: vertex indices into Vertices (path {piece, curve})", GH_ParamAccess.tree); + pManager.AddIntegerParameter("Face Indices", "FI", "Global surface face ids per piece", GH_ParamAccess.tree); + pManager.AddBrepParameter("Sphere", "S", "Bounding sphere of the decomposition as a closed Brep", GH_ParamAccess.item); + pManager.AddPointParameter("Sing Locations", "SL", "Nodal singularity locations (one per singularity)", GH_ParamAccess.list); + pManager.AddVectorParameter("Sing Directions", "SD", "Nodal singularity connector axis directions (one per singularity)", GH_ParamAccess.list); + pManager.AddIntegerParameter("Sing Parities", "SP", "Per singularity: parity (-1/0/1) on each piece", GH_ParamAccess.tree); + pManager.AddIntegerParameter("Sing On Pieces", "SOP", "Per piece: indices of the singularities on it", GH_ParamAccess.tree); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + string path = ""; + string mode = "auto"; + if (!DA.GetData(0, ref path)) return; + DA.GetData(1, ref mode); + mode = (mode ?? "auto").Trim().ToLowerInvariant(); + + double scale = 1.0; + DA.GetData(2, ref scale); + + GhExport content; + try + { + content = GhJsonIO.Load(path); + } + catch (Exception e) + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "Failed to read/parse JSON: " + e.Message); + return; + } + + if (content == null || content.decomposition_type != "surface") + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "Not a surface export (decomposition_type != 'surface')."); + return; + } + + var verts = GhJsonIO.ToVertices(content); + // scale the unified vertex set about the world origin; meshes and curves built from + // these points then come out scaled automatically (indices are unaffected) + if (scale != 1.0) + for (int i = 0; i < verts.Count; i++) + verts[i] = verts[i] * scale; + + var meshes = new DataTree(); + var meshFaces = new DataTree(); + var curves = new DataTree(); + var types = new DataTree(); + var curveIdx = new DataTree(); + var faceIdx = new DataTree(); + + if (content.pieces != null) + { + foreach (var piece in content.pieces) + { + var branch = new GH_Path(piece.piece_index); + + GhMesh chosen = PickMesh(piece, mode); + if (chosen?.triangles != null) + { + Mesh m = GhJsonIO.BuildMesh(chosen, verts); + if (m != null) meshes.Add(m, branch); + meshFaces.AddRange(chosen.triangles, branch); + } + else + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, $"Piece {piece.piece_index} has no usable mesh."); + } + + if (piece.face_indices != null) + faceIdx.AddRange(piece.face_indices, branch); + + if (piece.curves != null) + { + int ord = 0; + foreach (var c in piece.curves) + { + PolylineCurve pl = GhJsonIO.ToPolyline(c.vertex_indices, verts); + if (pl == null) continue; // <2 points (e.g. nodal singularity) + + curves.Add(pl, branch); + types.Add(c.type, branch); + curveIdx.AddRange(c.vertex_indices ?? Array.Empty(), new GH_Path(piece.piece_index, ord)); + ord++; + } + } + } + } + + DA.SetDataList(0, verts); + DA.SetDataTree(1, meshes); + DA.SetDataTree(2, meshFaces); + DA.SetDataTree(3, curves); + DA.SetDataTree(4, types); + DA.SetDataTree(5, curveIdx); + DA.SetDataTree(6, faceIdx); + + Brep sphere = GhJsonIO.ToSphereBrep(content.sphere); + if (sphere != null && scale != 1.0) + sphere.Transform(Transform.Scale(Point3d.Origin, scale)); + if (sphere != null) + DA.SetData(7, sphere); + else + AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, "No valid bounding sphere in the export."); + + // singularity / connector data (folded into the same JSON) + var singLocations = new List(); + var singDirections = new List(); + var singParities = new DataTree(); + var singOnPieces = new DataTree(); + + var sg = content.singularities; + if (sg != null) + { + if (sg.locations != null) + foreach (var p in sg.locations) + if (p != null && p.Length >= 3) + singLocations.Add(new Point3d(p[0] * scale, p[1] * scale, p[2] * scale)); + + if (sg.directions != null) + foreach (var d in sg.directions) + if (d != null && d.Length >= 3) + singDirections.Add(new Vector3d(d[0], d[1], d[2])); + + if (sg.parities != null) + for (int s = 0; s < sg.parities.Length; s++) + if (sg.parities[s] != null) + singParities.AddRange(sg.parities[s], new GH_Path(s)); + + if (sg.on_pieces != null) + for (int pc = 0; pc < sg.on_pieces.Length; pc++) + if (sg.on_pieces[pc] != null) + singOnPieces.AddRange(sg.on_pieces[pc], new GH_Path(pc)); + } + + DA.SetDataList(8, singLocations); + DA.SetDataList(9, singDirections); + DA.SetDataTree(10, singParities); + DA.SetDataTree(11, singOnPieces); + } + + private static GhMesh PickMesh(GhPiece piece, string mode) + { + if (mode == "raw") return piece.mesh_raw; + if (mode == "smooth") return piece.mesh_smooth ?? piece.mesh_raw; + return piece.mesh_smooth ?? piece.mesh_raw; // auto + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("lego.png"); + + public override Guid ComponentGuid + { + get { return new Guid("B5E3C7A2-1D4F-4E8A-9C6B-2F7A0D9E13A4"); } + } + } +} diff --git a/grasshopper/bertini_real/SurfaceSphereCaps.cs b/grasshopper/bertini_real/SurfaceSphereCaps.cs new file mode 100644 index 0000000..1b526c0 --- /dev/null +++ b/grasshopper/bertini_real/SurfaceSphereCaps.cs @@ -0,0 +1,252 @@ +using System; +using System.Collections.Generic; +using Grasshopper; +using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; +using Grasshopper.Kernel.Types; +using Rhino; +using Rhino.Geometry; + +namespace bertini_real +{ + /// + /// Builds the spherical "cap" mesh(es) that close a surface piece where it meets the + /// decomposition's bounding sphere. The cap is built on the piece mesh's OWN naked boundary + /// loop that lies on the sphere -- never on the separately-sampled sphere curve -- so a raw + /// piece gets a raw cap and a sampled piece a sampled cap (raw with raw, sampled with sampled), + /// and the boundary vertices are shared so a later join welds watertight. + /// + /// For each closed on-sphere boundary loop, two candidate fan caps (one to each pole of the + /// loop's mean direction) are built and the smaller-area one is kept. Works in mesh topology + /// space so coincident vertices (nodal singularities) merge, and skips degenerate edges/faces. + /// + public class SurfaceSphereCaps : GH_Component + { + public SurfaceSphereCaps() + : base("Sphere Caps", "SphereCaps", + "Cap a surface piece's on-sphere boundary loops with faceted spherical caps (smaller-area side)", + "bertini_real", "Surface") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddMeshParameter("Meshes", "M", "Surface piece meshes (from Surface Read GH JSON)", GH_ParamAccess.tree); + pManager.AddBrepParameter("Sphere", "S", "Bounding sphere Brep (from Surface Read GH JSON)", GH_ParamAccess.item); + pManager.AddNumberParameter("Tolerance", "T", "Distance tolerance for testing whether a boundary vertex lies on the sphere", GH_ParamAccess.item, 1e-3); + pManager.AddIntegerParameter("Resolution", "R", "Radial subdivisions of the cap (rings from boundary to pole); higher = smoother, follows the sphere", GH_ParamAccess.item, 4); + Params.Input[2].Optional = true; + Params.Input[3].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddMeshParameter("Caps", "C", "Spherical cap mesh(es) per piece (one per on-sphere loop)", GH_ParamAccess.tree); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + GH_Structure meshes; + if (!DA.GetDataTree(0, out meshes)) return; + + Brep sphereBrep = null; + if (!DA.GetData(1, ref sphereBrep) || sphereBrep == null) + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "No sphere supplied."); + return; + } + + double tol = 1e-3; + DA.GetData(2, ref tol); + + int resolution = 4; + DA.GetData(3, ref resolution); + if (resolution < 1) resolution = 1; + + // recover center/radius from the sphere Brep's bounding box (exact for a sphere) + BoundingBox bb = sphereBrep.GetBoundingBox(true); + Point3d center = bb.Center; + double radius = bb.Diagonal.X / 2.0; + if (radius <= 0) + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "Sphere has non-positive radius."); + return; + } + + var caps = new DataTree(); + + for (int b = 0; b < meshes.PathCount; b++) + { + GH_Path path = meshes.get_Path(b); + var branch = meshes.get_Branch(path); + + foreach (var goo in branch) + { + var gm = goo as GH_Mesh; + if (gm?.Value == null) continue; + Mesh mesh = gm.Value; + + List> loops = Capping.OnSphereBoundaryLoops(mesh, center, radius, tol, out bool unclean); + if (unclean) + AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, + $"Piece at path {path}: on-sphere boundary is not all clean closed loops (open arc or non-manifold junction); those parts were not capped."); + + foreach (var loop in loops) + { + Mesh cap = BuildSmallerCap(mesh, loop, center, radius, tol, resolution); + if (cap != null && cap.Faces.Count > 0) + caps.Add(cap, path); + } + } + } + + DA.SetDataTree(0, caps); + } + + /// + /// Builds the smaller-area cap over a closed loop of topology-vertex indices. The pole is + /// placed on the sphere along the loop's mean direction; the cap is subdivided into + /// radial rings, each slerped along the sphere from the + /// (fixed, shared) boundary toward the pole, so it follows the sphere instead of pinching + /// to a flat point. Ring 0 keeps the exact boundary vertices so a later join welds + /// watertight. Degenerate (near-zero-area) triangles are skipped. + /// + private static Mesh BuildSmallerCap(Mesh mesh, List loop, Point3d center, double radius, double tol, int resolution) + { + var topo = mesh.TopologyVertices; + int n = loop.Count; + + var pts = new Point3d[n]; // exact boundary points (ring 0) + var dirs = new Vector3d[n]; // unit directions from center + var mean = new Vector3d(0, 0, 0); + for (int k = 0; k < n; k++) + { + pts[k] = topo[loop[k]]; + Vector3d d = pts[k] - center; + dirs[k] = d; + dirs[k].Unitize(); + mean += d; + } + if (mean.IsTiny()) mean = new Vector3d(0, 0, 1); + mean.Unitize(); + + // pick the pole side that yields the smaller cap + int sign = SmallerCapSign(pts, center, radius, mean); + Vector3d poleDir = sign * mean; + Point3d pole = center + radius * poleDir; + + // The loop is traversed in an arbitrary direction, so the natural winding below may + // come out facing the sphere center (inward). The cap closes a solid that sits INSIDE + // the sphere, so its outward normal must point radially away from the center. Reverse + // the loop order when the winding would be inward, so the cap is outward by construction + // and agrees with the (outward) piece -- no seam fold for the later UnifyNormals to fix. + if (CapWindsInward(pts, pole, center)) + { + Array.Reverse(pts); + Array.Reverse(dirs); + } + + int R = Math.Max(1, resolution); + var cap = new Mesh(); + + // ring 0: exact boundary; rings 1..R-1: slerped toward the pole + for (int k = 0; k < n; k++) cap.Vertices.Add(pts[k]); + for (int r = 1; r < R; r++) + { + double t = (double)r / R; + for (int k = 0; k < n; k++) + { + Vector3d dir = Slerp(dirs[k], poleDir, t); + cap.Vertices.Add(center + radius * dir); + } + } + int poleIdx = cap.Vertices.Add(pole); + + int Idx(int r, int k) => r * n + k; + + // quad strips between consecutive full rings + for (int r = 0; r < R - 1; r++) + { + for (int k = 0; k < n; k++) + { + int k2 = (k + 1) % n; + Capping.AddTri(cap, Idx(r, k), Idx(r, k2), Idx(r + 1, k2)); + Capping.AddTri(cap, Idx(r, k), Idx(r + 1, k2), Idx(r + 1, k)); + } + } + // innermost ring fans to the pole + for (int k = 0; k < n; k++) + { + int k2 = (k + 1) % n; + Capping.AddTri(cap, Idx(R - 1, k), Idx(R - 1, k2), poleIdx); + } + + if (cap.Faces.Count == 0) return null; + cap.Normals.ComputeNormals(); + cap.Compact(); + return cap; + } + + /// + /// True when the cap built from (in their current order, fanned toward + /// ) would have its faces pointing toward the sphere center instead of + /// away from it. The actual cap triangles share the orientation of the simple cone fan + /// (pts[k] -> pts[k+1] -> pole), so we sum that fan's face normals dotted with the outward + /// radial direction; a negative total means the winding is inward and the loop should reverse. + /// + private static bool CapWindsInward(Point3d[] pts, Point3d pole, Point3d center) + { + int n = pts.Length; + double radialDot = 0.0; + for (int k = 0; k < n; k++) + { + int k2 = (k + 1) % n; + Vector3d nrm = Vector3d.CrossProduct(pts[k2] - pts[k], pole - pts[k]); + Point3d mid = 0.5 * (pts[k] + pts[k2]); // edge midpoint + Vector3d outward = mid - center; // radially outward from the sphere center + radialDot += nrm * outward; + } + return radialDot < 0.0; + } + + private static int SmallerCapSign(Point3d[] pts, Point3d center, double radius, Vector3d mean) + { + double Area(int sign) + { + Point3d apex = center + sign * radius * mean; + double area = 0.0; + int n = pts.Length; + for (int k = 0; k < n; k++) + { + int k2 = (k + 1) % n; + area += 0.5 * Vector3d.CrossProduct(pts[k] - apex, pts[k2] - apex).Length; + } + return area; + } + return Area(1) <= Area(-1) ? 1 : -1; + } + + /// Spherical interpolation of two unit vectors; linear fallback when (anti)parallel. + private static Vector3d Slerp(Vector3d v0, Vector3d v1, double t) + { + double dot = Math.Max(-1.0, Math.Min(1.0, v0 * v1)); + double omega = Math.Acos(dot); + double so = Math.Sin(omega); + if (omega < 1e-9 || so < 1e-9) + { + Vector3d lin = (1.0 - t) * v0 + t * v1; + if (lin.IsTiny()) return v0; + lin.Unitize(); + return lin; + } + return (Math.Sin((1.0 - t) * omega) / so) * v0 + (Math.Sin(t * omega) / so) * v1; + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("lego.png"); + + public override Guid ComponentGuid + { + get { return new Guid("A7D2E4F1-6C90-4B33-9E27-5C8F1A0B7E64"); } + } + } +} diff --git a/grasshopper/bertini_real/SurfaceSpreadByConnectors.cs b/grasshopper/bertini_real/SurfaceSpreadByConnectors.cs new file mode 100644 index 0000000..1701444 --- /dev/null +++ b/grasshopper/bertini_real/SurfaceSpreadByConnectors.cs @@ -0,0 +1,192 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using Grasshopper; +using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; +using Grasshopper.Kernel.Types; +using Rhino.Geometry; + +namespace bertini_real +{ + /// + /// Spreads pieces along their connector axes to show how they would assemble -- a directional + /// assembly explosion, as opposed to the radial Spread Pieces. + /// + /// The pieces and their singularities form a graph (each singularity is an edge joining exactly + /// two pieces, with a connector direction). This does a breadth-first traversal from a fixed + /// root: each piece's displacement is its parent's displacement PLUS a step of Distance along + /// the connecting axis, oriented so the child slides off the rod away from the parent. So a + /// chain telescopes outward and the root stays put -- robust for trees and chains, and it + /// spanning-trees any cycles. The root defaults to the most-connected piece (a natural hub); + /// set Root to fix a specific piece. Disconnected groups are each rooted independently. + /// + /// Wire Sing Locations / Sing Directions / Sing On Pieces straight from the reader. + /// + public class SurfaceSpreadByConnectors : GH_Component + { + public SurfaceSpreadByConnectors() + : base("Spread By Connectors", "SpreadConn", + "Assembly explosion: BFS the connector graph and slide each piece off its rods along the connector axes", + "bertini_real", "Surface") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddGeometryParameter("Geometry", "G", "Per-piece geometry, one branch per piece (meshes, or mesh + connectors)", GH_ParamAccess.tree); + pManager.AddPointParameter("Sing Locations", "SL", "Singularity locations (from Surface Read GH JSON)", GH_ParamAccess.list); + pManager.AddVectorParameter("Sing Directions", "SD", "Singularity connector directions (from Surface Read GH JSON)", GH_ParamAccess.list); + pManager.AddIntegerParameter("Sing On Pieces", "SOP", "Per piece: indices of the singularities on it (from Surface Read GH JSON)", GH_ParamAccess.tree); + pManager.AddNumberParameter("Distance", "D", "Separation distance per connector along the assembly chain (model units). 0 = no move.", GH_ParamAccess.item, 1.0); + pManager.AddIntegerParameter("Root", "R", "Piece index to hold fixed (-1 = auto: the most-connected piece)", GH_ParamAccess.item, -1); + Params.Input[4].Optional = true; + Params.Input[5].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddGeometryParameter("Geometry", "G", "Spread-apart geometry (same tree structure)", GH_ParamAccess.tree); + pManager.AddVectorParameter("Translations", "T", "Translation applied to each piece", GH_ParamAccess.tree); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + GH_Structure geometry; + if (!DA.GetDataTree(0, out geometry)) return; + + var locations = new List(); + var directions = new List(); + GH_Structure onPieces; + DA.GetDataList(1, locations); + DA.GetDataList(2, directions); + if (!DA.GetDataTree(3, out onPieces)) return; + + double distance = 1.0; + DA.GetData(4, ref distance); + int root = -1; + DA.GetData(5, ref root); + + // pieces we have geometry for, with their centers and branch paths + var pieces = new List(); + var centroid = new Dictionary(); + var branchPath = new Dictionary(); + for (int b = 0; b < geometry.PathCount; b++) + { + GH_Path path = geometry.get_Path(b); + int pi = path.Indices.Length > 0 ? path.Indices[path.Indices.Length - 1] : b; + BoundingBox bb = BoundingBox.Empty; + bool any = false; + foreach (var goo in geometry.get_Branch(path)) + if (goo is IGH_GeometricGoo gg && gg.IsValid) { bb.Union(gg.Boundingbox); any = true; } + if (!any || centroid.ContainsKey(pi)) continue; + pieces.Add(pi); + centroid[pi] = bb.Center; + branchPath[pi] = path; + } + if (pieces.Count == 0) return; + + // invert Sing On Pieces -> which pieces each singularity touches (restricted to our pieces) + var singToPieces = new Dictionary>(); + foreach (int pi in pieces) + { + var sopPath = new GH_Path(pi); + if (!onPieces.PathExists(sopPath)) continue; + foreach (var goo in onPieces.get_Branch(sopPath)) + if (goo is GH_Integer gi) + { + if (!singToPieces.TryGetValue(gi.Value, out var lst)) { lst = new List(); singToPieces[gi.Value] = lst; } + if (!lst.Contains(pi)) lst.Add(pi); + } + } + + // adjacency: a singularity joining exactly two pieces is an edge + var adj = new Dictionary>(); + void Link(int a, int bb2, int s) + { + if (!adj.TryGetValue(a, out var la)) { la = new List<(int, int)>(); adj[a] = la; } + la.Add((bb2, s)); + } + foreach (var kv in singToPieces) + if (kv.Value.Count == 2) + { + Link(kv.Value[0], kv.Value[1], kv.Key); + Link(kv.Value[1], kv.Value[0], kv.Key); + } + + // BFS from a root, accumulating displacement along the connector chain + var disp = new Dictionary(); + foreach (int pi in pieces) disp[pi] = Vector3d.Zero; + var visited = new HashSet(); + + // root order: an explicit Root first, then most-connected pieces (one root per component) + var order = pieces.OrderByDescending(pi => adj.TryGetValue(pi, out var l) ? l.Count : 0).ToList(); + if (root >= 0 && centroid.ContainsKey(root)) { order.Remove(root); order.Insert(0, root); } + + int components = 0; + foreach (int start in order) + { + if (visited.Contains(start)) continue; + components++; + var queue = new Queue(); + queue.Enqueue(start); + visited.Add(start); + disp[start] = Vector3d.Zero; + + while (queue.Count > 0) + { + int cur = queue.Dequeue(); + if (!adj.TryGetValue(cur, out var nbrs)) continue; + foreach (var (nbr, sing) in nbrs) + { + if (visited.Contains(nbr)) continue; + if (sing < 0 || sing >= directions.Count || sing >= locations.Count) continue; + + Vector3d axis = directions[sing]; + if (axis.IsTiny()) continue; + axis.Unitize(); + + // orient the axis so the child slides away from the singularity (and parent) + double along = (centroid[nbr] - locations[sing]) * axis; + Vector3d step = distance * (along >= 0 ? 1.0 : -1.0) * axis; + + disp[nbr] = disp[cur] + step; + visited.Add(nbr); + queue.Enqueue(nbr); + } + } + } + + if (components > 1) + AddRuntimeMessage(GH_RuntimeMessageLevel.Remark, + $"{components} disconnected connector groups; each is rooted independently and may overlap at the origin."); + + // emit + var outGeo = new DataTree(); + var outVec = new DataTree(); + foreach (int pi in pieces) + { + GH_Path path = branchPath[pi]; + Transform xf = Transform.Translation(disp[pi]); + foreach (var goo in geometry.get_Branch(path)) + { + if (!(goo is IGH_GeometricGoo gg) || !gg.IsValid) continue; + IGH_GeometricGoo moved = gg.DuplicateGeometry(); + moved = moved.Transform(xf); + outGeo.Add(moved, path); + } + outVec.Add(disp[pi], path); + } + + DA.SetDataTree(0, outGeo); + DA.SetDataTree(1, outVec); + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("transform.png"); + + public override Guid ComponentGuid + { + get { return new Guid("9F26C7B4-3A81-4D60-8E15-7C0B2F95E6A3"); } + } + } +} diff --git a/grasshopper/bertini_real/SurfaceSpreadPieces.cs b/grasshopper/bertini_real/SurfaceSpreadPieces.cs new file mode 100644 index 0000000..6a9c16b --- /dev/null +++ b/grasshopper/bertini_real/SurfaceSpreadPieces.cs @@ -0,0 +1,133 @@ +using System; +using System.Collections.Generic; +using Grasshopper; +using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; +using Grasshopper.Kernel.Types; +using Rhino.Geometry; + +namespace bertini_real +{ + /// + /// Spreads the pieces of a decomposition apart so they can be seen individually without + /// touching. (Named "Spread" rather than "Explode" to avoid Grasshopper's sense of explode = + /// decompose into constituents.) Each piece is one tree branch and is translated radially + /// outward from the common center by Distance model units (an absolute distance, matching + /// Spread By Connectors): Distance = 0 leaves everything in place, larger spreads them further. + /// + /// Operates on any per-piece GEOMETRY tree, so it accepts bare meshes (Surface Read GH JSON / + /// Close Piece) or a piece's mesh together with its connectors (Surface Group By Piece) -- all + /// items in a branch move together, so a piece and its connectors stay assembled. + /// + public class SurfaceSpreadPieces : GH_Component + { + public SurfaceSpreadPieces() + : base("Spread Pieces", "Spread", + "Move decomposition pieces apart (radially from their common center) so they can be seen separated", + "bertini_real", "Surface") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddGeometryParameter("Geometry", "G", "Per-piece geometry, one branch per piece: meshes, or a piece's mesh + connectors from Surface Group By Piece", GH_ParamAccess.tree); + pManager.AddNumberParameter("Distance", "D", "Distance each piece moves outward from the center (model units). 0 = no move.", GH_ParamAccess.item, 1.0); + pManager.AddPointParameter("Center", "C", "Center to spread away from (default: average of the piece centers)", GH_ParamAccess.item); + Params.Input[2].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddGeometryParameter("Geometry", "G", "Spread-apart geometry (same tree structure)", GH_ParamAccess.tree); + pManager.AddVectorParameter("Translations", "T", "Translation applied to each piece", GH_ParamAccess.tree); + pManager.AddPointParameter("Center", "C", "The center the pieces were spread from", GH_ParamAccess.item); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + GH_Structure geometry; + if (!DA.GetDataTree(0, out geometry)) return; + + double distance = 1.0; + DA.GetData(1, ref distance); + + Point3d center = Point3d.Unset; + bool hasCenter = DA.GetData(2, ref center); + + // per-piece center = center of the branch's combined bounding box (all geometry in it) + var paths = new List(); + var pieceCenter = new List(); + for (int b = 0; b < geometry.PathCount; b++) + { + GH_Path path = geometry.get_Path(b); + BoundingBox bb = BoundingBox.Empty; + bool any = false; + foreach (var goo in geometry.get_Branch(path)) + { + if (!(goo is IGH_GeometricGoo gg) || !gg.IsValid) continue; + bb.Union(gg.Boundingbox); + any = true; + } + if (!any) continue; + paths.Add(path); + pieceCenter.Add(bb.Center); + } + + if (paths.Count == 0) return; + + if (paths.Count == 1) + AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, + "Only one tree branch received, so there is nothing to spread (one piece per branch). " + + "The per-piece tree structure was probably flattened upstream (Scale, a wire, or a Flatten on this input) -- " + + "keep one branch per piece coming in."); + + // overall center: explicit input, else the average of the piece centers + Point3d origin; + if (hasCenter) + { + origin = center; + } + else + { + double x = 0, y = 0, z = 0; + foreach (var c in pieceCenter) { x += c.X; y += c.Y; z += c.Z; } + origin = new Point3d(x / pieceCenter.Count, y / pieceCenter.Count, z / pieceCenter.Count); + } + + var outGeo = new DataTree(); + var outVec = new DataTree(); + + for (int i = 0; i < paths.Count; i++) + { + GH_Path path = paths[i]; + // absolute distance: move each piece Distance units along its outward (unit) direction + Vector3d dir = pieceCenter[i] - origin; + if (!dir.IsTiny()) dir.Unitize(); + Vector3d t = distance * dir; + Transform xf = Transform.Translation(t); + + foreach (var goo in geometry.get_Branch(path)) + { + if (!(goo is IGH_GeometricGoo gg) || !gg.IsValid) continue; + // duplicate so the input geometry is left untouched, then translate + IGH_GeometricGoo moved = gg.DuplicateGeometry(); + moved = moved.Transform(xf); + outGeo.Add(moved, path); + } + + outVec.Add(t, path); + } + + DA.SetDataTree(0, outGeo); + DA.SetDataTree(1, outVec); + DA.SetData(2, origin); + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("transform.png"); + + public override Guid ComponentGuid + { + get { return new Guid("C2A6F39B-4E78-4D15-8B0A-3F9E1C7D6052"); } + } + } +} diff --git a/grasshopper/bertini_real/SurfaceUntangleCurves.cs b/grasshopper/bertini_real/SurfaceUntangleCurves.cs new file mode 100644 index 0000000..05c5761 --- /dev/null +++ b/grasshopper/bertini_real/SurfaceUntangleCurves.cs @@ -0,0 +1,104 @@ +using System; +using System.Collections; +using Grasshopper; +using Grasshopper.Kernel; +using Grasshopper.Kernel.Data; +using Grasshopper.Kernel.Types; + +namespace bertini_real +{ + /// + /// Splits the parallel Curves + Curve Types trees coming out of "Surface Read GH JSON" into + /// one output per curve type, so you can grab just the type you want without manual tree + /// filtering / get-item juggling. Per-piece tree paths are preserved on every output, so a + /// curve stays associated with the piece it came from. + /// + public class SurfaceUntangleCurves : GH_Component + { + public SurfaceUntangleCurves() + : base("Untangle Curves By Type", "UntangleCurves", + "Split the embedded curves from Surface Read GH JSON into one output per type", + "bertini_real", "Surface") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddCurveParameter("Curves", "C", "Embedded curves tree from Surface Read GH JSON", GH_ParamAccess.tree); + pManager.AddTextParameter("Curve Types", "T", "Curve Types tree from Surface Read GH JSON (parallel to Curves)", GH_ParamAccess.tree); + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddCurveParameter("Critical", "Cr", "Critical curve pieces", GH_ParamAccess.tree); + pManager.AddCurveParameter("Sphere", "Sp", "Sphere curve pieces", GH_ParamAccess.tree); + pManager.AddCurveParameter("Singular", "Si", "Singular curve pieces", GH_ParamAccess.tree); + pManager.AddCurveParameter("Midslice", "Mi", "Midslice curve pieces", GH_ParamAccess.tree); + pManager.AddCurveParameter("Critslice", "Cs", "Critslice curve pieces", GH_ParamAccess.tree); + pManager.AddCurveParameter("Other", "Ot", "Curves whose type is unknown / unrecognized", GH_ParamAccess.tree); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + GH_Structure curves; + GH_Structure types; + if (!DA.GetDataTree(0, out curves)) return; + if (!DA.GetDataTree(1, out types)) return; + + var critical = new GH_Structure(); + var sphere = new GH_Structure(); + var singular = new GH_Structure(); + var midslice = new GH_Structure(); + var critslice = new GH_Structure(); + var other = new GH_Structure(); + + for (int b = 0; b < curves.PathCount; b++) + { + GH_Path path = curves.get_Path(b); + IList curveBranch = curves.get_Branch(path); + IList typeBranch = types.PathExists(path) ? types.get_Branch(path) : null; + + if (typeBranch == null || typeBranch.Count != curveBranch.Count) + AddRuntimeMessage(GH_RuntimeMessageLevel.Warning, + $"Types do not line up with Curves at path {path}; missing tags treated as 'Other'."); + + for (int i = 0; i < curveBranch.Count; i++) + { + var gc = curveBranch[i] as GH_Curve; + if (gc == null) continue; + + string t = "unknown"; + if (typeBranch != null && i < typeBranch.Count && typeBranch[i] is GH_String gs && gs.Value != null) + t = gs.Value.Trim().ToLowerInvariant(); + + GH_Structure target; + switch (t) + { + case "critical": target = critical; break; + case "sphere": target = sphere; break; + case "singular": target = singular; break; + case "midslice": target = midslice; break; + case "critslice": target = critslice; break; + default: target = other; break; + } + + target.Append(gc, path); + } + } + + DA.SetDataTree(0, critical); + DA.SetDataTree(1, sphere); + DA.SetDataTree(2, singular); + DA.SetDataTree(3, midslice); + DA.SetDataTree(4, critslice); + DA.SetDataTree(5, other); + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("transform.png"); + + public override Guid ComponentGuid + { + get { return new Guid("3F1B8D6C-9A42-4C71-B8E5-71D0A6F4C982"); } + } + } +} diff --git a/grasshopper/bertini_real/UtilityCenteredCylinder.cs b/grasshopper/bertini_real/UtilityCenteredCylinder.cs new file mode 100644 index 0000000..a0eed32 --- /dev/null +++ b/grasshopper/bertini_real/UtilityCenteredCylinder.cs @@ -0,0 +1,102 @@ +using System; +using Grasshopper.Kernel; +using Rhino.Geometry; + +namespace bertini_real +{ + /// + /// A closed (capped) cylinder/prism centered on a plane's origin and running along its Z axis + /// from -Length/2 to +Length/2. Sides picks the cross-section: 1 = a true round cylinder, + /// 3 = triangular prism, 4 = square, 5 = pentagon, ... (Sides = 2 is degenerate.) Radius is + /// the circumradius (distance to the polygon vertices). Handy as a connector / boolean cutter + /// -- centered so it straddles its placement point, closed so mesh booleans bite. + /// + public class UtilityCenteredCylinder : GH_Component + { + public UtilityCenteredCylinder() + : base("Centered Closed Cylinder", "CenCyl", + "A capped cylinder/prism centered on the plane origin along its Z axis. Sides: 1 = round, 3+ = N-gon prism.", + "bertini_real", "Utility") + { + } + + protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager) + { + pManager.AddNumberParameter("Radius", "R", "Radius (circumradius for polygons -- distance to the vertices)", GH_ParamAccess.item, 5.0); + pManager.AddNumberParameter("Length", "L", "Length (centered on the plane origin)", GH_ParamAccess.item, 100.0); + pManager.AddIntegerParameter("Sides", "N", "Cross-section: 1 = round cylinder, 3 = triangle, 4 = square, 5 = pentagon, ...", GH_ParamAccess.item, 1); + pManager.AddPlaneParameter("Plane", "P", "Center plane; the cylinder runs along its Z axis", GH_ParamAccess.item, Plane.WorldXY); + Params.Input[0].Optional = true; + Params.Input[1].Optional = true; + Params.Input[2].Optional = true; + Params.Input[3].Optional = true; + } + + protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager) + { + pManager.AddBrepParameter("Cylinder", "C", "Closed cylinder/prism Brep", GH_ParamAccess.item); + } + + protected override void SolveInstance(IGH_DataAccess DA) + { + double radius = 5.0; + double length = 100.0; + int sides = 1; + Plane plane = Plane.WorldXY; + DA.GetData(0, ref radius); + DA.GetData(1, ref length); + DA.GetData(2, ref sides); + DA.GetData(3, ref plane); + + if (radius <= 0 || length <= 0) + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "Radius and Length must be positive."); + return; + } + if (sides < 1 || sides == 2) + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "Sides must be 1 (round) or >= 3 (polygon); 2 is degenerate."); + return; + } + + // base plane is half the length below the center, so the result is centered + var bottom = new Plane(plane.Origin - plane.ZAxis * (length / 2.0), plane.XAxis, plane.YAxis); + + Brep brep; + if (sides == 1) + { + var cylinder = new Cylinder(new Circle(bottom, radius), length); + brep = cylinder.ToBrep(true, true); // cap both ends -> closed solid + } + else + { + // regular n-gon profile in the base plane, extruded by Length and capped + var pts = new Point3d[sides + 1]; + for (int i = 0; i < sides; i++) + { + double a = 2.0 * Math.PI * i / sides; + pts[i] = bottom.PointAt(radius * Math.Cos(a), radius * Math.Sin(a)); + } + pts[sides] = pts[0]; + + var profile = new Polyline(pts).ToPolylineCurve(); + Extrusion extrusion = Extrusion.Create(profile, length, true); + if (extrusion == null) + { + AddRuntimeMessage(GH_RuntimeMessageLevel.Error, "Failed to build the prism."); + return; + } + brep = extrusion.ToBrep(); + } + + DA.SetData(0, brep); + } + + protected override System.Drawing.Bitmap Icon => IconLoader.GetIcon("lego.png"); + + public override Guid ComponentGuid + { + get { return new Guid("8A3F9C21-6D74-4E08-B5A2-1C7E0F934D6B"); } + } + } +} diff --git a/grasshopper/bertini_real/bertini_real.csproj b/grasshopper/bertini_real/bertini_real.csproj index 437d641..161bf10 100644 --- a/grasshopper/bertini_real/bertini_real.csproj +++ b/grasshopper/bertini_real/bertini_real.csproj @@ -28,6 +28,15 @@ + + + /Applications/Rhino 8.app/Contents/Frameworks/RhCore.framework/Versions/A/Resources/ref/net48/Grasshopper.dll + + + /Applications/Rhino 8.app/Contents/Frameworks/RhCore.framework/Versions/A/Resources/ref/net48/RhinoCommon.dll + + + diff --git a/grasshopper/docs/usage.puml b/grasshopper/docs/usage.puml new file mode 100644 index 0000000..90e9549 --- /dev/null +++ b/grasshopper/docs/usage.puml @@ -0,0 +1,128 @@ +@startuml +' Bertini_real Grasshopper components — single-source usage graph. +' Render with PlantUML: plantuml grasshopper/docs/usage.puml (-> usage.png/svg) + +title Bertini_real Grasshopper — single-source pipeline + +skinparam backgroundColor #FDFDFD +skinparam shadowing false +skinparam defaultFontName Helvetica +skinparam ArrowColor #5A6B7B +skinparam rectangle { + BorderColor #3A5A78 + BackgroundColor #EAF2FA + RoundCorner 12 +} +skinparam artifact { + BorderColor #7A6A3A + BackgroundColor #FBF4DD +} + +' ---- the one source ---- +artifact "br_gh_export.json\n(written by Python\nexport_gh_json / scripts)" as JSON + +note top of JSON + **No Rhino?** The same pipeline + (pieces -> caps -> close -> boolean) + is available in pure Python -- see + python/docs/tutorials/capping_and_joining.rst + (and scripts/close_pieces.py). +end note + +' ---- the one reader ---- +rectangle "**Surface Read GH JSON**\l--\lout:\l• Vertices (Point list)\l• Meshes (tree {piece})\l• Mesh Faces (tree {piece})\l• Curves (tree {piece})\l• Curve Types (tree {piece})\l• Curve Indices (tree {piece,curve})\l• Face Indices (tree {piece})\l• Sphere (Brep)\l• Sing Locations (list)\l• Sing Directions (list)\l• Sing Parities (tree {sing})\l• Sing On Pieces (tree {piece})\l" as READER + +JSON --> READER : File Path + +' ====================================================================== +' Branch 1 — curves +' ====================================================================== +rectangle "**Untangle Curves By Type**\l--\lout (each tree {piece}):\l• Critical • Sphere • Singular\l• Midslice • Critslice • Other\l" as UNTANGLE + +READER --> UNTANGLE : Curves\n+ Curve Types + +' ====================================================================== +' Branch 2 — geometry: cap -> close -> boolean -> (spread) +' ====================================================================== +rectangle "**Sphere Caps** / **Flat Caps**\l--\lin: Meshes, Sphere,\l Tolerance, Resolution\lout: Caps (tree {piece})\l(Flat Caps fans to the loop centroid;\l drop-in alternative to Sphere Caps)\l" as CAPS + +rectangle "**Close Piece**\l--\lin: Meshes, Caps, Tolerance\lout:\l• Closed (tree {piece})\l• Is Closed (tree {piece})\l• Naked Edges (tree {piece})\l" as CLOSE + +rectangle "**Boolean Piece**\l--\lin: Solid, Features, Operations(±1)\lout:\l• Result (tree {piece})\l• Is Closed (tree {piece})\l• Report (tree {piece})\l" as BOOL + +rectangle "**Color By Function**\l--\lin: Meshes, Function(x,y,z),\l Colours, Domain\lout: Meshes (vertex-colored),\l Values, Domain\l" as COLOR + +rectangle "**Spread Pieces** / **Spread By Connectors**\l--\lin: Geometry, Distance\l (By Connectors also takes Sing\l Locations/Directions/On Pieces, Root)\lout: Geometry, Translations\l(radial vs. assembly explosion\l along the connector axes)\l" as SPREAD + +READER --> CAPS : Meshes +READER --> CAPS : Sphere +READER --> CLOSE : Meshes +CAPS --> CLOSE : Caps +CLOSE --> BOOL : Closed (as Solid) + +' ====================================================================== +' Branch 3 — singularity connectors -> features -> boolean +' ====================================================================== +rectangle "**Surface Place Components**\l--\lin: Sing Locations, Sing Directions,\l Sing Parities, Sing On Pieces,\l Plug+/-, Socket+/- (any subset)\lout (each tree {piece}):\l• Plugs positive / negative\l• Sockets positive / negative\l" as PLACE + +rectangle "**Connectors To Features**\l--\lin: Plugs+/-, Sockets+/-\lout:\l• Features (tree {piece}, ordered)\l• Operations (tree {piece}, signs ±1)\l" as CONNFEAT + +rectangle "**Utility: Centered Closed Cylinder**\l--\lin: Radius, Length, Sides, Plane\lout: closed Brep (round, or N-gon prism)\l" as UTIL + +READER --> PLACE : Sing Locations\nSing Directions\nSing Parities\nSing On Pieces +UTIL --> PLACE : Plug+ / Socket+\n(connector geometry, true size) +PLACE --> CONNFEAT : Plugs+/-\nSockets+/- +CONNFEAT --> BOOL : Features\n+ Operations +UTIL ..> BOOL : Features (direct) + +' ====================================================================== +' Grouping / viewing +' ====================================================================== +rectangle "**Surface Group By Piece**\l--\lin: Meshes, Plugs+/-, Sockets+/-\lout:\l• Pieces with connectors\l (tree {piece}: mesh + connectors)\l• Piece Indices\l" as GROUP + +CLOSE ..> GROUP : Meshes\n(or Boolean Result) +PLACE --> GROUP : Plugs+/-\nSockets+/- + +BOOL --> COLOR : Result +COLOR --> SPREAD : Meshes (colored) +GROUP ..> SPREAD : Pieces with connectors +CLOSE ..> COLOR : Closed (alt.) +CLOSE ..> SPREAD : Closed (alt.) + +' ---- notes ---- +note right of READER + **One file, one reader.** + The whole graph is fed by + br_gh_export.json through this + single component (no second JSON). +end note + +note bottom of BOOL + **Boolean Piece is an ordered fold.** + Start from Solid, then apply each Feature + in order: +1 union, -1 subtract (default + subtract). Order matters -- each step acts + on the running result. (Python twin: + surface.mesh_boolean_fold, manifold3d.) +end note + +note bottom of COLOR + Colors meshes by a scalar f(x,y,z) at + each vertex. Mesh vertex colors survive + Spread, so color before or after spreading. +end note + +note as LEGEND + **Data-shape conventions** + • per-piece -> tree {piece} + • per-singularity -> flat list (index = sing id) + • per-curve idx -> tree {piece, curve} + • Vertices is one shared Point list; + meshes/curves index into it. + • per-piece geometry -> "Geometry"/"G" + (any type); meshes -> "Meshes"/"M". + solid arrow = typical wire, dashed = alternative +end note +LEGEND .. SPREAD + +@enduml diff --git a/python/bertini_real/curve/__init__.py b/python/bertini_real/curve/__init__.py index 8bbf100..8ee7be7 100644 --- a/python/bertini_real/curve/__init__.py +++ b/python/bertini_real/curve/__init__.py @@ -7,6 +7,22 @@ import numpy as np import copy +import json + + +def _points_to_xyz(points): + """ + convert an iterable of points (each of length 1, 2, or 3+) into a list of [x, y, z] + triples, padding missing coordinates with 0. used by the Grasshopper JSON export so + the unified vertex set is always 3d for Rhino. + """ + out = [] + for p in points: + x = float(p[0]) if len(p) > 0 else 0.0 + y = float(p[1]) if len(p) > 1 else 0.0 + z = float(p[2]) if len(p) > 2 else 0.0 + out.append([x, y, z]) + return out @@ -70,13 +86,55 @@ def __str__(self): + def to_vertices(self): + """ + a list of all the Vertices on this piece of a curve. See also to_points, which gives them as numpy points + + degenerate edges are skipped, and adjacent endpoints are unified into one point. loses all edge or connectivity structure. + """ + + + if not self.directed_edges: + raise NotImplementedError('insert code memoizing / computing the directed edges') + + # unpack a few things + vertices = self.curve.vertices # these have already been dehomogenized + c = self.curve + + the_vertices = [] # built up over time. difficult to pre-allocate. yagni. + prev_point_index = -1 + + for edge_index, direction in self.directed_edges: + + if is_edge_degenerate(c.edges[edge_index]): + continue + + + + if len(c.sampler_data)>0: + point_indices = c.sampler_data[edge_index] + else: + point_indices = c.edges[edge_index] + + if direction==EdgeDirection.backward: + point_indices = point_indices[::-1] + + vertices_this_edge = [] + for ii in point_indices: + if ii != prev_point_index: + vertices_this_edge.append(vertices[ii]) + prev_point_index = ii + + the_vertices.extend(vertices_this_edge) + + return the_vertices def to_points(self): """ computes a numpy array of points, in order, for this piece of a curve. - generate edges are skipped, and adjacent endpoints are unified into one point. + degenerate edges are skipped, and adjacent endpoints are unified into one point. """ @@ -97,7 +155,8 @@ def to_points(self): - if len(c.sampler_data)>0: + # an unsampled curve has sampler_data == None; fall back to the raw edge indices + if c.sampler_data and len(c.sampler_data)>0: point_indices = c.sampler_data[edge_index] else: point_indices = c.edges[edge_index] @@ -119,6 +178,46 @@ def to_points(self): return points + def to_point_indices(self): + """ + computes the ordered list of global vertex indices for this piece of a curve. + + this is the index-space twin of `to_points`: degenerate edges are skipped, and + adjacent duplicate endpoints are unified. the returned indices point into + `self.curve.vertices` (the unified, shared vertex set), so an embedded curve and the + surface mesh refer to the same points. used by the Grasshopper JSON export. + """ + + if not self.directed_edges: + raise NotImplementedError('insert code memoizing / computing the directed edges') + + c = self.curve + + indices = [] + prev_point_index = -1 + + for edge_index, direction in self.directed_edges: + + if is_edge_degenerate(c.edges[edge_index]): + continue + + # an unsampled curve has sampler_data == None; fall back to the raw edge indices + if c.sampler_data and len(c.sampler_data) > 0: + point_indices = c.sampler_data[edge_index] + else: + point_indices = c.edges[edge_index] + + if direction == EdgeDirection.backward: + point_indices = point_indices[::-1] + + for ii in point_indices: + if ii != prev_point_index: + indices.append(int(ii)) + prev_point_index = ii + + return indices + + class Curve(Decomposition): """ a Curve @@ -214,8 +313,80 @@ def break_into_pieces(self, edge_indices = None): - + def to_vertices(self): + """ + a list of all the Vertices on this piece of a curve. See also to_points, which gives them as numpy points + + degenerate edges are skipped, and adjacent endpoints are unified into one point. loses all edge or connectivity structure. + """ + + # unpack a few things + vertices = self.vertices # these have already been dehomogenized + c = self + + the_vertices = [] # built up over time. difficult to pre-allocate. yagni. + prev_point_index = -1 + + for edge_index in range(self.num_edges): + e = self.edges[edge_index] + + if is_edge_degenerate(e): + continue + + if len(c.sampler_data)>0: + point_indices = c.sampler_data[edge_index] + else: + point_indices = e + + vertices_this_edge = [] + for ii in point_indices: + if ii != prev_point_index: + vertices_this_edge.append(vertices[ii]) + prev_point_index = ii + + the_vertices.extend(vertices_this_edge) + + return the_vertices + + + + def export_gh_json(self, filename="br_gh_export.json"): + """ + write a self-contained JSON describing this curve for the Grasshopper plugin. + + the file holds one unified vertex set (`vertices`); each curve piece is an ordered + list of indices into that set. see `Surface.export_gh_json` for the surface analogue. + """ + points = self.extract_points() + + pieces = self.break_into_pieces(set(range(self.num_edges))) + + curve_pieces = [] + for ii, p in enumerate(pieces): + curve_pieces.append({ + "piece_index": ii, + "type": "standalone", + "curve_name": self.inputfilename, + "vertex_indices": p.to_point_indices(), + }) + + contents = { + "format_version": 2, + "decomposition_type": "curve", + "source_directory": self.directory, + "num_variables": self.num_variables, + "vertices": _points_to_xyz(points), + "vertex_count": len(points), + "sphere": self._sphere_dict(), + "curve_pieces": curve_pieces, + } + + with open(filename, "w") as f: + json.dump(contents, f, indent=2) + + print("wrote " + filename) + return filename def parse_edge(self, directory): diff --git a/python/bertini_real/data/__init__.py b/python/bertini_real/data/__init__.py index e843552..299f8ad 100644 --- a/python/bertini_real/data/__init__.py +++ b/python/bertini_real/data/__init__.py @@ -193,3 +193,13 @@ def gather_and_save(): import dill dill.dump(b, fileObject) fileObject.close() + + +def gather_and_export_gh(filename="br_gh_export.json"): + """ + gather a decomposition from the current directory and write a self-contained JSON + for the Grasshopper plugin. works for both curves and surfaces (both expose + `export_gh_json`). returns the filename written. + """ + d = gather() + return d.export_gh_json(filename) diff --git a/python/bertini_real/decomposition/__init__.py b/python/bertini_real/decomposition/__init__.py index 742a880..574f741 100644 --- a/python/bertini_real/decomposition/__init__.py +++ b/python/bertini_real/decomposition/__init__.py @@ -85,6 +85,18 @@ def read_input(self, directory): + def _sphere_dict(self): + """ + the bounding sphere of the decomposition (center + radius from the decomp file), + as a plain dict for the Grasshopper export. center is padded/truncated to [x, y, z]. + """ + c = list(self.center) + x = float(c[0]) if len(c) > 0 else 0.0 + y = float(c[1]) if len(c) > 1 else 0.0 + z = float(c[2]) if len(c) > 2 else 0.0 + return {"center": [x, y, z], "radius": float(self.radius)} + + def extract_points(self, indices=None): """ Helper method Extract points from vertices as a list diff --git a/python/bertini_real/surface/__init__.py b/python/bertini_real/surface/__init__.py index e5417a1..f619aad 100644 --- a/python/bertini_real/surface/__init__.py +++ b/python/bertini_real/surface/__init__.py @@ -16,7 +16,7 @@ import bertini_real.exception as br_except import numpy as np from bertini_real.decomposition import Decomposition -from bertini_real.curve import Curve, CurvePiece, is_edge_degenerate +from bertini_real.curve import Curve, CurvePiece, is_edge_degenerate, _points_to_xyz from bertini_real.vertex import Vertex from bertini_real.vertex import VertexType from bertini_real.util import ReversableList @@ -56,6 +56,24 @@ +def _mesh_triangles(mesh): + """ + flatten a `trimesh.Trimesh`'s faces into a dict for the Grasshopper JSON export. + + the triangle entries are indices into the surface's unified (global) vertex set, + because `as_mesh_raw`/`as_mesh_smooth` build the mesh from `extract_points()` with + `keep_all_vertices=True` (so trimesh does not reindex). returns None if mesh is None. + """ + if mesh is None: + return None + + triangles = np.asarray(mesh.faces, dtype=int).reshape(-1).tolist() + return { + "triangles": triangles, + "triangle_count": len(mesh.faces), + } + + def export_mesh(mesh, basename, autoname_using_folder=False, file_type=_default_file_type, verbose=True): """ Saves a mesh (generated elsewhere) to disk, @@ -153,9 +171,283 @@ def copy_all_scad_files_here(): +def _slerp(v0, v1, t): + """Spherical interpolation of two unit vectors; linear fallback when (anti)parallel.""" + dot = np.clip(np.dot(v0, v1), -1.0, 1.0) + omega = np.arccos(dot) + so = np.sin(omega) + if omega < 1e-9 or so < 1e-9: + lin = (1.0 - t) * v0 + t * v1 + n = np.linalg.norm(lin) + return v0 if n < 1e-12 else lin / n + return np.sin((1.0 - t) * omega) / so * v0 + np.sin(t * omega) / so * v1 + + +def _on_sphere_boundary_loops(mesh, center, radius, tol): + """ + Ordered loops (lists of vertex indices) of `mesh`'s naked boundary edges that lie on the + sphere of the given center/radius. Only clean degree-2 cycles are returned; open arcs or + non-manifold junctions are dropped. + """ + from collections import Counter, defaultdict + + V = np.asarray(mesh.vertices) + edge_count = Counter() + for tri in mesh.faces: + a, b, c = int(tri[0]), int(tri[1]), int(tri[2]) + for u, v in ((a, b), (b, c), (c, a)): + edge_count[frozenset((u, v))] += 1 + naked = [tuple(e) for e, cnt in edge_count.items() if cnt == 1 and len(e) == 2] + + def on_sphere(i): + return abs(np.linalg.norm(V[i][:3] - center) - radius) < tol + + adj = defaultdict(list) + for e in naked: + a, b = tuple(e) + if a != b and on_sphere(a) and on_sphere(b): + adj[a].append(b) + adj[b].append(a) + + loops = [] + seen = set() + for start in list(adj): + if start in seen: + continue + loop = [start] + seen.add(start) + prev, cur, ok = -1, start, True + while True: + nbrs = adj[cur] + if len(nbrs) != 2: + ok = False + break + nxt = nbrs[0] if nbrs[0] != prev else nbrs[1] + if nxt == start: + break + if nxt in seen: + ok = False + break + loop.append(nxt) + seen.add(nxt) + prev, cur = cur, nxt + if ok and len(loop) >= 3: + loops.append(loop) + return loops + + +def sphere_cap_meshes(mesh, center, radius, resolution=4, tol=1e-3): + """ + Faceted spherical caps that close `mesh`'s naked boundary loops lying on the sphere. + + The Python twin of the Grasshopper "Sphere Caps" component: it caps the mesh's OWN + boundary (so the caps share its vertices and weld watertight), keeps the smaller-area + side of each loop, and subdivides each cap into `resolution` radial rings slerped along + the sphere. Returns a list of `trimesh.Trimesh`. + """ + center = np.asarray(center, dtype=float)[:3] + V = np.asarray(mesh.vertices) + caps = [] + + for loop in _on_sphere_boundary_loops(mesh, center, radius, tol): + pts = np.array([V[i][:3] for i in loop]) + dirs = np.array([(p - center) / np.linalg.norm(p - center) for p in pts]) + mean = pts.mean(axis=0) - center + nrm = np.linalg.norm(mean) + mean = mean / nrm if nrm > 1e-12 else np.array([0.0, 0.0, 1.0]) + + def fan_area(sign): + apex = center + sign * radius * mean + return sum(np.linalg.norm(np.cross(pts[k] - apex, pts[(k + 1) % len(pts)] - apex)) / 2.0 + for k in range(len(pts))) + + sign = 1 if fan_area(1) <= fan_area(-1) else -1 + pole_dir = sign * mean + pole = center + radius * pole_dir + + R = max(1, int(resolution)) + n = len(loop) + verts = [p for p in pts] + for r in range(1, R): + t = r / R + for k in range(n): + verts.append(center + radius * _slerp(dirs[k], pole_dir, t)) + pole_i = len(verts) + verts.append(pole) + verts = np.array(verts) + + def vid(r, k): + return r * n + k + + faces = [] + + def add(i, j, k): + # skip only truly degenerate (collapsed-edge) triangles; keep thin ones + if (np.linalg.norm(verts[i] - verts[j]) < 1e-9 or + np.linalg.norm(verts[j] - verts[k]) < 1e-9 or + np.linalg.norm(verts[i] - verts[k]) < 1e-9): + return + faces.append([i, j, k]) + + for r in range(R - 1): + for k in range(n): + k2 = (k + 1) % n + add(vid(r, k), vid(r, k2), vid(r + 1, k2)) + add(vid(r, k), vid(r + 1, k2), vid(r + 1, k)) + for k in range(n): + add(vid(R - 1, k), vid(R - 1, (k + 1) % n), pole_i) + + if faces: + caps.append(trimesh.Trimesh(verts, np.array(faces), process=False)) + + return caps + + +def flat_cap_meshes(mesh, center, radius, resolution=1, tol=1e-3): + """ + Flat caps that close `mesh`'s on-sphere boundary loops with a fan to each loop's centroid -- + the flat alternative to sphere_cap_meshes, and the Python twin of the "Flat Caps" component. + + Same on-sphere boundary as sphere_cap_meshes, but the apex is the loop centroid (a flat fill) + rather than a point on the sphere, so e.g. a Bertini cylinder gets flat disk ends. Caps the + mesh's OWN boundary (welds watertight), with `resolution` concentric linearly-interpolated + rings (1 = a single flat fan). Returns a list of trimesh.Trimesh. + """ + center = np.asarray(center, dtype=float)[:3] + V = np.asarray(mesh.vertices) + caps = [] + + for loop in _on_sphere_boundary_loops(mesh, center, radius, tol): + pts = np.array([V[i][:3] for i in loop]) + centroid = pts.mean(axis=0) + + R = max(1, int(resolution)) + n = len(loop) + verts = [p for p in pts] + for r in range(1, R): + t = r / R + for k in range(n): + verts.append(pts[k] + t * (centroid - pts[k])) + apex_i = len(verts) + verts.append(centroid) + verts = np.array(verts) + + def vid(r, k): + return r * n + k + + faces = [] + + def add(i, j, k): + if (np.linalg.norm(verts[i] - verts[j]) < 1e-9 or + np.linalg.norm(verts[j] - verts[k]) < 1e-9 or + np.linalg.norm(verts[i] - verts[k]) < 1e-9): + return + faces.append([i, j, k]) + + for r in range(R - 1): + for k in range(n): + k2 = (k + 1) % n + add(vid(r, k), vid(r, k2), vid(r + 1, k2)) + add(vid(r, k), vid(r + 1, k2), vid(r + 1, k)) + for k in range(n): + add(vid(R - 1, k), vid(R - 1, (k + 1) % n), apex_i) + + if faces: + caps.append(trimesh.Trimesh(verts, np.array(faces), process=False)) + + return caps + + +def join_meshes(meshes): + """ + Concatenate meshes and merge coincident vertices into one (ideally watertight) trimesh. + The Python twin of the Grasshopper "Close Piece" weld. Returns None if nothing to join. + """ + meshes = [m for m in meshes if m is not None and len(m.faces) > 0] + if not meshes: + return None + + all_v = [] + all_f = [] + for m in meshes: + base = len(all_v) + all_v.extend(np.asarray(m.vertices).tolist()) + for f in np.asarray(m.faces): + all_f.append([int(f[0]) + base, int(f[1]) + base, int(f[2]) + base]) + + # process=True merges coincident vertices, welding the shared cap/piece boundary + joined = trimesh.Trimesh(np.array(all_v), np.array(all_f), process=True) + # make winding consistent / normals outward, so a watertight result is a proper "volume" + # (manifold3d booleans require this, and it fixes inverted/negative-volume pieces) + joined.fix_normals() + return joined + + +def spread_pieces(meshes, factor=0.5, center=None): + """ + Move each mesh radially away from the common center by factor*(its center - overall center), + for an exploded view (the Python twin of "Spread Pieces"). Returns new translated copies; + the inputs are left untouched. + """ + meshes = list(meshes) + centers = [np.asarray(m.bounds).mean(axis=0) for m in meshes] + if center is None: + center = np.mean(centers, axis=0) if centers else np.zeros(3) + else: + center = np.asarray(center, dtype=float)[:3] + + out = [] + for m, c in zip(meshes, centers): + moved = m.copy() + moved.apply_translation(factor * (c - center)) + out.append(moved) + return out + + +def mesh_boolean_fold(solid, features, signs=None): + """ + Fold an ordered sequence of boolean operations onto a solid mesh -- the Python twin of the + Grasshopper "Boolean Piece" component. + + solid: a trimesh.Trimesh (should be watertight; booleans on open meshes are unreliable). + features: list of trimesh.Trimesh to boolean in, IN ORDER. ("Feature" in the solid-modeling + sense -- an ordered additive/subtractive operation on a body.) + signs: list parallel to features; +1 = union, <=0 = subtract. Defaults to all subtract. + + Order matters: each step acts on the result of the previous, e.g. signs [+1, -1, +1, -1] + means union(f0), then subtract(f1), then union(f2), then subtract(f3). Uses trimesh's + exact 'manifold' backend (the manifold3d package), which is robust on clean manifolds. + Returns the resulting trimesh.Trimesh. + """ + try: + import manifold3d # noqa: F401 -- the exact boolean backend trimesh will use + except ImportError as e: + raise ImportError( + "mesh booleans need the 'manifold3d' package (pip install manifold3d)") from e + + import warnings + + features = list(features) + if signs is None: + signs = [-1] * len(features) + if len(signs) != len(features): + raise ValueError("signs must be parallel to features") + + if not solid.is_watertight: + warnings.warn("boolean solid is not watertight; the result may be wrong") + + result = solid.copy() + for feature, sign in zip(features, signs): + if sign > 0: + result = trimesh.boolean.union([result, feature], engine='manifold') + else: + result = trimesh.boolean.difference([result, feature], engine='manifold') + return result + + class SurfacePiece(): - """ - A "Piece" of an algebraic surface. Essentially, a union of Faces, with some additional interface. + """ + A "Piece" of an algebraic surface. Essentially, a union of Faces, with some additional interface. """ def __init__(self, indices, surface): @@ -247,7 +539,7 @@ def flatten_and_unique(list_nD): # type critical def point_singularities(self): - """ Compute singularity points from a SurfacePiece object + """ Compute the indices of the singularity points from a SurfacePiece object :rtype: A list of indices of point singularities """ @@ -405,6 +697,39 @@ def edge_pieces(self): + def to_gh_dict(self, piece_index, include_smooth=True): + """ + assemble this piece's data for the Grasshopper JSON export. + + meshes are expressed purely as triangle indices into the surface's unified vertex + set; embedded curves as ordered vertex-index lists into the same set. no vertex + coordinates live here -- they are shared at the top level of the export. + """ + + mesh_raw = _mesh_triangles(self.surface.as_mesh_raw(self.indices)) + + mesh_smooth = None + if include_smooth and self.surface.is_sampled(): + try: + mesh_smooth = _mesh_triangles(self.surface.as_mesh_smooth(self.indices)) + except br_except.SurfaceNotSampled: + mesh_smooth = None + + curves = [] + for cp in self.edge_pieces(): + curves.append({ + "type": self.surface._curve_type_for_name(cp.curve.inputfilename), + "curve_name": cp.curve.inputfilename, + "vertex_indices": cp.to_point_indices(), + }) + + return { + "piece_index": piece_index, + "face_indices": list(self.indices), + "mesh_smooth": mesh_smooth, + "mesh_raw": mesh_raw, + "curves": curves, + } @@ -436,6 +761,53 @@ def export_raw(self, basename=_default_piece_basename_raw,autoname_using_folder= self.surface.export_raw(self.indices,filename_no_ext,autoname_using_folder,file_type) + def as_mesh(self, smooth=None): + """ + The `trimesh.Trimesh` for this piece. smooth=None picks smooth when the surface is + sampled, else raw (raw with raw, sampled with sampled). + """ + if smooth is None: + smooth = self.surface.is_sampled() + if smooth: + return self.surface.as_mesh_smooth(self.indices) + return self.surface.as_mesh_raw(self.indices) + + + def sphere_caps(self, smooth=None, resolution=4, tol=1e-3): + """ + The faceted spherical cap mesh(es) closing this piece where it meets the bounding + sphere. See the module-level `sphere_cap_meshes`. Returns a list of trimesh. + """ + return sphere_cap_meshes(self.as_mesh(smooth), self.center, self.radius, resolution, tol) + + + def flat_caps(self, smooth=None, resolution=1, tol=1e-3): + """ + The faceted FLAT cap mesh(es) closing this piece where it meets the bounding sphere, with + the apex at each loop's centroid. See the module-level `flat_cap_meshes`. Returns a list + of trimesh. + """ + return flat_cap_meshes(self.as_mesh(smooth), self.center, self.radius, resolution, tol) + + + def as_closed_mesh(self, smooth=None, resolution=None, tol=1e-3, flat=False): + """ + This piece joined with its cap(s) into a single welded (ideally watertight) + `trimesh.Trimesh` -- the Rhino-free equivalent of (Sphere|Flat) Caps + Close Piece. A + piece bounded only by the sphere comes out watertight; one abutting a singular curve stays + open there (check `.is_watertight`). + + flat=False uses spherical caps (hugging the sphere); flat=True uses flat fans to the loop + centroid. resolution defaults to 4 for spherical, 1 for flat. + """ + if resolution is None: + resolution = 1 if flat else 4 + mesh = self.as_mesh(smooth) + if flat: + caps = flat_cap_meshes(mesh, self.center, self.radius, resolution, tol) + else: + caps = sphere_cap_meshes(mesh, self.center, self.radius, resolution, tol) + return join_meshes([mesh] + caps) def solidify_smooth(self, distance=_default_solidify_thickness, basename=_default_piece_basename_smooth, autoname_using_folder=False,file_type=_default_file_type): @@ -901,99 +1273,238 @@ def curve_with_name(self, curve_name): raise RuntimeError(f'unable to find a curve with name {curve_name} in this surface') - def write_piece_data(self): + + def _curve_type_for_name(self, curve_name): """ - Opens and edits current scad data to set the orientation and location of a plug and socket + classify an embedded curve by its `inputfilename` into one of the closed-vocabulary + type tags used by the Grasshopper export. match order mirrors `curve_with_name`. """ + if curve_name == self.critical_curve.inputfilename: + return "critical" + + if curve_name == self.sphere_curve.inputfilename: + return "sphere" + + for c in self.critical_point_slices: + if curve_name == c.inputfilename: + return "critslice" + + for c in self.midpoint_slices: + if curve_name == c.inputfilename: + return "midslice" + + if curve_name in self.singular_names: + return "singular" + + return "unknown" + + + def export_gh_json(self, filename="br_gh_export.json", include_smooth=True): + """ + write a self-contained JSON describing this surface for the Grasshopper plugin. + + the file holds one unified vertex set (`vertices`); each nonsingular piece carries + only triangle indices (raw and, when sampled, smooth) and the embedded curve pieces + as ordered vertex-index lists -- all indices into the shared `vertices`. this keeps + the surface mesh and its embedded curves referring to the same points in Rhino. + """ + + # prime the extract_points memo cache with the full (no-arg) point set first, so + # later per-piece mesh construction does not poison it with a partial set. + points = self.extract_points() + pieces = self.separate_into_nonsingular_pieces() - allPoints=[] - #create a list of the centroid coordinates of each piece - centroids = [] - for p in pieces: - centroids.append(p.centroid()) - + contents = { + "format_version": 2, + "decomposition_type": "surface", + "source_directory": self.directory, + "num_variables": self.num_variables, + "vertices": _points_to_xyz(points), + "vertex_count": len(points), + "sphere": self._sphere_dict(), + "is_sampled": self.is_sampled(), + "pieces": [p.to_gh_dict(ii, include_smooth) for ii, p in enumerate(pieces)], + } + + # fold the singularity / connector data (locations, tangent-cone directions, parities) + # into the same file, so Grasshopper has a single source instead of a second JSON. + try: + sing = self.singularity_connector_data() + except Exception as e: + print("WARNING: could not compute singularity connector data ({}); " + "exporting without singularities".format(e)) + sing = {"piece_names": [], "on_pieces": [], "locations": [], + "directions": [], "parities": []} + contents["singularities"] = { + "piece_names": sing["piece_names"], + "locations": sing["locations"], + "directions": sing["directions"], + "parities": sing["parities"], + "on_pieces": sing["on_pieces"], + } + + # verify each piece's sphere curves are closed loops (they always should be, barring a + # decomposition problem); warn loudly if not, so the issue is visible before Grasshopper. + for pc in contents["pieces"]: + for cv in pc["curves"]: + if cv["type"] == "sphere": + vi = cv["vertex_indices"] + if len(vi) < 2 or vi[0] != vi[-1]: + print("WARNING: piece {} has a non-closed sphere curve ({}); " + "the decomposition may be incomplete".format( + pc["piece_index"], cv["curve_name"])) + + with open(filename, "w") as f: + json.dump(contents, f, indent=2) + + print("wrote " + filename) + return filename + + + def all_curves(self): + + the_curves = [] + + the_curves.append(self.critical_curve) + + the_curves.append(self.sphere_curve) - # compute a list of nodal singularities, and which pieces they're connected to - pieces_connected_to_sing = defaultdict(list) - sings_on_pieces = {} #sings are in order of the piece index + for c in self.singular_curves: + the_curves.append(c) - for ii, p in enumerate(pieces): - sing_this_piece = p.point_singularities() - sings_on_pieces[ii] = sing_this_piece + for c in self.critical_point_slices: + the_curves.append(c) - # a dictionary keyed by the integer index of the singularity, with value a list of the pieces on which it is incident - for s in sing_this_piece: - pieces_connected_to_sing[s].append(ii) + for c in self.midpoint_slices: + the_curves.append(c) + return the_curves - # only put plug/socket at sing that's connected to two pieces - #then assign that sing to k(ey) and assign [piece1,piece2] to v(value) - wanted_sing_connections = {k:v for k,v in pieces_connected_to_sing.items() if len(v)==2} - def unit_vector(vector): - """Helper function to find a unit vector of a vector""" + def all_singular_points(self): + """ + get absolutely all of the singular points. + """ + + # there's a baked-in assumption that this surface is NOT contained in a higher-dimensional object. this is valid right now because the top-dimensional thing Bertini_real can decompose is a surface. - magnitude = np.linalg.norm(vector) - unit=[] - for i in range(len(vector)): - unit.append(vector[i]/magnitude) - return unit - directions = defaultdict(list) # explicitly keyed by the singularities - sing_directions = {} - sing_locations = {} - for sing_index,connected_pieces in wanted_sing_connections.items(): - ind_connected_piece_0 = connected_pieces[0] - ind_connected_piece_1 = connected_pieces[1] + the_singularites = [] + for v in self.vertices: + if v.is_of_type(VertexType.singular): + the_singularites.append(v) - # find the centroid of each piece by the index of the piece - centroid_0 = centroids[ind_connected_piece_0] - centroid_1 = centroids[ind_connected_piece_1] + return the_singularites - sing_coords =self.vertices[sing_index].point.real + def isolated_singularities(self): + VertexType = bertini_real.vertex.VertexType - #calculate the unit vectors by traveling from the centroid of the piece to the singularity - unit_0 = unit_vector(np.subtract(centroid_0, sing_coords)) - unit_1 = unit_vector(np.subtract(centroid_1, sing_coords)) + the_singularites = [] + for v in self.vertices: + if v.is_of_type(VertexType.singular) and v.is_of_type(VertexType.singular): + the_singularites.append(v) - #find unit vector resultant of unit_0 and flipped unit_1 - direction0 = unit_vector(np.add(unit_0, np.multiply(unit_1, -1))) - direction1 = np.multiply(direction0,-1) - directions[sing_index] = [direction0, direction1] + return the_singularites - sing_directions[sing_index] = (direction0) - sing_locations[sing_index] = (list(sing_coords)) - piece_names = [] - singularities_on_pieces = [] + def singularity_connector_data(self): + """ + Compute the data needed to place plug/socket connectors at nodal singularities. + + For each nodal singularity that joins exactly two nonsingular pieces, find its + location, the connector axis direction (the tangent-cone direction, from the Hessian + of the defining polynomial at the singularity), and the per-piece parity (which side + gets the plug vs the socket); also record which singularities lie on each piece. + + Needs the `bertini` parser and `sympy`, but only when there is at least one qualifying + singularity. Returns a dict of pure-Python (JSON-safe) values: + { "piece_names": [str, ...], # per piece + "on_pieces": [[int, ...], ...], # per piece: compact singularity indices + "locations": [[x, y, z], ...], # per singularity + "directions": [[x, y, z], ...], # per singularity + "parities": [[int, ...], ...] } # per singularity: a value per piece (-1/0/1) + All lists are empty when there are no qualifying singularities. + """ + pieces = self.separate_into_nonsingular_pieces() + piece_names = [p.generate_filename_smooth() for p in pieces] - singindex2int = {sing_index:ii for ii,sing_index in enumerate(wanted_sing_connections.keys())} - int2singindex = {ii:sing_index for ii,sing_index in enumerate(wanted_sing_connections.keys())} + # nodal singularities and which pieces each is incident to + sings_on_pieces = {} + pieces_connected_to_sing = defaultdict(list) + for ii, p in enumerate(pieces): + sings_this_piece = p.point_singularities() # indices into the vertex set + sings_on_pieces[ii] = sings_this_piece + for s in sings_this_piece: + pieces_connected_to_sing[s].append(ii) - # print(singindex2int) - # print(int2singindex) + # only singularities joining exactly two pieces receive a connector + wanted = {k: v for k, v in pieces_connected_to_sing.items() if len(v) == 2} + singindex2int = {s: i for i, s in enumerate(wanted.keys())} - #organize the data computed above to the scad files - for ii,p in enumerate(pieces): - sings_this_piece = [] + on_pieces = [] + for ii in range(len(pieces)): + on_pieces.append([singindex2int[s] for s in wanted if s in sings_on_pieces[ii]]) - for sing_index,connected_pieces in wanted_sing_connections.items(): - if sing_index in sings_on_pieces[ii]: - sings_this_piece.append(singindex2int[sing_index]) - singularities_on_pieces.append(sings_this_piece) - piece_names.append(pieces[ii].generate_filename_smooth()) + def unit_vector(vector): + return vector / np.linalg.norm(vector) - sing_directions_as_list = [sing_directions[sing_index] for sing_index in wanted_sing_connections.keys()] - sing_locations_as_list = [sing_locations[sing_index] for sing_index in wanted_sing_connections.keys()] + locations = [] + directions = [] + + if wanted: + import bertini as b2 + import sympy + + bsys = b2.parse.system(self.input.split('INPUT')[1]) + f = bsys.function(0) + F = sympy.S(str(f).replace('unnamed_function', '').replace('function', '').replace('f', '')) + variables = sorted(F.free_symbols, key=lambda s: s.name) + H = sympy.hessian(F, variables) + hessian_evalme = sympy.lambdify(variables, H, modules='numpy') + + for sing_index in wanted.keys(): + sing_coords = self.vertices[sing_index].point.real + + # tangent-cone direction: eigenvector of the Hessian belonging to the + # odd-one-out (smallest) eigenvalue + M = hessian_evalme(*sing_coords) + q = np.linalg.eig(M) + axis = np.real(q.eigenvectors[:, np.argmin(q.eigenvalues)]) + direction0 = unit_vector(np.asarray(axis, dtype=float)) + + directions.append([float(x) for x in direction0]) + locations.append([float(x) for x in sing_coords]) + + parities = [[0 for _ in range(len(pieces))] for _ in range(len(wanted))] + for s, ps in wanted.items(): + parities[singindex2int[s]][ps[0]] = -1 + parities[singindex2int[s]][ps[1]] = 1 + + return { + "piece_names": piece_names, + "on_pieces": on_pieces, + "locations": locations, + "directions": directions, + "parities": parities, + } + + + def write_piece_data(self): + """ + Opens and edits current scad data to set the orientation and location of a plug and socket + """ + + data = self.singularity_connector_data() + piece_names = data["piece_names"] + singularities_on_pieces = data["on_pieces"] + sing_directions_as_list = data["directions"] + sing_locations_as_list = data["locations"] + parity_of_sing_by_piece = data["parities"] + allPoints = [] - parity_of_sing_by_piece = [ [0 for jj in range(len(pieces))] for ii in range(len(wanted_sing_connections)) ] - for sing_index, ps in wanted_sing_connections.items(): - parity_of_sing_by_piece[singindex2int[sing_index]][ps[0]] = -1 - parity_of_sing_by_piece[singindex2int[sing_index]][ps[1]] = 1 - #open and auto write the data(piece file names (without extensions), all sings of pieces, sing directions in order of sing index, sing coords in order of sing index) of the piece with open("br_surf_piece_data.scad", "w") as f: f.write(f'piece_names = [') @@ -1008,20 +1519,33 @@ def unit_vector(vector): f.write(f'conn_size = 0.01;\n') #hard coded, but needs to be automatically computed print('br_surf_piece_data.scad') + # Option 2: custom JSON encoder + class NumpyEncoder(json.JSONEncoder): + def default(self, obj): + if isinstance(obj, np.ndarray): + return obj.tolist() + if isinstance(obj, np.integer): + return int(obj) + if isinstance(obj, np.floating): + return float(obj) + return super().default(obj) + + + #open and auto write piece data to a json file with open("br_surf_piece_data.json", "w") as j: j.write(json.dumps({"piece_names": piece_names, "singularities_on_pieces": singularities_on_pieces, "sing_directions": sing_directions_as_list, "sing_locations": sing_locations_as_list, - "parities" : parity_of_sing_by_piece},indent=4)) + "parities" : parity_of_sing_by_piece},indent=4,cls=NumpyEncoder)) print('wrote br_surf_piece_data.json') - with open("centroids.json", "w") as c: - for centroid in centroids: - c.write(str(centroid)+"\n") - print('wrote centroids.json') + # with open("centroids.json", "w") as c: + # for centroid in centroids: + # c.write(str(centroid)+"\n") + # print('wrote centroids.json') @@ -1029,6 +1553,7 @@ def unit_vector(vector): for point in allPoints: a.write("\n".join([str(s) for s in point]) + "\n") print('wrote allPoints.json') + def as_mesh_smooth(self, which_faces=None, keep_all_vertices=True): """ @@ -1216,19 +1741,15 @@ def as_mesh_raw(self, which_faces=None, keep_all_vertices=True): elif case == 5: break - t1 = [points[curr_edge[0]], points[curr_edge[1]], - points[face['midpoint']]] - t2 = [points[curr_edge[1]], points[curr_edge[2]], - points[face['midpoint']]] - - t3 = (curr_edge[0], curr_edge[1], face['midpoint']) - t4 = (curr_edge[1], curr_edge[2], face['midpoint']) - - T.append(t1) - T.append(t2) - - TT.append(t3) - TT.append(t4) + # fan the curve edge to the face midpoint; skip degenerate triangles (a repeated + # vertex index, e.g. from a degenerate curve edge). these are zero-area, and if + # kept they make the mesh non-manifold -- an (a, a, mid) face contributes the + # {a, mid} edge twice, which is what produced the 4-shared edges and duplicate + # faces in the raw mesh. + for tri in ((curr_edge[0], curr_edge[1], face['midpoint']), + (curr_edge[1], curr_edge[2], face['midpoint'])): + if len(set(tri)) == 3: + TT.append(tri) faces = [TT] vertex = [] @@ -1245,7 +1766,11 @@ def as_mesh_raw(self, which_faces=None, keep_all_vertices=True): face_np_array = np.array(face) - raw_mesh = trimesh.Trimesh(vertex_np_array, face_np_array) + # honor keep_all_vertices like as_mesh_smooth: process=False keeps the full global + # vertex set so the faces index into extract_points() (the unified set), and does not + # merge coincident-but-distinct vertices (which would fuse sheets at singularities). + should_trimesh_process = False if keep_all_vertices == True else True + raw_mesh = trimesh.Trimesh(vertex_np_array, face_np_array, process=should_trimesh_process) raw_mesh.fix_normals() return raw_mesh diff --git a/python/docs/index.rst b/python/docs/index.rst index 2ed6ea5..5bd9904 100644 --- a/python/docs/index.rst +++ b/python/docs/index.rst @@ -35,6 +35,7 @@ Tutorials ✏️ tutorials/matplotlib tutorials/glumpy tutorials/mesh_export + tutorials/capping_and_joining tutorials/anaglypy tutorials/snap_together .. tutorials/grasshopper_connection diff --git a/python/docs/tutorials/capping_and_joining.rst b/python/docs/tutorials/capping_and_joining.rst new file mode 100644 index 0000000..d6951cc --- /dev/null +++ b/python/docs/tutorials/capping_and_joining.rst @@ -0,0 +1,118 @@ +Capping pieces and joining them with connectors 🔩 +=========================================================================== + +After ``bertini_real`` decomposes a surface, you can split it into its nonsingular +**pieces**, **cap** each piece where it was cut by the bounding sphere, and then do solid +modeling on the resulting closed meshes -- for example, drilling a hole in one piece and +adding a matching rod to its neighbor so they snap together. + +Everything below is pure Python (using ``trimesh``); the same operations are available as +Grasshopper components (*Sphere Caps* / *Flat Caps*, *Close Piece*, *Boolean Piece*). + +.. note:: + The mesh booleans use ``trimesh``'s exact ``manifold`` backend, so you need the + ``manifold3d`` package installed (``pip install manifold3d``). Watertight solids require + a **sampled** decomposition. + +We use **Ding Dong** here 🔔 -- two pieces (a wide cone and a small bell) that meet at a +single nodal singularity. + +:: + + import bertini_real + surface = bertini_real.data.gather() # run from the decomposition's folder + pieces = surface.separate_into_nonsingular_pieces() + + +Spherical vs. flat caps +*********************************** + +Each piece's opening on the bounding sphere can be filled two ways. ``as_closed_mesh()`` +gives a **spherical** cap that continues the sphere; ``as_closed_mesh(flat=True)`` gives a +**flat** fan across the opening. The sphere still shows where the piece was cut. + +:: + + domed = pieces[0].as_closed_mesh() # spherical cap (default) + flat = pieces[0].as_closed_mesh(flat=True) # flat cap + +Spherical caps -- the underside bulges to follow the sphere: + +.. image:: capping_and_joining_pictures/dingdong_sphere_caps.png + :width: 350 + +Flat caps -- the same opening filled flat (interesting on blocky/raw decompositions): + +.. image:: capping_and_joining_pictures/dingdong_flat_caps.png + :width: 350 + + +Joining two pieces with a rod +*********************************** + +The two pieces meet at a singularity; ``singularity_connector_data()`` gives its location +and a direction (the tangent-cone axis) plus a per-piece parity that says which side gets +the **socket** (the hole) and which gets the **plug** (the rod). + +:: + + import numpy as np + import trimesh + from bertini_real.surface import mesh_boolean_fold, spread_pieces + + sing = surface.singularity_connector_data() + location = np.array(sing["locations"][0]) + direction = np.array(sing["directions"][0]); direction /= np.linalg.norm(direction) + parity = sing["parities"][0] + socket = pieces[parity.index(-1)].as_closed_mesh() # this one gets the hole + plug = pieces[parity.index(1)].as_closed_mesh() # this one gets the rod + +Build a square prism at the singularity, oriented along the direction. Make the rod a +little smaller than the hole (``clearance``) so it can slide into the socket: + +:: + + def square_prism(width, length): + box = trimesh.creation.box(extents=(width, width, length)) # centered, along Z + box.apply_transform(trimesh.geometry.align_vectors([0, 0, 1.0], direction)) + box.apply_translation(location) + return box + + width, length, clearance = 0.30, 1.40, 0.85 + hole = square_prism(width, length) + rod = square_prism(width * clearance, length) + +Fold the booleans onto each piece -- subtract the hole (``-1``), union the rod (``+1``): + +:: + + socket_holed = mesh_boolean_fold(socket, [hole], [-1]) + plug_rodded = mesh_boolean_fold(plug, [rod], [+1]) + +The plug, with a square rod grown out of it at the singularity: + +.. image:: capping_and_joining_pictures/dingdong_square_rod.png + :width: 350 + +Pulled apart (``spread_pieces``) you can see the rod on the cone and the receiving hole in +the underside of the bell: + +:: + + a, b = spread_pieces([socket_holed, plug_rodded], factor=0.8) + +.. image:: capping_and_joining_pictures/dingdong_exploded.png + :width: 350 + +And seated together -- the rod joins the two pieces: + +.. image:: capping_and_joining_pictures/dingdong_joined.png + :width: 350 + +Export anything for printing with ``trimesh``'s ``mesh.export("piece.stl")``, or drive the +whole flow from the command line with ``python/scripts/close_pieces.py`` (it takes +``--flat``, ``--spread``, and ``--smooth/--raw``). + +---- + +*The figures above were generated by* ``capping_and_joining_pictures/make_images.py``. diff --git a/python/docs/tutorials/capping_and_joining_pictures/dingdong_exploded.png b/python/docs/tutorials/capping_and_joining_pictures/dingdong_exploded.png new file mode 100644 index 0000000..0289156 Binary files /dev/null and b/python/docs/tutorials/capping_and_joining_pictures/dingdong_exploded.png differ diff --git a/python/docs/tutorials/capping_and_joining_pictures/dingdong_flat_caps.png b/python/docs/tutorials/capping_and_joining_pictures/dingdong_flat_caps.png new file mode 100644 index 0000000..c7f4c40 Binary files /dev/null and b/python/docs/tutorials/capping_and_joining_pictures/dingdong_flat_caps.png differ diff --git a/python/docs/tutorials/capping_and_joining_pictures/dingdong_joined.png b/python/docs/tutorials/capping_and_joining_pictures/dingdong_joined.png new file mode 100644 index 0000000..cf923e8 Binary files /dev/null and b/python/docs/tutorials/capping_and_joining_pictures/dingdong_joined.png differ diff --git a/python/docs/tutorials/capping_and_joining_pictures/dingdong_sphere_caps.png b/python/docs/tutorials/capping_and_joining_pictures/dingdong_sphere_caps.png new file mode 100644 index 0000000..be390db Binary files /dev/null and b/python/docs/tutorials/capping_and_joining_pictures/dingdong_sphere_caps.png differ diff --git a/python/docs/tutorials/capping_and_joining_pictures/dingdong_square_hole.png b/python/docs/tutorials/capping_and_joining_pictures/dingdong_square_hole.png new file mode 100644 index 0000000..bc43db8 Binary files /dev/null and b/python/docs/tutorials/capping_and_joining_pictures/dingdong_square_hole.png differ diff --git a/python/docs/tutorials/capping_and_joining_pictures/dingdong_square_rod.png b/python/docs/tutorials/capping_and_joining_pictures/dingdong_square_rod.png new file mode 100644 index 0000000..1b51842 Binary files /dev/null and b/python/docs/tutorials/capping_and_joining_pictures/dingdong_square_rod.png differ diff --git a/python/docs/tutorials/capping_and_joining_pictures/make_images.py b/python/docs/tutorials/capping_and_joining_pictures/make_images.py new file mode 100644 index 0000000..f98d11d --- /dev/null +++ b/python/docs/tutorials/capping_and_joining_pictures/make_images.py @@ -0,0 +1,114 @@ +""" +Regenerate the figures for tutorials/capping_and_joining.rst from the dingdong decomposition. + +Run from the repo with the package importable, e.g.: + PYTHONPATH=python python python/docs/tutorials/capping_and_joining_pictures/make_images.py \ + --dingdong test/surface/dingdong/output_dim_2_comp_0 +Needs: trimesh, manifold3d (mesh booleans), bertini (tangent-cone directions), matplotlib. +""" +import argparse +import os + +import numpy as np +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d.art3d import Poly3DCollection + +import trimesh +from bertini_real.surface import Surface, mesh_boolean_fold, spread_pieces + +HERE = os.path.dirname(os.path.abspath(__file__)) + + +def render(items, outpath, elev=22, azim=-60): + """items: list of (trimesh, rgb-tuple).""" + fig = plt.figure(figsize=(5, 5)) + ax = fig.add_subplot(111, projection="3d") + allpts = [] + for mesh, color in items: + V = np.asarray(mesh.vertices) + F = np.asarray(mesh.faces) + tris = V[F] + pc = Poly3DCollection(tris) + shade = 0.45 + 0.55 * np.clip(mesh.face_normals[:, 2], 0, 1) + cols = np.clip(np.array(color)[None, :] * shade[:, None], 0, 1) + pc.set_facecolor(cols) + pc.set_edgecolor((0, 0, 0, 0.10)) + pc.set_linewidth(0.1) + ax.add_collection3d(pc) + allpts.append(V) + P = np.vstack(allpts) + mins, maxs = P.min(0), P.max(0) + c, r = (mins + maxs) / 2.0, (maxs - mins).max() / 2.0 + ax.set_xlim(c[0] - r, c[0] + r) + ax.set_ylim(c[1] - r, c[1] + r) + ax.set_zlim(c[2] - r, c[2] + r) + ax.set_axis_off() + ax.view_init(elev=elev, azim=azim) + fig.tight_layout() + fig.savefig(os.path.join(HERE, outpath), dpi=130) + plt.close(fig) + print("wrote", outpath) + + +def square_prism(width, length, loc, direction): + box = trimesh.creation.box(extents=(width, width, length)) # centered at origin, along Z + box.apply_transform(trimesh.geometry.align_vectors([0, 0, 1.0], direction)) + box.apply_translation(loc) + return box + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--dingdong", default="test/surface/dingdong/output_dim_2_comp_0") + args = ap.parse_args() + + s = Surface(args.dingdong) + pieces = s.separate_into_nonsingular_pieces() + blue = (0.36, 0.56, 0.85) + gold = (0.85, 0.66, 0.30) + + # --- flat vs spherical caps, on the piece that meets the sphere --- + # the big "cone" piece has a large base cap; a low side view shows the dome-vs-flat silhouette + cap_piece = max(range(len(pieces)), key=lambda i: pieces[i].as_mesh().bounding_box.volume) + render([(pieces[cap_piece].as_closed_mesh(flat=False), blue)], "dingdong_sphere_caps.png", elev=8) + render([(pieces[cap_piece].as_closed_mesh(flat=True), blue)], "dingdong_flat_caps.png", elev=8) + + # --- square hole in one piece, matching (smaller) rod in the other --- + sd = s.singularity_connector_data() + loc = np.array(sd["locations"][0], dtype=float) + direction = np.array(sd["directions"][0], dtype=float) + direction = direction / np.linalg.norm(direction) + parity = sd["parities"][0] + socket_idx = parity.index(-1) # gets the hole + plug_idx = parity.index(1) # gets the rod + + # size the rod/hole to the SMALL (socket) piece so the hole leaves a ring of material + width = 0.30 + length = 1.4 # centered on the singularity; reaches into both pieces + clearance = 0.85 # rod a bit smaller than the hole, so it slides in + + socket = pieces[socket_idx].as_closed_mesh() + plug = pieces[plug_idx].as_closed_mesh() + + hole = square_prism(width, length, loc, direction) + rod = square_prism(width * clearance, length, loc, direction) + + socket_holed = mesh_boolean_fold(socket, [hole], [-1]) + plug_rodded = mesh_boolean_fold(plug, [rod], [+1]) + + render([(socket_holed, blue)], "dingdong_square_hole.png") + render([(plug_rodded, gold)], "dingdong_square_rod.png") + + # assembled (rod seated in the hole) and exploded + render([(socket_holed, blue), (plug_rodded, gold)], "dingdong_joined.png") + exploded = spread_pieces([socket_holed, plug_rodded], factor=0.8) + render([(exploded[0], blue), (exploded[1], gold)], "dingdong_exploded.png") + + print("watertight: socket_holed={}, plug_rodded={}".format( + socket_holed.is_watertight, plug_rodded.is_watertight)) + + +if __name__ == "__main__": + main() diff --git a/python/scripts/close_pieces.py b/python/scripts/close_pieces.py new file mode 100644 index 0000000..8f23b42 --- /dev/null +++ b/python/scripts/close_pieces.py @@ -0,0 +1,78 @@ +# Rhino-free pipeline: split a surface decomposition into its nonsingular pieces, cap each +# where it meets the bounding sphere, join into (ideally watertight) solids, optionally spread +# them apart for viewing, and export each as an STL. This mirrors the Grasshopper components +# (Sphere Caps + Close Piece + Spread Pieces) for people who don't have a Rhino license. +# +# Watertight solids require a *sampled* decomposition: the raw (blocky) mesh is non-manifold, +# so raw pieces close visually but are not watertight. The script reports this honestly. +# +# Usage (from a decomposition's working directory, or pass the output_dim_2_comp_0 folder): +# python close_pieces.py [folder] [--resolution N] [--spread F] [--smooth | --raw] +# [--basename NAME] [--combined] + +import argparse +import sys + +import bertini_real as br +from bertini_real.surface import Surface, spread_pieces, join_meshes + + +def main(argv=None): + ap = argparse.ArgumentParser(description="Close surface pieces with sphere caps and export STLs (no Rhino needed).") + ap.add_argument("folder", nargs="?", default=None, + help="path to an output_dim_2_comp_0 folder; omit to gather() from the current directory") + ap.add_argument("--resolution", type=int, default=None, help="cap subdivisions (default 4 spherical, 1 flat)") + ap.add_argument("--flat", action="store_true", help="flat caps (fan to the loop centroid) instead of spherical") + ap.add_argument("--spread", type=float, default=0.0, help="exploded-view factor; 0 = pieces in place (default 0)") + ap.add_argument("--smooth", dest="smooth", action="store_true", default=None, help="force sampled (smooth) meshes") + ap.add_argument("--raw", dest="smooth", action="store_false", help="force raw (blocky) meshes") + ap.add_argument("--basename", default="br_closed_piece", help="output STL basename (default br_closed_piece)") + ap.add_argument("--combined", action="store_true", help="also export one combined STL of all pieces") + args = ap.parse_args(argv) + + if args.folder: + surface = Surface(args.folder) + else: + decomposition = br.data.gather() + if not isinstance(decomposition, Surface): + sys.exit("this script is for surface (dimension 2) decompositions") + surface = decomposition + + pieces = surface.separate_into_nonsingular_pieces() + print("{} nonsingular piece(s)".format(len(pieces))) + + closed = [] + for i, piece in enumerate(pieces): + mesh = piece.as_closed_mesh(smooth=args.smooth, resolution=args.resolution, flat=args.flat) + closed.append(mesh) + if mesh is None: + print(" piece {}: empty".format(i)) + else: + print(" piece {}: {} faces, watertight={}".format(i, len(mesh.faces), mesh.is_watertight)) + + # report the singularity / connector data too (locations + tangent-cone directions) + try: + sing = surface.singularity_connector_data() + print("{} nodal singularity connector(s)".format(len(sing["locations"]))) + except Exception as e: + print("(singularity connector data unavailable: {})".format(e)) + + meshes = [m for m in closed if m is not None] + if args.spread: + meshes = spread_pieces(meshes, factor=args.spread) + print("spread pieces apart by factor {}".format(args.spread)) + + for i, mesh in enumerate(meshes): + outname = "{}_{}.stl".format(args.basename, i) + mesh.export(outname) + print("wrote {}".format(outname)) + + if args.combined and meshes: + combined = join_meshes(meshes) + if combined is not None: + combined.export(args.basename + "_all.stl") + print("wrote {}_all.stl".format(args.basename)) + + +if __name__ == "__main__": + main() diff --git a/python/scripts/export_for_grasshopper.py b/python/scripts/export_for_grasshopper.py new file mode 100644 index 0000000..fee5804 --- /dev/null +++ b/python/scripts/export_for_grasshopper.py @@ -0,0 +1,70 @@ +# Generate the intermediary JSON consumed by the bertini_real Grasshopper components. +# +# This writes one self-contained file (default: br_gh_export.json) holding a single unified +# vertex set; for a surface, each nonsingular piece is a mesh (triangle indices) plus the +# curves embedded on it (vertex-index polylines); for a curve, each piece is a polyline. +# Point the "Surface Read GH JSON" (or "Curve Read GH JSON") Grasshopper component at it. +# +# Usage: +# # from a decomposition's working directory (the folder containing `Dir_Name`): +# python export_for_grasshopper.py [output.json] +# +# # or point it explicitly at an output_dim_X_comp_Y folder from anywhere: +# python export_for_grasshopper.py path/to/output_dim_2_comp_0 [output.json] +# +# Works for both curves (dim 1) and surfaces (dim 2). + +import os +import re +import sys + +import bertini_real as br +from bertini_real.curve import Curve +from bertini_real.surface import Surface + + +def export_from_directory(directory, filename): + """Build the decomposition directly from an output_dim_X_comp_Y folder and export it.""" + basename = os.path.basename(os.path.normpath(directory)) + match = re.search(r"dim_(\d+)", basename) + if not match: + raise SystemExit( + f"could not determine the dimension from folder name {basename!r}; " + "expected something like output_dim_2_comp_0" + ) + + dimension = int(match.group(1)) + if dimension == 1: + decomposition = Curve(directory) + elif dimension == 2: + decomposition = Surface(directory) + else: + raise SystemExit(f"dimension {dimension} not supported (only curves=1 and surfaces=2)") + + return decomposition.export_gh_json(filename) + + +def main(argv): + directory = None + filename = "br_gh_export.json" + + # a positional arg that is an existing directory is the decomposition folder; + # any other positional arg is the output filename. + for arg in argv[1:]: + if os.path.isdir(arg): + directory = arg + else: + filename = arg + + if directory is not None: + out = export_from_directory(directory, filename) + else: + # idiomatic path: gather from the current working directory (needs `Dir_Name`) + out = br.data.gather_and_export_gh(filename) + + # export_gh_json already prints the filename it wrote; point at the next step. + print(f"open {os.path.abspath(out)} with the 'Surface Read GH JSON' / 'Curve Read GH JSON' component") + + +if __name__ == "__main__": + main(sys.argv) diff --git a/python/scripts/prep_surf_for_grasshopper.py b/python/scripts/prep_surf_for_grasshopper.py new file mode 100644 index 0000000..8753540 --- /dev/null +++ b/python/scripts/prep_surf_for_grasshopper.py @@ -0,0 +1,45 @@ + + +import bertini_real as br +import json +import os +import numpy as np + + +class NumpyEncoder(json.JSONEncoder): + """ Special json encoder for numpy types """ + def default(self, obj): + if isinstance(obj, (np.int_, np.intc, np.intp, np.int8, + np.int16, np.int32, np.int64, np.uint8, + np.uint16, np.uint32, np.uint64)): + return int(obj) + elif isinstance(obj, (np.float16, np.float32, + np.float64)): + return float(obj) + elif isinstance(obj, (np.ndarray,)): + return obj.tolist() + return json.JSONEncoder.default(self, obj) + +surf = br.data.read_most_recent() + +surface_name = 'br_piece' # os.getcwd().split('/')[-1] +surf.write_piece_data() + +pieces = surf.separate_into_nonsingular_pieces() + +for p in pieces: + p.export_smooth() # basename=f"{surface_name}_smooth" + p.export_raw() # basename=f"{surface_name}_raw" + p.solidify_smooth(0.02, basename=f"{surface_name}_solidified_smooth_0.02") # basename=f"{surface_name}_solidified_smooth_0.02" + p.solidify_raw(0.02, basename=f"{surface_name}_solidified_smooth_0.02") # basename=f"{surface_name}_solidified_raw_0.02" + + edge_pieces = p.edge_pieces() + + edge_pieces_as_points = [] + for ii,edge_piece in enumerate(edge_pieces): + edge_pieces_as_points.append( (edge_piece.inputfilename()+"_"+str(ii),edge_piece.to_points()) ) + + with open(p.generate_filename_no_ext(basename="touching_curves")+'.json','w') as f: + + # print(edge_pieces_as_points) + json.dump(edge_pieces_as_points, f,cls=NumpyEncoder, indent=4) \ No newline at end of file diff --git a/python/setup.py b/python/setup.py index 85dde02..ee572cc 100644 --- a/python/setup.py +++ b/python/setup.py @@ -4,6 +4,9 @@ extras = { 'optional': [ + ], + 'test': [ + 'pytest', ] } @@ -28,7 +31,9 @@ 'algopy', 'sympy', 'scipy', - 'networkx'], + 'networkx', + 'bertini2>=3.0.0', + 'manifold3d'], extras_require=extras, package_dir={'bertini_real': 'bertini_real'}, package_data={'bertini_real': ['surface/scad/*.scad']}, # for plugs and sockets on pieces of surfaces diff --git a/python/tests/conftest.py b/python/tests/conftest.py new file mode 100644 index 0000000..707c94a --- /dev/null +++ b/python/tests/conftest.py @@ -0,0 +1,10 @@ +""" +Make the tests import the in-repo `bertini_real` package even when a copy is also installed +in site-packages, by putting the repo `python/` directory first on sys.path. +""" +import os +import sys + +_PYTHON_DIR = os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) +if _PYTHON_DIR not in sys.path: + sys.path.insert(0, _PYTHON_DIR) diff --git a/python/tests/test_booleans.py b/python/tests/test_booleans.py new file mode 100644 index 0000000..38d9e12 --- /dev/null +++ b/python/tests/test_booleans.py @@ -0,0 +1,52 @@ +""" +Tests for mesh_boolean_fold (the ordered boolean fold; Python twin of the Grasshopper +Boolean Piece component). Uses synthetic primitives, so no decomposition data is needed. +""" +import pytest + +trimesh = pytest.importorskip("trimesh") +pytest.importorskip("manifold3d") # the exact boolean backend + + +def test_default_is_subtract(): + from bertini_real.surface import mesh_boolean_fold + box = trimesh.creation.box(extents=(2, 2, 2)) # volume 8 + drill = trimesh.creation.cylinder(radius=0.3, height=4) + out = mesh_boolean_fold(box, [drill]) # no signs -> subtract + assert out.is_watertight + assert out.volume < box.volume + + +def test_union_then_subtract_order(): + from bertini_real.surface import mesh_boolean_fold + box = trimesh.creation.box(extents=(2, 2, 2)) + add = trimesh.creation.box(extents=(2, 2, 2)) + add.apply_translation((1.5, 0, 0)) # overlaps box -> grows it + drill = trimesh.creation.cylinder(radius=0.3, height=10) + + out = mesh_boolean_fold(box, [add, drill], signs=[+1, -1]) + assert out.is_watertight + # the drill pierces the unioned body, so the hole is present in the combined solid + assert out.volume < trimesh.boolean.union([box, add], engine='manifold').volume + + +def test_order_matters(): + """subtract-then-union differs from union-then-subtract when the tools overlap.""" + from bertini_real.surface import mesh_boolean_fold + box = trimesh.creation.box(extents=(2, 2, 2)) + add = trimesh.creation.box(extents=(2, 2, 2)) + add.apply_translation((1.0, 0, 0)) + drill = trimesh.creation.cylinder(radius=0.4, height=10) + drill.apply_translation((1.0, 0, 0)) # sits where `add` will be + + union_first = mesh_boolean_fold(box, [add, drill], signs=[+1, -1]) + subtract_first = mesh_boolean_fold(box, [drill, add], signs=[-1, +1]) + # union-then-subtract leaves the hole; subtract-then-union back-fills it -> larger volume + assert subtract_first.volume > union_first.volume + + +def test_signs_length_mismatch_raises(): + from bertini_real.surface import mesh_boolean_fold + box = trimesh.creation.box(extents=(1, 1, 1)) + with pytest.raises(ValueError): + mesh_boolean_fold(box, [box], signs=[+1, -1]) diff --git a/python/tests/test_gh_export.py b/python/tests/test_gh_export.py new file mode 100644 index 0000000..754dbb3 --- /dev/null +++ b/python/tests/test_gh_export.py @@ -0,0 +1,221 @@ +""" +Tests for the Grasshopper JSON export (Surface.export_gh_json / Curve.export_gh_json and the +helpers they rely on). + +Two tiers: + * fast pure-logic tests with no decomposition data; and + * fixture-based invariant tests that run against real decompositions under the repo `test/` + tree when present (skipped otherwise, so the suite still passes without the data). +""" +import json +import os + +import numpy as np +import pytest + + +# the closed vocabulary of curve-type tags the C# side branches on +VOCAB = {"critical", "sphere", "singular", "midslice", "critslice", "unknown", "standalone"} + +_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) + + +def _decomp(rel): + full = os.path.join(_REPO, rel) + return full if os.path.isdir(full) else None + + +WHITNEY = _decomp("test/surface/whitney/output_dim_2_comp_0") # sampled, has a singular curve +SPHERE = _decomp("test/surface/sphere/output_dim_2_comp_0") # unsampled +EISTUTE = _decomp("test/curve/intersections_of_surfaces/eistute_sphere/output_dim_1_comp_0") +NORDSTRAND = _decomp("test/surface/nordstrands_weird/output_dim_2_comp_0") # nodal singularities + + +# --------------------------------------------------------------------------- # +# fast, pure-logic tests (no decomposition data required) +# --------------------------------------------------------------------------- # + +def test_mesh_triangles_flattens(): + from bertini_real.surface import _mesh_triangles + + class FakeMesh: + faces = np.array([[0, 1, 2], [2, 3, 4]]) + + out = _mesh_triangles(FakeMesh()) + assert out["triangles"] == [0, 1, 2, 2, 3, 4] + assert out["triangle_count"] == 2 + + +def test_mesh_triangles_none(): + from bertini_real.surface import _mesh_triangles + assert _mesh_triangles(None) is None + + +def test_points_to_xyz_pads_and_truncates(): + from bertini_real.curve import _points_to_xyz + out = _points_to_xyz([[1.0], [1.0, 2.0], [1.0, 2.0, 3.0, 4.0]]) + assert out == [[1.0, 0.0, 0.0], [1.0, 2.0, 0.0], [1.0, 2.0, 3.0]] + + +def test_curve_type_for_name(): + from bertini_real.surface import Surface + + s = Surface.__new__(Surface) # bypass __init__/file IO + + class C: + def __init__(self, name): + self.inputfilename = name + + s.critical_curve = C("input_critical_curve") + s.sphere_curve = C("input_surf_sphere") + s.critical_point_slices = [C("crit0")] + s.midpoint_slices = [C("mid0")] + s.singular_names = ["sing0"] + + assert s._curve_type_for_name("input_critical_curve") == "critical" + assert s._curve_type_for_name("input_surf_sphere") == "sphere" + assert s._curve_type_for_name("crit0") == "critslice" + assert s._curve_type_for_name("mid0") == "midslice" + assert s._curve_type_for_name("sing0") == "singular" + assert s._curve_type_for_name("nope") == "unknown" + + +# --------------------------------------------------------------------------- # +# fixture-based invariant tests (need real decomposition data + trimesh) +# --------------------------------------------------------------------------- # + +trimesh = pytest.importorskip("trimesh") + + +def _check_surface_invariants(contents, surface): + assert contents["decomposition_type"] == "surface" + assert contents["vertex_count"] == len(contents["vertices"]) + assert all(len(v) == 3 for v in contents["vertices"]) + assert len(contents["sphere"]["center"]) == 3 + assert contents["sphere"]["radius"] > 0 + assert len(contents["pieces"]) == len(surface.separate_into_nonsingular_pieces()) + + vc = contents["vertex_count"] + for p in contents["pieces"]: + assert p["mesh_raw"] is not None + if not contents["is_sampled"]: + assert p["mesh_smooth"] is None + for key in ("mesh_raw", "mesh_smooth"): + m = p[key] + if m is None: + continue + assert len(m["triangles"]) == 3 * m["triangle_count"] + if m["triangles"]: + assert 0 <= min(m["triangles"]) and max(m["triangles"]) < vc + for cv in p["curves"]: + assert cv["type"] in VOCAB + if cv["vertex_indices"]: + assert max(cv["vertex_indices"]) < vc + + +@pytest.mark.skipif(not WHITNEY, reason="whitney decomposition not present") +def test_surface_export_invariants_whitney(tmp_path): + from bertini_real.surface import Surface + s = Surface(WHITNEY) + contents = json.load(open(s.export_gh_json(str(tmp_path / "w.json")))) + _check_surface_invariants(contents, s) + # whitney has a singular curve -- make sure that type is actually emitted + types = {cv["type"] for p in contents["pieces"] for cv in p["curves"]} + assert "singular" in types + + +@pytest.mark.skipif(not WHITNEY, reason="whitney decomposition not present") +def test_index_polyline_matches_points_whitney(): + """to_point_indices mapped through the unified vertices must equal the proven to_points().""" + from bertini_real.surface import Surface + s = Surface(WHITNEY) + pts = s.extract_points() + for piece in s.separate_into_nonsingular_pieces(): + for cp in piece.edge_pieces(): + idx = cp.to_point_indices() + via_idx = np.array([pts[i] for i in idx])[:, :3] if idx else np.empty((0, 3)) + via_pts = cp.to_points() + assert via_idx.shape == via_pts.shape + if via_pts.size: + assert np.allclose(via_idx, via_pts) + + +@pytest.mark.skipif(not WHITNEY, reason="whitney decomposition not present") +def test_piece_as_closed_mesh_and_spread(): + """Rhino-free pipeline: sampled pieces cap+join to watertight solids, and spread moves them.""" + from bertini_real.surface import Surface, spread_pieces + s = Surface(WHITNEY) + pieces = s.separate_into_nonsingular_pieces() + closed = [p.as_closed_mesh(resolution=4) for p in pieces] + assert all(m is not None and m.is_watertight for m in closed) + + moved = spread_pieces(closed, factor=0.5) + assert len(moved) == len(closed) + # spreading must actually displace at least one piece, and leave the originals untouched + assert any(not np.allclose(a.bounds.mean(axis=0), b.bounds.mean(axis=0)) + for a, b in zip(closed, moved)) + + +@pytest.mark.skipif(not WHITNEY, reason="whitney decomposition not present") +def test_piece_flat_caps_watertight(): + """The flat-cap alternative also welds to a watertight solid.""" + from bertini_real.surface import Surface + s = Surface(WHITNEY) + for p in s.separate_into_nonsingular_pieces(): + flat = p.as_closed_mesh(flat=True) + assert flat is not None and flat.is_watertight + + +@pytest.mark.skipif(not SPHERE, reason="unsampled sphere decomposition not present") +def test_unsampled_surface_has_null_smooth(tmp_path): + from bertini_real.surface import Surface + s = Surface(SPHERE) + contents = json.load(open(s.export_gh_json(str(tmp_path / "s.json")))) + assert contents["is_sampled"] is False + _check_surface_invariants(contents, s) + for p in contents["pieces"]: + assert p["mesh_smooth"] is None + + +@pytest.mark.skipif(not NORDSTRAND, reason="nordstrand decomposition not present") +def test_surface_export_singularities(tmp_path): + pytest.importorskip("bertini") # tangent-cone directions need the bertini parser + from bertini_real.surface import Surface + s = Surface(NORDSTRAND) + contents = json.load(open(s.export_gh_json(str(tmp_path / "n.json")))) + sg = contents["singularities"] + n_pieces = len(contents["pieces"]) + + # one direction per singularity, both xyz + assert len(sg["locations"]) == len(sg["directions"]) + assert all(len(p) == 3 for p in sg["locations"]) + assert all(len(d) == 3 for d in sg["directions"]) + + # one parity row per singularity, one entry per piece; a connector joins exactly two pieces + assert len(sg["parities"]) == len(sg["locations"]) + for row in sg["parities"]: + assert len(row) == n_pieces + assert row.count(-1) == 1 and row.count(1) == 1 + + # on_pieces is per piece; its singularity ids are in range + assert len(sg["on_pieces"]) == n_pieces + n_sing = len(sg["locations"]) + for ids in sg["on_pieces"]: + for i in ids: + assert 0 <= i < n_sing + + +@pytest.mark.skipif(not EISTUTE, reason="eistute_sphere curve decomposition not present") +def test_curve_export_invariants(tmp_path): + from bertini_real.curve import Curve + c = Curve(EISTUTE) + contents = json.load(open(c.export_gh_json(str(tmp_path / "c.json")))) + assert contents["decomposition_type"] == "curve" + assert contents["vertex_count"] == len(contents["vertices"]) + assert len(contents["sphere"]["center"]) == 3 + assert contents["sphere"]["radius"] > 0 + vc = contents["vertex_count"] + for p in contents["curve_pieces"]: + assert p["type"] in VOCAB + if p["vertex_indices"]: + assert max(p["vertex_indices"]) < vc diff --git a/python/tests/test_raw_mesh.py b/python/tests/test_raw_mesh.py new file mode 100644 index 0000000..8b0f577 --- /dev/null +++ b/python/tests/test_raw_mesh.py @@ -0,0 +1,78 @@ +""" +Tests for as_mesh_raw producing a clean, manifold mesh. + +Regression coverage for the fix where as_mesh_raw emitted degenerate triangles (from +degenerate curve edges) -- a triangle with a repeated vertex index contributes a self-edge +that gets counted twice, which made the raw mesh non-manifold (edges shared by 4 faces) and +introduced duplicate faces, so raw pieces could not be closed into watertight solids. Also +covers honoring keep_all_vertices so the raw faces index the global (unified) vertex set. +""" +import os +from collections import Counter + +import numpy as np +import pytest + +_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) + + +def _decomp(rel): + full = os.path.join(_REPO, rel) + return full if os.path.isdir(full) else None + + +WHITNEY = _decomp("test/surface/whitney/output_dim_2_comp_0") # sampled (raw still derivable) +NORDSTRAND = _decomp("test/surface/nordstrands_weird/output_dim_2_comp_0") # raw, nodal singularities + +trimesh = pytest.importorskip("trimesh") + + +def _max_edge_face_count(mesh): + ec = Counter() + for tri in mesh.faces: + a, b, c = int(tri[0]), int(tri[1]), int(tri[2]) + for u, v in ((a, b), (b, c), (c, a)): + ec[frozenset((u, v))] += 1 + return max(ec.values(), default=0) + + +def _num_duplicate_faces(mesh): + F = np.asarray(mesh.faces) + return len(F) - len(np.unique(np.sort(F, axis=1), axis=0)) + + +def _num_degenerate_faces(mesh): + F = np.asarray(mesh.faces) + return sum(1 for t in F if len({int(t[0]), int(t[1]), int(t[2])}) < 3) + + +def _assert_clean_manifold(mesh, global_point_count): + assert _num_degenerate_faces(mesh) == 0, "raw mesh has degenerate triangles" + assert _num_duplicate_faces(mesh) == 0, "raw mesh has duplicate faces" + # every edge borders 1 (boundary) or 2 (interior) faces -- never more + assert _max_edge_face_count(mesh) <= 2, "raw mesh is non-manifold (edge shared by >2 faces)" + # keep_all_vertices=True: faces index the global/unified vertex set + assert len(mesh.vertices) == global_point_count, "raw mesh vertices are not the global set" + + +@pytest.mark.skipif(not WHITNEY, reason="whitney decomposition not present") +def test_raw_mesh_manifold_whitney(): + from bertini_real.surface import Surface + s = Surface(WHITNEY) + gp = len(s.extract_points()) + for p in s.separate_into_nonsingular_pieces(): + _assert_clean_manifold(p.as_mesh(smooth=False), gp) + + +@pytest.mark.skipif(not NORDSTRAND, reason="nordstrand decomposition not present") +def test_raw_mesh_manifold_and_closes_nordstrand(): + from bertini_real.surface import Surface + s = Surface(NORDSTRAND) + gp = len(s.extract_points()) + pieces = s.separate_into_nonsingular_pieces() + for p in pieces: + _assert_clean_manifold(p.as_mesh(smooth=False), gp) + # the headline: every raw piece now caps + joins into a watertight solid + for p in pieces: + closed = p.as_closed_mesh(smooth=False, resolution=4) + assert closed is not None and closed.is_watertight