Added meshing to the block structures and working on non-linear mass

Pavel Govyadinov
1 parent 68f91f03
Showing 1 changed file with 122 additions and 14 deletions Show diff stats
voronoi_test.py
@@ -7,6 +7,10 @@ Created on Tue Sep  3 13:15:54 2019
 """
 from scipy.spatial import Voronoi, voronoi_plot_2d
+from scipy.spatial import Delaunay, delaunay_plot_2d
+
+from triangle import triangulate
+from triangle import plot as triangle_plot
 from scipy.interpolate import interp1d
 #import matplotlib._cntr as cntr
 from shapely.geometry import Point
@@ -30,13 +34,59 @@ class Polygon_mass:
         self.get_aabb()
         self.gen_polygon()
         self.torque = []
+        self.vel = 0.0
+        self.aa = 0.0
+      
+        
+    def clear_torques(self):
+        self.torque = []        
     def add_torque(self, p, f):
+        #direction of the torque = cross of (r, f)
+        #magnitude = ||r||*||f||*sin(theta)
+        #r = level arm vector
         d = self.CoM - p
         r = np.linalg.norm(self.CoM - p)
         theta = math.acos(np.dot(d, f)/np.dot(d, d)/np.dot(f, f))
-        torque = r * math.sin(theta) * f
-        self.torque.append(torque)
+        torque = math.sin(theta) * np.sqrt(f[0]*f[0]+f[1]*f[1]) * np.sqrt(r[0]*r[0]+r[1]*r[1])
+        #if < 0 then clockwise, else counter
+        direction = np.cross(r, f)
+        if direction < 0:
+            self.torque.append([torque, "clockwise"])
+        else:
+            self.torque.append([torque, "counterclock"])
+            
+    def calculate_moment(self):
+        
+        #returns the area of a triangle defined by two points
+        def area(t):
+            output = 1.0/2.0*abs(t[0,0]*(t[1,1]-t[2,1]) + t[1,0]*(t[2,1]-t[0,1]) + t[2,0]*(t[0,1]-t[1,1]))
+            return output
+        
+        def center(t):
+            output = np.asarray([(t[0,0]+t[1,0]+t[2,0])/3.0, (t[0,1]+t[1,1]+t[2,1])/3.0])
+        
+        segs = []
+        pts = np.asarray(self.polygon.exterior.xy).T
+        pts = pts[:-1, :]
+        for k in range(pts.shape[0]-1):
+            segs.append([k, k+1])
+        segs.append([pts.shape[0]-1, 0])
+        mesh = dict(vertices=pts, segments=segs)
+        #print(self.polygon.area())
+        tri = triangulate(mesh, 'pq20D')
+        moment = 0.0
+        #NEED TO ADD MASS maybe?
+        for i in range(tri['triangles'].shape[0]):
+            t = tri['vertices'][tri['triangles'][i]]
+            A = area(t)
+            C = center(t)
+            Moi = A+A*np.linalg.norm(C-self.CoM)**2
+            moment += Moi
+            
+        self.MoI = moment        
+#        triangle_plot(plt.gca(), **tri)
+#        plt.gca().set_title(str(self.polygon.area))
     def rotate(self, phi, direction = "counterclock"):
         if("counterclock"):
@@ -62,6 +112,35 @@ class Polygon_mass:
         plt.scatter(p[:,0], p[:,1], color='r')
         plt.scatter(self.CoM[0], self.CoM[1], marker='*')
+        
+        segs = []
+        pts = np.asarray(self.polygon.exterior.xy).T
+        pts = pts[:-1, :]
+        for k in range(pts.shape[0]-1):
+            segs.append([k, k+1])
+        segs.append([pts.shape[0]-1, 0])
+        points = self.G.vertex_properties["pos"].get_2d_array(range(2)).T
+        segs2 = []
+        n_pts = pts.shape[0]
+        for e in self.G.edges():
+            segs2.append([int(e.source())+n_pts, int(e.target())+n_pts])
+        
+        
+        
+        pts = np.concatenate((pts, points))
+        segs = segs + segs2
+        mesh = dict(vertices=pts, segments=segs)
+        #print(self.polygon.area())
+ #       tri = triangulate(mesh, 'pq20D')
+        triangle_plot(plt.gca(), **mesh)
+        for i in range(len(segs)):
+            plt.plot((pts[segs[i][0]][0], pts[segs[i][0]][1]), (pts[segs[i][1]][0], pts[segs[i][1]][1]), color='b')
+        plt.gca().set_title(str(self.polygon.area))
+        
+        #tri = Delaunay(np.asarray(self.polygon.exterior.coords.xy).T)
+        #tri = triangulate(mesh, 'pq20a' + str(self.polygon.area/100.0)+'D')
+        #delaunay_plot_2d(tri)
+        
 #        for n, contour in enumerate(self.cn):
 #            X = interp1d(np.arange(0, x.shape[0]), x)
 #            Y = interp1d(np.arange(0, y.shape[0]), y)
@@ -73,14 +152,14 @@ class Polygon_mass:
 #        level = (mx-mn)/5.5
 #        cn = plt.contour(x, y, D, levels = [level])
-        for e in self.G.edges():
-            coord = self.G.vertex_properties["pos"][e.source()]
-            coord2 = self.G.vertex_properties["pos"][e.target()]
-            X = [coord[0], coord2[0]]
-            Y = [coord[1], coord2[1]]
-            #all_plots.plot(x, y, 'go--', linewidth=1, markersize=1)
-            plt.plot(X, Y, 'go--', linewidth=1, markersize=1)
-            
+#        for e in self.G.edges():
+#            coord = self.G.vertex_properties["pos"][e.source()]
+#            coord2 = self.G.vertex_properties["pos"][e.target()]
+#            X = [coord[0], coord2[0]]
+#            Y = [coord[1], coord2[1]]
+#            #all_plots.plot(x, y, 'go--', linewidth=1, markersize=1)
+#            plt.plot(X, Y, 'go--', linewidth=1, markersize=1)
+#            
         plt.plot(*self.polygon.exterior.xy, color = 'r')
         plt.show()
@@ -211,8 +290,33 @@ class Polygon_mass:
         C_x = (1.0 / (6.0 * A)) * C_x
         C_y = (1.0 / (6.0 * A)) * C_y
-        return np.array([C_x, C_y])    
-        
+        return np.array([C_x, C_y])  
+    
+#Graph G
+#List of Polygonmass objects with the same cluster index.
+def get_torques(G, masses):
+    for i in masses:
+        i.clear_torques()
+    for e in G.edges():
+        #if the source and target cluster is not equal to each other
+        #add an inter subgraph edge.
+        if(G.vertex_properties["clusters"][e.source()] != G.vertex_properties["clusters"][e.target()]):
+            pf0 = G.vertex_properties["pos"][e.target()].get_2D_array(range(2))
+            pf1 = G.vertex_properties["pos"][e.source()].get_2D_array(range(2))
+            
+            t0 = G.vertex_properties["clusters"][e.source()]
+            t1 = G.vertex_properties["clusters"][e.target()]
+            
+            v0_index = G.vertex_properties["idx"][e.source()]
+            v1_index = G.vertex_properties["idx"][e.target()]
+            
+            p0 = masses[t0].G.vertex_properties["pos"][masses[t0].G.vertex(v0_index)].get_2D_array(range(2))
+            p1 = masses[t1].G.vertex_properties["pos"][masses[t1].G.vertex(v1_index)].get_2D_array(range(2))
+            
+            f0 = pf0-p1
+            f1 = pf1-p0
+            masses[t0].add_torque(p0, f0)
+            masses[t1].add_torque(p1, f1)
 def voronoi_polygons(voronoi, diameter):
@@ -517,15 +621,19 @@ def gen_image(G, G_c, itr, bb_flag = False, bb = None, reposition = False):
         no_links.plot(x, y, 'go--', linewidth=1, markersize=1)
         voronoi.plot(x, y, 'go--', linewidth=1, markersize=1)
-    #for every subgraph generate a layout and plot the result
+    #for every subgraph generate a layout and plot the resulting Polygon_mass
+    masses = []
     for i in range(num_clusters):
         g, center = gen_subclusters(G, G_c, i, reposition)
         d = G_c.vertex_properties["pos"][i] - center
         t = Polygon_mass(g)
+        masses.append(g)
         #t.distancefield()
         for v in g.vertices():
             G.vertex_properties["pos"][g.vertex_properties["idx"][v]] = g.vertex_properties["pos"][v] + d
-            #g.vertex_properties["pos"][g.vertex_properties["idx"][v]] = g.vertex_properties["pos"][v] + d
+            
+            
+        #g.vertex_properties["pos"][g.vertex_properties["idx"][v]] = g.vertex_properties["pos"][v] + d
         #sub_pts = g.vertex_properties["pos"].get_2d_array(range(2)).T
     #all_plots.scatter(pts[:,0], pts[:, 1], marker="*")
     #    for e in g.edges():