Add VTK Matrix/Basis Generators, Multilinear Forms, and High-Order Material Matrices

.gitignore

@@ -1,3 +1,5 @@
 __pycache__
 
-*.lock
+# don't commit the lock; it may get outdated to quickly
+poetry.lock
+

README.md

@@ -23,3 +23,16 @@ If we use basis functions for order $\mathcal{O}$, we have $\frac{1}{6} \times \
 | $M_{kl}$           | $\int_T \Psi_k \Psi_l dx$                     | `M3`               |
 |                    | $\int_T \Phi_k \Phi_l dx$                     | `M2`               |
 | $F_{kl}^{-,j}$     | $\int_T \Psi_k \Psi_l dx$                     | `rT`               |
+
+## Usage
+
+Use poetry to install or use the package. In particular, consider the commands
+
+```
+# install
+poetry install
+# format files
+poetry run black src
+# run unit tests
+poetry run coverage run -m nose2
+```

src/seissol_matrices/base.py

@@ -0,0 +1,22 @@
+import numpy as np
+
+
+def collocate(basis, points):
+    # takes a basis object, and a points array
+    # of the form: npoints × dim
+    # outputs a collocation matrix of the form
+    # npoints × nbasis
+
+    # as far as the author of these lines is aware,
+    # this module has no broadcasing functionality yet
+    # for basis functions. Which is sad.
+
+    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)
+    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
@@ -85,19 +85,25 @@ def eval_diff_basis(self, x, i, k):
     def number_of_basis_functions(self, o):
         pass
 
+    def dim(self):
+        pass
+
 
 ################################################################################
 
 
 class BasisFunctionGenerator1D(BasisFunctionGenerator):
+    def dim(self):
+        return 1
+
     def eval_basis(self, x, i):
         r_num = 2 * x - 1.0
         r_den = 1.0
         return singularity_free_jacobi_polynomial(i, 0, 0, r_num, r_den)
 
     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
@@ -113,6 +119,9 @@ def number_of_basis_functions(self):
 
 
 class BasisFunctionGenerator2D(BasisFunctionGenerator):
+    def dim(self):
+        return 2
+
     def unroll_index(self, i):
         n = i[0] + i[1]
         tri = 0.5 * n * (n + 1)
@@ -171,6 +180,9 @@ def number_of_basis_functions(self):
 
 
 class BasisFunctionGenerator3D(BasisFunctionGenerator):
+    def dim(self):
+        return 3
+
     def unroll_index(self, i):
         n = i[0] + i[1] + i[2]
         tet = (n * (n + 1) * (n + 2)) / 6.0