Use matplotlib tri interpolators (#25)

* Use matplotlib tri interpolators. This is much faster than scipy interpolators because we can use the existing Delaunay triangulation

.gitignore

@@ -1,3 +1,4 @@
+.vscode
 .DS_Store
 
 # Byte-compiled / optimized / DLL files

docs/notebooks/logo.ipynb

@@ -232,7 +232,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Constructing Voronoi polygons: 100%|█████| 4347/4347 [00:00<00:00, 22193.19it/s]\n"
+      "Constructing Voronoi polygons: 100%|██████████████████████████████████████████████████████████████████████████| 4347/4347 [00:00<00:00, 22861.61it/s]\n"
      ]
     },
     {
@@ -275,7 +275,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Simulating: 100%|██████████████████████████▉| 800/800 [00:49<00:00, 16.26tau/s ]\n"
+      "Simulating: 100%|███████████████████████████████████████████████████████████████████████████████████████████████▉| 800/800 [00:49<00:00, 16.12tau/s ]\n"
      ]
     }
    ],
@@ -8744,7 +8744,7 @@
     {
      "data": {
       "text/html": [
-       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>tdgl</td><td>0.3.1</td></tr><tr><td>Numpy</td><td>1.24.3</td></tr><tr><td>SciPy</td><td>1.10.1</td></tr><tr><td>matplotlib</td><td>3.7.1</td></tr><tr><td>jax</td><td>None</td></tr><tr><td>numba</td><td>0.57.0</td></tr><tr><td>IPython</td><td>8.14.0</td></tr><tr><td>Python</td><td>3.10.11 | packaged by conda-forge | (main, May 10 2023, 19:01:19) [Clang 14.0.6 ]</td></tr><tr><td>OS</td><td>posix [darwin]</td></tr><tr><td>Number of CPUs</td><td>Physical: 10, Logical: 10</td></tr><tr><td>BLAS Info</td><td>OPENBLAS</td></tr><tr><td colspan='2'>Thu Jul 20 12:33:35 2023 PDT</td></tr></table>"
+       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>tdgl</td><td>0.4.0</td></tr><tr><td>Numpy</td><td>1.24.3</td></tr><tr><td>SciPy</td><td>1.10.1</td></tr><tr><td>matplotlib</td><td>3.7.1</td></tr><tr><td>jax</td><td>None</td></tr><tr><td>numba</td><td>0.57.0</td></tr><tr><td>IPython</td><td>8.14.0</td></tr><tr><td>Python</td><td>3.10.11 | packaged by conda-forge | (main, May 10 2023, 19:01:19) [Clang 14.0.6 ]</td></tr><tr><td>OS</td><td>posix [darwin]</td></tr><tr><td>Number of CPUs</td><td>Physical: 10, Logical: 10</td></tr><tr><td>BLAS Info</td><td>OPENBLAS</td></tr><tr><td colspan='2'>Thu Aug 03 14:51:31 2023 PDT</td></tr></table>"
       ],
       "text/plain": [
        "<IPython.core.display.HTML object>"

docs/notebooks/mesh.ipynb

@@ -29,7 +29,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Constructing Voronoi polygons: 100%|████████████████████████████████████████████████████████████████████████████| 525/525 [00:00<00:00, 22006.23it/s]\n"
+      "Constructing Voronoi polygons: 100%|████████████████████████████████████████████████████████████████████████████| 525/525 [00:00<00:00, 21276.69it/s]\n"
      ]
     }
    ],

docs/notebooks/polygons.ipynb

@@ -469,9 +469,9 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Constructing Voronoi polygons: 100%|███████| 459/459 [00:00<00:00, 16784.24it/s]\n",
-      "Constructing Voronoi polygons: 100%|███████| 459/459 [00:00<00:00, 16851.97it/s]\n",
-      "Constructing Voronoi polygons: 100%|███████| 459/459 [00:00<00:00, 16564.87it/s]\n"
+      "Constructing Voronoi polygons: 100%|████████████████████████████████████████████████████████████████████████████| 459/459 [00:00<00:00, 16241.63it/s]\n",
+      "Constructing Voronoi polygons: 100%|████████████████████████████████████████████████████████████████████████████| 459/459 [00:00<00:00, 17675.39it/s]\n",
+      "Constructing Voronoi polygons: 100%|████████████████████████████████████████████████████████████████████████████| 459/459 [00:00<00:00, 17257.42it/s]\n"
      ]
     },
     {
@@ -520,7 +520,7 @@
     {
      "data": {
       "text/html": [
-       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>tdgl</td><td>0.3.1</td></tr><tr><td>Numpy</td><td>1.24.3</td></tr><tr><td>SciPy</td><td>1.10.1</td></tr><tr><td>matplotlib</td><td>3.7.1</td></tr><tr><td>jax</td><td>None</td></tr><tr><td>numba</td><td>0.57.0</td></tr><tr><td>IPython</td><td>8.14.0</td></tr><tr><td>Python</td><td>3.10.11 | packaged by conda-forge | (main, May 10 2023, 19:01:19) [Clang 14.0.6 ]</td></tr><tr><td>OS</td><td>posix [darwin]</td></tr><tr><td>Number of CPUs</td><td>Physical: 10, Logical: 10</td></tr><tr><td>BLAS Info</td><td>OPENBLAS</td></tr><tr><td colspan='2'>Thu Jul 20 12:31:45 2023 PDT</td></tr></table>"
+       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>tdgl</td><td>0.4.0</td></tr><tr><td>Numpy</td><td>1.24.3</td></tr><tr><td>SciPy</td><td>1.10.1</td></tr><tr><td>matplotlib</td><td>3.7.1</td></tr><tr><td>jax</td><td>None</td></tr><tr><td>numba</td><td>0.57.0</td></tr><tr><td>IPython</td><td>8.14.0</td></tr><tr><td>Python</td><td>3.10.11 | packaged by conda-forge | (main, May 10 2023, 19:01:19) [Clang 14.0.6 ]</td></tr><tr><td>OS</td><td>posix [darwin]</td></tr><tr><td>Number of CPUs</td><td>Physical: 10, Logical: 10</td></tr><tr><td>BLAS Info</td><td>OPENBLAS</td></tr><tr><td colspan='2'>Thu Aug 03 14:38:40 2023 PDT</td></tr></table>"
       ],
       "text/plain": [
        "<IPython.core.display.HTML object>"

docs/notebooks/py-mesh.ipynb

@@ -184,7 +184,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Constructing Voronoi polygons: 100%|██████████████████████████████████████████████████████████████████████████| 4190/4190 [00:00<00:00, 22017.73it/s]\n"
+      "Constructing Voronoi polygons: 100%|██████████████████████████████████████████████████████████████████████████| 4190/4190 [00:00<00:00, 22386.66it/s]\n"
      ]
     },
     {
@@ -361,7 +361,7 @@
     {
      "data": {
       "text/html": [
-       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>tdgl</td><td>0.3.0</td></tr><tr><td>Numpy</td><td>1.24.3</td></tr><tr><td>SciPy</td><td>1.10.1</td></tr><tr><td>matplotlib</td><td>3.7.1</td></tr><tr><td>jax</td><td>None</td></tr><tr><td>numba</td><td>0.57.0</td></tr><tr><td>IPython</td><td>8.14.0</td></tr><tr><td>Python</td><td>3.10.11 | packaged by conda-forge | (main, May 10 2023, 19:01:19) [Clang 14.0.6 ]</td></tr><tr><td>OS</td><td>posix [darwin]</td></tr><tr><td>Number of CPUs</td><td>Physical: 10, Logical: 10</td></tr><tr><td>BLAS Info</td><td>OPENBLAS</td></tr><tr><td colspan='2'>Thu Jun 08 13:56:32 2023 PDT</td></tr></table>"
+       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>tdgl</td><td>0.4.0</td></tr><tr><td>Numpy</td><td>1.24.3</td></tr><tr><td>SciPy</td><td>1.10.1</td></tr><tr><td>matplotlib</td><td>3.7.1</td></tr><tr><td>jax</td><td>None</td></tr><tr><td>numba</td><td>0.57.0</td></tr><tr><td>IPython</td><td>8.14.0</td></tr><tr><td>Python</td><td>3.10.11 | packaged by conda-forge | (main, May 10 2023, 19:01:19) [Clang 14.0.6 ]</td></tr><tr><td>OS</td><td>posix [darwin]</td></tr><tr><td>Number of CPUs</td><td>Physical: 10, Logical: 10</td></tr><tr><td>BLAS Info</td><td>OPENBLAS</td></tr><tr><td colspan='2'>Thu Aug 03 14:38:28 2023 PDT</td></tr></table>"
       ],
       "text/plain": [
        "<IPython.core.display.HTML object>"

docs/notebooks/quickstart.ipynb

@@ -277,7 +277,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Constructing Voronoi polygons: 100%|█████| 4633/4633 [00:00<00:00, 22649.16it/s]\n"
+      "Constructing Voronoi polygons: 100%|██████████████████████████████████████████████████████████████████████████| 4633/4633 [00:00<00:00, 22269.09it/s]\n"
      ]
     }
    ],
@@ -368,8 +368,8 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Thermalizing: 100%|████████████████████████▉| 100/100 [00:09<00:00, 10.46tau/s ]\n",
-      "Simulating: 100%|██████████████████████████▉| 150/150 [00:25<00:00,  5.78tau/s ]\n"
+      "Thermalizing: 100%|█████████████████████████████████████████████████████████████████████████████████████████████▉| 100/100 [00:10<00:00,  9.80tau/s ]\n",
+      "Simulating: 100%|███████████████████████████████████████████████████████████████████████████████████████████████▉| 150/150 [00:26<00:00,  5.65tau/s ]\n"
      ]
     }
    ],
@@ -8097,7 +8097,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Simulating: 100%|██████████████████████████▉| 200/200 [00:16<00:00, 11.93tau/s ]\n"
+      "Simulating: 100%|███████████████████████████████████████████████████████████████████████████████████████████████▉| 200/200 [00:17<00:00, 11.49tau/s ]\n"
      ]
     }
    ],
@@ -8199,15 +8199,15 @@
       "\tTotal fluxoid: 1.01 Phi_0\n",
       "\n",
       "Round hole:\n",
-      "\tFluxoid(flux_part=1.3657295446325866, supercurrent_part=0.6250145703049892) Phi_0\n",
+      "\tFluxoid(flux_part=1.3657295446325866, supercurrent_part=0.6250145703049893) Phi_0\n",
       "\tTotal fluxoid: 1.99 Phi_0\n",
       "\n",
       "Square hole:\n",
       "\tFluxoid(flux_part=1.3654288409190132, supercurrent_part=0.6221730478571399) Phi_0\n",
       "\tTotal fluxoid: 1.99 Phi_0\n",
       "\n",
       "Bottom vortex:\n",
-      "\tFluxoid(flux_part=0.6065407514631895, supercurrent_part=0.4053668212615612) Phi_0\n",
+      "\tFluxoid(flux_part=0.6065407514631895, supercurrent_part=0.40536682126156137) Phi_0\n",
       "\tTotal fluxoid: 1.01 Phi_0\n",
       "\n"
      ]
@@ -12372,7 +12372,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "Simulating: 100%|██████████████████████████▉| 200/200 [00:38<00:00,  5.14tau/s ]\n"
+      "Simulating: 100%|███████████████████████████████████████████████████████████████████████████████████████████████▉| 200/200 [00:41<00:00,  4.84tau/s ]\n"
      ]
     }
    ],
@@ -27402,7 +27402,7 @@
     {
      "data": {
       "text/html": [
-       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>tdgl</td><td>0.3.1</td></tr><tr><td>Numpy</td><td>1.24.3</td></tr><tr><td>SciPy</td><td>1.10.1</td></tr><tr><td>matplotlib</td><td>3.7.1</td></tr><tr><td>jax</td><td>None</td></tr><tr><td>numba</td><td>0.57.0</td></tr><tr><td>IPython</td><td>8.14.0</td></tr><tr><td>Python</td><td>3.10.11 | packaged by conda-forge | (main, May 10 2023, 19:01:19) [Clang 14.0.6 ]</td></tr><tr><td>OS</td><td>posix [darwin]</td></tr><tr><td>Number of CPUs</td><td>Physical: 10, Logical: 10</td></tr><tr><td>BLAS Info</td><td>OPENBLAS</td></tr><tr><td colspan='2'>Thu Jul 20 12:40:22 2023 PDT</td></tr></table>"
+       "<table><tr><th>Software</th><th>Version</th></tr><tr><td>tdgl</td><td>0.4.0</td></tr><tr><td>Numpy</td><td>1.24.3</td></tr><tr><td>SciPy</td><td>1.10.1</td></tr><tr><td>matplotlib</td><td>3.7.1</td></tr><tr><td>jax</td><td>None</td></tr><tr><td>numba</td><td>0.57.0</td></tr><tr><td>IPython</td><td>8.14.0</td></tr><tr><td>Python</td><td>3.10.11 | packaged by conda-forge | (main, May 10 2023, 19:01:19) [Clang 14.0.6 ]</td></tr><tr><td>OS</td><td>posix [darwin]</td></tr><tr><td>Number of CPUs</td><td>Physical: 10, Logical: 10</td></tr><tr><td>BLAS Info</td><td>OPENBLAS</td></tr><tr><td colspan='2'>Thu Aug 03 14:59:02 2023 PDT</td></tr></table>"
       ],
       "text/plain": [
        "<IPython.core.display.HTML object>"

tdgl/device/device.py

@@ -3,7 +3,7 @@
 import time
 from contextlib import contextmanager, nullcontext
 from operator import attrgetter, itemgetter
-from typing import Any, Dict, List, NamedTuple, Sequence, Tuple, Union
+from typing import Any, Dict, List, NamedTuple, Optional, Sequence, Tuple, Union
 
 import h5py
 import matplotlib.pyplot as plt
@@ -12,6 +12,7 @@
 from IPython.display import HTML
 from matplotlib.patches import PathPatch
 from matplotlib.path import Path
+from matplotlib.tri import Triangulation
 from shapely import affinity
 from shapely.geometry import Point
 
@@ -104,7 +105,8 @@ def __init__(
         # It should never be changed after instantiation.
         self._length_units = length_units
 
-        self.mesh = None
+        self.mesh: Optional[Mesh] = None
+        self._triangulation: Optional[Triangulation] = None
 
     @property
     def length_units(self) -> str:
@@ -202,6 +204,18 @@ def V0(self, conductivity: Union[pint.Quantity, None] = None) -> pint.Quantity:
         J0 = self.K0 / self.thickness
         return (self.coherence_length * J0 / conductivity).to("volts")
 
+    @property
+    def triangulation(self) -> Optional[Triangulation]:
+        """Matplotlib triangulation of the mesh."""
+        if self.mesh is None:
+            return None
+        if self._triangulation is None:
+            xi = self.layer.coherence_length
+            sites = xi * self.mesh.sites
+            triangles = self.mesh.elements
+            self._triangulation = Triangulation(sites[:, 0], sites[:, 1], triangles)
+        return self._triangulation
+
     def terminal_info(self) -> Tuple[TerminalInfo, ...]:
         """Returns a tuple of ``TerminalInfo`` objects,
         one for each current terminal in the device.

tdgl/finite_volume/mesh.py

@@ -1,8 +1,9 @@
-from typing import List, Sequence, Tuple, Union
+from typing import List, Optional, Sequence, Tuple, Union
 
 import h5py
 import matplotlib.pyplot as plt
 import numpy as np
+from matplotlib.tri import Triangulation
 
 from ..geometry import close_curve
 from .edge_mesh import EdgeMesh
@@ -60,6 +61,7 @@ def __init__(
         self.dual_sites = dual_sites
         self.edge_mesh = edge_mesh
         self.voronoi_polygons = voronoi_polygons
+        self._triangulation: Optional[Triangulation] = None
 
     @property
     def x(self) -> np.ndarray:

tdgl/solution/data.py

@@ -1,11 +1,11 @@
 import dataclasses
 import os
-from typing import Any, Dict, List, Optional, Sequence, Tuple, Union
+from typing import Any, Dict, List, Literal, Optional, Sequence, Tuple, Union
 
 import h5py
 import matplotlib.pyplot as plt
+import matplotlib.tri as mtri
 import numpy as np
-from scipy import interpolate
 from tqdm import tqdm
 
 from ..finite_volume.mesh import Mesh
@@ -447,7 +447,7 @@ def get_current_through_paths(
     solution_path: os.PathLike,
     paths: Union[np.ndarray, List[np.ndarray]],
     dataset: Optional[str] = None,
-    interp_method: str = "linear",
+    interp_method: Literal["linear", "cubic"] = "linear",
     units: Optional[str] = None,
     with_units: bool = True,
     progress_bar: bool = True,
@@ -474,20 +474,18 @@ def get_current_through_paths(
 
     solution = Solution.from_hdf5(solution_path)
     device = solution.device
-    mesh = device.mesh
+    tri = device.triangulation
     ureg = device.ureg
 
     valid_methods = ("linear", "cubic")
     if interp_method not in valid_methods:
         raise ValueError(
             f"Interpolation method must be one of {valid_methods} (got {interp_method})."
         )
-    if interp_method == "linear":
-        interpolator = interpolate.LinearNDInterpolator
-        interp_kwargs = dict(fill_value=0)
-    else:  # "cubic"
-        interpolator = interpolate.CloughTocher2DInterpolator
-        interp_kwargs = dict(fill_value=0)
+    interp_type = {
+        "linear": mtri.LinearTriInterpolator,
+        "cubic": mtri.CubicTriInterpolator,
+    }[interp_method]
 
     valid_datasets = ("supercurrent", "normal_current", None)
     if dataset not in valid_datasets:
@@ -516,7 +514,6 @@ def get_current_through_paths(
 
     step_min, step_max = solution.data_range
     times = solution.times
-    sites = device.points
     raw_currents = [np.zeros_like(times) for _ in paths]
     with h5py.File(solution_path, "r") as h5file:
         for i in tqdm(
@@ -527,12 +524,16 @@ def get_current_through_paths(
                 K = np.array(grp["normal_current"]) + np.array(grp["supercurrent"])
             else:
                 K = np.array(grp[dataset])
-            K = mesh.get_quantity_on_site(K)
-            K_interp = interpolator(sites, K, **interp_kwargs)
+            K = device.mesh.get_quantity_on_site(K)
+            Kx_interp = interp_type(tri, K[:, 0])
+            Ky_interp = interp_type(tri, K[:, 1])
             for j, (path, lengths, normals, ix) in enumerate(
                 zip(paths, edge_lengths, unit_normals, in_device)
             ):
-                K_path = K_interp(path)
+                Kx_path = Kx_interp(path[:, 0], path[:, 1]).data
+                Ky_path = Ky_interp(path[:, 0], path[:, 1]).data
+                K_path = np.array([Kx_path, Ky_path]).T
+                K_path[~np.isfinite(K_path).all(axis=1)] = 0
                 # Evaluate the sheet current at the edge centers
                 K_edge = (K_path[:-1] + K_path[1:]) / 2
                 K_dot_n = (K_edge * normals).sum(axis=1)