Formatting and bugfixing

src/seissol_matrices/base.py

@@ -1,5 +1,6 @@
 import numpy as np
 
+
 def collocate(basis, points):
     # takes a basis object, and a points array
     # of the form: npoints × dim
@@ -13,9 +14,9 @@ def collocate(basis, points):
     assert basis.dim() == points.shape[1]
 
     nbasis = basis.number_of_basis_functions()
-    
+
     coll = np.empty((nbasis, points.shape[0]))
     for i in range(nbasis):
         for j in range(points.shape[0]):
-            coll[i,j] = basis.eval_basis(points[j,:], i)
+            coll[i, j] = basis.eval_basis(points[j, :], i)
     return coll

src/seissol_matrices/basis_functions.py

@@ -40,7 +40,7 @@ def singularity_free_jacobi_polynomial(n, a, b, x, y):
 
 def singularity_free_jacobi_polynomial_and_derivatives(n, a, b, x, y):
     if n == 0:
-        return 1.0, 0.0, 0.0
+        return np.ones(np.shape(x)), np.zeros(np.shape(x)), np.zeros(np.shape(x))
     pm_1 = 1.0
     ddx_pm_1 = 0.0
     ddy_pm_1 = 0.0
@@ -103,7 +103,7 @@ def eval_basis(self, x, i):
 
     def eval_diff_basis(self, x, i, k):
         if i == 0:
-            return 0
+            return np.zeros(x.shape)
         else:
             r_num = 2 * x - 1.0
             r_den = 1.0

src/seissol_matrices/dg_matrices.py

@@ -5,7 +5,9 @@
 import quad_rules.GaussJacobi
 import quad_rules.Quadrature
 
-import base
+import itertools
+
+from seissol_matrices import base
 
 from seissol_matrices import basis_functions
 
@@ -32,29 +34,35 @@ def __init__(self, o, d):
                 ]
             )
             self.face_generator = dg_generator(o, 2)
-            self.quadrule_finder = quad_rules.JaskowiecSukumar.JaskowiecSukumar().find_best_rule
+            self.quadrule_finder = (
+                quad_rules.JaskowiecSukumar.JaskowiecSukumar().find_best_rule
+            )
             self.generator_finder = basis_functions.BasisFunctionGenerator3D
         elif self.dim == 2:
             self.generator = basis_functions.BasisFunctionGenerator2D(self.order)
             self.geometry = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
-            self.quadrule_finder = quad_rules.WitherdenVincentTri.WitherdenVincentTri().find_best_rule
+            self.quadrule_finder = (
+                quad_rules.WitherdenVincentTri.WitherdenVincentTri().find_best_rule
+            )
             self.generator_finder = basis_functions.BasisFunctionGenerator2D
         elif self.dim == 1:
             self.generator = basis_functions.BasisFunctionGenerator1D(self.order)
             self.geometry = np.array([[0.0], [1.0]])
-            self.quadrule_finder = quad_rules.GaussJacobi.GaussJacobi(0, 0).find_best_rule
+            self.quadrule_finder = quad_rules.GaussJacobi.GaussJacobi(
+                0, 0
+            ).find_best_rule
             self.generator_finder = basis_functions.BasisFunctionGenerator1D
         else:
             raise Execption("Can only generate 1D, 2D or 2D basis functions")
 
         self.nodes, self.weights = self.get_quadrature_rule(2 * self.order)
 
     def get_quadrature_rule(self, order):
-        n, w = self.quadrule_finder(order)        
+        n, w = self.quadrule_finder(order)
 
         return quad_rules.Quadrature.transform(n, w, self.geometry)
 
-    def multilinear_form(self, face, der, order = None, side = None):
+    def multilinear_form(self, face, der, order=None, side=None):
         """
         Given an array der, we compute a tensor V[i...] of dimension
         len(der) defined by (given in C++ parameter pack notation)
@@ -78,14 +86,14 @@ def multilinear_form(self, face, der, order = None, side = None):
         # instead of the two current ones.
 
         # we only support a first derivative at the moment
-        assert np.all([d in (None,*(i for i in range(self.dim))) for d in der])
+        assert np.all([d in (None, *(i for i in range(self.dim))) for d in der])
 
         if order is None:
             order = [self.order] * len(der)
-        
+
         if side is None:
-            side = [-1] * len(der)
-        
+            side = [None] * len(der)
+
         assert len(order) == len(der)
         assert len(side) == len(der)
 
@@ -94,47 +102,58 @@ def get_basis_function(o, d, s):
                 basis = self.generator_finder(o)
             else:
                 basis = self.face_generator.generator_finder(o)
-            basiseval_pre = basis.eval_basis if d is None else lambda x,i: basis.eval_diff_basis(x,i,d)
+            basiseval_pre = (
+                basis.eval_basis
+                if d is None
+                else lambda x, i: basis.eval_diff_basis(x, i, d)
+            )
             if face:
                 if s is None:
-                    basiseval = lambda x, i: basiseval_pre(self.volume_to_face_parametrisation(x, face), i)
+                    basiseval = lambda x, i: basiseval_pre(
+                        self.volume_to_face_parametrisation(x, face), i
+                    )
                 elif s < 0:
                     basiseval = basiseval_pre
                 else:
-                    basiseval = lambda x, i: basiseval_pre(self.face_to_face_parametrisation(x, face), i)
+                    basiseval = lambda x, i: basiseval_pre(
+                        self.face_to_face_parametrisation(x, face), i
+                    )
             else:
                 # volume * face does not make sense when evaluating on the whole volume
                 assert s is None
                 basiseval = basiseval_pre
             return basiseval, basis.number_of_basis_functions()
 
-        basis_functions = [get_basis_function(o,d,s) for o,d,s in zip(order, der, side)]
-
+        basis_functions = [
+            get_basis_function(o, d, s) for o, d, s in zip(order, der, side)
+        ]
+
         if face:
             nodes, weights = self.face_generator.get_quadrature_rule(np.prod(order))
         else:
             nodes, weights = self.get_quadrature_rule(np.prod(order))
-        sizes = [bn for _,bn in basis_functions]
+
+        sizes = [bn for _, bn in basis_functions]
         tensor = np.empty(sizes)
         generic_eval_basis = lambda x, index: np.prod(
-            basis_functions[k][0](x, i) for k,i in enumerate(index)
+            [basis_functions[k][0](x, i) for k, i in enumerate(index)], axis=0
         )
-        for index in itertools.product(range(size) for size in sizes):
+        for index in itertools.product(*[range(size) for size in sizes]):
             eval_basis = lambda x: generic_eval_basis(x, index)
-            tensor[*index] = quad_rules.Quadrature.quad(nodes, weights, eval_basis)
+            tensor[index] = quad_rules.Quadrature.quad(nodes, weights, eval_basis)
         return tensor
 
     def mass_matrix(self):
         if not np.any(self.M == None):
             return self.M
-        
+
         self.M = self.multilinear_form(False, [None, None])
         return self.M
 
     def stiffness_matrix(self, dim):
         if not np.any(self.K[dim] == None):
             return self.K[dim]
-        
+
         self.K[dim] = self.multilinear_form(False, [dim, None])
         return self.K[dim]
 
@@ -215,17 +234,18 @@ def rDivM(self, side):
         matrix = self.rT(side)
         mass = self.mass_matrix()
         return np.linalg.solve(mass, matrix.T)
-    
+
     def collocate_volume(self, points):
         return base.collocate(self.generator, points)
-    
+
     def collocate_face(self, points, side):
         # points are meant to be 2D here
 
         # this method wants dim × npoints; but points is given the other way. So, we transpose it twice
         projected = self.volume_to_face_parametrisation(points.T, side).T
         return base.collocate(self.generator, projected)
 
+
 if __name__ == "__main__":
     from seissol_matrices import json_io
 

src/seissol_matrices/generate.py

@@ -3,25 +3,34 @@
 from dg_matrices import dg_generator
 from vtk_points import vtk_lagrange_2d, vtk_lagrange_3d
 
+
 def main():
     assert len(sys.argv) > 1
-    if sys.argv[1] == 'vtk':
-        import dg_matrices 
+    if sys.argv[1] == "vtk":
+        import dg_matrices
+
         vtkpoints2 = {}
         vtkpoints3 = {}
         for deg in range(0, 8):
             vtkpoints2[deg] = vtk_lagrange_2d(deg)
             vtkpoints3[deg] = vtk_lagrange_3d(deg)
-            json_io.write_matrix(vtkpoints2[deg], f'vtk2d{deg}', f'vtkbase.json')
-            json_io.write_matrix(vtkpoints3[deg], f'vtk3d{deg}', f'vtkbase.json')
+            json_io.write_matrix(vtkpoints2[deg], f"vtk2d{deg}", f"vtkbase.json")
+            json_io.write_matrix(vtkpoints3[deg], f"vtk3d{deg}", f"vtkbase.json")
         for basisorder in range(2, 8):
             dggen = dg_generator(basisorder, 3)
             for deg in range(0, 8):
-                json_io.write_matrix(dggen.collocate_volume(vtkpoints3[deg]),
-                    f'coll{basisorder}vd{deg}', f'vtko{basisorder}.json')
+                json_io.write_matrix(
+                    dggen.collocate_volume(vtkpoints3[deg]),
+                    f"coll{basisorder}vd{deg}",
+                    f"vtko{basisorder}.json",
+                )
                 for f in range(4):
-                    json_io.write_matrix(dggen.collocate_face(vtkpoints2[deg], f),
-                        f'coll{basisorder}f{f}d{deg}', f'vtko{basisorder}.json')
+                    json_io.write_matrix(
+                        dggen.collocate_face(vtkpoints2[deg], f),
+                        f"coll{basisorder}f{f}d{deg}",
+                        f"vtko{basisorder}.json",
+                    )
+
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()

src/seissol_matrices/vtk_points.py

@@ -1,79 +1,109 @@
 import numpy as np
 
+
 def vtk_lagrange_2d(d, start=0, divd=None):
     if d < 0:
         return np.empty((0, 2))
     if d == 0:
-        return np.array([[1/3, 1/3]])
+        return np.array([[1 / 3, 1 / 3]])
     else:
         if divd is None:
             divd = d
         s = 1 / divd
         base = 0
         obase = 0
         allrings = []
-        for r in range(start, d//3 + 1):
-            if d == 3*r:
-                allrings += [[s*r, s*r]]
+        for r in range(start, d // 3 + 1):
+            if d == 3 * r:
+                allrings += [[s * r, s * r]]
             else:
-                ringcorner = [[s*r, s*r], [s*(d - 2*r), s*r], [s*r, s*(d - 2*r)]]
-                ringborder = [[s*(r + i), s*r] for i in range(1, d - 3*r)]
-                ringborder += [[s*(d - 2*r - i), s*(r + i)] for i in range(1, d - 3*r)]
-                ringborder += [[s*r, s*(d - 2*r - i)] for i in range(1, d - 3*r)]
+                ringcorner = [
+                    [s * r, s * r],
+                    [s * (d - 2 * r), s * r],
+                    [s * r, s * (d - 2 * r)],
+                ]
+                ringborder = [[s * (r + i), s * r] for i in range(1, d - 3 * r)]
+                ringborder += [
+                    [s * (d - 2 * r - i), s * (r + i)] for i in range(1, d - 3 * r)
+                ]
+                ringborder += [
+                    [s * r, s * (d - 2 * r - i)] for i in range(1, d - 3 * r)
+                ]
                 allrings += ringcorner + ringborder
         return np.array(allrings)
 
+
 def vtk_lagrange_3d(d):
     if d == 0:
-        return np.array([[1/4, 1/4, 1/4]])
+        return np.array([[1 / 4, 1 / 4, 1 / 4]])
     else:
         s = 1 / d
         base = 0
         obase = 0
         allrings = []
-        for r in range(d//4 + 1):
-            if d == 4*r:
-                allrings += [[s*r, s*r, s*r]]
+        for r in range(d // 4 + 1):
+            if d == 4 * r:
+                allrings += [[s * r, s * r, s * r]]
             else:
-                ringcorner = [[s*r, s*r, s*r], [s*(d - 3*r), s*r, s*r], [s*r, s*(d - 3*r), s*r], [s*r, s*r, s*(d - 3*r)]]
-                ringborder = [[s*(r + i), s*r, s*r] for i in range(1, d - 4*r)]
-                ringborder += [[s*(d - 3*r - i), s*(r + i), s*r] for i in range(1, d - 4*r)]
-                ringborder += [[s*r, s*(d - 3*r - i), s*r] for i in range(1, d - 4*r)]
-                ringborder += [[s*r, s*r, s*(r + i)] for i in range(1, d - 4*r)]
-                ringborder += [[s*(d - 3*r - i), s*r, s*(r + i)] for i in range(1, d - 4*r)]
-                ringborder += [[s*r, s*(d - 3*r - i), s*(r + i)] for i in range(1, d - 4*r)]
+                ringcorner = [
+                    [s * r, s * r, s * r],
+                    [s * (d - 3 * r), s * r, s * r],
+                    [s * r, s * (d - 3 * r), s * r],
+                    [s * r, s * r, s * (d - 3 * r)],
+                ]
+                ringborder = [[s * (r + i), s * r, s * r] for i in range(1, d - 4 * r)]
+                ringborder += [
+                    [s * (d - 3 * r - i), s * (r + i), s * r]
+                    for i in range(1, d - 4 * r)
+                ]
+                ringborder += [
+                    [s * r, s * (d - 3 * r - i), s * r] for i in range(1, d - 4 * r)
+                ]
+                ringborder += [[s * r, s * r, s * (r + i)] for i in range(1, d - 4 * r)]
+                ringborder += [
+                    [s * (d - 3 * r - i), s * r, s * (r + i)]
+                    for i in range(1, d - 4 * r)
+                ]
+                ringborder += [
+                    [s * r, s * (d - 3 * r - i), s * (r + i)]
+                    for i in range(1, d - 4 * r)
+                ]
                 subface = vtk_lagrange_2d(d - r, r + 1, d)
-                ringfaces = [(i, s*r, j) for i,j in subface]
-                ringfaces += [(j, (1 - s*r) - i - j, i) for i,j in subface]
-                ringfaces += [(s*r, j, i) for i,j in subface]
-                ringfaces += [(j, i, s*r) for i,j in subface]
+                ringfaces = [(i, s * r, j) for i, j in subface]
+                ringfaces += [(j, (1 - s * r) - i - j, i) for i, j in subface]
+                ringfaces += [(s * r, j, i) for i, j in subface]
+                ringfaces += [(j, i, s * r) for i, j in subface]
                 allrings += ringcorner + ringborder + ringfaces
         return np.array(allrings)
 
+
 def vtk_lagrange_2d_from_vtk(d):
     import vtk
+
     tet = vtk.vtkLagrangeTriangle()
-    points = ((d+1)*(d+2))//2
+    points = ((d + 1) * (d + 2)) // 2
     tet.GetPointIds().SetNumberOfIds(points)
     tet.GetPoints().SetNumberOfPoints(points)
     tet.Initialize()
     pointlist = []
     for i in range(points):
-        barycentric = [0,0,0]
+        barycentric = [0, 0, 0]
         tet.ToBarycentricIndex(i, barycentric)
         pointlist += [[b / d for b in barycentric[:2]]]
     return np.array(pointlist)
 
+
 def vtk_lagrange_3d_from_vtk(d):
     import vtk
+
     tet = vtk.vtkLagrangeTetra()
-    points = ((d+1)*(d+2)*(d+3))//6
+    points = ((d + 1) * (d + 2) * (d + 3)) // 6
     tet.GetPointIds().SetNumberOfIds(points)
     tet.GetPoints().SetNumberOfPoints(points)
     tet.Initialize()
     pointlist = []
     for i in range(points):
-        barycentric = [0,0,0,0]
+        barycentric = [0, 0, 0, 0]
         tet.ToBarycentricIndex(i, barycentric)
         pointlist += [[b / d for b in barycentric[:3]]]
     return np.array(pointlist)