D. Mesh, Topology and Pattern

Form finding a cable-net or generating a compression-only vault usually starts with a 2D Mesh, which represents the horizontal projection of the eventual structure. The line segments of this Mesh is called Pattern. For example, a pattern could be created from an existing structure.

Nave Vault of the Sherborne Abbey, Dorset, UK (photo by Lawrence Lew)

1. Interactive Mesh Constructor in Rhino

Pattern generation usually starts with a base Mesh, and then optimize or modify the mesh topologically later. Here are a few ways to create a base mesh with objects in Rhino.

1.a: From Mesh

A Rhino mesh object can be converted to a COMPAS mesh object.

guid = compas_rhino.select_mesh(message='Select Mesh')
compas_rhino.rs.HideObjects(guid)
mesh = RhinoMesh.from_guid(guid).to_compas(cls=Mesh)

1.b: From Surface

A non-trimmed Rhino NURBS surface or a Rhino polysurface can be used to create a Mesh. Note that the control point numbers in u-, v-directions of the surface won't influence the mesh.

guid = compas_rhino.select_surface(message='select an untrimmed surface or a polysurface')
compas_rhino.rs.HideObject(guid)
mesh = RhinoSurface.from_guid(guid).to_compas(cls=Mesh)

1.c: From Lines

The most simple and manual way is to draw the edges of the Mesh as Rhino lines. Each edge should be an individual line in the XY plane; all lines should be broken at all line intersections. In other words, these lines may not be overlapping.

A mesh data structure is the network of faces. Thus, if the input set of lines must consist of closed loops of lines representing the faces. If the lines are not closed, such as a "leaf" edge, which contains a vertex with degree 1 will be omitted.

guids = compas_rhino.select_lines(message='Select Lines')
compas_rhino.rs.HideObjects(guids)
lines = compas_rhino.get_line_coordinates(guids)
basegrid = Mesh.from_lines(lines, delete_boundary_face=True)

input lines

mesh from lines

1.d: From Points

For strictly two-dimensional points in the XY plane, they can also be used to create a base mesh.

guids = compas_rhino.select_points(message='Select Points')
compas_rhino.rs.HideObjects(guids)
points = compas_rhino.get_point_coordinates(guids)
basegrid = Mesh.from_points(points, delete_boundary_face=True)

input points

mesh from points

Mesh.from_lines() and Mesh.from_points() creates two different types of two-dimensional face shapes that are commonly used. The former one is quadrilateral, quad mesh and the latter one is triangle. Quadrilateral and triangle faces have 4 and 3 vertices on every face. A quadrilateral mesh could be turned into a triangle mesh by splitting each set of four into two sets of three. A quadrilateral face could be non-planar. By moving the vertices along z-axis in the earlier part of this session has led to 4 non-planar quad faces.

Mesh.quads_to_triangles() could convert a quad mesh to a triangle mesh.

triangle and quad mesh

2. Ex. Delaunay triangulation

Triangulation is a method that is fast and robust for create a base mesh. Here a set of random vertices in the boundary is used to compute a Delaunay triangulation Delaunay triangulations maximize the minimum angle of all the angles of the triangles in the triangulation; they tend to avoid sliver triangles. We can use function compas.geometry.delaunay_from_points.

import random as r

from compas.geometry import delaunay_from_points
from compas.datastructures import Mesh
from compas_rhino.artists import MeshArtist

# ==============================================================================
# Construct a mesh from random points and triangulation
# ==============================================================================
width = 100
height = 100

vertices = []

# create 100 random points
for i in range(100):
    x = r.randint(0, width)
    y = r.randint(0, height)

    vertices.append([x, y, 0])

# points on the boundary
boundary = [[0, 0, 0], [0, 100, 0], [100, 100, 0], [100, 0, 0]]
vertices.extend(boundary)

faces = delaunay_from_points(vertices, boundary=boundary)

mesh = Mesh().from_vertices_and_faces(vertices, faces)

# ==============================================================================
# Visualization
# ==============================================================================
artist = MeshArtist(mesh, layer="CSD2::Triangulation")
artist.clear_layer()
artist.draw_faces(join_faces=True)
artist.draw_edges()
artist.redraw()

3. Remote Procedure Call and Proxy

Python has a default implementation in C: CPython. However, Rhino uses another implementation of the language, IronPython, based on .NET. As such, it implements the default core language and most of the standard library. Some Python packages available through PIP are pure Python and rely on elements supported by IronPython. Some other modules are, however, specific to the implementation. In particular, NumPy uses native C code, which is supported by CPython but not by IronPython.

To work around this limitation, COMPAS provides a mechanism to access the functionality of a CPython environment seamlessly from any other Python environment through a Remote Procedure Call(RPC).

Continuing the example of Delaunay triangulation in the last session, there are two functions in COMPAS, delaunay_from_points() and delaunay_from_points_numpy(), the former could be called in Rhino with pure IronPython, whereas the latter needs a Proxy.

from compas.rpc import Proxy
proxy = Proxy()
proxy.package = 'compas.geometry'
numpy_faces = proxy.delaunay_from_points_numpy(vertices)
mesh = Mesh().from_vertices_and_faces(vertices, numpy_faces)

import time
start = time.time()
faces = delaunay_from_points(vertices, boundary=boundary)
end = time.time()
print('computing time:', end - start)

Comparing the computing time, delaunay_from_points() takes 0.0419921875 s and delaunay_from_points_numpy() takes 0.01100921630859375s, which is around four times faster. The comparative advantage of using numpy and Proxy would expand along with the increase in the input scale.

4. Topological Modification

4.a: Delete Vertices

Mesh.delete_vertex(key) would delete a vertex from the mesh and everything that is attached to it. For example, holes can be created by deleting interior vertices is a base mesh, or a pattern.

while True:
    key = cablenet_mesh.get_any_vertex()
    # the vertex to be deleted is not on the boundary
    if cablenet_mesh.is_vertex_on_boundary(key) is True:
        continue
    else:
        # the faces connected to the vertex are not on the boundary
        for fkey in cablenet_mesh.vertex_faces(key):
            if cablenet_mesh.is_face_on_boundary(fkey) is True:
                break
        else:
            break
        continue
        
cablenet_mesh.delete_vertex(key=key)

4.b: Delete Faces

while True:
    fkey = cablenet_mesh.get_any_face()
    # the vertex to be deleted is not on the boundary
    if cablenet_mesh.is_face_on_boundary(fkey) is False:
        break

cablenet_mesh.delete_face(fkey)

5. Subdivision

Mesh Class contains different subdivision methods to densify the original mesh.

5.a: subdivision tracking attributes

Here we create a function to track the hierarchy of each subdivision using the mesh attributes.

def subd_tracked(mesh, k):
    levels = []
    for _ in range(k):
        mesh.update_default_face_attributes({'children': None})
        mesh.update_default_edge_attributes({'child': None})
        subd = mesh_subdivide(mesh, scheme='quad', k=1)
        subd.update_default_face_attributes({'children': None})
        subd.update_default_edge_attributes({'child': None})
        for fkey in mesh.faces():
            children = []
            for u, v in mesh.face_halfedges(fkey):
                u_nbrs = subd.vertex_neighbors(u)
                v_nbrs = subd.vertex_neighbors(v)
                shared = list(set(u_nbrs) & set(v_nbrs))
                children += shared
                mesh.set_edge_attribute((u, v), 'child', shared[0])
            keys = list(set.intersection(*[set(subd.vertex_neighbors(key)) for key in children]))
            mesh.set_face_attribute(fkey, 'children', children + keys)
        levels.append(subd)
        mesh = subd
    return levels
    
S1, S2 = subd_tracked(cablenet_mesh, k=2)

5.b: visualize vertices hierarchy

Here original nodes are white dots, children are red and grandchildren are black.

S1, S2 = subd_tracked(cablenet_mesh, k=2)
for fkey in cablenet_mesh.faces():
    children = cablenet_mesh.face_attribute(fkey, 'children')
    for child in children:
        vertexcolor[child] = (255, 0, 0)
    m_vkey = children[-1]
    child_fkeys = S1.vertex_faces(m_vkey)
    for child_fkey in child_fkeys:
        grandchildren = S1.face_attribute(child_fkey, 'children')
        for grandchild in grandchildren:
            vertexcolor[grandchild] = (0, 0, 0)

5.c: visualize edge hierarchy

S1, S2 = subd_tracked(cablenet_mesh, k=2)
for u, v in cablenet_mesh.edges():
    child = cablenet_mesh.edge_attribute((u, v), 'child')
    u_grandchild = S1.edge_attribute((u, child), 'child')
    v_grandchild = S1.edge_attribute((v, child), 'child')

    vertexcolor[child] = (255, 0, 0)
    vertexcolor[u_grandchild] = (0, 0, 0)
    vertexcolor[v_grandchild] = (0, 0, 0)

5.d: visualize face hierarchy

The following image shows you the change of face keys. The smallest fkeys are shown in green, and the largest ones are shown in red. In each subdivision iteration, old fkey would be deleted and new fkey would be created. Thus, if you try to show the fkeys after 2 times of subdivisoon, you will not be able to the former fkeys.

Note that adding and removing elements will not cause identifiers to be renumbered. Therefore, after certain topological operations (e.g. subdivision), vertex and face identifiers no longer necessarily form contiguous sequences. This needs to be taken into account when converting sequences of vertices, faces, and edges to lists, for example for numerical calculation. To transparently convert non-contiguous sequences of identifiers to contiguous list indices, use “key/index maps”.

key_index = {key: index for index, key in enumerate(mesh.vertices())}

key_index = mesh.key_index()
vertices = list(mesh.vertices())
edges = [(key_index[u], key_index[v]) for u, v in mesh.edges()]
faces = [[key_index[key] for key in mesh.face_vertices(face)] for face in mesh.faces()]

5.e: mesh subdivision methods

An input mesh can be subdivided under different subdivision algorithms.

mesh_subdivide_quad() could subdivide a mesh such that all faces are quads. The former example is a quad scheme subdivision.

sub_mesh = cablenet_mesh.subdivide(scheme='quad', k=1)

mesh_subdivide_tri() could subdivide a mesh using simple insertion of vertices.

sub_mesh = cablenet_mesh.subdivide(scheme='tri', k=1)

mesh_subdivide_catmullclark() subdivides a mesh using the Catmull-Clark algorithm. Catmull-Clark subdivision is similar to the quad subdivision, but with smoothing after every level of further subdivision.

sub_mesh = cablenet_mesh.subdivide(scheme='catmullclark', k=1)

mesh_subdivide_corner() subdivides a mesh by cutting corners.

sub_mesh = cablenet_mesh.subdivide(scheme='corner', k=1)