Dynamic relaxation
In previous posts, we described how to read network data from a file and convert it into matrices relevant for structural calculations. Here, we will use the method of dynamic relaxation to compute an equilibrium shape of the provided network data.
01.from numpy import array, zeros, float6402.from scipy.sparse import coo_matrix, diags03. 04.# Input data:05.# -----------06.# vertices: a list of unique xyz coordinates07.# edges: a list of pairs of vertex indices08.# points: index list of fixed vertices09. 10.m = len(edges)11.n = len(vertices)12. 13.fixed = points14.free = list(set(range(n)) - set(fixed))15. 16.data = array([-1]*m + [1]*m)17.rows = array(range(m) + range(m))18.cols = array(edges).reshape(2*m, order='F')19.C = coo_matrix((data, (rows, cols))).tocsr()20.Ci = C[:, free]21.Cit = Ci.transpose()22.Cf = C[:, fixed] 23.xyz = array(vertices)24.p = zeros((n, 3), dtype=float64)25.q = array([1.0 for i in range(m)]).reshape((-1, 1))26.Q = diags([q.flatten()], [0])27. 28.v = zeros((n, 3), dtype=float64)29.r = zeros((n, 3), dtype=float64)30.kmax = 100031.dt = 1.032.mass = 1.033. 34.for k in xrange(kmax):35. xyz0 = xyz.copy()36. v0 = v.copy()37. dx = v0 * dt38. xyz[free] = xyz0[free] + dx[free] 39. r[free] = p[free] - Cit.dot(Q.dot(C.dot(xyz)))40. a = 0.1 * r[free] / mass41. dv = dt * a42. v[free] = v0[free] + dv43. dx = v * dt44. xyz[free] = xyz0[free] + dx[free]With this file (‘http://block.arch.ethz.ch/labs/samples/saddle.obj’) as input, you should get the following figure of equilibrium. For details on plotting network data, read this post.
