Skip to content

Commit

Permalink
Changed the way faults are represented.
Browse files Browse the repository at this point in the history
  • Loading branch information
Ricardo Serrano committed Apr 26, 2017
1 parent 981760e commit 3f7c095
Show file tree
Hide file tree
Showing 2 changed files with 54 additions and 36 deletions.
64 changes: 43 additions & 21 deletions geomodelr/cpp/match.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -207,57 +207,79 @@ vector<triangle> test_start( const vector<point3>& pa, const vector<point3>& pb,
auto next_a = [&] ( int i ) -> int {
return size_t(i+1) < pa.size() ? i+1:-1;
};

auto next_b = [&] ( int i ) -> int {
if ( binv ) {
return i-1;
} else {
return size_t(i+1) < pb.size() ? i+1:-1;
}
};

auto b_idx = [&] ( int i ) -> int {
return pa.size() + i;
};

auto norm = [] ( const point3& v ) {
return std::sqrt( geometry::dot_product( v, v ) );
};

auto normalize = [norm] ( point3& v ) {
geometry::divide_value( v, norm( v ) );
};

auto vector_between = [normalize] ( const point3& pn, const point3& p ) {
point3 v = pn;
geometry::subtract_point( v, p );
normalize(v);
return v;
};

double length_a = 0.0;
double length_b = 0.0;

for ( size_t i = 0; i < pa.size()-1; i++ ) {
length_a += geometry::distance(pa[i], pa[i+1]);
}

for ( size_t i = 0; i < pb.size()-1; i++ ) {
length_b += geometry::distance(pb[i], pb[i+1]);
}

vector<triangle> result;

int n_a, n_b;
int c_a = 0;
int c_b = binv ? pb.size()-1 : 0;
n_a = next_a(c_a);
n_b = next_b(c_b);

double prop_a = 0.0;
double prop_b = 0.0;

while ( n_a != -1 || n_b != -1 ) { // Get next points if advance a, and if advance b, and check if they are not -1 both.
// Get triangle if advance a, centered in b.
point3 v = pa[c_a];
geometry::subtract_point( v, pb[c_b] );
geometry::divide_value( v, std::sqrt( geometry::dot_product( v, v ) ) );
double next_prop_a = std::numeric_limits<double>::infinity();
double next_prop_b = std::numeric_limits<double>::infinity();

double angle_a = 0;
double angle_b = 0;
if ( n_a != -1 ) {
// Get vector if advance a.
point3 vn = pa[n_a];
geometry::subtract_point( vn, pb[c_b] );
geometry::divide_value( vn, std::sqrt( geometry::dot_product( vn, vn ) ) );

// Get angle if advance a.
angle_a = std::acos(geometry::dot_product( v, vn ));
next_prop_a = prop_a + (geometry::distance( pa[n_a], pa[c_a] )/length_a);
}

if ( n_b != -1 ) {
// Get vector if advance b.
point3 vn = pb[n_b];
geometry::subtract_point( vn, pa[c_a] );
geometry::divide_value( vn, std::sqrt( geometry::dot_product( vn, vn ) ) );
// Get angle if advance b.
angle_b = std::acos(-geometry::dot_product( v, vn ));
next_prop_b = prop_b + (geometry::distance( pb[n_b], pb[c_b] )/length_b);
}
std::cerr << "angles " << angle_a << " " << angle_b << "\n";

// Check angle if advance b.
if ( angle_a < angle_b ) {
if ( next_prop_b < next_prop_a ) {
result.push_back( triangle( c_a, b_idx(c_b), b_idx(n_b) ) );
c_b = n_b;
prop_b = next_prop_b;
} else {
result.push_back( triangle( c_a, n_a, b_idx(c_b) ) );
c_a = n_a;
prop_a = next_prop_a;
}

n_a = next_a(c_a);
n_b = next_b(c_b);
}
Expand Down
26 changes: 11 additions & 15 deletions geomodelr/test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ def test_faultplane_for_lines(self):
la = [(0, 0, 0), (1,1,0), (2,2,0), (3,3,0)]
lb = [(0, 3, 1), (1,2,1), (2,1,1), (3,0,1)]
# This is a weird testcase that should work. Star in two directions.
self.assertEqual(cpp.faultplane_for_lines(la, lb),[(0, 1, 4), (1, 2, 4), (2, 4, 5), (2, 5, 6), (2, 3, 6), (3, 6, 7)])
self.assertEqual(cpp.faultplane_for_lines(la, lb),[(0, 1, 4), (1, 4, 5), (1, 2, 5), (2, 5, 6), (2, 3, 6), (3, 6, 7)])

# This is a weird testcase.
# Open circles in different directions.
Expand All @@ -69,15 +69,12 @@ def test_faultplane_for_lines(self):
# One starts when the other finishes.
la = [(-1,-2,0), (-2,-1,0), (-2,1,0), (-1, 2,0), (1, 2,0), (2, 1,0), (2,-1,0), (1, -2, 0)]
lb = [(-1, 2, 1), (-2, 3, 1), (-2, 5, 1), (-1, 6, 1), (1, 6, 1), (2, 5, 1), (2, 3, 1), (1, 2, 1)]
faults = cpp.faultplane_for_lines(la, lb)
expected = [(0,1, 8), (1, 8, 9), (1, 2, 9),
(2, 9,10), (2, 3, 10), (3, 10, 11),
(3,11,12), (3, 4, 12), (4, 12, 13),
(4, 5,13), (5, 13, 14), (5, 6, 14),
(6,14,15), (6, 7, 15)]
print faults
print expected
self.assertEqual(faults, expected)

self.assertEqual(cpp.faultplane_for_lines(la, lb), [(0, 1, 8), (1, 8, 9), (1, 2, 9),
(2, 9, 10), (2, 3, 10), (3, 10, 11),
(3, 4, 11), (4, 11, 12), (4, 5, 12),
(5, 12, 13), (5, 6, 13), (6, 13, 14),
(6, 7, 14), (7, 14, 15)])
# Case from Aburra's Valley
la = [(827675.59327569162, 1165500.0, 1628.0922045583072),
(827765.51647690905, 1165500.0, 1378.975507981105),
Expand All @@ -97,11 +94,10 @@ def test_faultplane_for_lines(self):
(827315.46193332784, 1169500.0, -284.6906919181347),
(840876.00000000093, 1169500.0, -30000)]

self.assertEqual(cpp.faultplane_for_lines(la, lb), [(0, 8, 9), (0, 9, 10), (0, 1, 10),
(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, 15, 16), (6, 7, 16)])
self.assertEqual(cpp.faultplane_for_lines(la, lb), [(0, 8, 9), (0, 1, 9), (1, 9, 10), (1, 2, 10),
(2, 10, 11), (2, 3, 11), (3, 11, 12), (3, 4, 12),
(4, 12, 13), (4, 5, 13), (5, 13, 14), (5, 6, 14),
(6, 14, 15), (6, 7, 15), (7, 15, 16)])

# Test that you can create cross sections and that bugs throw something in python (and not segfault).
def test_sections(self):
Expand Down

0 comments on commit 3f7c095

Please sign in to comment.