Index: DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/GeometryGenerator.cs
===================================================================
diff -u -r4052 -r4134
--- DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/GeometryGenerator.cs	(.../GeometryGenerator.cs)	(revision 4052)
+++ DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/GeometryGenerator.cs	(.../GeometryGenerator.cs)	(revision 4134)
@@ -19,6 +19,7 @@
 // Stichting Deltares and remain full property of Stichting Deltares at all times.
 // All rights reserved.
 
+using Deltares.DamEngine.Data.Standard;
 using System;
 using System.Collections.Generic;
 using System.Diagnostics;
@@ -451,7 +452,7 @@
                             point4.X = attachedCurveList[index2].GetEndPoint().X;
                             point4.Z = attachedCurveList[index2].GetEndPoint().Z;
 
-                            double angle = Routines2D.FindAngle(point1, point2, point3, point4);
+                            double angle = FindAngle(point1, point2, point3, point4);
 
                             if (angle < minAngle)
                             {
@@ -848,7 +849,7 @@
             // check if it is an inner loop
             for (var innerIndex = 0; innerIndex < newPointCount; innerIndex++)
             {
-                PointInPolygon location = Routines2D.CheckIfPointIsInPolygon(loop, newPolygon[innerIndex].X, newPolygon[innerIndex].Z);
+                PointInPolygon location = CheckIfPointIsInPolygon(loop, newPolygon[innerIndex].X, newPolygon[innerIndex].Z);
 
                 if (location == PointInPolygon.InsidePolygon)
                 {
@@ -863,7 +864,7 @@
             // check if it has an inner loop
             for (var innerIndex1 = 0; innerIndex1 < existingLoopPointCount; innerIndex1++)
             {
-                PointInPolygon location = Routines2D.CheckIfPointIsInPolygon(newLoop, polygon[innerIndex1].X, polygon[innerIndex1].Z);
+                PointInPolygon location = CheckIfPointIsInPolygon(newLoop, polygon[innerIndex1].X, polygon[innerIndex1].Z);
 
                 if (location == PointInPolygon.InsidePolygon)
                 {
@@ -898,6 +899,142 @@
         aNewSurfaceList.AddRange(newlyDetectedSurfaceList);
     }
 
+    /// <summary>
+    /// Find if points Ax,ay are between a boundary polygon.
+    /// </summary>
+    /// <param name="polygon">The polygon.</param>
+    /// <param name="x">The x.</param>
+    /// <param name="z">The z.</param>
+    /// <returns></returns>
+    private PointInPolygon CheckIfPointIsInPolygon(GeometryLoop polygon, double x, double z)
+    {
+        // This is done by calculating and adding the angles from the point to all points
+        // of the polygon (using the ArcTan2 function). If the sum of these angles is
+        // aproximate zero the point lies outside the polygon.
+        // If the sum in aproximate + or - 2Pi the point lies in the polygon;
+        // Performance tests have proven that this algorithm performs better than the 
+        // polynomial winding number algorithm described at:
+        // http://geomalgorithms.com/a03-_inclusion.html
+
+        const double epsilon = 1e-10;
+        var result = PointInPolygon.OutsidePolygon;
+
+        // add the last point at the end, cos the polygon must be closed, so the last node equals the first
+        var pointAdded = false;
+        if (polygon.CalcPoints.Count > 0)
+        {
+            if (!polygon.CalcPoints[0].LocationEquals(polygon.CalcPoints[polygon.CalcPoints.Count - 1]))
+            {
+                polygon.CalcPoints.Add(polygon.CalcPoints[0]);
+                pointAdded = true;
+            }
+        }
+
+        if (polygon.CalcPoints.Count > 2)
+        {
+            double deltaX = polygon[0].X - x;
+            double deltaZ = polygon[0].Z - z;
+            double absX = Math.Abs(deltaX);
+            double absZ = Math.Abs(deltaZ);
+            if ((absX < epsilon) && (absZ < epsilon))
+            {
+                UndoAddIfNeeded(polygon, pointAdded);
+                return PointInPolygon.OnPolygonEdge;
+            }
+
+            double al1 = Math.Atan2(deltaZ, deltaX);
+            double som = 0;
+            var index = 1;
+
+            while (index < polygon.CalcPoints.Count)
+            {
+                deltaX = polygon[index].X - x;
+                deltaZ = polygon[index].Z - z;
+                absX = Math.Abs(deltaX);
+                absZ = Math.Abs(deltaZ);
+                if ((absX < epsilon) && (absZ < epsilon))
+                {
+                    UndoAddIfNeeded(polygon, pointAdded);
+                    return PointInPolygon.OnPolygonEdge;
+                }
+
+                double al2 = Math.Atan2(deltaZ, deltaX);
+                double al3 = al2 - al1;
+
+                if (al3 < -Math.PI)
+                {
+                    al3 = al3 + (2.0 * Math.PI);
+                }
+
+                if (al3 > Math.PI)
+                {
+                    al3 = al3 - (2.0 * Math.PI);
+                }
+
+                if (((Math.PI - al3) < epsilon) || ((Math.PI + al3) < epsilon))
+                {
+                    UndoAddIfNeeded(polygon, pointAdded);
+                    return PointInPolygon.OnPolygonEdge;
+                }
+
+                som = som + al3;
+                al1 = al2;
+                index++;
+            }
+
+            if ((som > (1.9 * Math.PI)) || (som < (-1.9 * Math.PI)))
+            {
+                result = PointInPolygon.InsidePolygon;
+            }
+            else
+            {
+                result = PointInPolygon.OutsidePolygon;
+            }
+        }
+
+        UndoAddIfNeeded(polygon, pointAdded);
+        return result;
+    }
+
+    private void UndoAddIfNeeded(GeometryLoop polygon, bool needed)
+    {
+        if (needed)
+        {
+            polygon.CalcPoints.Remove(polygon.CalcPoints[polygon.CalcPoints.Count - 1]);
+        }
+    }
+
+    /// <summary>
+    ///  Determines angle between 2 lines.
+    ///  Clockwise Input: point1, commonpoint, point2, commonpoint, normalvect(ex [0,0,-1]), validate  
+    /// </summary>
+    /// <param name="line1Point1"></param>
+    /// <param name="line1Point2"></param>
+    /// <param name="line2Point1"></param>
+    /// <param name="line2Point2"></param>
+    /// <returns></returns>
+    private double FindAngle(Point2D line1Point1, Point2D line1Point2, Point2D line2Point1,
+                                   Point2D line2Point2)
+    {
+        double x1 = line1Point2.X - line1Point1.X;
+        double z1 = line1Point2.Z - line1Point1.Z;
+
+        double x2 = line2Point2.X - line2Point1.X;
+        double z2 = line2Point2.Z - line2Point1.Z;
+
+        double angle1 = GeneralMathRoutines.RadianToDegree(Math.Atan2(z1, x1));
+        double angle2 = GeneralMathRoutines.RadianToDegree(Math.Atan2(z2, x2));
+
+        double angle = angle2 - angle1;
+
+        if (angle < 0)
+        {
+            angle += 360;
+        }
+
+        return angle;
+    }
+
     #region Nested type: DirectionCurve
 
     internal struct DirectionCurve