C1.4 Under the Hood

For those curious about the implementation details of the dynamic relaxation algorithm in the COMPAS framework, we show the general structure of code and their overall logic here. Looking through the code of compas.numerical.dr.dr.py: Firstly, general helper variables, functions, and classes are defined at the start:

# Butcher Tableau for Runge-Kutta integration scheme
K = [
    [0.0],
    [0.5, 0.5],
    [0.5, 0.0, 0.5],
    [1.0, 0.0, 0.0, 1.0],
]

class Coeff():
    """Damping coefficients."""
    def __init__(self, c):
        self.c = c
        self.a = (1 - c * 0.5) / (1 + c * 0.5)
        self.b = 0.5 * (1 + self.a)


def norm_vector(vector): ...


def norm_vectors(vectors): ...


def adjacency_from_edges(edges): ...

Next, inside the body of the dr function, the input arguments are pre-processed and default values are assigned to the remaining iteration parameters:

    # --------------------------------------------------------------------------
    # attribute lists
    # --------------------------------------------------------------------------
    X = vertices                                                           # m
    P = loads                                                              # kN
    Q_pre = qpre or [0.0 for _ in range(e)]                                # kN/m
    F_pre = fpre or [0.0 for _ in range(e)]                                # kN
    L_pre = lpre or [0.0 for _ in range(e)]                                # m
    L_init = linit or [0.0 for _ in range(e)]                              # m
    E = E or [0.0 for _ in range(e)]                                       # kN/mm2 => GPa
    radius = radius or [0.0 for _ in range(e)]                             # mm

    # --------------------------------------------------------------------------
    # sectional properties
    # --------------------------------------------------------------------------
    A = [3.14159 * r ** 2 for r in radius]                                 # mm2
    EA = [_e * _a for _e, _a in zip(E, A)]                                 # kN

    # --------------------------------------------------------------------------
    # initial values
    # --------------------------------------------------------------------------
    Q = [1.0 for _ in range(e)]
    L = [sum((X[i][axis] - X[j][axis]) ** 2
        for axis in (0, 1, 2)) ** 0.5 for i, j in iter(edges)]
    # if none of the initial lengths are set,
    # set the initial lengths to the current lengths    if all(l == 0 for l in L_init):
        L_init = L[:]
    F = [q * l for q, l in zip(Q, L)]
    M = [sum(0.5 * dt ** 2 * Q[ij_e[(i, j)]] for j in i_nbrs[i]) for i in range(n)]
    V = [[0.0, 0.0, 0.0] for _ in range(n)]
    R = [[0.0, 0.0, 0.0] for _ in range(n)]
    dX = [[0.0, 0.0, 0.0] for _ in range(n)]

Next, an update_residuals and Runge-Kutta integration function are declared. The latter is nothing more than a refined method for deriving the nodal residuals at each step. Runge-Kutta integration evaluates multiple forward steps, and returns a weighted average of the residuals for each forward step. Depending on the number of forward steps (2 or 4), different weighting factors are used. This gives a faster and more robust convergence. For a concise mathematical background, those interested are referred here.

The logic core of the algorithm is defined in the iteration section. It can be seen that the different edge target inputs (force density, force, length and strain) all compete through their summation into one force density per edge. Hence, the algorithm allows to input these various edge target variables, but providing competing values for them all at once may give unexpected and/or uncontrollable results.

# --------------------------------------------------------------------------
# iteration
# --------------------------------------------------------------------------
for k in range(kmax):
    # get force density contribution form each source
    Q_fpre = [f / l if l else 0 for f, l in zip(F_pre, L)]
    Q_lpre = [f / l if l else 0 for f, l in zip(F, L_pre)]
    Q_ea = [ea * (l - linit) / (linit * l) if (l * linit)
            else 0 for ea, l, linit in zip(EA, L, L_init)]

    # set total force density as sum from each source
    Q = [qq + qf + ql + qea for qq, qf, ql, qea
        in zip(Q_pre, Q_fpre, Q_lpre, Q_ea)]
    # set fictitious nodal mass
    M = [sum(0.5 * dt ** 2 * Q[ij_e[(i, j)]] for j in i_nbrs[i]) for i in range(n)]

    # store geometry of previous iteration
    X0 = deepcopy(X)
    # get damped nodal velocities
    V0 = [[ca * V[i][axis] for axis in (0, 1, 2)] for i in range(n)]
    # get damped nodal accelerations from Runge-Kutta integration
    dV = runge_kutta(X0, V0, steps=2)
    # update nodal coordinates
    for i in free:
        V[i] = [V0[i][axis] + dV[i][axis] for axis in (0, 1, 2)]
        dX[i] = [V[i][axis] * dt for axis in (0, 1, 2)]
        X[i] = [X0[i][axis] + dX[i][axis] for axis in (0, 1, 2)]
        
    # update edge lengths, forces, nodal residuals
    L = [sum((X[i][axis] - X[j][axis]) ** 2 for axis in (0, 1, 2)) ** 0.5 for i, j in iter(edges)]
    F = [q * l for q, l in zip(Q, L)]
    update_R()

    # stopping critertia
    crit1 = max(norm_vectors([R[i] for i in free]))
    crit2 = max(norm_vectors([dX[i] for i in free]))

    # callback
    if callback:
        callback(k, X, (crit1, crit2), callback_args)

    # convergence
    if crit1 < tol1:
        break
    if crit2 < tol2:
        break

# --------------------------------------------------------------------------
# update
# --------------------------------------------------------------------------
update_R()

return X, Q, F, L, R

One of the stopping criteria is defined as before: the maximum residual reaching under a certain threshold. An alternative criterion is provided as a kind of energy convergence: the maximum nodal displacement is a measure of how much kinematic energy there is in the dynamic system. Finally, the updated nested lists of coordinates, force densities, edge lengths and residuals are returned.