diff --git a/docs/examples/example_spokewheel.py b/docs/examples/example_spokewheel.py new file mode 100644 index 0000000..8fd66f9 --- /dev/null +++ b/docs/examples/example_spokewheel.py @@ -0,0 +1,317 @@ +from compas.colors import Color +from compas.datastructures import Graph +from compas.geometry import Cylinder +from compas.geometry import Ellipse +from compas.geometry import Line +from compas.geometry import NurbsCurve +from compas.geometry import Vector +from compas.itertools import pairwise +from compas.tolerance import TOL +from compas_dr.constraints import Constraint +from compas_dr.numdata import InputData +from compas_dr.solvers import dr_constrained_numpy +from compas_viewer import Viewer +from compas_viewer.config import Config + +OUTER = 0 +INNER = 40 +SPOKE = 1.5 +STRUT = -2 +LOAD = -1 + +SLIDE = True + +# ============================================================================= +# Geometry +# ============================================================================= + +outer: NurbsCurve = NurbsCurve.from_ellipse(Ellipse(major=20, minor=0.75 * 20)) +inner: NurbsCurve = NurbsCurve.from_ellipse(Ellipse(major=12, minor=0.75 * 12)) + +outer_ring = outer.to_polyline(45) +top_ring = inner.translated([0, 0, 1]).to_polyline(45) +bottom_ring = inner.translated([0, 0, -1]).to_polyline(45) + +# ============================================================================= +# Graph +# ============================================================================= + +graph = Graph() + +graph.update_default_edge_attributes( + outer_ring=False, + top_ring=False, + bottom_ring=False, + top_spoke=False, + bottom_spoke=False, + strut=False, + force=0, + qpre=0, + fpre=0, + lpre=0, +) +graph.update_default_node_attributes( + outer_ring=False, + top_ring=False, + bottom_ring=False, + panel=None, + anchor=False, + residual=None, + load=None, + constraint=None, +) + +# ============================================================================= +# Graph components +# ============================================================================= + +gkey_node = {} + +for i, point in enumerate(outer_ring[:-1]): + gkey = TOL.geometric_key(point) + if gkey not in gkey_node: + node = graph.add_node(x=point[0], y=point[1], z=point[2], outer_ring=True, panel=i) + gkey_node[gkey] = node + +for i, point in enumerate(bottom_ring[:-1]): + gkey = TOL.geometric_key(point) + if gkey not in gkey_node: + node = graph.add_node(x=point[0], y=point[1], z=point[2], bottom_ring=True, panel=i) + gkey_node[gkey] = node + +for i, point in enumerate(top_ring[:-1]): + gkey = TOL.geometric_key(point) + if gkey not in gkey_node: + node = graph.add_node(x=point[0], y=point[1], z=point[2], top_ring=True, panel=i) + gkey_node[gkey] = node + +for (a, b), (aa, bb) in zip(pairwise(outer_ring), pairwise(bottom_ring)): + u = gkey_node[TOL.geometric_key(a)] + v = gkey_node[TOL.geometric_key(b)] + uu = gkey_node[TOL.geometric_key(aa)] + vv = gkey_node[TOL.geometric_key(bb)] + + graph.add_edge(u, v, outer_ring=True) + graph.add_edge(uu, vv, bottom_ring=True) + graph.add_edge(u, uu, bottom_spoke=True) + +for (a, b), (aa, bb) in zip(pairwise(outer_ring), pairwise(top_ring)): + u = gkey_node[TOL.geometric_key(a)] + v = gkey_node[TOL.geometric_key(b)] + uu = gkey_node[TOL.geometric_key(aa)] + vv = gkey_node[TOL.geometric_key(bb)] + + graph.add_edge(uu, vv, top_ring=True) + graph.add_edge(u, uu, top_spoke=True) + +for a, b in zip(bottom_ring[:-1], top_ring[:-1]): + u = gkey_node[TOL.geometric_key(a)] + v = gkey_node[TOL.geometric_key(b)] + + graph.add_edge(u, v, strut=True) + +# ============================================================================= +# Constraints +# ============================================================================= + +constraint = Constraint(outer) + +graph.nodes_attribute( + name="anchor", + value=True, + keys=graph.nodes_where(outer_ring=True), +) + +graph.nodes_attribute( + name="constraint", + value=constraint, + keys=graph.nodes_where(outer_ring=True), +) + +# for node in [0, 11, 21, 22, 32]: +# graph.unset_node_attribute(node, "constraint") + +for node in [11, 33]: + graph.unset_node_attribute(node, "constraint") + +# ============================================================================= +# Force densities +# ============================================================================= + +graph.edges_attribute(name="qpre", value=OUTER, keys=graph.edges_where(outer_ring=True)) + +graph.edges_attribute(name="qpre", value=INNER, keys=graph.edges_where(bottom_ring=True)) +graph.edges_attribute(name="qpre", value=INNER, keys=graph.edges_where(top_ring=True)) + +graph.edges_attribute(name="qpre", value=SPOKE, keys=graph.edges_where(bottom_spoke=True)) +graph.edges_attribute(name="qpre", value=SPOKE, keys=graph.edges_where(top_spoke=True)) + +graph.edges_attribute(name="fpre", value=STRUT, keys=graph.edges_where(strut=True)) + +graph.nodes_attribute(name="load", value=Vector(0, 0, LOAD), keys=graph.nodes_where(bottom_ring=True)) + +# ============================================================================= +# FormFinding +# ============================================================================= + +points = graph.nodes_attributes(names="xyz") +edges = list(graph.edges()) +fixed = list(graph.nodes_where(anchor=True)) +qpre = list(graph.edges_attribute(name="qpre")) +fpre = list(graph.edges_attribute(name="fpre")) +loads = [graph.node_attribute(node, name="load") or Vector(0, 0, 0) for node in graph.nodes()] + +if not SLIDE: + constraints = [None] * graph.number_of_nodes() +else: + constraints = graph.nodes_attribute(name="constraint") + +indata = InputData( + vertices=points, + edges=edges, + fixed=fixed, + loads=loads, + qpre=qpre, + fpre=fpre, +) + +result = dr_constrained_numpy( + indata=indata, + constraints=constraints, + tol1=1e-6, + tol2=1e-9, + rk_steps=4, + kmax=100, + callback=lambda k, x, crit1, crit2, callback_args: print(k), +) + +for node in graph.nodes(): + graph.node_attributes(node, "xyz", result.xyz[node]) + graph.node_attribute(node, "residual", Vector(*result.residuals[node])) + +for index, edge in enumerate(graph.edges()): + graph.edge_attribute(edge, "force", result.forces[index, 0]) + +# ============================================================================= +# Pre-process visualisation +# ============================================================================= + +scale_tension = 0.002 +scale_compression = 0.02 + +compression = [] +tension = [] + +for edge in graph.edges(): + force = graph.edge_attribute(edge, name="force") + + if force > 0: + radius = +force * scale_tension + tension.append(Cylinder.from_line_and_radius(graph.edge_line(edge), radius)) + + if force < 0: + radius = -force * scale_compression + compression.append(Cylinder.from_line_and_radius(graph.edge_line(edge), radius)) + +tol = 0.01 + +reactions = [] +residuals = [] +tangents = [] + +for node in graph.nodes(): + if graph.node_attribute(node, name="anchor"): + location = graph.node_point(node) + residual = graph.node_attribute(node, name="residual") + constraint = graph.node_attribute(node, name="constraint") + if residual: + reactions.append(Line.from_point_and_vector(location, residual * -0.3)) + if constraint: + constraint.location = location + constraint.residual = residual + constraint.compute_tangent() + tangents.append(Line.from_point_and_vector(location, constraint.tangent)) + + elif graph.node_attribute(node, name="constraint"): + residual: Vector = graph.node_attribute(node, name="residual") + if residual and residual.length > tol: + residuals.append(Line.from_point_and_vector(graph.node_point(node), residual)) + +loads = [] + +for node in graph.nodes(): + load: Vector = graph.node_attribute(node, name="load") + if load and load.length: + loads.append(Line.from_point_and_vector(graph.node_point(node), load)) + +supports = [graph.node_point(node) for node in graph.nodes_where(anchor=True)] + +# ============================================================================= +# Visualisation +# ============================================================================= + +loadcolor = Color.green().darkened(50) +reactioncolor = Color.green() +residualcolor = Color.cyan() +supportcolor = Color(0.2, 0.2, 0.2) +tensioncolor = Color.red() +compressioncolor = Color.blue() + +config = Config() +config.camera.target = [0, 0, 0] +config.camera.position = [25, -40, 30] + +viewer = Viewer(config=config) + +viewer.scene.add(graph, show_points=True) + +viewer.scene.add( + supports[:10], + name="Supports", + pointcolor=supportcolor, + pointsize=10, +) +viewer.scene.add( + tension, + name="Tension", + facecolor=tensioncolor, + linecolor=tensioncolor.contrast, +) +viewer.scene.add( + compression, + name="Compression", + facecolor=Color.blue(), + show_lines=False, +) +viewer.scene.add( + reactions, + name="Reactions", + linecolor=reactioncolor, + linewidth=3, +) +viewer.scene.add( + residuals, + name="Residuals", + linecolor=residualcolor, + linewidth=3, +) +viewer.scene.add( + tangents, + name="Tangents", + linecolor=Color.magenta(), + linewidth=3, +) +viewer.scene.add( + loads, + name="Loads", + linecolor=loadcolor, + linewidth=3, +) +viewer.scene.add( + constraint.geometry, + name="Constraints", + linecolor=Color.cyan(), + linewidth=2, +) + +viewer.show()