moving the intersection calculation to seperate point2D class.

under the namespae Opm::detail

CMakeLists_files.cmake

@@ -85,6 +85,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/polymer/SimulatorPolymer.hpp
 	opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp
 	opm/polymer/TransportSolverTwophaseCompressiblePolymer.hpp
+	opm/polymer/Point2D.hpp
     opm/polymer/TransportSolverTwophasePolymer.hpp
     opm/polymer/fullyimplicit/PolymerPropsAd.hpp
     opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp

opm/polymer/Point2D.hpp

@@ -0,0 +1,102 @@
+/*
+  Copyright 2015 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_POINT2DINTERSECTION_INCLUDED
+#define OPM_POINT2DINTERSECTION_INCLUDED
+
+namespace Opm {
+
+    namespace detail {
+        class Point2D
+        {
+            public:
+                Point2D(const double xi, const double yi)
+                    : x_(xi),
+                      y_(yi) {
+
+                      }
+
+                Point2D()
+                {
+                }
+
+                const double getX() const
+                {
+                    return x_;
+                }
+
+                const double getY() const
+                {
+                    return y_;
+                }
+
+                void setX(const double x)
+                {
+                    x_ = x;
+                }
+
+                void setY(const double y){
+                    y_ = y;
+                }
+
+                /// Finding the intersection point of a line segment and a line.
+                /// return true, if found.
+                static bool findIntersection(Point2D line_segment1[2], Point2D line2[2], Point2D& intersection_point)
+                {
+
+                    const double x1 = line_segment1[0].getX();
+                    const double y1 = line_segment1[0].getY();
+                    const double x2 = line_segment1[1].getX();
+                    const double y2 = line_segment1[1].getY();
+
+                    const double x3 = line2[0].getX();
+                    const double y3 = line2[0].getY();
+                    const double x4 = line2[1].getX();
+                    const double y4 = line2[1].getY();
+
+                    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.setX(x);
+                        intersection_point.setY(y);
+                        return true;
+                    } else {
+                       return false;
+                    }
+                }
+
+            private:
+                double x_;
+                double y_;
+
+        }; // class Point2D
+
+    } // namespace detail
+
+} // namespace Opm
+
+#endif // OPM_POINT2DINTERSECTION_INCLUDED
+

opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp

@@ -285,16 +285,6 @@ namespace Opm {
                              int nc,
                              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);
-
         /// Computing the shear multiplier based on the water velocity/shear rate with PLYSHLOG keyword
         bool computeShearMultLog(std::vector<double>& water_vel, std::vector<double>& visc_mult, std::vector<double>& shear_mult);
 

opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp

@@ -25,6 +25,7 @@
 #define OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED
 
 #include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
+#include <opm/polymer/Point2D.hpp>
 
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
@@ -892,40 +893,6 @@ namespace Opm {
         extraAddWellEq(state, xw, cq_ps, cmix_s, cqt_is, well_cells);
     }
 
-
-    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;
-        }
-    }
-
     template<class Grid>
     bool
     BlackoilPolymerModel<Grid>::computeShearMultLog(std::vector<double>& water_vel, std::vector<double>& visc_mult, std::vector<double>& shear_mult)
@@ -988,20 +955,20 @@ namespace Opm {
             }
 
             if (foundSegment == true) {
-                Point2D lineSegment[2];
-                lineSegment[0] = Point2D{logShearWaterVel[iIntersection], logShearVRF[iIntersection]};
-                lineSegment[1] = Point2D{logShearWaterVel[iIntersection + 1], logShearVRF[iIntersection + 1]};
+                detail::Point2D lineSegment[2];
+                lineSegment[0] = detail::Point2D{logShearWaterVel[iIntersection], logShearVRF[iIntersection]};
+                lineSegment[1] = detail::Point2D{logShearWaterVel[iIntersection + 1], logShearVRF[iIntersection + 1]};
 
-                Point2D line[2];
-                line[0] = Point2D{0, logWaterVelO};
-                line[1] = Point2D{logWaterVelO, 0};
+                detail::Point2D line[2];
+                line[0] = detail::Point2D{0, logWaterVelO};
+                line[1] = detail::Point2D{logWaterVelO, 0};
 
-                Point2D intersectionPoint;
+                detail::Point2D intersectionPoint;
 
-                bool foundIntersection = findIntersection(lineSegment, line, intersectionPoint);
+                bool foundIntersection = detail::Point2D::findIntersection(lineSegment, line, intersectionPoint);
 
                 if (foundIntersection) {
-                    shear_mult[i] = std::exp(intersectionPoint.y);
+                    shear_mult[i] = std::exp(intersectionPoint.getY());
                 } else {
                     std::cerr << " failed in finding the solution for shear-thinning multiplier " << std::endl;
                     return false; // failed in finding the solution.