Add some collocation matrices

src/seissol_matrices/base.py

@@ -1,7 +1,7 @@
 import numpy as np
 
 
-def collocate(basis, points):
+def collocate(basis, points, deriv=None):
     # takes a basis object, and a points array
     # of the form: npoints × dim
     # outputs a collocation matrix of the form
@@ -15,8 +15,12 @@ def collocate(basis, points):
 
     nbasis = basis.number_of_basis_functions()
 
+    evalfun = lambda j, i: basis.eval_basis(points[j, :], i)
+    if deriv is not None:
+        evalfun = lambda j, i: basis.eval_diff_basis(points[j, :], i, deriv)
+
     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] = evalfun(j, i)
     return coll

src/seissol_matrices/dg_matrices.py

@@ -264,15 +264,15 @@ def rDivM(self, side, matorder=None):
             matrix = np.linalg.solve(facemass, prematrix).transpose((0, 2, 1))
         return np.linalg.solve(mass, matrix)
 
-    def collocate_volume(self, points):
-        return base.collocate(self.generator, points)
+    def collocate_volume(self, points, deriv=None):
+        return base.collocate(self.generator, points, deriv)
 
-    def collocate_face(self, points, side):
+    def collocate_face(self, points, side, deriv=None):
         # 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)
+        return base.collocate(self.generator, projected, deriv)
 
 
 if __name__ == "__main__":

src/seissol_matrices/generate.py

@@ -5,6 +5,7 @@
 from vtk_points import vtk_lagrange_2d, vtk_lagrange_3d
 from seissol_matrices import quad_points
 from plasticity_matrices import PlasticityGenerator
+import numpy as np
 
 
 def main():
@@ -136,6 +137,69 @@ def main():
                             json_io.write_tensor(
                                 V3mTo2n, f"homV3mTo2n({a},{b})", filename
                             )
+    elif sys.argv[1] == "elementwise":
+        for basisorder in range(2, 9):
+            dggen = dg_generator(basisorder, 3)
+            json_io.write_tensor(
+                dggen.nodes,
+                f"ew_quad_nodes_vv",
+                f"elemwise-collocate-p{basisorder}.json",
+            )
+            json_io.write_tensor(
+                dggen.weights,
+                f"ew_quad_weights_v",
+                f"elemwise-collocate-p{basisorder}.json",
+            )
+            json_io.write_tensor(
+                dggen.face_generator.nodes,
+                f"ew_quad_nodes_ff",
+                f"elemwise-collocate-p{basisorder}.json",
+            )
+            json_io.write_tensor(
+                dggen.face_generator.weights,
+                f"ew_quad_weights_f",
+                f"elemwise-collocate-p{basisorder}.json",
+            )
+            json_io.write_tensor(
+                dggen.collocate_volume(dggen.nodes).T,
+                f"ew_collocate_f_vv",
+                f"elemwise-collocate-p{basisorder}.json",
+            )
+            json_io.write_tensor(
+                np.stack(
+                    [dggen.collocate_volume(dggen.nodes, d).T for d in range(3)],
+                    axis=-1,
+                ),
+                f"ew_collocate_df_vv",
+                f"elemwise-collocate-p{basisorder}.json",
+            )
+            json_io.write_tensor(
+                dggen.face_generator.collocate_volume(dggen.face_generator.nodes).T,
+                f"ew_collocate_f_ff",
+                f"elemwise-collocate-p{basisorder}.json",
+            )
+            json_io.write_tensor(
+                np.stack(
+                    [
+                        dggen.volume_to_face_parametrisation(dggen.nodes.T, f).T
+                        for f in range(4)
+                    ],
+                    axis=-1,
+                ),
+                f"ew_quad_nodes_vf",
+                f"elemwise-collocate-p{basisorder}.json",
+            )
+            json_io.write_tensor(
+                np.stack(
+                    [
+                        dggen.collocate_face(dggen.face_generator.nodes, f).T
+                        for f in range(4)
+                    ],
+                    axis=-1,
+                ),
+                f"ew_collocate_f_vf",
+                f"elemwise-collocate-p{basisorder}.json",
+            )
 
 
 if __name__ == "__main__":