Add IP matrices

src/seissol_matrices/plasticity/stroud.py

@@ -0,0 +1,41 @@
+import sys
+
+sys.path.append("../..")
+from quad_rules import GaussJacobi
+import numpy as np
+
+np.set_printoptions(precision=15, floatmode="fixed")
+
+
+def stroud(n):
+    points_0, weights_0 = GaussJacobi.GaussJacobi(0, 0).compute_nodes_and_weights(n)
+    points_1, weights_1 = GaussJacobi.GaussJacobi(1, 0).compute_nodes_and_weights(n)
+    points_2, weights_2 = GaussJacobi.GaussJacobi(2, 0).compute_nodes_and_weights(n)
+
+    x = np.zeros((n**3, 3))
+    w = np.zeros(n**3)
+    for i, index in enumerate(np.ndindex(n, n, n)):
+        x[i, :] = np.array([points_2[index[0]], points_1[index[1]], points_0[index[2]]])
+        w[i] = (
+            weights_2[index[0]] * weights_1[index[1]] * weights_0[index[2]] * 0.5**6
+        )
+    x = 0.5 * x + 0.5
+
+    y = np.zeros((n**3, 3))
+    y[:, 0] = x[:, 0]
+    y[:, 1] = x[:, 1] * (1 - y[:, 0])
+    y[:, 2] = x[:, 2] * (1 - y[:, 0] - y[:, 1])
+
+    return y, w
+
+
+if __name__ == "__main__":
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("n", type=int)
+    args = parser.parse_args()
+
+    y, w = stroud(args.n)
+    for p in y:
+        print(f"{p[0]} {p[1]} {p[2]}")

src/seissol_matrices/plasticity_matrices.py

@@ -0,0 +1,66 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import os
+from seissol_matrices import basis_functions, dg_matrices, xml_io
+from seissol_matrices.plasticity import stroud
+
+
+def parse_nodes(filename):
+    with open(filename) as f:
+        lines = [[float(e) for e in l.rstrip().split()] for l in f]
+        return np.array(lines)
+
+
+class PlasticityGenerator:
+    def __init__(self, order):
+        self.bf3_generator = basis_functions.BasisFunctionGenerator3D(order)
+        self.dg3_generator = dg_matrices.dg_generator(order, d=3)
+        self.order = order
+
+    def generate_Vandermonde(self, mode):
+        assert mode in ["ip", "nb"]
+        if mode == "nb":
+            nodes = parse_nodes(
+                f"{os.path.dirname(__file__)}/plasticity/nb_{self.order}.txt"
+            )
+        else:
+            nodes, _ = stroud.stroud(self.order + 1)
+        m = self.bf3_generator.number_of_basis_functions()
+        n = nodes.shape[0]
+
+        vandermonde = np.zeros((n, m))
+
+        for i in range(n):
+            for j in range(m):
+                node = nodes[i]
+                vandermonde[i, j] = self.bf3_generator.eval_basis(
+                    [node[0], node[1], node[2]], j
+                )
+
+        return vandermonde
+
+    def generate_Vandermonde_inv(self, mode):
+        vandermonde = self.generate_Vandermonde(mode)
+        if mode == "nb":
+            vandermonde_inv = np.linalg.solve(vandermonde, np.eye(vandermonde.shape[0]))
+        else:
+            mass = self.dg3_generator.mass_matrix()
+            nodes, weights = stroud.stroud(self.order + 1)
+            vandermonde_inv = np.zeros(vandermonde.T.shape)
+            for l in range(vandermonde_inv.shape[0]):
+                for i in range(vandermonde_inv.shape[1]):
+                    vandermonde_inv[l, i] = vandermonde[i, l] * weights[i] / mass[l, l]
+
+        return vandermonde_inv
+
+
+if __name__ == "__main__":
+    for mode in ["nb", "ip"]:
+        for order in range(2, 8):
+            generator = PlasticityGenerator(order)
+            vandermonde = generator.generate_Vandermonde(mode)
+            vandermonde_inv = generator.generate_Vandermonde_inv(mode)
+            filename = f"output/plasticity_{mode}_matrices_{order}.xml"
+            xml_io.write_matrix(vandermonde, "v", filename)
+            xml_io.write_matrix(vandermonde_inv, "vInv", filename)

tests/test_plasticity.py

@@ -5,67 +5,113 @@
 import os
 import requests
 
-from seissol_matrices import plasticity
+from seissol_matrices import plasticity_matrices
 from seissol_matrices import xml_io
 from seissol_matrices import quad_points
 from tests import helper
 
 
+def filename(mode, order):
+    return f"plasticity_{mode}_{order}.xml"
+
+
 class abstract_tester(object):
-    def test_nb(self):
-        vandermonde = self.generator.generate_Vandermonde("nb")
+    def test(self):
+        vandermonde = self.generator.generate_Vandermonde(self.mode)
         vandermonde_from_file = xml_io.read_matrix("v", self.filename)
         self.compare_matrices(vandermonde, vandermonde_from_file)
 
-        vandermonde_inv = self.generator.generate_Vandermonde_inv("nb")
+        vandermonde_inv = self.generator.generate_Vandermonde_inv(self.mode)
         vandermonde_inv_from_file = xml_io.read_matrix("vInv", self.filename)
         self.compare_matrices(vandermonde_inv, vandermonde_inv_from_file)
 
 
-def setUpClassFromOrder(cls, order):
-    cls.filename = f"plasticity_nb_{order}.xml"
+def downlaod_matrix(mode, order):
     # point to some old commit
-    url = f"https://raw.githubusercontent.com/SeisSol/SeisSol/c99ca8f814b7184eb435272f9d4f63b03b8b6cf4/generated_code/matrices/plasticity_nb_matrices_{order}.xml"
+    url = f"https://raw.githubusercontent.com/SeisSol/SeisSol/c99ca8f814b7184eb435272f9d4f63b03b8b6cf4/generated_code/matrices/plasticity_{mode}_matrices_{order}.xml"
     r = requests.get(url, allow_redirects=True)
-    open(cls.filename, "wb").write(r.content)
+    fn = filename(mode, order)
+    open(fn, "wb").write(r.content)
+
+
+def setUpClassFromOrder(cls, mode, order):
+    downlaod_matrix(mode, order)
+    cls.filename = filename(mode, order)
+    cls.mode = mode
     cls.order = order
-    cls.generator = plasticity.PlasticityGenerator(cls.order)
+    cls.generator = plasticity_matrices.PlasticityGenerator(cls.order)
+
+
+class test_nb_2(abstract_tester, helper.helper):
+    @classmethod
+    def setUpClass(cls):
+        setUpClassFromOrder(cls, "nb", 2)
+
+
+class test_nb_3(abstract_tester, helper.helper):
+    @classmethod
+    def setUpClass(cls):
+        setUpClassFromOrder(cls, "nb", 3)
+
+
+class test_nb_4(abstract_tester, helper.helper):
+    @classmethod
+    def setUpClass(cls):
+        setUpClassFromOrder(cls, "nb", 4)
+
+
+class test_nb_5(abstract_tester, helper.helper):
+    @classmethod
+    def setUpClass(cls):
+        setUpClassFromOrder(cls, "nb", 5)
+
+
+class test_nb_6(abstract_tester, helper.helper):
+    @classmethod
+    def setUpClass(cls):
+        setUpClassFromOrder(cls, "nb", 6)
+
+
+class test_nb_7(abstract_tester, helper.helper):
+    @classmethod
+    def setUpClass(cls):
+        setUpClassFromOrder(cls, "nb", 7)
 
 
-class test_plasticity_2(abstract_tester, helper.helper):
+class test_ip_2(abstract_tester, helper.helper):
     @classmethod
     def setUpClass(cls):
-        setUpClassFromOrder(cls, 2)
+        setUpClassFromOrder(cls, "ip", 2)
 
 
-class test_plasticity_3(abstract_tester, helper.helper):
+class test_ip_3(abstract_tester, helper.helper):
     @classmethod
     def setUpClass(cls):
-        setUpClassFromOrder(cls, 3)
+        setUpClassFromOrder(cls, "ip", 3)
 
 
-class test_plasticity_4(abstract_tester, helper.helper):
+class test_ip_4(abstract_tester, helper.helper):
     @classmethod
     def setUpClass(cls):
-        setUpClassFromOrder(cls, 4)
+        setUpClassFromOrder(cls, "ip", 4)
 
 
-class test_plasticity_5(abstract_tester, helper.helper):
+class test_ip_5(abstract_tester, helper.helper):
     @classmethod
     def setUpClass(cls):
-        setUpClassFromOrder(cls, 5)
+        setUpClassFromOrder(cls, "ip", 5)
 
 
-class test_plasticity_6(abstract_tester, helper.helper):
+class test_ip_6(abstract_tester, helper.helper):
     @classmethod
     def setUpClass(cls):
-        setUpClassFromOrder(cls, 6)
+        setUpClassFromOrder(cls, "ip", 6)
 
 
-class test_plasticity_7(abstract_tester, helper.helper):
+class test_ip_7(abstract_tester, helper.helper):
     @classmethod
     def setUpClass(cls):
-        setUpClassFromOrder(cls, 7)
+        setUpClassFromOrder(cls, "ip", 7)
 
 
 if __name__ == "__main__":