Merge branch 'parametric-mesher'

.github/workflows/build-and-test.yaml

@@ -20,7 +20,7 @@ jobs:
       - name: Install dependencies
         run: |
           sudo apt-get update
-          pip install --upgrade pip setuptools wheel auditwheel patchelf
+          pip install --upgrade setuptools wheel auditwheel patchelf
       - name: Build the wheel
         run: python setup.py bdist_wheel
       - name: Build the source distribution
@@ -53,7 +53,7 @@ jobs:
       - name: Install dependencies
         shell: powershell
         run: |
-          pip install --upgrade pip setuptools wheel
+          pip install --upgrade setuptools wheel
       - name: Build the wheel
         run: python setup.py bdist_wheel
       - name: Upload wheel artifact
@@ -127,7 +127,6 @@ jobs:
           path: ./
       - name: Install Twine
         run: |
-          pip install --upgrade pip
           pip install twine
       - name: Upload to PyPI
         env:
@@ -152,7 +151,6 @@ jobs:
           path: ./
       - name: Install Twine
         run: |
-          pip install --upgrade pip
           pip install twine
       - name: Upload to PyPI
         env:

.gitignore

@@ -7,3 +7,4 @@ __pycache__
 *.o
 build/*
 dist/*
+tests/world_out.stl

setup.py

@@ -42,7 +42,7 @@
     },
     long_description = open('README.md').read(),
     long_description_content_type='text/markdown',
-    install_requires=['matplotlib', 'vedo', 'numpy', 'gmsh>=4.9', 'pygmsh>=7.1.13', 'scipy'],
+    install_requires=['matplotlib', 'vedo', 'numpy', 'scipy'],
     project_urls = {
         'Documentation': "https://leon.science/traceon",
         'Code': "https://github.com/leon-vv/traceon",

tests/geometric_tests.py

@@ -0,0 +1,194 @@
+import os.path as path
+
+import unittest
+from math import *
+
+import traceon.mesher as M
+from traceon.geometry import *
+import traceon.plotting as P
+
+class MeshTests(unittest.TestCase):
+
+    def test_loading_mesh(self):
+        input_file = path.join(path.dirname(__file__), 'world.stl')
+        output_file = path.join(path.dirname(__file__), 'world_out.stl')
+
+        m = M.Mesh.read_file(input_file)
+        m.write_file(output_file)
+        m2 = M.Mesh.read_file(output_file)
+
+        assert np.allclose(m.points, m2.points)
+        assert np.allclose(m.triangles, m2.triangles)
+
+class PathTests(unittest.TestCase):
+    def test_from_irregular_function(self):
+        f = lambda x: [x, x**(3/2), 0.]
+        y = Path.from_irregular_function(f)
+        assert np.isclose(y.path_length, 1.43970987337155), y.path_length
+
+        # Half the path length is reached at the special point x=0.5662942656295281
+        a = 0.5662942656295281
+        yhalf = Path.from_irregular_function(lambda x: [a*x, (a*x)**(3/2), 0.])
+        assert np.isclose(yhalf.path_length, y.path_length/2)
+
+        assert np.allclose(f(0.), y(0.))
+        assert np.allclose(f(1.), y(y.path_length))
+        assert np.allclose(f(a), y(0.5*y.path_length))
+        assert np.allclose(yhalf.endpoint(), y.middle_point())
+        assert np.allclose(f(1.), y.endpoint())
+
+    def test_spline_through_points(self):
+        points = [
+            [0, 0, 0],
+            [0, 1, 1],
+            [1, 1, 1],
+            [1, 0, 1],
+            [0, 0, 2]]
+
+        path = Path.spline_through_points(points)
+        assert np.allclose(path.starting_point(), points[0])
+        assert np.allclose(path.endpoint(), points[-1])
+
+    def test_irregular_path_length(self):
+        y = Path.from_irregular_function(lambda x: np.array([x, x**2, 0.]))
+        assert np.isclose(y.path_length, 1.478942857544597), y.path_length
+        y = Path.from_irregular_function(lambda x: np.array([5*x, (5*x)**2, 0.]))
+        assert np.isclose(y.path_length, 25.87424479037671), y.path_length
+        y = Path.from_irregular_function(lambda x: np.array([0., 5*x, (5*x)**2]))
+        assert np.isclose(y.path_length, 25.87424479037671), y.path_length
+        y = Path.from_irregular_function(lambda x: np.array([(5*x)**2, 0., 5*x]))
+
+    def test_move(self):
+        y = Path.from_irregular_function(lambda x: np.array([x, x, 0.]))
+        y = y.move(dz=1.)
+        assert np.allclose(y(sqrt(2)), [1., 1., 1.])
+
+    def test_rotate(self):
+        y = Path.from_irregular_function(lambda x: np.array([x, x, 0.]))
+        y = y.rotate(Rz=pi/2)
+        assert np.allclose(y(sqrt(0)), [0., 0., 0.])
+        assert np.allclose(y(sqrt(2)), [-1., 1., 0.])
+        y = y.rotate(Rz=pi/2)
+        assert np.allclose(y(sqrt(0)), [0., 0., 0.])
+        assert np.allclose(y(sqrt(2)), [-1., -1., 0.])
+        y = y.rotate(Ry=pi/2)
+        assert np.allclose(y(sqrt(0)), [0., 0., 0.])
+        assert np.allclose(y(sqrt(2)), [0., -1., 1.])
+        y = y.rotate(Rx=-pi/2)
+        assert np.allclose(y(sqrt(0)), [0., 0., 0.])
+        assert np.allclose(y(sqrt(2)), [0., 1., 1.])
+
+    def test_rotate_around_point(self):
+        origin = np.array([2., 1., -1.])
+        y = Path.from_irregular_function(lambda x: origin + np.array([x, x, 0.]))
+        y = y.rotate(Rz=pi/2, origin=origin)
+        assert np.allclose(y(sqrt(0)), origin)
+        assert np.allclose(y(sqrt(2)), origin + np.array([-1., 1., 0.]))
+        y = y.rotate(Rz=pi/2, origin=origin)
+        assert np.allclose(y(sqrt(0)), origin)
+        assert np.allclose(y(sqrt(2)), origin + np.array([-1., -1., 0.]))
+        y = y.rotate(Ry=pi/2, origin=origin)
+        assert np.allclose(y(sqrt(0)), origin)
+        assert np.allclose(y(sqrt(2)), origin + np.array([0., -1., 1.]))
+        y = y.rotate(Rx=-pi/2, origin=origin)
+        assert np.allclose(y(sqrt(0)), origin)
+        assert np.allclose(y(sqrt(2)), origin + np.array([0., 1., 1.]))
+
+    def test_discretize_path(self):
+        path_length = 10 
+        breakpoints = [3.33, 5., 9.]
+
+        u = discretize_path(path_length, breakpoints, 1.)
+
+        assert 0. in u
+        assert 10 in u
+
+        for b in breakpoints:
+            assert b in u
+
+    def test_arc(self):
+        r = 2.
+        p = Path.arc([0., 0., 0.], [r, 0., 0.], [0., r, 0.])
+        assert np.isclose(p.path_length, 1/4 * 2*pi*r)
+        assert np.allclose(p(1/8 * 2*pi*r), [r/sqrt(2), r/sqrt(2), 0.])
+
+        center = np.array([1., -1, -1])
+        p = Path.arc(center, center+np.array([r,0.,0.]), center+np.array([0.,r, 0.]))
+        assert np.isclose(p.path_length, 1/4 * 2*pi*r)
+        assert np.allclose(p(1/8 * 2*pi*r), center + np.array([r/sqrt(2), r/sqrt(2), 0.]))
+
+        r = 3
+        center = np.array([0, r, 0.])
+        p = Path.arc(center, [0., 0., 0.], [0., r, r])
+        assert np.isclose(p.path_length, 1/4* 2*pi*r)
+        assert np.allclose(p(1/8 * 2*pi*r), center + np.array([0., -r/sqrt(2), r/sqrt(2)]))
+
+        r = 3
+        center = np.array([0, r, 0.])
+        p = Path.arc(center, [0., 0., 0.], [0., r, r], reverse=True)
+        assert np.isclose(p.path_length, 3/4* 2*pi*r)
+        assert np.allclose(p(1/2 * 3/4*2*pi*r), center + np.array([0., +r/sqrt(2), -r/sqrt(2)]))
+
+        r = 3
+        center = np.array([0, r, 0.])
+        p = Path.arc(center, [0., 0., 0.], [0., r, -r])
+        assert np.isclose(p.path_length, 1/4* 2*pi*r)
+        assert np.allclose(p(1/2 * 1/4*2*pi*r), center + np.array([0., -r/sqrt(2), -r/sqrt(2)]))
+
+        r = 3
+        center = np.array([0, r, 0.])
+        p = Path.arc(center, [0., 0., 0.], [0., r, -r], reverse=True)
+        assert np.isclose(p.path_length, 3/4* 2*pi*r)
+        assert np.allclose(p(1/2 * 3/4*2*pi*r), center + np.array([0., r/sqrt(2), r/sqrt(2)]))
+
+
+
+
+
+class SurfaceTests(unittest.TestCase):
+
+    def test_spanned_by_paths(self):
+        y1 = Path.from_irregular_function(lambda x: [x, -1. - x, x]) 
+        y2 = Path.from_irregular_function(lambda x: [x, 1. + x, x]) 
+        surf = Surface.spanned_by_paths(y1, y2)
+
+        p1 = surf.path_length1
+        p2 = surf.path_length2
+
+        assert np.allclose(surf(0., 0.), [0., -1., 0.])
+        assert np.allclose(surf(0., p2), [0., 1., 0.])
+        assert np.allclose(surf(p1, 0.), [1., -2., 1.])
+        assert np.allclose(surf(p1, p2), [1., 2., 1.])
+
+    def test_non_degenerate_triangles(self):
+        # Found bug where a triangle had two common
+        # points, making it have zero area.
+        points = [[0, 0, 5],
+            [12, 0, 5],
+            [12, 0, 15],
+            [0, 0, 15]]
+
+        inner = Path.line(points[0], points[1])\
+            .line_to(points[2]).line_to(points[3]).revolve_z()
+
+        inner.name = 'test'
+        mesh = inner.mesh(mesh_size=10/2)
+
+        for i, t in enumerate(mesh.triangles):
+            p1, p2, p3 = mesh.points[t]
+            assert not np.all(p1 == p2), (i, t, p1, p2, p3)
+            assert not np.all(p2 == p3)
+            assert not np.all(p3 == p1)
+
+        assert len(mesh.physical_to_triangles['test']) == len(mesh.triangles)
+
+
+
+
+
+
+
+
+
+
+

tests/test_radial.py

@@ -27,9 +27,9 @@ class TestRadial(unittest.TestCase):
     def test_charge_radial_vertical(self):
         vertices = np.array([
             [1.0, 0.0, 0.0],
-            [1.0, 1.0, 0.0],
-            [1.0, 1/3, 0.0],
-            [1.0, 2/3, 0.0]])
+            [1.0, 0.0, 1.0],
+            [1.0, 0.0, 1/3],
+            [1.0, 0.0, 2/3]])
 
         correct = 2*pi
         approx = B.charge_radial(vertices, 1.0);
@@ -51,9 +51,9 @@ def test_charge_radial_horizontal(self):
     def test_charge_radial_skewed(self):
         vertices = np.array([
             [0.0, 0.0, 0.0],
-            [1.0, 1.0, 0.0],
-            [1/3, 1/3, 0.0],
-            [2/3, 2/3, 0.0]])
+            [1.0, 0.0, 1.0],
+            [1/3, 0.0, 1/3],
+            [2/3, 0.0, 2/3]])
 
         correct = pi*sqrt(2)
         approx = B.charge_radial(vertices, 1.0);
@@ -62,10 +62,10 @@ def test_charge_radial_skewed(self):
 
     def test_field_radial(self):
         vertices = np.array([
-            [1.0, 1.0, 0.0],
-            [2.0, 2.0, 0.0],
-            [1.0 + 1/3, 1.0 + 1/3, 0.0],
-            [1.0 + 2/3, 1.0 + 2/3, 0.0]])
+            [1.0, 0.0, 1.0],
+            [2.0, 0.0, 2.0],
+            [1.0 + 1/3, 0.0, 1.0 + 1/3],
+            [1.0 + 2/3, 0.0, 1.0 + 2/3]])
 
         r0, z0 = 2.0, -2.5
 
@@ -90,14 +90,12 @@ def Ez(r, z):
         assert np.allclose(B.field_radial(np.array([r0, z0]), charges, jac, pos)/epsilon_0, [Er, Ez], atol=0.0, rtol=1e-10)
 
     def test_rectangular_current_loop(self):
-        with G.Geometry(G.Symmetry.RADIAL) as geom:
-            points = [[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0]]
-            poly = geom.add_polygon(points)
-            geom.add_physical(poly, 'coil')
-            geom.set_mesh_size_factor(50)
-            mesh = geom.generate_triangle_mesh(False)
-
-        exc = E.Excitation(mesh)
+        coil = G.Surface.rectangle_xz(1.0, 2.0, 1.0, 2.0)
+        coil.name = 'coil'
+
+        mesh = coil.mesh(mesh_size=0.1)
+
+        exc = E.Excitation(mesh, E.Symmetry.RADIAL)
         exc.add_current(coil=5)
 
         field = S.solve_bem(exc)
@@ -180,20 +178,20 @@ def test_interpolation_current_loop(self):
             assert np.allclose(field, field_interp, atol=1e-3, rtol=5e-3)
 
     def test_mag_pot_derivatives(self):
-        with G.Geometry(G.Symmetry.RADIAL) as geom:
-            points = [[0, 5], [5, 5], [5, -5], [0, -5]]
-            lines = [geom.add_line(geom.add_point(p1), geom.add_point(p2)) for p1, p2 in zip(points, points[1:])]
-            geom.add_physical(lines, 'boundary')
-            
-            r1 = geom.add_rectangle(1, 2, 2, 3, 0)
-            r2 = geom.add_rectangle(1, 2, -3, -2, 0)
-            
-            geom.add_physical(r1.curves, 'r1')
-            geom.add_physical(r2.curves, 'r2')
-            geom.set_mesh_size_factor(10)
-            mesh = geom.generate_line_mesh(False)
-
-        e = E.Excitation(mesh)
+        boundary = G.Path.line([0., 0., 5.], [5., 0., 5.])\
+            .line_to([5., 0., -5.])\
+            .line_to([0., 0., -5.])
+
+        r1 = G.Path.rectangle_xz(1, 2, 2, 3)
+        r2 = G.Path.rectangle_xz(1, 2, -3, -2)
+
+        boundary.name = 'boundary'
+        r1.name = 'r1'
+        r2.name = 'r2'
+
+        mesh = (boundary + r1 + r2).mesh(mesh_size=0.1)
+         
+        e = E.Excitation(mesh, E.Symmetry.RADIAL)
         e.add_magnetostatic_potential(r1 = 10)
         e.add_magnetostatic_potential(r2 = -10)
 
@@ -221,25 +219,24 @@ def test_mag_pot_derivatives(self):
     def test_rectangular_coil(self):
         # Field produced by a 1mm x 1mm coil, with inner radius 2mm, 1ampere total current
         # What is the field produced at (2.5mm, 4mm)
-        with G.Geometry(G.Symmetry.RADIAL) as geom:
-            rect = geom.add_rectangle(2, 3, 2, 3, 0)
-            geom.add_physical(rect.surface, 'coil')
-            geom.set_mesh_size_factor(5)
-            mesh = geom.generate_triangle_mesh(False)
+        coil = G.Surface.rectangle_xz(2, 3, 2, 3)
+        coil.name = 'coil'
 
-        exc = E.Excitation(mesh)
+        mesh = coil.mesh(mesh_size=0.1)
+
+        exc = E.Excitation(mesh, E.Symmetry.RADIAL)
         exc.add_current(coil=1)
         field = S.solve_bem(exc)
 
         assert np.isclose(np.sum(field.current_point_charges.jacobians), 1.0) # Area is 1.0
         assert np.isclose(np.sum(field.current_point_charges.charges[:, np.newaxis]*field.current_point_charges.jacobians), 1.0) # Total current is 1.0
 
         target = np.array([2.5, 0., 4.0])
-        correct_r = dblquad(lambda x, y: magnetic_field_of_loop(1.0, x, np.array([target[0], 0.,target[2]-y]))[0], 2, 3, 2, 3)[0]
-        correct_z = dblquad(lambda x, y: magnetic_field_of_loop(1.0, x, np.array([target[0], 0., target[2]-y]))[2], 2, 3, 2, 3)[0]
-
         computed = mu_0*field.current_field_at_point(np.array([target[0], target[2]]))
-        correct = np.array([correct_r, correct_z])
 
-        assert np.allclose(computed, correct, atol=1e-11) 
+        correct_r = dblquad(lambda x, y: mu_0*B.current_field_radial_ring(target[0], target[2], x, y)[0], 2, 3, 2, 3, epsrel=1e-4)[0]
+        correct_z = dblquad(lambda x, y: mu_0*B.current_field_radial_ring(target[0], target[2], x, y)[1], 2, 3, 2, 3, epsrel=1e-4)[0]
+        correct = np.array([correct_r, correct_z])
+
+        assert np.allclose(computed, correct, atol=0.0, rtol=1e-9)