Implement color cycler, bugfixes

manuvarkey · Jun 19, 2022 · 0012249 · 0012249
1 parent 9c4f2b9
commit 0012249
Show file tree

Hide file tree

Showing 7 changed files with 142 additions and 125 deletions.
diff --git a/Earth design simulation.ipynb b/Earth design simulation.ipynb
@@ -39,7 +39,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/home/eee/Projects/Earthing/earthing/__init__.py:692: RuntimeWarning: divide by zero encountered in true_divide\n",
+      "/home/eee/Projects/Earthing/earthing/__init__.py:698: RuntimeWarning: divide by zero encountered in true_divide\n",
       "  COEF = self.rho/(4*pi*a)*np.arctan(a/ALPHA)\n"
      ]
     }
@@ -98,14 +98,14 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Value of earth resistance =  [4.901   inf]  Ohm\n"
+      "Value of earth resistance =  [4.902   inf]  Ohm\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/home/eee/Projects/Earthing/earthing/__init__.py:730: RuntimeWarning: divide by zero encountered in true_divide\n",
+      "/home/eee/Projects/Earthing/earthing/__init__.py:736: RuntimeWarning: divide by zero encountered in true_divide\n",
       "  res = self.V / self.Ig\n"
      ]
     }
@@ -141,7 +141,7 @@
     {
      "data": {
       "text/plain": [
-       "array([4901.104, 2511.734])"
+       "array([4902.068, 2512.788])"
       ]
      },
      "execution_count": 5,
@@ -165,7 +165,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "1618.714  @  [5.  2.5]\n"
+      "1619.56  @  [5.  2.5]\n"
      ]
     }
    ],
@@ -186,7 +186,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "817.389  @  [9.961 5.043]\n"
+      "835.448  @  [9.961 4.973]\n"
      ]
     }
    ],
@@ -206,7 +206,7 @@
     {
      "data": {
       "text/plain": [
-       "445.0587225606989"
+       "447.86775373126284"
       ]
      },
      "execution_count": 8,

diff --git a/earthing/__init__.py b/earthing/__init__.py
@@ -29,6 +29,7 @@
 from matplotlib import path
 from mpl_toolkits.mplot3d import axes3d
 from matplotlib import cm
+from cycler import cycler
 
 __all__ = ["Network",
            "resistance_plate", "resistance_pipe", "resistance_strip",
@@ -40,7 +41,11 @@
 # Limit visibility of modules to __all__
 def __dir__():
     return __all__
-
+
+# Global variables
+def plot_cycler():
+    return (cycler(color=list('grcy')))
+eplot_linewidth = 3
 
 ## Resistance functions
 
@@ -508,7 +513,7 @@ def add_mesh(self, loc, Lx, Ly, Nx, Ny, w):
             self.add_strip(loc_start, loc_end, w)
         # Add second row offset by 2*w to avoid computation errors
         for i in range(Ny):
-            loc_start = np.array(loc) + np.array([i*Lx/(Ny-1), 0, -2*w])
+            loc_start = np.array(loc) + np.array([i*Lx/(Ny-1), 0, 0])
             loc_end = loc_start + np.array([0, Ly, 0])
             self.add_strip(loc_start, loc_end, w)
 
@@ -526,13 +531,14 @@ def add_plate(self, loc, w, h, n_cap=[1,0,0], h_cap=[0,0,1]):
 
     # Solution functions
 
-    def mesh_geometry(self, desc_size=0.25):
+    def mesh_geometry(self, desc_size=0.2):
         """Descretise geometry"""
 
         if not self.elements[-1]:
             raise Exception('No problem geometry found')
 
         self.descrete_elements = []
+        self.subnet_sizes = []
 
         # Descretise Geometry
         for subnet in self.elements:
@@ -987,15 +993,16 @@ def plot_surface_potential(self, xlim=(-20, 20), ylim=(-20, 20),
 
 
         # Plot problem geometry as lines
-        for subnet in self.elements:
+        styler = plot_cycler()
+        for subnet, style in zip(self.elements, styler):
             for element in subnet:
                 if isinstance(element, NetworkElementStrip) \
                    or isinstance(element, NetworkElementPipe):
                     start = element.loc
                     end = element.loc_end
                     X = [start[0], end[0]]
                     Y = [start[1], end[1]]
-                    ax.plot(X, Y, color='green')  # plot line
+                    ax.plot(X, Y, **style, linewidth=eplot_linewidth)  # plot line
                 if isinstance(element, NetworkElementPlate):
                     c1 = element.loc - element.w_cap*element.w/2 \
                          - element.h_cap*element.h/2
@@ -1004,7 +1011,7 @@ def plot_surface_potential(self, xlim=(-20, 20), ylim=(-20, 20),
                     c4 = c3 - element.w_cap * element.w
                     X = [c1[0], c2[0], c3[0], c4[0], c1[0]]
                     Y = [c1[1], c2[1], c3[1], c4[1], c1[1]]
-                    ax.plot(X, Y, color='green')  # plot line
+                    ax.plot(X, Y, **style, linewidth=eplot_linewidth)  # plot line
 
         ax.set_xticks(np.arange(*xlim, grid_spacing))
         ax.set_yticks(np.arange(*ylim, grid_spacing))
@@ -1049,6 +1056,32 @@ def plot_geometry_3d(self, xlim=(-20, 20), ylim=(-20, 20), zlim=(-5, 5),
                                        offset=0, cmap="plasma", alpha=0.7)
             cbar = plt.colorbar(contour_plot, pad=0.1)
             # cbar.set_label('Volts', loc='bottom')
+
+        # Plot direction vectors for descrete elements
+        if normal:
+            X = []
+            Y = []
+            Z = []
+            U = []
+            V = []
+            W = []    
+            for desc_element in self.descrete_elements:
+                if isinstance(desc_element, DescreteElementPlate):
+                    x, y, z = desc_element.loc
+                    u, v, w = desc_element.normal
+                    X.append(x)
+                    Y.append(y)
+                    Z.append(z)
+                    U.append(u*normal_scale)
+                    V.append(v*normal_scale)
+                    W.append(w*normal_scale)
+            ax.quiver(X, Y, Z, U, V, W)  # quiver
+
+        # Set plot general properties    
+        ax.set(xlim=xlim, ylim=ylim, zlim=zlim,
+               xlabel='X', ylabel='Y', zlabel='Z')
+        plt.title(title)
+        plt.show()
 
         if current_distribution:
             # Plot problem geometry with current as weight
@@ -1065,7 +1098,8 @@ def plot_geometry_3d(self, xlim=(-20, 20), ylim=(-20, 20), zlim=(-5, 5),
             # cbar.set_label('Amps', loc='bottom')
         else:
             # Plot problem geometry as lines
-            for subnet in self.elements:
+            styler = plot_cycler()
+            for subnet, style in zip(self.elements, styler):
                 for element in subnet:
                     if isinstance(element, NetworkElementStrip) \
                        or isinstance(element, NetworkElementPipe):
@@ -1074,7 +1108,7 @@ def plot_geometry_3d(self, xlim=(-20, 20), ylim=(-20, 20), zlim=(-5, 5),
                         X = [start[0], end[0]]
                         Y = [start[1], end[1]]
                         Z = [start[2], end[2]]
-                        ax.plot(X, Y, Z, color='green')  # plot line
+                        ax.plot(X, Y, Z, **style, linewidth=eplot_linewidth)  # plot line
                     if isinstance(element, NetworkElementPlate):
                         c1 = element.loc - element.w_cap*element.w/2 \
                              - element.h_cap*element.h/2
@@ -1084,33 +1118,6 @@ def plot_geometry_3d(self, xlim=(-20, 20), ylim=(-20, 20), zlim=(-5, 5),
                         X = [c1[0], c2[0], c3[0], c4[0], c1[0]]
                         Y = [c1[1], c2[1], c3[1], c4[1], c1[1]]
                         Z = [c1[2], c2[2], c3[2], c4[2], c1[2]]
-                        ax.plot(X, Y, Z, color='green')  # plot line
+                        ax.plot(X, Y, Z, **style, linewidth=eplot_linewidth)  # plot line
 
-        # Plot direction vectors for descrete elements
-        if normal:
-            X = []
-            Y = []
-            Z = []
-            U = []
-            V = []
-            W = []    
-            for desc_element in self.descrete_elements:
-                if isinstance(desc_element, DescreteElementPlate):
-                    x, y, z = desc_element.loc
-                    u, v, w = desc_element.normal
-                    X.append(x)
-                    Y.append(y)
-                    Z.append(z)
-                    U.append(u*normal_scale)
-                    V.append(v*normal_scale)
-                    W.append(w*normal_scale)
-            ax.quiver(X, Y, Z, U, V, W)  # quiver
 
-        # Set plot general properties    
-        ax.set(xlim=xlim, ylim=ylim, zlim=zlim,
-               xlabel='X', ylabel='Y', zlabel='Z')
-        plt.title(title)
-        plt.show()
-
-
-
diff --git a/setup.cfg b/setup.cfg
@@ -1,6 +1,6 @@
 [metadata]
 name = earthing
-version = 1.0.2
+version = 1.1.0
 author = Manu Varkey
 author_email = [email protected]
 description = A python library for design of earthing networks in electrical substations.

diff --git a/tests/Earth design - performance test.ipynb b/tests/Earth design - performance test.ipynb
@@ -18,6 +18,9 @@
     "from math import sqrt, log, pi, ceil, floor\n",
     "import time\n",
     "\n",
+    "import numpy as np\n",
+    "from numpy.linalg import norm\n",
+    "\n",
     "import matplotlib.pyplot as plt\n",
     "plt.rcParams[\"font.family\"] = \"Ubuntu Mono\"\n",
     "plt.rcParams['font.size'] = 9\n",
@@ -59,11 +62,11 @@
     "ylim=(-5, Ly+5)\n",
     "zlim=(-5, 2)\n",
     "grid=(25,25)\n",
-    "delta=0.2\n",
+    "delta=0.1\n",
     "Ig = 1000\n",
     "\n",
     "# Define network\n",
-    "network = Network(rho)\n",
+    "network = Network(rho, Ig)\n",
     "# network.add_strip([0, 0, -h], [L, 0, -h], strip_width)\n",
     "# network.add_mesh([0.1,0.1,-h], Lx, Ly, Nx, Ny, strip_width)\n",
     "network.add_rod([0,0,-h], radius_rod, Lr)\n",
@@ -83,16 +86,16 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "--- 0.001352548599243164 seconds ---\n",
-      "--- 0.002187013626098633 seconds ---\n"
+      "--- 0.003772735595703125 seconds ---\n",
+      "--- 0.0013759136199951172 seconds ---\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/home/eee/Projects/Earthing/tests/earthing/__init__.py:594: RuntimeWarning: divide by zero encountered in true_divide\n",
-      "  COEF = self.rho/(4*pi*AA)*np.arctan(AA/ALPHA)\n"
+      "/home/eee/Projects/Earthing/earthing/__init__.py:693: RuntimeWarning: divide by zero encountered in true_divide\n",
+      "  COEF = self.rho/(4*pi*a)*np.arctan(a/ALPHA)\n"
      ]
     }
    ],
@@ -121,7 +124,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Value of earth resistance =  56.084  Ohm\n"
+      "Value of earth resistance =  [58.521]  Ohm\n"
      ]
     }
    ],
@@ -171,18 +174,18 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "--- 0.12018370628356934 seconds ---\n",
-      "--- 0.0004189014434814453 seconds ---\n"
+      "--- 0.22726035118103027 seconds ---\n",
+      "--- 0.0021457672119140625 seconds ---\n"
      ]
     }
    ],
    "source": [
     "start_time = time.time()\n",
-    "network.solve_surface_potential(Ig, grid=grid, xlim=xlim, ylim=ylim)\n",
+    "network.solve_surface_potential(grid=grid, xlim=xlim, ylim=ylim)\n",
     "print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
     "\n",
     "start_time = time.time()\n",
-    "network.solve_surface_potential_fast(Ig, grid=grid, xlim=xlim, ylim=ylim)\n",
+    "network.solve_surface_potential_fast(grid=grid, xlim=xlim, ylim=ylim)\n",
     "print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
     "network.plot_geometry_3d(xlim=xlim, ylim=ylim, zlim=zlim, ground=False, ground_pot=True)"
    ]
@@ -208,7 +211,7 @@
     {
      "data": {
       "text/plain": [
-       "56084.0"
+       "array([58521.296])"
       ]
      },
      "execution_count": 8,
@@ -217,7 +220,7 @@
     }
    ],
    "source": [
-    "network.gpr(Ig)"
+    "network.gpr()"
    ]
   },
   {
@@ -229,7 +232,7 @@
     {
      "data": {
       "text/plain": [
-       "2367.3769208768726"
+       "2364.5471493780733"
       ]
      },
      "execution_count": 9,
@@ -239,7 +242,7 @@
    ],
    "source": [
     "# Touch voltage at fence 2m away\n",
-    "network.get_point_potential([0,-2,0], Ig) - network.get_point_potential([0,-3,0], Ig)"
+    "network.get_point_potential([0,-2,0]) - network.get_point_potential([0,-3,0])"
    ]
   },
   {
@@ -251,7 +254,7 @@
     {
      "data": {
       "text/plain": [
-       "53258.06549580934"
+       "array([55695.80818397])"
       ]
      },
      "execution_count": 10,
@@ -261,7 +264,7 @@
    ],
    "source": [
     "# Mesh voltage\n",
-    "network.gpr(Ig) - network.get_point_potential([Lx/Nx,Ly/Ny,0], Ig)"
+    "network.gpr() - network.get_point_potential([Lx/Nx,Ly/Ny,0])"
    ]
   },
   {
@@ -273,7 +276,7 @@
     {
      "data": {
       "text/plain": [
-       "23146.33896766483"
+       "25265.52798750884"
       ]
      },
      "execution_count": 11,
@@ -285,7 +288,7 @@
     "# Step voltage\n",
     "loc = np.array([-Lx, -Ly, 0])\n",
     "loc = loc/norm(loc)\n",
-    "network.get_point_potential([0,0,0], Ig) - network.get_point_potential(loc, Ig)"
+    "network.get_point_potential([0,0,0]) - network.get_point_potential(loc)"
    ]
   },
   {