Merge branch 'modflow'

geomodelr/cpp/faults.cpp

@@ -23,6 +23,7 @@
 // ====================== AUXILIAR FUNCTIONS =================================
 // Converts the python dictionary of faults into a C++ map.
 //      faults_python: python dictionary of faults.
+
 map<wstring, vector<triangle_pt> > pydict_to_map(const pydict& faults_python){
 
     pylist dict_keys = faults_python.keys();
@@ -764,4 +765,4 @@ pydict find_faults_topography_intersection_python(const pydict& fplanes, const p
 
     return map_to_pydict(find_faults_topography_intersection(faults_cpp,topography_array, z_max, z_min, x_inf, y_inf, dx, dy,
         rows, cols));
-}
+}

geomodelr/cpp/geomodel.cpp

@@ -220,9 +220,9 @@ BOOST_PYTHON_MODULE(cpp)
 
 	// Main exported class, Model.
 	python::class_<ModelPython>("Model", python::init<const pylist&, const pyobject&, const pyobject&, const pyobject&,
-					         const pyobject&, pylist&>())
+					         const pyobject&, pylist&, pydict&>())
 					    .def(python::init<const pylist&, const pyobject&,
-					    	 const pyobject&, pylist&>())
+					    	 const pyobject&, pylist&, pydict&>())
 					    .def("make_matches", &ModelPython::make_matches)
 					    .def("possible_closest", &ModelPython::possible_closest, python::args("point"), doc_possible_closest)
 					    .def("model_point", &ModelPython::model_point, python::args("point"), doc_model_point)
@@ -236,11 +236,18 @@ BOOST_PYTHON_MODULE(cpp)
 					    .def("intersect_plane", &ModelPython::intersect_plane, doc_intersect_plane)
 					    .def("intersect_planes", &ModelPython::intersect_planes, doc_intersect_planes)
 					    .def("intersect_topography", &ModelPython::intersect_topography)
+					    .def("find_unit_limits", &ModelPython::find_unit_limits)
 					    .def("info", &ModelPython::info)
 					    .add_property("bbox", &ModelPython::pybbox)
 					    .add_property("matches", &ModelPython::get_matches, &ModelPython::set_matches)
+					    .add_property("lines", &ModelPython::get_lines)
+					    .add_property("not_extended_lines", &ModelPython::get_not_extended_lines)
 					    .add_property("faults", &ModelPython::get_faults)
-					    .add_property("not_extended_faults", &ModelPython::get_not_extended);
+					    .add_property("not_extended_faults", &ModelPython::get_not_extended_faults)
+					    .add_property("fracts", &ModelPython::get_fracts)
+					    .add_property("not_extended_fracts", &ModelPython::get_not_extended_fracts)
+					    .add_property("veins", &ModelPython::get_veins)
+					    .add_property("not_extended_veins", &ModelPython::get_not_extended_veins);
 
 }
 

geomodelr/cpp/geomodel.hpp

@@ -23,5 +23,6 @@
 #include "match.hpp"
 #include "model.hpp"
 #include "faults.hpp"
+#include "speed_up.hpp"
 
 #endif

geomodelr/cpp/match.cpp

@@ -265,7 +265,7 @@ vector<triangle> faultplane_for_lines(const vector<point3>& l_a, const vector<po
 	return test_start( l_a, l_b, true  );
 }
 
-std::tuple<map<wstring, vector<triangle_pt>>, map<wstring, vector<size_t>>> Match::match_lines()
+std::tuple<map<wstring, vector<triangle_pt>>, map<wstring, vector<size_t>>> Match::match_lines( const map<wstring, wstring>& feature_types )
 {
 	/*
 	Creates the fault planes given the cross sections with faults with the same name.
@@ -373,6 +373,11 @@ std::tuple<map<wstring, vector<triangle_pt>>, map<wstring, vector<size_t>>> Matc
 	}
 	vector<value_f> envelopes;
 	for ( auto it = this->faults.begin(); it != this->faults.end(); it++ ) {
+		wstring ft = feature_types.find(it->first)->second;
+		if ( ft != L"FAULT" ) {
+			// Don't check other types that are not faults.
+			continue;
+		}
 		const auto& triangles = it->second;
 		for ( size_t i = 0; i < triangles.size(); i++ ) {
 			const auto& tr = triangles[i];

geomodelr/cpp/match.hpp

@@ -49,7 +49,7 @@ class Match {
 	Match( const Section * a, const Section * b );
 	virtual ~Match();
 	void match_polygons();
-	std::tuple<map<wstring, vector<triangle_pt>>, map<wstring, vector<size_t>>> match_lines();
+	std::tuple<map<wstring, vector<triangle_pt>>, map<wstring, vector<size_t>>> match_lines( const map<wstring, wstring>& feature_types );
 	std::tuple<int, int, int> crosses_triangles(const point2& pt, double cut) const;
 	void set( const vector<std::pair<int, int>>& match );
 };

geomodelr/cpp/model.cpp

@@ -111,7 +111,7 @@ void Model::make_matches() {
 		// Match the polygons.
 		this->match.back()->match_polygons();
 		// Get the matching faults.
-		auto m = this->match.back()->match_lines();
+		auto m = this->match.back()->match_lines( this->feature_types );
 		add_to_faults(m);
 	}
 
@@ -425,6 +425,10 @@ map<wstring,vector<line>> Model::intersect_plane(const line_3d& plane) const{
 	return find_faults_multiple_planes_intersection(this->global_faults, vector<line_3d>(1, plane));
 }
 
+std::pair<double, bool> Model::find_unit_limits(double xp, double yp,double z_max, double z_min, double eps) const{
+	return  find_unit_limits_cpp(this, xp, yp, z_max, z_min, eps);
+}
+
 vector<std::tuple<wstring, double, double>> Model::possible_closest( const point3& pt ) const {
 	if ( not this->sections.size() ) {
 		return vector<std::tuple<wstring, double, double>>();
@@ -597,6 +601,13 @@ pydict ModelPython::intersect_plane(const pylist& plane) const{
 	return (map_to_pydict(((Model *)this)->intersect_plane(plane_cpp)));
 }
 
+pytuple ModelPython::find_unit_limits(double xp, double yp, double z_max, double z_min, double eps) const{
+
+	std::pair<double, bool> output = ((Model *)this)->find_unit_limits(xp, yp, z_max, z_min, eps);
+	return python::make_tuple(output.first,output.second);
+}
+
+
 double ModelPython::height( const pyobject& pt ) const {
 	double d0 = python::extract<double>(pt[0]);
 	double d1 = python::extract<double>(pt[1]);
@@ -670,7 +681,7 @@ double TopographyPython::height( const pyobject& pypt ) const {
 	return ((Topography *)this)->height( point2(python::extract<double>(pypt[0]), python::extract<double>(pypt[1])) );
 }
 
-void ModelPython::fill_model( const pyobject& topography, const pylist& sections ) {
+void ModelPython::fill_model( const pyobject& topography, const pylist& sections, const pydict& feature_types ) {
 	if ( python::len(topography) > 0 ) {// Get all the information sent from python to build the topography.
 
 		const pyobject& point = python::extract<pyobject>(topography["point"]);
@@ -703,11 +714,16 @@ void ModelPython::fill_model( const pyobject& topography, const pylist& sections
 	for ( size_t i = 0; i < this->sections.size(); i++ ) {
 		this->cuts.push_back(this->sections[i]->cut);
 	}
+	pylist keys = feature_types.keys();
+	for ( int i = 0; i < python::len( keys ); i++ ) {
+		this->feature_types[python::extract<wstring>(keys[i])] = python::extract<wstring>(feature_types[keys[i]]);
+	}
 }
 
 ModelPython::ModelPython( const pyobject& bbox, const pyobject& base_point, 
 		    	 const pyobject& direction, const pyobject& map, 
-		    	 const pyobject& topography, const pylist& sections ): 
+		    	 const pyobject& topography, const pylist& sections, 
+			 const pydict& feature_types ): 
 	Model(make_tuple(std::tuple<double, double, double>(python::extract<double>(bbox[0]),python::extract<double>(bbox[1]),python::extract<double>(bbox[2])), 
 			 std::tuple<double, double, double>(python::extract<double>(bbox[3]),python::extract<double>(bbox[4]),python::extract<double>(bbox[5]))),
 		point2( python::extract<double>(base_point[0]), python::extract<double>(base_point[1]) ),
@@ -719,12 +735,13 @@ ModelPython::ModelPython( const pyobject& bbox, const pyobject& base_point,
 		pylist pl = python::extract<pylist>( sections[i] );
 		pl.insert( 2, sbbox );
 	}
-	this->fill_model( topography, sections );
+	this->fill_model( topography, sections, feature_types );
 }
 
 
 ModelPython::ModelPython( const pyobject& bbox, const pyobject& map, 
-			  const pyobject& topography, const pylist& sections ): 
+			  const pyobject& topography, const pylist& sections,
+			  const pydict& feature_types ): 
 	Model(make_tuple(std::tuple<double, double, double>(python::extract<double>(bbox[0]),python::extract<double>(bbox[1]),python::extract<double>(bbox[2])), 
 			 std::tuple<double, double, double>(python::extract<double>(bbox[3]),python::extract<double>(bbox[4]),python::extract<double>(bbox[5]))))
 {
@@ -734,18 +751,18 @@ ModelPython::ModelPython( const pyobject& bbox, const pyobject& map,
 		pylist pl = python::extract<pylist>( sections[i] );
 		pl.insert( 2, sbbox );
 	}
-	this->fill_model( topography, sections );
+
+	this->fill_model( topography, sections, feature_types );
 }
 
 
 void ModelPython::make_matches() {
 	// map<wstring, vector<triangle_pt>> faults = ((Model *)this)->make_matches();
-
 	((Model *)this)->make_matches();
 
 }
 
-pydict ModelPython::get_faults() const {
+pydict ModelPython::filter_lines( bool ext, const wstring& ft ) const {
 	auto python_point = [&]( const point3& pt ) {
 		return python::make_tuple(gx(pt), gy(pt), gz(pt));
 	};
@@ -756,31 +773,18 @@ pydict ModelPython::get_faults() const {
 
 	// Now convert faults to python and return.
 	pydict ret;
-	for ( auto it = this->global_faults.begin(); it != this->global_faults.end(); it++ ) {
+	auto add_line_triangles = [&]( map<wstring, vector<triangle_pt>>::const_iterator it ) {
 		pylist tris;
 		for ( const triangle_pt& t: it->second ) {
 			pytuple tr = python_triangle( t );
 			tris.append( tr );
 		}
 		ret[it->first] = tris;
-	}
-
-	return ret;
-}
-
-pydict ModelPython::get_not_extended() const {
-	auto python_point = [&]( const point3& pt ) {
-		return python::make_tuple(gx(pt), gy(pt), gz(pt));
-	};
-
-	auto python_triangle = [python_point]( const triangle_pt& tr ) {
-		return python::make_tuple(python_point(g0(tr)), python_point(g1(tr)), python_point(g2(tr)));
 	};
 
-	// Now convert faults to python and return.
-	pydict ret;
-	for ( auto it = this->global_faults.begin(); it != this->global_faults.end(); it++ ) {
+	auto add_line_triangles_not_ext = [&]( map<wstring, vector<triangle_pt>>::const_iterator it ) {
 		pylist tris;
+		// Find the extended lines of the given plane.
 		const auto ef = this->extended_faults.find(it->first);
 		for ( size_t i = 0; i < it->second.size(); i++ ) {
 			if ( ef != this->extended_faults.end() ) {
@@ -794,11 +798,38 @@ pydict ModelPython::get_not_extended() const {
 			tris.append( tr );
 		}
 		ret[it->first] = tris;
+	};
+
+	for ( auto it = this->global_faults.begin(); it != this->global_faults.end(); it++ ) {
+		if ( ft != L"ALL" ) 
+		{
+			wstring lft = (this->feature_types.find(it->first))->second;
+			if ( lft != ft ) {
+				continue;
+			}
+		}
+		if ( not ext ) {
+			add_line_triangles_not_ext( it );
+		} else {
+			add_line_triangles( it );
+		}
 	}
 
 	return ret;
 }
 
+// Direct methods to get different kinds of extruded lines.
+pydict ModelPython::get_lines() const { return this->filter_lines( true, L"ALL" ); }
+pydict ModelPython::get_not_extended_lines() const { return this->filter_lines( false, L"ALL" ); }
+
+pydict ModelPython::get_fracts( ) const { return this->filter_lines( true, L"FRACT" ); }
+pydict ModelPython::get_not_extended_fracts() const { return this->filter_lines( false, L"FRACT" ); }
+
+pydict ModelPython::get_faults( ) const { return this->filter_lines( true, L"FAULT" ); }
+pydict ModelPython::get_not_extended_faults() const { return this->filter_lines( false, L"FAULT" ); }
+
+pydict ModelPython::get_veins( ) const { return this->filter_lines( false, L"VEIN" ); }
+pydict ModelPython::get_not_extended_veins() const { return this->filter_lines( false, L"VEIN" ); }
 
 void ModelPython::set_matches( const pylist& matching ) {
 	this->clear_matches();
@@ -808,16 +839,13 @@ void ModelPython::set_matches( const pylist& matching ) {
 		const wstring& name1 = python::extract<wstring>(matching[i][0][0]);
 		const wstring& name2 = python::extract<wstring>(matching[i][0][1]);
 		vector<std::pair<int, int>> vmatch;
-
 		size_t mmatch = python::len(matching[i][1]);
-
 		for ( size_t j = 0; j < mmatch; j++ ) {
 			int a = python::extract<int>(matching[i][1][j][0]);
 			int b = python::extract<int>(matching[i][1][j][1]);
 
 			vmatch.push_back(std::make_pair(a, b));
 		}
-
 		cppmatching.push_back(std::make_tuple( std::make_tuple(name1, name2), vmatch ));
 	}
 	((Model *)this)->set_matches( cppmatching );