Hi @mike_002 and @inno,
Here’s a revision of my Python script from the weekend.
It mostly works now, after figuring out some of the bugs.
The script only produces polylines that then can be lofted with vanilla Grasshopper components, and thus made into three-dimensional geometry.
I’ve included a workflow with meshes and one with polysurfaces in the Grasshopper file. I’d use the meshes for prototyping and finding a desired design and the breps for manufacturing or really sharp renderings, since it’s more precise.
The top faces of the cells or cassettes are all planar, like in the reference.
The domes on their undersides aren’t totally spherical or round. I couldn’t figure that out, because sweeping somehow doesn’t seem to work with a planar and a non-planar rail, as well as planar section curve?
In this example, I’ve left them faceted to illustrate this, but the resolution can be dialled up much higher to make them much smoother, which in turn makes the model quite heavy.
The radius or size of the holes is determined randomly, but there are lots of parameters that you change to vary things.
Here’s the script:
from ghpythonlib import treehelpers
import Rhino.Geometry as rg
import random
import math
random.seed(Seed)
def normal_at(mesh, test_pt, max_distance=0.01):
"""Gets the normal at a mesh point that is close the the test point.
Args:
mesh (Rhino.Geometry.NgonMesh): A mesh with computed face normals
test_pt (Rhino.Geometry.Point3d): A test point to search from
max_distance (float): Optional upper bound on the distance
from the test point to the mesh, by default 0.01
Returns:
The normal at the closest mesh point.
"""
closest_mesh_pt = mesh.ClosestMeshPoint(test_pt, max_distance)
if closest_mesh_pt is None:
return
return mesh.NormalAt(closest_mesh_pt)
def average_vector(vectors):
"""Returns the arithmetic average vector of a list of vectors."""
a = rg.Vector3d.Zero
for v in vectors:
a += v
a /= len(vectors)
return a
def average_point(points):
"""Returns the arithmetic average point of a list of points."""
a = rg.Point3d.Origin
for pt in points:
a += pt
a /= len(points)
return a
def create_polyline(vertices, closed=True):
"""Returns a polyline connecting a number of vertices."""
pline = rg.Polyline()
for v in vertices:
pline.Add(v)
if closed:
pline.Add(vertices[0])
return pline
def tween_between(a, b, count):
"""Returns an x-amount of tween curves between the polylines a and b.
Both polylines must have the same number of vertices."""
div_pts = []
if count < 1:
return []
for j in xrange(len(a)):
line = rg.LineCurve(a[j], b[j])
params = line.DivideByCount(count + 1, False)
div_pts.append([line.PointAt(t) for t in params])
tween_crvs = []
for j in xrange(count):
pline = rg.Polyline()
for pts in div_pts:
pline.Add(pts[j])
tween_crvs.append(rg.PolylineCurve(pline))
return tween_crvs
a = []
if __name__ == "__main__":
NPB = [] # non-planar boundaries
BDS = [] # boundary-defining side subdivisions
OPB = [] # offset, planar boundaries
OPS = [] # planar, offset face-defining subdivisions
CHB = [] # circular, planar hole boundaries
DDA = [] # dome-defining arcs
DDR = [] # dome-defining rings
if NgonMesh is not None:
NgonMesh.FaceNormals.ComputeFaceNormals()
section_arcs = []
section_rings = []
offset_subdivs = []
boundary_subdivs = []
for i in xrange(NgonMesh.Ngons.Count): # loop each ngon face
vertex_indices = NgonMesh.Ngons.GetNgon(i).BoundaryVertexIndexList()
vertices = [
rg.Point3d(NgonMesh.Vertices[j]) for j in vertex_indices
]
boundary_pline = create_polyline(vertices) # face boundary
boundary_pline.MergeColinearSegments(AngleTol, True)
# Divide the initial, non-planar boundary
div_pts = [] # face boundary division points
ia = NumArcs + 1
for j in xrange(boundary_pline.SegmentCount):
segment = boundary_pline.SegmentAt(j)
div_pts.append(segment.From)
crv = rg.LineCurve(segment)
params = crv.DivideByCount(ia, False)
for t in params:
div_pts.append(crv.PointAt(t))
div_pts_normals = [
normal_at(NgonMesh, p, float("inf")) for p in div_pts
]
ngon_pline = create_polyline(div_pts)
NPB.append(ngon_pline)
ngon_center = NgonMesh.ClosestPoint(average_point(div_pts))
ngon_normal = average_vector(div_pts_normals)
ngon_normal.Unitize()
if Flip:
ngon_normal.Reverse()
# Define the offset, planar boundary
offset_center = ngon_center + ngon_normal * OutHeight
plane = rg.Plane(offset_center, ngon_normal)
ll = 1.0 if OutHeight == 0.0 else OutHeight + OutHeight * 0.5
offset_pts = []
for j in xrange(len(div_pts)):
start = div_pts[j]
end = start + (div_pts_normals[j] * ll)
ray = rg.Line(start, end)
_, t = rg.Intersect.Intersection.LinePlane(ray, plane)
offset_pts.append(ray.PointAt(t))
offset_pline = create_polyline(offset_pts) # is planar
OPB.append(offset_pline)
offset_center = average_point(offset_pts)
# Evaluate the random center hole radius
max_ngon_radius = float("inf")
for j in xrange(offset_pline.Count - 1):
dist2_vertex = offset_center.DistanceToSquared(offset_pline[j])
midpt = (offset_pline[j] + offset_pline[j + 1]) * 0.5
dist2_midpt = offset_center.DistanceToSquared(midpt)
min_dist = math.sqrt(min(dist2_vertex, dist2_midpt))
if min_dist < max_ngon_radius:
max_ngon_radius = min_dist
max_radius = max_ngon_radius - max_ngon_radius * MinOffset
min_radius = min(MinRadius, max_radius)
radius = random.uniform(min_radius, max_radius)
rest_radius = max_ngon_radius - radius
# Construct the backset circle which can define an inner thickness
backset_height = max(0.0, min(InHeight, OutHeight))
backset_center = offset_center - (plane.ZAxis * backset_height)
backset_plane = rg.Plane(backset_center, plane.ZAxis)
backset_circle = rg.Circle(backset_plane, backset_center, radius)
# Construct the dome-defining arcs
arcs = []
for j in xrange(len(offset_pts)):
closest_pt = backset_circle.ClosestPoint(offset_pts[j])
end_pts = [closest_pt, div_pts[j]] # arc start and end points
end_indices = [0, 1] # indices of the end points
lengths = [p.DistanceToSquared(offset_pts[j]) for p in end_pts]
sorted_indices = [
p for _, p in sorted(zip(lengths, end_indices))
]
fi = sorted_indices[-1]
li = sorted_indices[0]
tangent = offset_pts[j] - end_pts[fi]
arc = rg.Arc(end_pts[fi], tangent, end_pts[li])
if fi == end_indices[-1]:
arc.Reverse()
arcs.append(arc)
section_arcs.append(arcs)
# Divide the arcs
r = NumRings + 1 # number of rings
arc_pts = []
min_arc_div_dist2 = float("inf")
for arc in arcs:
crv = rg.ArcCurve(arc)
params = crv.DivideByCount(r, True)
div_pts = [crv.PointAt(t) for t in params]
dist2 = div_pts[0].DistanceToSquared(div_pts[1])
if dist2 < min_arc_div_dist2:
min_arc_div_dist2 = dist2
arc_pts.append(div_pts)
# Construct the dome-defining rings
rings = []
for j in xrange(r):
ring_pts = []
for points in arc_pts:
ring_pts.append(points[j])
ring_pline = create_polyline(ring_pts)
rings.append(ring_pline)
# Construct the offset center hole boundary
backset_ring = rings[0]
offset_ring = rg.Polyline()
for vtx in backset_ring:
offset_vtx = vtx + backset_plane.ZAxis * backset_height
offset_ring.Add(offset_vtx)
CHB.append(offset_ring)
if backset_height == 0.0: # backset_ring == rings[0]
rings.pop(0) # avoid duplicate polylines
section_rings.append([r.ToPolylineCurve() for r in rings])
# Subdivide the planar, offset face in relation to the arc divisions
num = math.floor((rest_radius) / math.sqrt(min_arc_div_dist2)) - 1
face_subdivs = tween_between(offset_pline, offset_ring, num)
offset_subdivs.append(face_subdivs)
# Subdivide the side faces in relation to the arc divisions
sides_subdivs = tween_between(ngon_pline, offset_pline, SideDivs)
boundary_subdivs.append(sides_subdivs)
BDS = treehelpers.list_to_tree(boundary_subdivs)
OPS = treehelpers.list_to_tree(offset_subdivs)
DDA = treehelpers.list_to_tree(section_arcs)
DDR = treehelpers.list_to_tree(section_rings)
voute_de_lefrevre.gh (88.6 KB)