Bugfix in calculate interval

lbmweights/lattice.py

@@ -130,13 +130,12 @@ class Lattice:
                 # positive weights
                 if i == 0:
                     cur_interval = sp.Union(cur_interval, sp.Interval(-oo, roots[i+1]))
-                elif i == (len(roots) - 1):
+                elif i == (len(roots) - 2):
                     cur_interval = sp.Union(cur_interval, sp.Interval(roots[i], oo))
                 else:
                     cur_interval = sp.Union(cur_interval, sp.Interval(roots[i], roots[i+1]))
             interval = interval.intersect(cur_interval)
-            pass
-        self._interval = interval
+        self._interval = sp.Complement(interval, sp.FiniteSet(0))
         return interval
 
     @timer
@@ -289,7 +288,8 @@ class Lattice:
         logger.info(self._interval.boundary)
         if self._interval == sp.EmptySet():
             return []
-        weights = [x.subs(c_s_sq, self._interval.boundary.inf) for x in weight_polynomials]
+        weights = [sp.N(x.eval(self._interval.boundary.inf)) for x in weight_polynomials]
+        weights = [x if abs(x) >= 1e-10 else 0 for x in weights]
         return weights
 
 

test/test_lattice.py

 import unittest
+import sympy as sp
 from lbmweights import Lattice
 from lbmweights.utils.mylog import logger
 
@@ -41,6 +42,9 @@ class TestV37(unittest.TestCase):
     def testOutput(self):
         weights = self.lattice.calculate_weights()
         self.assertEqual(len(weights), 9)
+        interval = self.lattice._interval
+        self.assertAlmostEqual(sp.N(interval.inf), 0.6979533, places=6)
+        self.assertAlmostEqual(sp.N(interval.sup), 0.8704738, places=6)
 
 
 class TestReducibleShell(unittest.TestCase):
@@ -52,6 +56,9 @@ class TestReducibleShell(unittest.TestCase):
     def testOutput(self):
         weights = self.lattice.calculate_weights()
         self.assertEqual(len(weights), 12)
+        interval = self.lattice._interval
+        self.assertAlmostEqual(sp.N(interval.inf), 0.7592510, places=6)
+        self.assertAlmostEqual(sp.N(interval.sup), 0.9054850, places=6)
 
 
 if __name__ == '__main__':