Fixing bug with faults and covering. Also added Makefile to compile f…

…aster.

.gitignore

@@ -55,6 +55,7 @@ target/
 
 # Swap files.
 *.swp
+*.o
 geomodelr/test_files/Version_06_07.json
 
 geomodelr/test_files/Version_06_07.json

geomodelr/cpp/Makefile

@@ -0,0 +1,25 @@
+
+cpp: cpp.so
+	cp cpp.so ..
+
+cpp.so: basic.o section.o match.o model.o geomodel.o faults.o speed_up.o polygon.o isosurfaces_vdb.o
+	g++ -pthread -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -Wl,-z,relro -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wdate-time -D_FORTIFY_SOURCE=2 -g -fstack-protector-strong -Wformat -Werror=format-security -Wl,-Bsymbolic-functions -Wl,-z,relro -Wdate-time -D_FORTIFY_SOURCE=2 -g -fstack-protector-strong -Wformat -Werror=format-security  basic.o section.o match.o model.o geomodel.o faults.o speed_up.o polygon.o isosurfaces_vdb.o -L/usr/local/lib -L/usr/lib/x86_64-linux-gnu/ -lboost_python -lopenvdb -o cpp.so
+
+polygon.o: polygon.cpp polygon.hpp
+	g++ -pthread -DNDEBUG -g -fwrapv -O2 -Wall -fno-strict-aliasing -Wdate-time -D_FORTIFY_SOURCE=2 -g -fstack-protector-strong -Wformat -Werror=format-security -fPIC -I/usr/local/include -I/usr/include/python2.7 -c polygon.cpp -std=c++11
+basic.o: basic.cpp basic.hpp
+	g++ -pthread -DNDEBUG -g -fwrapv -O2 -Wall -fno-strict-aliasing -Wdate-time -D_FORTIFY_SOURCE=2 -g -fstack-protector-strong -Wformat -Werror=format-security -fPIC -I/usr/local/include -I/usr/include/python2.7 -c basic.cpp -std=c++11
+section.o: section.cpp section.hpp
+	g++ -pthread -DNDEBUG -g -fwrapv -O2 -Wall -fno-strict-aliasing -Wdate-time -D_FORTIFY_SOURCE=2 -g -fstack-protector-strong -Wformat -Werror=format-security -fPIC -I/usr/local/include -I/usr/include/python2.7 -c section.cpp -std=c++11
+match.o: match.cpp match.hpp
+	g++ -pthread -DNDEBUG -g -fwrapv -O2 -Wall -fno-strict-aliasing -Wdate-time -D_FORTIFY_SOURCE=2 -g -fstack-protector-strong -Wformat -Werror=format-security -fPIC -I/usr/local/include -I/usr/include/python2.7 -c match.cpp -std=c++11
+model.o: model.cpp model.hpp
+	g++ -pthread -DNDEBUG -g -fwrapv -O2 -Wall -fno-strict-aliasing -Wdate-time -D_FORTIFY_SOURCE=2 -g -fstack-protector-strong -Wformat -Werror=format-security -fPIC -I/usr/local/include -I/usr/include/python2.7 -c model.cpp -std=c++11
+geomodel.o: geomodel.cpp geomodel.hpp
+	g++ -pthread -DNDEBUG -g -fwrapv -O2 -Wall -fno-strict-aliasing -Wdate-time -D_FORTIFY_SOURCE=2 -g -fstack-protector-strong -Wformat -Werror=format-security -fPIC -I/usr/local/include -I/usr/include/python2.7 -c geomodel.cpp -std=c++11
+faults.o: faults.cpp faults.hpp
+	g++ -pthread -DNDEBUG -g -fwrapv -O2 -Wall -fno-strict-aliasing -Wdate-time -D_FORTIFY_SOURCE=2 -g -fstack-protector-strong -Wformat -Werror=format-security -fPIC -I/usr/local/include -I/usr/include/python2.7 -c faults.cpp -std=c++11
+speed_up.o: speed_up.cpp speed_up.hpp
+	g++ -pthread -DNDEBUG -g -fwrapv -O2 -Wall -fno-strict-aliasing -Wdate-time -D_FORTIFY_SOURCE=2 -g -fstack-protector-strong -Wformat -Werror=format-security -fPIC -I/usr/local/include -I/usr/include/python2.7 -c speed_up.cpp -std=c++11
+isosurfaces_vdb.o: isosurfaces_vdb.cpp isosurfaces_vdb.hpp
+	g++ -pthread -DNDEBUG -g -fwrapv -O2 -Wall -fno-strict-aliasing -Wdate-time -D_FORTIFY_SOURCE=2 -g -fstack-protector-strong -Wformat -Werror=format-security -fPIC -I/usr/local/include -I/usr/include/python2.7 -c isosurfaces_vdb.cpp -std=c++11

geomodelr/cpp/match.cpp

@@ -66,13 +66,14 @@ void Match::match_polygons() {
 	// Check if any of the polygons is covered by a fault.
 	auto covered_by_fault = [&] ( const multi_polygon& intersect ) -> bool {
 		for ( const multi_polygon& flt: this->excluded_area | boost::adaptors::map_values ) {
+			// std::cerr << "fault" << geometry::wkt(flt) << " \n";
 			if ( geometry::covered_by( intersect, flt ) ) {
 				return true;
 			}
 		}
 		return false;
 	};
-	
+
 	map<wstring, vector<int>> units_a;
 	map<wstring, vector<int>> units_b;
 
@@ -113,6 +114,7 @@ void Match::match_polygons() {
 								}
 							}
 						} catch ( geometry::exception& e ) {
+							// std::cerr << "error intersecting\n"; 
 							if ( geometry::intersects( this->a->poly_trees[pols_a[i]]->boost_poly, 
 										   this->b->poly_trees[pols_b[j]]->boost_poly ) ) 
 							{
@@ -315,7 +317,7 @@ std::pair<vector<triangle>, bool> faultplane_for_lines(const vector<point3>& l_a
 	if ( angle <= M_PI/2.0 ) {
 		return std::make_pair(test_start( l_a, l_b, false ), false);
 	}
-	return std::make_pair(test_start( l_a, l_b, true  ), false);
+	return std::make_pair(test_start( l_a, l_b, true  ), true);
 }
 
 multi_polygon fix_polygon( polygon& pol ) {
@@ -325,12 +327,13 @@ multi_polygon fix_polygon( polygon& pol ) {
 	geometry::correct(pol);
 	geometry::remove_spikes(pol);
 
-	/*
 	string reason;
+	/*
 	if ( not geometry::is_valid( pol, reason ) ) {
+		std::cerr << geometry::wkt(pol) << "\n";
 		std::cerr << "invlid reason " << reason << "\n";
-	}
-	*/
+	}*/
+
 	// Check for infinite buckle with max_iters.
 	size_t max_iters = pol.outer().size();
 	while( split.size() ) {
@@ -348,10 +351,15 @@ multi_polygon fix_polygon( polygon& pol ) {
 		size_t n = curr.size();
 		for ( size_t i = 0; i < n-1; i++ ) {
 			segment s1(curr[i], curr[(i+1)%n]);
-			for ( size_t j = i+2; j < n-1; j++ ) {
+			for ( size_t j = i+2; j < n; j++ ) {
+				if ( (j+1)%n == i ) {
+					continue;
+				}
 				segment s2(curr[j], curr[(j+1)%n]);
 				std::vector<point2> output;
 				bool inter = boost::geometry::intersection(s1, s2, output);
+				// std::cerr << "output " << i << " " << j << " size " << output.size() << "\n";
+
 				if ( inter && output.size() == 1 ) {
 					polygon p1, p2;
 					ring& o1 = p1.outer();

geomodelr/cpp/model.cpp

@@ -1039,16 +1039,16 @@ double Topography::height(const point2& pt) const {
 	double B = this->heights[(i+1)*dims[1] + j];
 	double C = this->heights[(i+1)*dims[1] + j+1];
 	double D = this->heights[i*dims[1] + j+1];
-
+	
 	double x = 2.0*(gx(pos)- double(i)) - 1.0;
 	double y = 2.0*(gy(pos)- double(j)) - 1.0;
-
-	x = std::max(std::min(x,1.0),-1.0);
-	y = std::max(std::min(y,1.0),-1.0);
+	
+	x = std::max( std::min( x, 1.0), -1.0 );
+	y = std::max( std::min( y, 1.0), -1.0 );
 
 	double v1 = B + A + x*(B - A);
 	double v2 = C + D + x*(C - D);
-		
+
 	return 0.25*(v2 + v1 + y*(v2 - v1));
 }
 

geomodelr/cpp/polygon.cpp

@@ -44,8 +44,10 @@ double Polygon::distance_point_cover_faults(const point2& pt) const {
 		return 0.0;
 	} else {
 		std::pair<line_segment,double> ray_dist_pair = this->ray_distance(pt);
+		// std::cerr << geometry::wkt(ray_dist_pair.first) << " " << ray_dist_pair.second << "\n";
 		for ( auto it = this->section->fault_lines->qbegin( geometry::index::intersects(ray_dist_pair.first));
 			   it != this->section->fault_lines->qend(); it++ ) {
+
 			const point2& ps = this->section->line_ends[g1(*it)].first;
 			const point2& pe = this->section->line_ends[g1(*it)].second;
 			return std::min( this->ray_crossing( pt, ps ), this->ray_crossing(pt, pe) );
@@ -163,7 +165,7 @@ double Polygon::ray_crossing ( const point2& pt, const point2& nd ) const {
 	return mind;
 }
 
-double Polygon::set_distance_function( const wstring& s ) {
+void Polygon::set_distance_function( const wstring& s ) {
 	if ( s == L"cover" ) {
 		this->distance_point = std::bind(&Polygon::distance_point_cover_faults, this, std::placeholders::_1);
 	} else {

geomodelr/cpp/polygon.hpp

@@ -49,7 +49,7 @@ class Polygon {
 	double ray_crossing ( const point2& pt, const point2& nd ) const;
 	double distance_point_basic_faults(const point2& pt) const;
 	double distance_point_cover_faults(const point2& pt) const;
-	double set_distance_function( const wstring& s );
+	void set_distance_function( const wstring& s );
 };
 
 #endif

geomodelr/test.py

@@ -36,7 +36,7 @@
 import math
 from shapely.geometry import Polygon, Point
 
-class TestGeoModelR(unittest.TestCase):
+class Tests(unittest.TestCase):
     def setUp(self):
         pass
         #Profile GeoModelR
@@ -800,23 +800,23 @@ def query_func(p):
 
         # Calculate bounded.
         verts, triangs = isosurfaces.calculate_isosurface(m, "Anfibolitas", 50, method="OPENVDB" )
-        self.assertEqual(len(verts), 34748)
-        self.assertEqual(len(triangs), 69472)
+        self.assertEqual(len(verts), 22974)
+        self.assertEqual(len(triangs), 45924)
 
         # Calculate unbounded
         verts, triangs = isosurfaces.calculate_isosurface(m, "Anfibolitas", 50, False, method="OPENVDB" )
-        self.assertEqual(len(verts), 38660)
-        self.assertEqual(len(triangs), 68407)
+        self.assertEqual(len(verts), 23845)
+        self.assertEqual(len(triangs), 42130)
 
         # Filter by normal.
         verts, triangs = isosurfaces.calculate_isosurface(m, "Anfibolitas", 50, False, True, True, method="OPENVDB" )
-        self.assertEqual(len(verts), 13092)
-        self.assertEqual(len(triangs), 24029)
+        self.assertEqual(len(verts), 8813)
+        self.assertEqual(len(triangs), 15696)
 
         # Filter by normal, negative.
         verts, triangs = isosurfaces.calculate_isosurface(m, "Anfibolitas", 50, False, True, False, method="OPENVDB" )
-        self.assertEqual(len(verts), 23631)
-        self.assertEqual(len(triangs), 45203)
+        self.assertEqual(len(verts), 14528)
+        self.assertEqual(len(triangs), 27184)
 
     def test_modflow(self):
 
@@ -1429,7 +1429,22 @@ def test_soils(self):
         self.assertEqual(model.signed_distance( "unit1", ( 1.5, 0.0, 2.25 ) ), -0.25)
         self.assertEqual(model.signed_distance( "unit1", ( 1.5, 0.0, 2.0 ) ),   0.0 )
         self.assertEqual(model.signed_distance( "unit1", ( 1.5, 0.0, 1.5 ) ),   0.5 )
-
+
+    def test_v1(self):
+        this_dir, this_filename = os.path.split(__file__)
+        fn = os.path.join(this_dir, 'test_files', 'v1.json')
+        geo_model = model.model_from_file(fn)
+        geo_model.params = { "faults": "cover" }
+        self.assertEqual(geo_model.matches, [((u'name2', u'name1'), [(0, 0), (1, 1), (2, 2)])])
+        cls = geo_model.closest((834590.52,1174732.73,1450.80))
+        self.assertEqual(cls[0], "U2")
+        self.assertAlmostEqual(cls[1], 51.996853528295254)
+        cls = geo_model.closest((834588,1174732.73,1450.80))
+        self.assertEqual(cls[0], "U1")
+        self.assertAlmostEqual(cls[1], 0)
+
+
+
 def main(args=None):
     unittest.main()