From a9f55128d4159af445ec48c32351942eccd59e49 Mon Sep 17 00:00:00 2001
From: Kai Bao <kai.bao@sintef.no>
Date: Tue, 2 Jun 2015 17:05:55 +0200
Subject: [PATCH] adding the function to find the intersection point.

This apply to find the intersection point of the line and a line
segmention. Will be used in the shear multipler calculation of with
PLYSHLOG.

Solver is not the best place to put this function, while need suggestion
and fixed later.
---
 .../fullyimplicit/BlackoilPolymerModel.hpp    | 10 +++++++
 .../BlackoilPolymerModel_impl.hpp             | 30 +++++++++++++++++++
 2 files changed, 40 insertions(+)

diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
index ccdc3a04f..6a5eda68f 100644
--- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
+++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
@@ -261,6 +261,16 @@ namespace Opm {
                              int nw) const;
 
 
+        struct Point2D {
+            double x;
+            double y;
+        };
+
+        /// Finding the intersection point of a line segment and a line.
+        /// return true, if found.
+        bool findIntersection (Point2D line_segment1[2], Point2D line2[2], Point2D& intersection_point);
+
+
 
     };
 
diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
index 71d97a68c..dea3eeee8 100644
--- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
+++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
@@ -587,6 +587,36 @@ namespace Opm {
         return polymer_props_ad_.polymerWaterVelocityRatio(state.concentration);
     }
 
+    template<class Grid>
+    bool
+    BlackoilPolymerModel<Grid>::findIntersection (Point2D line_segment1[2], Point2D line2[2], Point2D& intersection_point){
+
+        const double x1 = line_segment1[0].x;
+        const double y1 = line_segment1[0].y;
+        const double x2 = line_segment1[1].x;
+        const double y2 = line_segment1[1].y;
+
+        const double x3 = line2[0].x;
+        const double y3 = line2[0].y;
+        const double x4 = line2[1].x;
+        const double y4 = line2[1].y;
+
+        const double d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4);
+
+        if (d == 0.) { return false; }
+
+        const double x = ((x3 - x4) * (x1 * y2 - y1 * x2) - (x1 - x2) * (x3 * y4 - y3 * x4)) / d;
+        const double y = ((y3 - y4) * (x1 * y2 - y1 * x2) - (y1 - y2) * (x3 * y4 - y3 * x4)) / d;
+
+        if( x >= std::min(x1,x2) && x <= std::max(x1,x2) ){
+            intersection_point.x = x;
+            intersection_point.y = y;
+            return true;
+        } else {
+            return false;
+        }
+    }
+
 } // namespace Opm
 
 #endif // OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED