Merge pull request #375 from jcapriot/bug/simplex_caching

Fix caching of internal projection matrices

Author: jcapriot

discretize/unstructured_mesh.py

@@ -405,7 +405,8 @@ def __get_inner_product_projection_matrices(
     ):
         if getattr(self, "_proj_stash", None) is None:
             self._proj_stash = {}
-        if i_type not in self._proj_stash:
+        key = (i_type, with_volume)
+        if key not in self._proj_stash:
             dim = self.dim
             n_cells = self.n_cells
             if i_type == "F":
@@ -456,8 +457,8 @@ def __get_inner_product_projection_matrices(
                 )
                 Ps.append(T @ P)
 
-            self._proj_stash[i_type] = (Ps, T_col_inds, T_ind_ptr)
-        Ps, T_col_inds, T_ind_ptr = self._proj_stash[i_type]
+            self._proj_stash[key] = (Ps, T_col_inds, T_ind_ptr)
+        Ps, T_col_inds, T_ind_ptr = self._proj_stash[key]
         if return_pointers:
             return Ps, (T_col_inds, T_ind_ptr)
         else:

docs/content/additional_resources.rst

@@ -14,14 +14,14 @@ found on the official websites of the packages
 Python for scientific computing
 -------------------------------
 
-* `Python for Scientists <https://sites.google.com/site/pythonforscientists/>`_ Links to commonly used packages, Matlab to Python comparison
+* `Learn Python <https://pyzo.org/learn.html>`_ Links to commonly used packages, Matlab to Python comparison
 * `Python Wiki <http://wiki.python.org/moin/NumericAndScientific>`_ Lists packages and resources for scientific computing in Python
 
 Numpy and Matlab
 ----------------
 
 * `NumPy for Matlab Users <https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html>`_
-* `Python vs Matlab <https://sites.google.com/site/pythonforscientists/python-vs-matlab>`_
+* `Python vs Matlab <https://pyzo.org/python_vs_matlab.html>`_
 
 Lessons in Python
 -----------------

tests/simplex/test_operators.py

@@ -1,5 +1,8 @@
 import numpy as np
+import pytest
+
 import discretize
+from discretize import SimplexMesh
 from discretize.utils import example_simplex_mesh
 
 
@@ -179,3 +182,16 @@ def test_grad_order(self):
         self.name = "SimplexMesh grad order test"
         self._test_type = "Grad"
         self.orderTest()
+
+
+@pytest.mark.parametrize("i_type", ["E", "F"])
+def test_simplex_projection_caching(i_type):
+    n = 5
+    mesh = SimplexMesh(*example_simplex_mesh((n, n)))
+    P1 = mesh._SimplexMesh__get_inner_product_projection_matrices(
+        i_type, with_volume=False, return_pointers=False
+    )
+    P2 = mesh._SimplexMesh__get_inner_product_projection_matrices(
+        i_type, with_volume=True, return_pointers=False
+    )
+    assert P1 is not P2