add old compgeom repo to new common repository

.gitignore

@@ -9,3 +9,5 @@
 *.pos 
 *.vtu
 *.vrml
+*~
+*.swo

src/porepy/utils/geometry.py

@@ -0,0 +1,53 @@
+import numpy as np
+
+#------------------------------------------------------------------------------#
+
+def project_plane( pts, normal = None ):
+    """ Project the points on a plane using local coordinates.
+
+    Parameters:
+    pts: np.ndarray, 3xn, the points.
+    normal: (optional) the normal of the plane, otherwise three points are
+    required.
+
+    Returns:
+    pts: np.array, 2xn, projected points on the plane in the local coordinates.
+
+    """
+
+    if normal is None:
+    	normal = compute_normal( pts )
+    else:
+    	normal = normal / np.linalg.norm( normal )
+
+    # Projection matrix onto tangential plane
+    T = np.identity(3) - np.outer( normal, normal )
+    # Projected points
+    pts = np.array( [ np.dot( T, p ) for p in pts.T ] )
+    # Disregard points on the origin??
+    index = np.where( np.sum( np.abs( pts ), axis = 0 ) != 0 )[0]
+    return pts[:,index]
+
+#------------------------------------------------------------------------------#
+
+def compute_normal( pts ):
+    """ Compute the normal of a set of points.
+
+    The algorithm computes the normal of the plane defined by the first three
+    points in the set. 
+    TODO: Introduce optional check that all points lay in the same plane
+    (should be separate subroutine).
+
+    Parameters:
+    pts: np.ndarray, 3xn, the points. Need n > 2.
+
+    Returns:
+    normal: np.array, 1x3, the normal.
+
+    """
+
+    assert( pts.shape[1] > 2 )
+    normal = np.cross( pts[:,0] - pts[:,1], pts[:,0] - pts[:,2] )
+    return normal / np.linalg.norm( normal )
+
+#------------------------------------------------------------------------------#

src/porepy/utils/sort_points.py

@@ -0,0 +1,93 @@
+import numpy as np
+from compgeom.basics import project_plane_matrix
+
+#------------------------------------------------------------------------------#
+
+def sort_point_pairs(lines, check_circular=True, ordering=False):
+    """ Sort pairs of numbers to form a chain.
+
+    The target application is to sort lines, defined by their
+    start end endpoints, so that they form a continuous polyline.
+
+    The algorithm is brute-force, using a double for-loop. This can
+    surely be imporved.
+
+    Parameters:
+    lines: np.ndarray, 2xn, the line pairs.
+    check_circular: Verify that the sorted polyline form a circle.
+                    Defaluts to true.
+    ordering: np.array, return in the original order if a line is flipped or not
+
+    Returns:
+    sorted_lines: np.ndarray, 2xn, sorted line pairs.
+
+    """
+
+    num_lines = lines.shape[1]
+    sorted_lines = -np.ones( (2, num_lines), dtype = lines.dtype )
+
+    # Start with the first line in input
+    sorted_lines[:, 0] = lines[:, 0]
+
+    # The starting point for the next line
+    prev = sorted_lines[1, 0]
+
+    # Keep track of which lines have been found, which are still candidates
+    found = np.zeros(num_lines, dtype=np.bool)
+    found[0] = True
+
+    # Order of the origin line list, store if they're flip or not to form the chain
+    is_ordered = np.zeros(num_lines, dtype=np.bool)
+    is_ordered[0] = True
+
+    # The sorting algorithm: Loop over all places in sorted_line to be filled,
+    # for each of these, loop over all members in lines, check if the line is still
+    # a candidate, and if one of its points equals the current starting point.
+    # More efficient and more elegant approaches can surely be found, but this
+    # will do for now.
+    for i in range(1, num_lines):  # The first line has already been found
+        for j in range(0, num_lines):
+            if not found[j] and lines[0, j] == prev:
+                sorted_lines[:, i] = lines[:, j]
+                found[j] = True
+                prev = lines[1, j]
+                is_ordered[j] = True
+                break
+            elif not found[j] and lines[1, j] == prev:
+                sorted_lines[:, i] = lines[::-1, j]
+                found[j] = True
+                prev = lines[0, j]
+                break
+    # By now, we should have used all lines
+    assert(np.all(found))
+    if check_circular:
+        assert sorted_lines[0, 0] == sorted_lines[1, -1]
+    if ordering: return sorted_lines, is_ordered
+    else:        return sorted_lines
+
+#------------------------------------------------------------------------------#
+
+def sort_point_plane( pts, centre, normal = None ):
+    """ Sort the points which lie on a plane.
+
+    The algorithm assumes a star-shaped disposition of the points with respect
+    the centre.
+
+    Parameters:
+    pts: np.ndarray, 3xn, the points.
+    centre: np.ndarray, 3x1, the face centre.
+    normal: (optional) the normal of the plane, otherwise three points are
+    required.
+
+    Returns:
+    map_pts: np.array, 1xn, sorted point ids.
+
+    """
+    R = project_plane_matrix( pts, normal )
+    pts = np.array( [ np.dot(R, p) for p in pts.T ] ).T
+    centre = np.dot(R, centre)
+    delta = np.array( [ p - centre for p in pts.T] ).T[0:2,:]
+    delta = np.array( [ d / np.linalg.norm(d) for d in delta.T] ).T
+    return np.argsort( np.arctan2( *delta ) )
+
+#------------------------------------------------------------------------------#

src/porepy/utils/tests/test_basics.py

@@ -0,0 +1,71 @@
+import numpy as np
+import unittest
+
+from compgeom import basics
+
+#------------------------------------------------------------------------------#
+
+class BasicsTest( unittest.TestCase ):
+
+#------------------------------------------------------------------------------#
+
+    def test_compute_normal_2d( self ):
+        pts = np.array( [ [ 0., 2., -1. ],
+                          [ 0., 4.,  2. ],
+                          [ 0., 0.,  0. ] ] )
+        normal = basics.compute_normal( pts )
+        normal_test = np.array([0.,0.,1.])
+        pt = pts[:,0]
+
+        assert np.allclose( np.linalg.norm( normal ), 1. )
+        assert np.allclose( [ np.dot( normal, p - pt ) \
+                              for p in pts[:,1:].T ],
+                            np.zeros( pts.shape[1] - 1 )  )
+        assert np.allclose( normal, normal_test )
+
+#------------------------------------------------------------------------------#
+
+    def test_compute_normal_3d( self ):
+        pts = np.array( [ [  2.,  0.,  1.,  1. ],
+                          [  1., -2., -1.,  1. ],
+                          [ -1.,  0.,  2., -8. ] ] )
+        normal_test = np.array( [7., -5., -1.] )
+        normal_test = normal_test / np.linalg.norm( normal_test )
+        normal = basics.compute_normal( pts )
+        pt = pts[:,0]
+
+        assert np.allclose( np.linalg.norm( normal ), 1. )
+        assert np.allclose( [ np.dot( normal, p - pt ) \
+                              for p in pts[:,1:].T ],
+                            np.zeros( pts.shape[1] - 1 )  )
+        assert np.allclose( normal, normal_test ) or \
+               np.allclose( normal, -1. * normal_test )
+
+#------------------------------------------------------------------------------#
+
+    def test_is_planar_2d( self ):
+        pts = np.array( [ [ 0., 2., -1. ],
+                          [ 0., 4.,  2. ],
+                          [ 2., 2.,  2. ] ] )
+        assert basics.is_planar( pts )
+
+#------------------------------------------------------------------------------#
+
+    def test_is_planar_3d( self ):
+        pts = np.array( [ [    0.,    1.,    0., 4./7. ],
+                          [    0.,    1.,    1.,    0. ],
+                          [ 5./8., 7./8., 7./4., 1./8. ] ] )
+        assert basics.is_planar( pts )
+
+#------------------------------------------------------------------------------#
+
+    def test_project_plane( self ):
+        pts = np.array( [ [  2.,  0.,  1.,  1. ],
+                          [  1., -2., -1.,  1. ],
+                          [ -1.,  0.,  2., -8. ] ] )
+        R = basics.project_plane_matrix( pts )
+        P_pts = np.dot( R, pts )
+
+        assert np.allclose( P_pts[2,:], 1.15470054 * np.ones(4) )
+
+#------------------------------------------------------------------------------#

src/porepy/utils/tests/test_ccw.py

@@ -0,0 +1,56 @@
+import numpy as np
+import unittest
+
+from compgeom import basics as cg
+
+
+class TestCCW(unittest.TestCase):
+
+    def setup(self):
+        p1 = np.array([0, 0])
+        p2 = np.array([1, 0])
+        p3 = np.array([1, 1])
+        return p1, p2, p3
+
+    def setup_close(self, y):
+        p1 = np.array([0, 0])
+        p2 = np.array([2, 0])
+        p3 = np.array([1, y])
+        return p1, p2, p3
+
+    def test_is_ccw(self):
+        p1, p2, p3 = self.setup()
+        assert cg.is_ccw_polyline(p1, p2, p3)
+
+    def test_not_ccw(self):
+        p1, p2, p3 = self.setup()
+        assert not cg.is_ccw_polyline(p1, p3, p2)
+
+    def test_tolerance(self):
+        p1, p2, p3 = self.setup_close(1e-6)
+
+        # Safety margin saves ut
+        assert cg.is_ccw_polyline(p1, p2, p3, tol=1e-4, default=True)
+
+        # Default kills us, even though we're inside safety margin
+        assert not cg.is_ccw_polyline(p1, p2, p3, tol=1e-4, default=False)
+
+        # Outside safety margin, and on the ccw side
+        assert cg.is_ccw_polyline(p1, p2, p3, tol=1e-8, default=False)
+
+    def test_tolerance_outside(self):
+        p1, p2, p3 = self.setup_close(-1e-6)
+
+        # Safety margin saves ut
+        assert cg.is_ccw_polyline(p1, p2, p3, tol=1e-4, default=True)
+
+        # Default kills us, even though we're inside safety margin
+        assert not cg.is_ccw_polyline(p1, p2, p3, tol=1e-4, default=False)
+
+        # Outside safety margin, and not on the ccw side
+        assert not cg.is_ccw_polyline(p1, p2, p3, tol=1e-8, default=False)
+
+
+    if __name__ == '__main__':
+        unittest.main()
+

src/porepy/utils/tests/test_distance_segments.py

@@ -0,0 +1,68 @@
+import unittest
+import numpy as np
+
+from compgeom import basics
+
+class TestSegmentDistance(unittest.TestCase):
+    def setup_2d_unit_square(self):
+        p00 = np.array([0, 0])
+        p10 = np.array([1, 0])
+        p01 = np.array([0, 1])
+        p11 = np.array([1, 1])
+        return p00, p10, p11, p01
+
+    def test_segment_no_intersect_2d(self):
+        p00, p10, p11, p01 = self.setup_2d_unit_square()
+        d = basics.distance_segment_segment(p00, p01, p11, p10)
+        assert d == 1
+
+    def test_segment_intersect_2d(self):
+        p00, p10, p11, p01 = self.setup_2d_unit_square()
+        d = basics.distance_segment_segment(p00, p11, p10, p01)
+        assert d == 0
+
+    def test_line_passing(self):
+        # Lines not crossing
+        p1 = np.array([0, 0])
+        p2 = np.array([1, 0])
+        p3 = np.array([2, -1])
+        p4 = np.array([2, 1])
+        d = basics.distance_segment_segment(p1, p2, p3, p4)
+        assert d == 1
+
+    def test_share_point(self):
+        # Two lines share a point
+        p1 = np.array([0, 0])
+        p2 = np.array([0, 1])
+        p3 = np.array([1, 1])
+        d = basics.distance_segment_segment(p1, p2, p2, p3)
+        assert d == 0
+
+    def test_intersection_3d(self):
+        p000 = np.array([0, 0, 0])
+        p111 = np.array([1, 1, 1])
+        p100 = np.array([1, 0, 0])
+        p011 = np.array([0, 1, 1])
+        d = basics.distance_segment_segment(p000, p111, p100, p011)
+        assert d == 0
+
+    def test_changed_order_3d(self):
+        # The order of the start and endpoints of the segments should not matter
+        dim = 3
+        p1 = np.random.rand(1, 3)[0]
+        p2 = np.random.rand(1, 3)[0]
+        p3 = np.random.rand(1, 3)[0]
+        p4 = np.random.rand(1, 3)[0]
+        d1 = basics.distance_segment_segment(p1, p2, p3, p4)
+        d2 = basics.distance_segment_segment(p2, p1, p3, p4)
+        d3 = basics.distance_segment_segment(p1, p2, p4, p3)
+        d4 = basics.distance_segment_segment(p2, p1, p4, p3)
+        d5 = basics.distance_segment_segment(p4, p3, p2, p1)
+        assert np.allclose(d1, d2)
+        assert np.allclose(d1, d3)
+        assert np.allclose(d1, d4)
+        assert np.allclose(d1, d5)
+
+    if __name__ == '__main__':
+        unittest.main()
+

src/porepy/utils/tests/test_doctests.py

@@ -0,0 +1,8 @@
+import unittest
+import doctest
+from compgeom import basics
+
+test_suite = unittest.TestSuite()
+test_suite.addTest(doctest.DocTestSuite(basics))
+
+unittest.TextTestRunner().run(test_suite)