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)