cadquery

cadquery copied to clipboard

Attached code

import cadquery as cq
import numpy as np

# Define a cubic lattice
n = 4 # Dimension of the lattice (n x n x n)
spacing = 1  # Spacing between points

# Generate lattice points
points = np.array([(x, y, z) for x in range(n) for y in range(n) for z in range(n)]) * spacing

# User-defined radius for the spheres and cylinders
sphere_radius = float(input("Enter sphere radius: "))
cylinder_radius = float(input("Enter cylinder radius: "))

# Initialize a new CadQuery model for the lattice
model = cq.Workplane("XY")

# Generate spheres at each lattice point
for point in points:
    # Convert each point from numpy array to a tuple
    model = model.union(cq.Workplane("XY").sphere(sphere_radius).translate(tuple(point)))

# Function to compute connections for a cubic lattice
def cubic_connections(points, n, spacing):
    connections = []
    for x in range(n):
        for y in range(n):
            for z in range(n):
                idx = x * n * n + y * n + z
                # Connect to the six neighboring points within the bounds
                if x < n - 1: connections.append((idx, idx + n * n))  # +x direction
                if y < n - 1: connections.append((idx, idx + n))      # +y direction
                if z < n - 1: connections.append((idx, idx + 1))      # +z direction
    return connections

# Compute connections
connections = cubic_connections(points, n, spacing)

# Add cylinders based on the connections
for start, end in connections:
    start_point = points[start]
    end_point = points[end]
    
    # Calculate the vector between start and end points
    vec = end_point - start_point
    height = np.linalg.norm(vec)  # The height of the cylinder is the distance between points

    if height == 0:  # Check if the connection is valid
        continue  # Skip creating a cylinder for zero-length connections

    # Calculate the axis and angle to rotate the cylinder between points
    z_axis = np.array([0, 0, 1])
    axis = np.cross(z_axis, vec)
    angle = np.degrees(np.arccos(np.clip(np.dot(z_axis, vec) / height, -1.0, 1.0)))

    if np.linalg.norm(axis) == 0:  # Handle case when the axis vector is zero
        # If the start and end points are perfectly aligned with z-axis, just align the cylinder along the z-axis
        axis = z_axis
        angle = 0  # No rotation needed

    # Create a cylinder and align it between the two points
    cylinder = (cq.Workplane("XY")
                .cylinder(height, cylinder_radius)
                .rotate((0, 0, 0), tuple(axis), angle)  # Rotate the cylinder along the calculated axis
                .translate(tuple((start_point + end_point) / 2)))  # Place at the midpoint of the connection
    
    # Union the cylinder with the existing model
    model = model.union(cylinder)

# Create a bounding box from the points without any expansion
min_point = np.min(points, axis=0)  # Minimum coordinates
max_point = np.max(points, axis=0)  # Maximum coordinates

# Create a box using the min and max points
bounding_box = (
    cq.Workplane("XZ")
    .box(max_point[0] - min_point[0], max_point[1] - min_point[1], max_point[2] - min_point[2])
    .translate((min_point[0] + (max_point[0] - min_point[0]) / 2,
                 min_point[1] + (max_point[1] - min_point[1]) / 2,
                 min_point[2] + (max_point[2] - min_point[2]) / 2))
)

# Perform boolean difference: subtract lattice from the bounding box
final_model = bounding_box.cut(model)

# Export the resulting model to STL format
cq.exporters.export(final_model, "tortous.stl")

# Optional: display the model if running in Jupyter Notebook or CQ-editor
try:
    show_object(final_model)
except NameError:
    print("Run in CQ-editor or Jupyter Notebook with CadQuery integration to view the model.")