C2.2 Form Finding

Simple form finding of a mesh with default attributes. (by Lotte)

Objective

We will learn how to form find the cable net with the dynamic relaxation method using the compas function dr. For doing so we must create the data structure with the required attributes and compile it to the numerical data. Afterwards, we must update our datastructure and visualise it.

Procedure

o. Paths and Imports

We import the comps mesh datastructure and the dynamic relaxation compas function dr.

from compas.datastructures import Mesh
from compas.numerical import dr

We will build upon the mesh data structure that we imported and saved in the JSON file in the data folder in C2.1 Geometry Import section.

# ==============================================================================
# Paths
# ==============================================================================

HERE = os.path.dirname(__file__)
DATA = os.path.abspath(os.path.join(HERE, '../..', 'data'))
FILE_I = os.path.join(DATA, 'cablenmesh_import.json')

a. Datastructure

We want to create the COMPAS mesh data structure, then we want to extend it with attributes and methods related to form finding and analysis of flexible cable net formwork for concrete shells, and we want to set its boundary conditions.

a1. Create the Mesh

We start by creating a COMPAS mesh from the JSON file that we generated by the import from Rhino.

# ==============================================================================
# Cablenet mesh datastructure
# ==============================================================================

# a1. create the mesh from imported geometry
mesh = Mesh.from_json(FILE_I)

We use a mesh instead of a network because it provides functionalities for the concrete weight and fabric shuttering etc.

a2. Update Default Vertex and Edge Attributes

The necessary default attributes for the vertices are (similar to nodes in section C0.):

'rx': 0.0, 'ry': 0.0, 'rz': 0.0 – component of a residual force.
'px': 0.0, 'py': 0.0, 'pz': 0.0 – component of an externally applied load.
'is_anchor': False, – to indicate that a vertex is anchored and can take reaction forces in XYZ.

The necessary attributes for the edges are:

'q': 1.0 – force densities of an edge.
'f': 0.0 – force in an edge.
'l': 0.0 – stressed length of an edge.
'l0': 0.0 –unstressed length of an edge.

The meaning of the unstressed length l0 will be explained later in section D1. Cables.

# a2. set default vertex and edge attributes
dva = {
    'rx': 0.0,            # X-component of an residual force.
    'ry': 0.0,            # Y-component of an residual force.
    'rz': 0.0,            # Z-component of an residual force.
    'px': 0.0,            # X-component of an externally applied load.
    'py': 0.0,            # Y-component of an externally applied load.
    'pz': 0.0,            # Z-component of an externally applied load.
    'is_anchor': False,   # Indicate that a vertex is anchored and can take reaction forces in XYZ.
}
mesh.update_default_vertex_attributes(dva)

dea = {
    'q': 1.0,             # Force densities of an edge.
    'f': 0.0,             # Force in an edge.
    'l': 0.0,             # Stressed Length of an edge.
    'l0': 0.0,            # Unstressed Length of an edge.
}
mesh.update_default_edge_attributes(dea)

a3. Set Boundary Conditions

We identify the anchors on the boundary and update the 'is_anchor' attribute of these vertices:

# a3. Boundary conditions
boundary = mesh.vertices_on_boundaries()[0]
mesh.vertices_attribute('is_anchor', True, keys=boundary)

We visualise the cable net with anchored vertices in red:

# ==============================================================================
# Visualize
# ==============================================================================

baselayer = "DF21_C2::02_Basic Form Finding"

artist = MeshArtist(mesh, layer=baselayer+"::InputMesh")
artist.clear_layer()

artist.draw_vertices(color={vertex: (255, 0, 0) for vertex in mesh.vertices_where({'is_anchor': True})})
artist.draw_edges()
artist.draw_faces()

b. Form Finding

Now we want to form find the cablemesh into its equilibrium state under pure prestress. For this, we will use the formerly introduced COMPAS function dr. As shown in C0. Principles, the function does not 'understand' the datastructure, it requires the input data in lists. Thus we provide it with lists, it form finds the equilibrium numerically and it outputs lists, but it has no idea of the mesh itself, so we must always update the mesh attributes ourselves afterwards.

to see what this function does exactly check out section C1.4 Under the Hood.

b1. Compile to Numeric Data

Check here to see what the function requires as input parameters and returns. Note that only the ones not marked with 'optional' are required.

Again, the numerical function inputs are lists that are ordered in the index and not the keys (vertex ids), so we must convert for the input the keys to indices and for the returns the indices to keys!

For this, we use the Mesh.key_index() method that returns a dictionary that maps vertex dictionary keys to the corresponding index in a vertex list or array.

# ==============================================================================
# b1. Compile numerical data as fofin input > NEW
# ==============================================================================

# to map vertex dictionary keys to the corresponding index in a vertex list
vertex_index = mesh.key_index()

# complile lists (ordered in indices) of geometry, loads and force densities
X = mesh.vertices_attributes('xyz')
P = mesh.vertices_attributes(['px', 'py', 'pz'])
Q = mesh.edges_attribute('q')

# translate anchored vertex keys to index list of fixed vertices
fixed = [vertex_index[vertex] for vertex in mesh.vertices_where({'is_anchor': True})]  # noqa: E501
# translate edges from vertex keys to indexes
edges = [(vertex_index[u], vertex_index[v]) for u, v in mesh.edges()]

b2. Compute Equilibrium

Now we hand this over to the dr function:

# ==============================================================================
# b2. Compute equilibrium
# ==============================================================================

X, Q, F, L, R = dr(X, edges, fixed, P, Q)

If you were to run it like this with the visualisation you would see that the geometry actually hasn't changed at all yet.

b3. Update the Cablemesh Datastructure

After the numerical calculation, we must update the mesh datastructure with its new equilibrium xyz-coordinates and the remaining residual forces by setting its vertex attributes and the new force density, force and length by setting its edge attributes.

# ==============================================================================
# b3. Update cablemesh
# ==============================================================================

for vertex in mesh.vertices():
    # translate back from indices to keys
    index = vertex_index[vertex]
    mesh.vertex_attributes(vertex, 'xyz', X[index])
    mesh.vertex_attributes(vertex, ['rx', 'ry', 'rz'], R[index])

for index, edge in enumerate(mesh.edges()):
    # translate back from indices to keys
    mesh.edge_attribute(edge, 'q', Q[index])
    mesh.edge_attribute(edge, 'f', F[index])
    mesh.edge_attribute(edge, 'l', L[index])

The simple visualisation remains the same as in a. Datastructures. but changing the output layer to layer=baselayer+"::Mesh

c. Into fofin function

Since we will be using the three steps from b. Form Finding over and over again, it is convenient to wrap them into a function that we can call later.

# ==============================================================================
# Helpers
# ==============================================================================

def fofin(mesh):
    """Compute equilibrium and update the geometry.
    """
    # ==============================================================================
    # Compile numerical data as fofin input
    # ==============================================================================

    # to map vertex dictionary keys to the corresponding index in a vertex list
    vertex_index = mesh.key_index()

    # complile lists (ordered in indices) of geometry, loads, force densities
    X = mesh.vertices_attributes('xyz')
    P = mesh.vertices_attributes(['px', 'py', 'pz'])
    Q = mesh.edges_attribute('q')

    # translate anchored vertex keys to index list of fixed vertices
    fixed = [vertex_index[vertex] for vertex in mesh.vertices_where({'is_anchor': True})]  # noqa: E501
    # translate edges from vertex keys to indexes
    edges = [(vertex_index[u], vertex_index[v]) for u, v in mesh.edges()]

    # ==============================================================================
    # Compute equilibrium
    # ==============================================================================

    X, Q, F, L, R = dr(X, edges, fixed, P, Q)

    # ==============================================================================
    # Update cablenet mesh
    # ==============================================================================

    for vertex in mesh.vertices():
        # translate back from indices to keys
        index = vertex_index[vertex]
        mesh.vertex_attributes(vertex, 'xyz', X[index])
        mesh.vertex_attributes(vertex, ['rx', 'ry', 'rz'], R[index])

    for index, edge in enumerate(mesh.edges()):
        # translate back from indices to keys
        mesh.edge_attribute(edge, 'q', Q[index])
        mesh.edge_attribute(edge, 'f', F[index])
        mesh.edge_attribute(edge, 'l', L[index])

# ==============================================================================
# c. Compute equilibrium and update the geometry > NEW
# ==============================================================================

fofin(mesh)

The cablemesh is very flat and thus has very low geometric stiffness.

d. Visualization

If we want to visualise the reaction forces, residual forces and edge forces to better understand what's going on, we can do so as follows:

Import the functions from the C2_3_visualisation script:

from C2_3_visualisation import draw_reactions
from C2_3_visualisation import draw_residuals
from C2_3_visualisation import draw_forces
from C2_3_visualisation import draw_loads

And call them after the MeshArtist:

# ==============================================================================
# Visualize
# ==============================================================================

baselayer = "DF21_C2::02_Basic Form Finding"

artist = MeshArtist(mesh, layer=baselayer+"::Mesh")
artist.clear_layer()

artist.draw_vertices(color={vertex: (255, 0, 0) for vertex in mesh.vertices_where({'is_anchor': True})})  # noqa: E501
artist.draw_edges()
artist.draw_faces()

draw_reactions(mesh, baselayer=baselayer)
draw_residuals(mesh, baselayer=baselayer)
draw_forces(mesh, baselayer=baselayer)
draw_loads(mesh, baselayer=baselayer)

You will not be able to see any external loads or residuals yet, because these are 0.

So now, let's find out how these are implemented in C2.3 Visualisation.

PreviousC2.1 Geometry Import NextC2.3 Visualisation

Last updated 2 years ago