B1. Selfweight

Selfweight of the Concrete Shell

Objectives

You will learn how to form find the cable net under self-weight that adapts to each form-finding iteration:

In module A. we form found the cable net for its pure prestress state without external loading. Now, if the concrete weight is added as an external load 'pz', the geometry will find a new configuration in equilibrium for this loading state. However, in the new geometric configuration, the loads from the concrete's self-weight have changed with the tributary areas of the vertices. Thus the cable-net must be again form found for the new loading stage. This must be repeated until the difference in loads/residuals is small enough to consider convergence.

Procedure

• a. update the default attributes with the concrete thickness and density

• b. compute the self-weight per vertex and add as a load attribute

• c. find the new form under self-weight of the previous form

• d. update the residuals from the difference of self-weight

• e. repeat recomputing the equilibrium until convergence

o. Paths

We will build upon the mesh data structure from the form-finding section. We will output the resulting mesh in a new JSON file.

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

HERE = os.path.dirname(__file__)
DATA = os.path.abspath(os.path.join(HERE, '..', 'data'))
FILE_I = os.path.join(DATA, 'cablenet_fofin.json')
FILE_O = os.path.join(DATA, 'cablenet_fofin_sw.json')

a. Update Default Attributes

Add to the cable net mesh data structure the thickness (0.05m) of the concrete shell as default vertex attribute and the concrete weight density (14 kN/m3) as general mesh attribute.

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

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

# set new default vertex attributes
dva = {
    't': 0.05,          # Thickness of the concrete shell.
}
mesh.update_default_vertex_attributes(dva)

# set mesh attributes
mesh.attributes['density'] = 14.0  # Density of the lightweight concrete.

b. Compute the Selfweight

We will compute the self-weight per vertex based on its tributary area, shell thickness and concrete density, and add this load as the vertex attribute 'pz'.

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

def selfweight(mesh):

    density = mesh.attributes['density']

    # for all vertices that are free
    for key, attr in mesh.vertices_where({'is_anchor': False}, True):

        # get the thickness of the vertex
        thickness = attr['t']

        # compute the tributary area of the vertex
        area = mesh.vertex_area(key)

        # compute the concrete weight
        selfweight = area * thickness * density

        # apply selfweight by updating vertex attribute (acting in negative z)
        attr['pz'] = -selfweight

We will visualise the loads per vertex.

To access the visualisation functions from A3. Visualisation, the file must be in the same folder. So just copy it over.

# ==============================================================================
# Fofin for selfweight
# ==============================================================================

# compute selfweight for the current geometry
selfweight(mesh)

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

baselayer = "CSD2::B1::FormFindingSelfweight"

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})})
artist.draw_edges()
artist.draw_faces()

draw_loads(mesh, baselayer=baselayer, scale=10)

c. Form Finding under Selfweight

First, wrap the form-finding procedure into a function. Then, find the new form under self-weight.

The form-finding procedure with the compilation of numerical data and the geometric update of the data structure must be repeated for every iteration, thus it is more efficient and clean to put it all into a function that can be called in every iteration.

It's exactly the same code as we used in A2. Form-Finding, so just copy paste it :)

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

...

def fofin(mesh):
    # ==============================================================================
    # 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 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})]
    # 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 = fd(X, edges, fixed, Q, P)

    # ==============================================================================
    # 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][0])
        mesh.edge_attribute(edge, 'f', F[index][0])
        mesh.edge_attribute(edge, 'l', L[index][0])

Also, visualise the reaction forces and residuals (there are none calculated yet).

# ==============================================================================
# Fofin for selfweight
# ==============================================================================

# compute selfweight for the current geometry
selfweight(mesh)

# form finding with selfweight loads and geometry update
fofin(mesh)

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

baselayer = "CSD2::B1::FormFindingSelfweight"

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})})
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, scale=10)

d. Update Residuals with Difference of Selfweights

The residuals that are retrieved after the fofin update do not represent the change in self-weights between the geometry before and after the form-finding step (the tributary area changes in size). Thus, the difference in self-weight must be added to the residuals after the self-weight is recalculated for the updated geometry.

use Mesh.key_index to translate the list sorted by indices to the respective key!

def update_residual(mesh, loads_previous):

    key_index = mesh.key_index()

    R = []
    for key in mesh.vertices_where({'is_anchor': False}):

        index = key_index[key]

        p_previous = loads_previous[index]
        p = mesh.vertex_attributes(key, ['px', 'py', 'pz'])
        r = mesh.vertex_attributes(key, ['rx', 'ry', 'rz'])

        r = add_vectors(r, subtract_vectors(p, p_previous))

        mesh.vertex_attributes(key, ['rx', 'ry', 'rz'], r)

        r_length = length_vector(r)
        R.append(r_length)

    residuals = sum(R)

    return residuals

Visualise the resulting residuals.

Store the previous loads before computing the new self-weight!

# ==============================================================================
# Fofin for selfweight
# ==============================================================================

# compute selfweight for the current geometry
selfweight(mesh)

# store previous selfweight loads (zero)
loads_previous = mesh.vertices_attributes(('px', 'py', 'pz'))

# form finding with selfweight loads and geomtry update
fofin(mesh)

# recompute selfweight for the updated geometry
selfweight(mesh)

# compute the residuals with difference of selfweight from updated to previous geometry
residuals = update_residual(mesh, loads_previous)
print(residuals)

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

baselayer = "CSD2::B1::FormFindingSelfweight"

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})})
artist.draw_edges()
artist.draw_faces()

draw_reactions(mesh, baselayer=baselayer)
draw_residuals(mesh, baselayer=baselayer, scale=10)
# draw_loads(mesh, baselayer=baselayer, scale=10)

e. Iterative Procedure

Now, all steps would need to be repeated until the residuals converge below a certain admissible threshold. Instead of copying the code over and over, we wrap the procedure into a for-loop.

Set a maximum number of iterations so that in case you don't find convergence because of a bug, your computer does not crash.

# ==============================================================================
# Fofin for selfweight
# ==============================================================================

# define maximum iterations and tolerance for residuals
tol = 0.001
kmax = 10

# store previous selfweight loads (zero)
loads_previous = mesh.vertices_attributes(('px', 'py', 'pz'))

# compute selfweight for the current geometry
selfweight(mesh)

# for all k smaller than kmax
for k in range(kmax):

    # recompute the residuals with difference of selfweight from updated to previous geometry
    residuals = update_residual(mesh, loads_previous)

    # stopping criteria if updated residual is smaller than tolerance
    print('k', k, 'residuals', residuals)
    if residuals < tol:
        print('Convergence!')
        break

    # form finding with selfweight loads and geomtry update
    fofin(mesh)

    # store previous selfweight loads (computed before the updated geometry)
    loads_previous = mesh.vertices_attributes(('px', 'py', 'pz'))

    # recompute selfweight for the updated geometry
    selfweight(mesh)

Visualise the converged cable net under self-weight with its cable forces. You can also observe the shape change of the previously form-found cable net and the cable net with the concrete load. Export the file for the next step.

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

baselayer = "CSD2::B1::FormFindingSelfweight0"

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})})
artist.draw_edges()
artist.draw_faces()

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

# ==============================================================================
# Export
# ==============================================================================

mesh.to_json(FILE_O)

to compute the self-weight, compas_fofin already provides a class compas_fofin.loads.SelfweightCalculator that is called in the compas_fofin update_xyz_numpy function and that provides more features (e.g. checks if the face is loaded or not).