Merge pull request #5 from eisenforschung/import_tests

Import tests

tests/unit/elasticity/test_elasticity.py

@@ -0,0 +1,109 @@
+import numpy as np
+import unittest
+from elaston.linear_elasticity.elasticity.linear_elasticity import LinearElasticity
+from elaston.linear_elasticity.elasticity import tools
+
+
+def create_random_C(isotropic=False):
+    C11_range = np.array([0.7120697386322292, 1.5435656086034886])
+    coeff_C12 = np.array([0.65797601, -0.0199679])
+    coeff_C44 = np.array([0.72753844, -0.30418746])
+    C = np.zeros((6, 6))
+    C11 = C11_range[0] + np.random.random() * C11_range.ptp()
+    C12 = np.polyval(coeff_C12, C11) + 0.2 * (np.random.random() - 0.5)
+    C44 = np.polyval(coeff_C44, C11) + 0.2 * (np.random.random() - 0.5)
+    C[:3, :3] = C12
+    C[:3, :3] += (C11 - C12) * np.eye(3)
+    if isotropic:
+        C[3:, 3:] = np.eye(3) * (C[0, 0] - C[0, 1]) / 2
+    else:
+        C[3:, 3:] = C44 * np.eye(3)
+    return tools.C_from_voigt(C)
+
+
+class TestElasticity(unittest.TestCase):
+    def test_frame(self):
+        medium = LinearElasticity(np.random.random((6, 6)))
+        self.assertAlmostEqual(np.linalg.det(medium.orientation), 1)
+        medium.orientation = np.random.random((3, 3))
+        self.assertAlmostEqual(np.linalg.det(medium.orientation), 1)
+
+    def test_orientation(self):
+        elastic_tensor = create_random_C()
+        epsilon = np.random.random((3, 3))
+        epsilon += epsilon.T
+        sigma = np.einsum("ijkl,kl->ij", elastic_tensor, epsilon)
+        medium = LinearElasticity(elastic_tensor)
+        medium.orientation = np.array([[1, 1, 1], [1, 0, -1]])
+        sigma = np.einsum("iI,jJ,IJ->ij", medium.orientation, medium.orientation, sigma)
+        sigma_calc = np.einsum(
+            "ijkl,kK,lL,KL->ij",
+            medium.elastic_tensor,
+            medium.orientation,
+            medium.orientation,
+            epsilon,
+        )
+        self.assertTrue(np.allclose(sigma - sigma_calc, 0))
+
+    def test_youngs_modulus(self):
+        medium = LinearElasticity(np.eye(6))
+        self.assertTrue(np.allclose(medium.youngs_modulus, 1))
+
+    def test_poissons_ratio(self):
+        medium = LinearElasticity(np.eye(6))
+        self.assertTrue(np.allclose(medium.poissons_ratio, 0))
+
+    def test_isotropic(self):
+        medium = LinearElasticity(create_random_C(isotropic=True))
+        self.assertTrue(medium._is_isotropic)
+
+    def test_energy(self):
+        elastic_tensor = create_random_C()
+        medium = LinearElasticity(elastic_tensor)
+        r_max = 1e6 * np.random.random() + 10
+        r_min_one = 10 * np.random.random()
+        r_min_two = 10 * np.random.random()
+        E_one = medium.get_dislocation_energy([0, 0, 1], r_min_one, r_max)
+        E_two = medium.get_dislocation_energy([0, 0, 1], r_min_two, r_max)
+        self.assertGreater(E_one, 0)
+        self.assertGreater(E_two, 0)
+        self.assertAlmostEqual(
+            E_one / np.log(r_max / r_min_one), E_two / np.log(r_max / r_min_two)
+        )
+
+    def test_force(self):
+        elastic_tensor = create_random_C()
+        medium = LinearElasticity(elastic_tensor)
+        medium.orientation = [[1, -2, 1], [1, 1, 1], [-1, 0, 1]]
+        lattice_constant = 3.52
+        partial_one = np.array([-0.5, 0, np.sqrt(3) / 2]) * lattice_constant
+        partial_two = np.array([0.5, 0, np.sqrt(3) / 2]) * lattice_constant
+        stress = medium.get_dislocation_stress([0, 10, 0], partial_one)
+        force = medium.get_dislocation_force(stress, [0, 1, 0], partial_two)
+        self.assertAlmostEqual(force[1], 0)
+        self.assertAlmostEqual(force[2], 0)
+
+    def test_elastic_tensor_input(self):
+        C = create_random_C()
+        medium = LinearElasticity([C[0, 0, 0, 0], C[0, 0, 1, 1], C[0, 1, 0, 1]])
+        self.assertTrue(np.allclose(C, medium.elastic_tensor))
+
+    def test_isotropy_tolerance(self):
+        medium = LinearElasticity(create_random_C())
+        with self.assertRaises(ValueError):
+            medium.isotropy_tolerance = -1
+
+    def test_bulk_modulus(self):
+        medium = LinearElasticity(create_random_C())
+        self.assertGreater(medium.bulk_modulus, 0)
+
+    def test_greens_function(self):
+        medium = LinearElasticity(create_random_C(isotropic=True))
+        self.assertAlmostEqual(
+            medium.get_greens_function([1, 1, 1], isotropic=True)[0, 0],
+            medium.get_greens_function([1, 1, 1], isotropic=False)[0, 0]
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/unit/elasticity/test_eshelby.py

@@ -0,0 +1,75 @@
+import numpy as np
+import unittest
+from elaston.linear_elasticity.eshelby import Eshelby
+from elaston.linear_elasticity import tools
+
+
+def create_random_HL(b=None):
+    C = np.zeros((6, 6))
+    C[:3, :3] = np.random.random()
+    C[:3, :3] += np.random.random()*np.eye(3)
+    C[3:, 3:] = np.random.random()*np.eye(3)
+    C = tools.C_from_voigt(C)
+    if b is None:
+        b = np.random.random(3)
+    return Eshelby(C, b)
+
+
+class TestEschelby(unittest.TestCase):
+    def test_p(self):
+        hl = create_random_HL()
+        self.assertTrue(
+            np.allclose(np.absolute(np.linalg.det(hl._get_pmat(hl.p))), 0),
+            'p-matrix has a full dimension'
+        )
+
+    def test_Ak(self):
+        hl = create_random_HL()
+        self.assertTrue(
+            np.allclose(np.absolute(np.einsum('nk,nik->ni', hl.Ak, hl._get_pmat(hl.p))), 0),
+            'Ak is not the kernel of the p-matrix'
+        )
+
+    def test_DAk(self):
+        hl = create_random_HL()
+        self.assertTrue(
+            np.isclose(np.real(np.einsum('nk,n->k', hl.Ak, hl.D)), hl.burgers_vector).all(),
+            'Magnitude not corresponding to the Burgers vector'
+        )
+
+    def test_Dq(self):
+        hl = create_random_HL()
+        F = np.einsum('n,ij->nij', hl.p, hl.elastic_tensor[:, 1, :, 1])
+        F += hl.elastic_tensor[:, 1, :, 0]
+        F = np.einsum('nik,nk->ni', F, hl.Ak)
+        self.assertTrue(
+            np.allclose(np.real(np.einsum('nk,n->k', F, hl.D)), 0),
+            'Equilibrium condition not satisfied'
+        )
+
+    def test_displacement(self):
+        hl = create_random_HL(b=[0, 0, 1])
+        positions = (np.random.random((100, 2))-0.5)*10
+        d_analytical = np.arctan2(*positions.T[:2][::-1])/2/np.pi
+        self.assertTrue(
+            np.all(np.absolute(hl.get_displacement(positions)[:, -1]-d_analytical) < 1.0e-4),
+            'Screw dislocation displacement field not reproduced'
+        )
+
+    def test_strain(self):
+        hl = create_random_HL(b=[0, 0, 1])
+        positions = (np.random.random((100, 2))-0.5)*10
+        strain_analytical = positions[:, 0]/np.sum(positions**2, axis=-1)/4/np.pi
+        self.assertTrue(
+            np.all(np.absolute(hl.get_strain(positions)[:, 1, 2]-strain_analytical) < 1.0e-4),
+            'Screw dislocation strain field (yz-component) not reproduced'
+        )
+        strain_analytical = -positions[:,1]/np.sum(positions**2, axis=-1)/4/np.pi
+        self.assertTrue(
+            np.all(np.absolute(hl.get_strain(positions)[:, 0, 2]-strain_analytical) < 1.0e-4),
+            'Screw dislocation strain field (xz-component) not reproduced'
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/unit/elasticity/test_green.py

@@ -0,0 +1,73 @@
+import numpy as np
+import unittest
+from elaston.linear_elasticity.green import Anisotropic, Isotropic
+from elaston.linear_elasticity import tools
+
+
+def create_random_C(isotropic=False):
+    C11_range = np.array([0.7120697386322292, 1.5435656086034886])
+    coeff_C12 = np.array([0.65797601, -0.0199679])
+    coeff_C44 = np.array([0.72753844, -0.30418746])
+    C = np.zeros((6, 6))
+    C11 = C11_range[0]+np.random.random()*C11_range.ptp()
+    C12 = np.polyval(coeff_C12, C11)+0.2*(np.random.random()-0.5)
+    C44 = np.polyval(coeff_C44, C11)+0.2*(np.random.random()-0.5)
+    C[:3, :3] = C12
+    C[:3, :3] += (C11-C12)*np.eye(3)
+    if isotropic:
+        C[3:, 3:] = np.eye(3)*(C[0, 0]-C[0, 1])/2
+    else:
+        C[3:, 3:] = C44*np.eye(3)
+    return tools.C_from_voigt(C)
+
+
+class TestGreen(unittest.TestCase):
+    def test_derivative(self):
+        aniso = Anisotropic(create_random_C())
+        positions = np.tile(np.random.random(3), 2).reshape(-1, 3)
+        dz = 1.0e-6
+        index = np.random.randint(3)
+        positions[1, index] += dz
+        G_an = aniso.get_greens_function(positions.mean(axis=0), derivative=1)[:, index, :]
+        G_num = np.diff(aniso.get_greens_function(positions), axis=0)/dz
+        self.assertTrue(np.isclose(G_num-G_an, 0).all())
+        G_an = aniso.get_greens_function(positions.mean(axis=0), derivative=2)[:, :, :, index]
+        G_num = np.diff(aniso.get_greens_function(positions, derivative=1), axis=0)/dz
+        self.assertTrue(np.isclose(G_num-G_an, 0).all())
+
+    def test_comp_iso_aniso(self):
+        shear_modulus = 52.5
+        lame_parameter = 101.3
+        poissons_ratio = 0.5 / (1 + shear_modulus / lame_parameter)
+        C_11 = lame_parameter + 2 * shear_modulus
+        C_12 = lame_parameter
+        C_44 = shear_modulus
+        iso = Isotropic(poissons_ratio, shear_modulus)
+        aniso = Anisotropic(tools.C_from_voigt(tools.coeff_to_voigt([C_11, C_12, C_44])))
+        x = np.random.randn(100, 3) * 10
+        for i in range(3):
+            self.assertLess(
+                np.ptp(
+                    iso.get_greens_function(x, derivative=i)
+                    - aniso.get_greens_function(x, derivative=i)
+                ),
+                1e-8,
+                msg=f"Aniso- and isotropic Green's F's give different results for derivative={i}"
+            )
+
+    def test_memory(self):
+        aniso = Anisotropic(create_random_C())
+        positions = np.tile(np.random.random(3), 2).reshape(-1, 3)
+        G_normal = aniso.get_greens_function(positions)
+        G_save = aniso.get_greens_function(positions, save_memory=True)
+        self.assertTrue(np.all(np.isclose(G_normal, G_save)))
+        G_unique = aniso.get_greens_function(positions, check_unique=True)
+        self.assertTrue(np.all(np.isclose(G_normal, G_unique)))
+        G = aniso.get_greens_function(
+            positions, check_unique=True, save_memory=True
+        )
+        self.assertTrue(np.all(np.isclose(G_normal, G)))
+
+
+if __name__ == "__main__":
+    unittest.main()