Put the material mode index last

src/seissol_matrices/base.py

@@ -1,6 +1,19 @@
 import numpy as np
 
 
+def material_index_last(tensor):
+    """
+    Move the material mode index of a tensor from the leading to the
+    trailing axis.
+
+    The multilinear forms yield the material basis as the leading axis.
+    yateto assigns stride one to the leading index and aligns it to the
+    vector width, so the material index belongs at the trailing, slowest
+    varying end instead.
+    """
+    return np.moveaxis(tensor, 0, -1)
+
+
 def collocate(basis, points):
     # takes a basis object, and a points array
     # of the form: npoints × dim

src/seissol_matrices/dg_matrices.py

@@ -165,27 +165,29 @@ def kDivM(self, dim, matorder=None):
         mass = self.mass_matrix()
         if matorder is None:
             stiffness = self.stiffness_matrix(dim)
-        else:
-            stiffness = self.multilinear_form(
-                False, [None, dim, None], order=[matorder, self.order, self.order]
-            )
-        return np.linalg.solve(mass, stiffness)
+            return np.linalg.solve(mass, stiffness)
+
+        stiffness = self.multilinear_form(
+            False, [None, dim, None], order=[matorder, self.order, self.order]
+        )
+        return base.material_index_last(np.linalg.solve(mass, stiffness))
 
     def kDivMT(self, dim, matorder=None):
         mass = self.mass_matrix()
         if matorder is None:
             stiffness = self.stiffness_matrix(dim)
-            sT = stiffness.T
-        else:
-            correct1 = self.multilinear_form(
-                False, [dim, None, None], order=[matorder, self.order, self.order]
-            )
-            correct2 = self.multilinear_form(
-                False, [None, None, dim], order=[matorder, self.order, self.order]
-            )
-            sT = correct1 + correct2
+            return np.linalg.solve(mass, stiffness.T)
 
-        return np.linalg.solve(mass, sT)
+        # the material sits inside the derivative, so the weak form transpose
+        # picks up the term in which the material basis is differentiated
+        differentiated_material = self.multilinear_form(
+            False, [dim, None, None], order=[matorder, self.order, self.order]
+        )
+        differentiated_basis = self.multilinear_form(
+            False, [None, None, dim], order=[matorder, self.order, self.order]
+        )
+        sT = differentiated_material + differentiated_basis
+        return base.material_index_last(np.linalg.solve(mass, sT))
 
     def face_to_face_parametrisation(self, x, side):
         # implement Dumbser & Käser, 2006 Table 2 b)
@@ -252,17 +254,18 @@ def rDivM(self, side, matorder=None):
         mass = self.mass_matrix()
         if matorder is None:
             matrix = self.rT(side).T
-        else:
-            prematrix = self.multilinear_form(
-                True,
-                [None, None, None],
-                order=[matorder, self.order, self.order],
-                side=[None, None, side],
-                dim=[3, 2, 3],
-            )
-            facemass = self.face_generator.mass_matrix()
-            matrix = np.linalg.solve(facemass, prematrix).transpose((0, 2, 1))
-        return np.linalg.solve(mass, matrix)
+            return np.linalg.solve(mass, matrix)
+
+        prematrix = self.multilinear_form(
+            True,
+            [None, None, None],
+            order=[matorder, self.order, self.order],
+            side=[None, None, side],
+            dim=[3, 2, 3],
+        )
+        facemass = self.face_generator.mass_matrix()
+        matrix = np.linalg.solve(facemass, prematrix).transpose((0, 2, 1))
+        return base.material_index_last(np.linalg.solve(mass, matrix))
 
     def collocate_volume(self, points):
         return base.collocate(self.generator, points)

src/seissol_matrices/dr_matrices.py

@@ -1,6 +1,7 @@
 #!/usr/bin/env python
 
 import numpy as np
+from seissol_matrices import base
 from seissol_matrices import basis_functions
 from seissol_matrices import dg_matrices
 from seissol_matrices import quad_points
@@ -79,9 +80,11 @@ def V3mTo2nTWDivM(self, a, b, matorder=None):
         for i in range(n):
             W[i, i] = weights[i]
 
-        matrixT = matrix.T if matorder is None else matrix.transpose((0, 2, 1))
+        if matorder is None:
+            return np.linalg.solve(mass, np.dot(matrix.T, W))
 
-        return np.linalg.solve(mass, np.dot(matrixT, W))
+        matrixT = matrix.transpose((0, 2, 1))
+        return base.material_index_last(np.linalg.solve(mass, np.dot(matrixT, W)))
 
     def quadpoints(self):
         points = self.quadrule.points()

tests/test_material_order.py

@@ -0,0 +1,74 @@
+#!/usr/bin/env python3
+
+import unittest
+import numpy as np
+
+from seissol_matrices import dg_matrices
+from seissol_matrices import dr_matrices
+from seissol_matrices import quad_points
+
+
+class abstract_tester(object):
+    def compare(self, a, b):
+        self.assertEqual(a.shape, b.shape)
+        np.testing.assert_allclose(a, b, rtol=0.0, atol=1.2e-11)
+
+    def test_kDivM_material_index_last(self):
+        for dim in range(3):
+            matrix = self.generator.kDivM(dim, 2)
+            self.assertEqual(matrix.shape, (self.nbf, self.nbf, 4))
+
+    def test_kDivM_constant_material(self):
+        for dim in range(3):
+            self.compare(
+                self.generator.kDivM(dim, 1)[:, :, 0], self.generator.kDivM(dim)
+            )
+
+    def test_kDivMT_constant_material(self):
+        for dim in range(3):
+            self.compare(
+                self.generator.kDivMT(dim, 1)[:, :, 0], self.generator.kDivMT(dim)
+            )
+
+    def test_rDivM_constant_material(self):
+        for side in range(4):
+            self.compare(
+                self.generator.rDivM(side, 1)[:, :, 0], self.generator.rDivM(side)
+            )
+
+    def test_V3mTo2nTWDivM_constant_material(self):
+        for a in range(4):
+            for b in range(4):
+                self.compare(
+                    self.dr_generator.V3mTo2nTWDivM(a, b, 1)[:, :, 0],
+                    self.dr_generator.V3mTo2nTWDivM(a, b),
+                )
+
+
+def setUpClassFromOrder(cls, order):
+    cls.order = order
+    cls.nbf = order * (order + 1) * (order + 2) // 6
+    cls.generator = dg_matrices.dg_generator(order, 3)
+    cls.dr_generator = dr_matrices.dr_generator(order, quad_points.stroud(order + 1))
+
+
+class test_material_order_2(abstract_tester, unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        setUpClassFromOrder(cls, 2)
+
+
+class test_material_order_3(abstract_tester, unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        setUpClassFromOrder(cls, 3)
+
+
+class test_material_order_4(abstract_tester, unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        setUpClassFromOrder(cls, 4)
+
+
+if __name__ == "__main__":
+    unittest.main()