Fixed problems with faults, that did not triangulate when flat.

geomodelr/cpp/match.cpp

@@ -461,7 +461,23 @@ vector<triangle> Match::triangulate( const vector<point3>& l_a, const vector<poi
 	for ( const auto& p: l_b ) {
 		points.append(python::make_tuple(gx(p), gy(p), gz(p)));
 	}
-	pyobject pytris = Match::pytriangulate(points, na);
+	pyobject pytris;
+	try {
+		pytris = Match::pytriangulate(points, na);
+	} catch ( const python::error_already_set& ex ) {
+		PyObject *e, *v, *t;
+		PyErr_Fetch(&e, &v, &t);
+		if (!e) {
+			throw GeomodelrException("There was an error triangulating the faults");
+		}
+
+		pyobject e_obj(python::handle<>(python::allow_null(e)));
+		pyobject v_obj(python::handle<>(python::allow_null(v)));
+		pyobject t_obj(python::handle<>(python::allow_null(t)));
+
+
+		throw GeomodelrException("There was an error triangulating the faults");
+	}
 	vector<triangle> tris;
 	int nt = python::len(pytris);
 	for ( int i = 0; i < nt; i++ ) {

geomodelr/shared.py

@@ -20,6 +20,7 @@
 from numpy import linalg as la
 from scipy.spatial import Delaunay
 import itertools
+import math
 
 class GeometryException(Exception):
     pass
@@ -40,28 +41,76 @@ def shape_list(shape, sh_type):
            return [shape]
        return []
 
+def over_point( points ):
+    """
+    Returns a point perpendicular to the plane and above all the others. 
+    This is to avoid errors with the triangulation.
+    """
+    mxx = -float('inf')
+    mxy = -float('inf')
+    mxz = -float('inf')
+    mnx =  float('inf')
+    mny =  float('inf')
+    mnz =  float('inf')
+
+    for p in points:
+        mxx = max(p[0], mxx)
+        mxy = max(p[1], mxy)
+        mxz = max(p[2], mxz)
+        mnx = min(p[0], mnx)
+        mny = min(p[1], mny)
+        mnz = min(p[2], mnz)
+
+    maxd = -1
+    maxp = None
+    for p in points:
+        dx = max(mxx-p[0], p[0]-mnx)
+        dy = max(mxy-p[0], p[0]-mny)
+        dz = max(mxz-p[0], p[0]-mnz)
+        sq = math.sqrt(dx*dx+dy*dy+dz*dz)
+        if sq > maxd:
+            maxd = sq
+            maxp = p
+
+    p0 = np.array(maxp)
+    A = np.array( map( lambda p: [p[0], p[1], 1], points ) )
+    y = np.array( map( lambda p: p[2], points ) )
+    x = la.lstsq(A, y)[0]
+    d = math.sqrt( 1 + x[0]*x[0] + x[1]*x[1] )
+    x /= d
+    n = np.array([-x[0], -x[1], 1.0/d])
+    sqd = math.sqrt((mxx-mnx)*(mxx-mnx)+(mxy-mny)*(mxy-mny)+(mxz-mnz)*(mxz-mnz))
+
+    return p0 + n*sqd
+
 # Triangulates a set of points and returns the triangles filtered by the given filter.
 def triangulate(points, fltr=lambda x: True):
+
     tetras = Delaunay(points)
     tris = set()
     for tet in tetras.simplices.copy():
         for tri in itertools.combinations(tet,3):
             srt = tuple(sorted(tri))
-            if not tri in tris and fltr(srt):
+            if fltr(srt):
                 tris.add(srt)
     return list(tris)
 
 # Triangulates a set of points filtering them because they are on both sides of the line.
 def opposite_triangles(points, na):
+    n = len(points)
+    pt = over_point(points)
     def is_opposite(tri):
+        if n in tri:
+            return False
         if not tri[0] < na or not tri[2] >= na:
             return False
         if tri[0]+1 == tri[1] and tri[1] < na:
             return True
         if tri[1]+1 == tri[2] and tri[1] >= na:
             return True
         return False
-    return map( lambda tri: map(int, tri), triangulate(points, is_opposite) )
+
+    return map( lambda tri: map(int, tri), triangulate(points+[pt], is_opposite) )
 
 # Calculates the distance between two parallel lines.
 def line_side(surface_line, cs_line):

geomodelr/test.py

@@ -97,7 +97,7 @@ def test_faultplane_for_lines(self):
                                                                (1, 10, 11), (1, 11, 12), (1, 12, 13), 
                                                                (1, 2, 13), (2, 3, 13), (3, 13, 14), 
                                                                (3, 4, 14), (4, 5, 14), (5, 14, 15), 
-                                                               (5, 6, 15), (6, 7, 15), (7, 15, 16)])
+                                                               (5, 6, 15), (6, 15, 16), (6, 7, 16)])
 
         # Test that you can create cross sections and that bugs throw something in python (and not segfault).
     def test_sections(self):

geomodelr/utils.py

@@ -205,10 +205,12 @@ def should_refine_cell(idx):
 
     checked = np.zeros((elems.shape[0],), dtype=bool)
     for r in xrange(oct_refine):
+
         elems_ti = []
         points_ti = []
         units_ti = []
         creamids = {}
+
         for idx in xrange(elems.shape[0]):
             if not checked[idx] and should_refine_cell(idx):
                 newelems = []
@@ -232,7 +234,7 @@ def should_refine_cell(idx):
                 points_ti += newpoints
             else:
                 checked[idx] = True
-
+        
         units = np.insert(units, units.shape[0], units_ti)
         points = np.insert(points, points.shape[0], np.array(points_ti), axis=0)
         elems = np.insert(elems, elems.shape[0], elems_ti, axis=0)