using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Runtime.Serialization;
using Deltares.Geometry;
using Deltares.Geotechnics;
using Deltares.Geotechnics.GeotechnicalGeometry;
using Deltares.Geotechnics.Soils;
using Deltares.Mathematics;
using Deltares.Stability;
using Deltares.Standard.Extensions;

//MSTab next generation namespaces

namespace Deltares.MStab.IO.Classic
{
    [Serializable]
    public class GeometryException : Exception
    {
        //
        // For guidelines regarding the creation of new exception types, see
        //    http://msdn.microsoft.com/library/default.asp?url=/library/en-us/cpgenref/html/cpconerrorraisinghandlingguidelines.asp
        // and
        //    http://msdn.microsoft.com/library/default.asp?url=/library/en-us/dncscol/html/csharp07192001.asp
        //

        public GeometryException() {}

        public GeometryException(string message)
            : base(message) {}

        public GeometryException(string message, Exception inner)
            : base(message, inner) {}

        protected GeometryException(
            SerializationInfo info,
            StreamingContext context)
            : base(info, context) {}
    }

    internal class Geometry
    {
        private const double Epsilon = 0.005;
        private Dictionary<long, List<long>> curveDictionary = null; //Dictionary<oldcurve, replaced by new curves>
        private Dictionary<long, long> surfaceDictionary = null; //Dictionary<new surface, points at old layer index>

        public Dictionary<long, List<long>> CurveDictionary
        {
            get
            {
                return curveDictionary;
            }
        }

        public Dictionary<long, long> SurfaceDictionary
        {
            get
            {
                return surfaceDictionary;
            }
        }

        public void GetGeometryData(int handle,
                                    StabilityModel stabilityModel,
                                    ref Dictionary<GeometrySurface, int> aSurfaceDictionary,
                                    ref Dictionary<long, List<GeometryCurve>> aBoundaryDictionary,
                                    ref Dictionary<GeometryCurve, List<GeometryCurve>> aCurveDictionary,
                                    ref List<GeometryPoint> waternetGeometryPoints)
        {
            try
            {
                if (stabilityModel == null)
                {
                    throw new GeometryException("SoilGeometryModel Object not instantiated");
                }

                if (stabilityModel.Geometry == null)
                {
                    throw new GeometryException("GeometryData Object not instantiated");
                }

                if (stabilityModel.SoilModel == null)
                {
                    throw new GeometryException("SoilMaterials Object not instantiated");
                }

                //Step 1: Load Classic Geometry
                List<ClassicGeometryFunctionDefinitions.GeometryPoint> classicGeometryPoints = null;
                List<ClassicGeometryFunctionDefinitions.GeometryCurve> classicGeometryCurves = null;
                ClassicGeometryFunctionDefinitions.GeometryBoundaries classicGeometryBoundaries = null;
                ClassicGeometryFunctionDefinitions.GeometryLayers classicGeometryLayers = null;
                ClassicGeometryFunctionDefinitions.GeometryParam classicGeometryConstants = null;

                try
                {
                    GetClassicGeometry.GetClassicGeometryPoints(handle, ref classicGeometryPoints);
                    GetClassicGeometry.GetClassicGeometryCurves(handle, ref classicGeometryCurves);
                    GetClassicGeometry.GetClassicGeometryBoundaries(handle, ref classicGeometryBoundaries);
                    GetClassicGeometry.GetClassicGeometryLayers(handle, ref classicGeometryLayers);
                    GetClassicGeometry.GetClassicGeometryConstants(handle, ref classicGeometryConstants);
                }
                catch (Exception e)
                {
                    Debug.Assert(false, "Could not receive classic dataobjects" + e.ToString());
                }

                // Step 1a: Set boundaries of Geometry such that everything fits in and no points on boundary
                stabilityModel.Geometry.Left = Double.MinValue;
                stabilityModel.Geometry.Right = Double.MaxValue;
                stabilityModel.Geometry.Bottom = Double.MinValue;

                //Step 2 Add geometry points into the new geometry
                waternetGeometryPoints = new List<GeometryPoint>();
                var pointDictionary = new Dictionary<long, long>();
                int counter = 0;
                try
                {
                    for (int index = 0; index < classicGeometryPoints.Count; index++)
                    {
                        //Make sure all geometry points are unique
                        List<GeometryPoint> pointList = stabilityModel.Geometry.Points.FindAll(x => ((x.X == classicGeometryPoints[index].x) &&
                                                                                              (x.Y == classicGeometryPoints[index].z) &&
                                                                                              (x.Z == classicGeometryPoints[index].y)));

                        if ((pointList.Count == 0) || (classicGeometryPoints[index].pointType == (int) PointType.ptWaterPoint))
                        {
                            var geoPoint = new GeometryPoint();
                            geoPoint.Name = String.Format("{0}", counter++);

                            geoPoint.X = classicGeometryPoints[index].x;
                            geoPoint.Y = classicGeometryPoints[index].z;
                            geoPoint.Z = classicGeometryPoints[index].y;

                            if (classicGeometryPoints[index].pointType == (int) PointType.ptWaterPoint)
                            {
                                waternetGeometryPoints.Add(geoPoint);
                            }
                            else
                            {
                                stabilityModel.Geometry.Create(geoPoint);
                            }
                        }
                    }
                }
                catch (Exception e)
                {
                    Debug.Assert(false, "Error while adding classic geometrypoints in new geometry" + e.ToString());
                }

                //Step 3 Add Geometry curves in the new geometry
                var classic2NewCurves = new Dictionary<long, GeometryCurve>();
                try
                {
                    for (int index = 0; index < classicGeometryCurves.Count; index++)
                    {
                        var curve = new GeometryCurve();
                        curve.Name = String.Format("Curve {0}", index + 1);

                        if (classicGeometryPoints[classicGeometryCurves[index].pointIndices[0] - 1].pointType != (int) PointType.ptWaterPoint)
                        {
                            curve.HeadPoint = stabilityModel.Geometry.Points.Find(x => ((x.X == classicGeometryPoints[classicGeometryCurves[index].pointIndices[0] - 1].x) &&
                                                                                 (x.Y == classicGeometryPoints[classicGeometryCurves[index].pointIndices[0] - 1].z) &&
                                                                                 (x.Z == classicGeometryPoints[classicGeometryCurves[index].pointIndices[0] - 1].y)));
                            curve.EndPoint = stabilityModel.Geometry.Points.Find(x => ((x.X == classicGeometryPoints[classicGeometryCurves[index].pointIndices[1] - 1].x) &&
                                                                                (x.Y == classicGeometryPoints[classicGeometryCurves[index].pointIndices[1] - 1].z) &&
                                                                                (x.Z == classicGeometryPoints[classicGeometryCurves[index].pointIndices[1] - 1].y)));
                            stabilityModel.Geometry.Create(curve);

                            classic2NewCurves.Add(index, curve); //Synchronize ID's (hypothetically it is possible that ref id's are not incremental)
                        }
                    }
                }
                catch (Exception e)
                {
                    Debug.Assert(false, "Error in conversion classic2new curves" + e.ToString());
                }

                try
                {
                    CreateSoilSurface(classicGeometryLayers, classicGeometryBoundaries, classic2NewCurves,
                                      ref stabilityModel, ref aSurfaceDictionary, ref aBoundaryDictionary, ref aCurveDictionary);
                }
                catch (Exception e)
                {
                    Debug.Assert(false, "Error while creating the soil surface" + e.ToString());
                }

                //Step 4: Create complete boundary loop                
                var geometryBoundaryLoop = new GeometryLoop();
                geometryBoundaryLoop.Name = "Geometry boundary loop";

                var tempLoop = new GeometryLoop();
                try
                {
                    foreach (var curve in stabilityModel.Geometry.Curves)
                    {
                        int curveCounter = 0;
                        foreach (var x in stabilityModel.Geometry.Surfaces)
                        {
                            bool loopbelongstocurve = x.OuterLoop.CurveList.Find(y => y == curve) != null;

                            if (loopbelongstocurve)
                            {
                                curveCounter += 1;
                            }
                        }
                        if (curveCounter == 1)
                        {
                            tempLoop.CurveList.Add(curve);
                        }
                    }
                }
                catch (Exception e)
                {
                    Debug.Assert(false, "Could not create the boundary loop" + e.ToString());
                }

                //Curves in geometryloop should be counter clockwise                    
                geometryBoundaryLoop.CurveList.Add(tempLoop.CurveList[0]);
                var curveref = tempLoop.CurveList[0]; //Initial curve
                tempLoop.CurveList.Remove(tempLoop.CurveList[0]);

                GeometryPoint refBeginPoint = null;

                if (curveref.HeadPoint.X < curveref.EndPoint.X)
                {
                    refBeginPoint = curveref.HeadPoint;
                }
                else if ((curveref.HeadPoint.X == curveref.EndPoint.X) && (curveref.HeadPoint.Z < curveref.EndPoint.Z))
                {
                    refBeginPoint = curveref.HeadPoint;
                }
                else
                {
                    refBeginPoint = curveref.EndPoint;
                }

                try
                {
                    int clockwise = 0;
                    bool curveAdded = true;

                    while (curveAdded && tempLoop.CurveList.Count > 0)
                    {
                        for (int j = 0; j < tempLoop.CurveList.Count; j++)
                        {
                            curveAdded = false;
                            var curvecmp = (GeometryCurve) tempLoop.CurveList[j];

                            if (curvecmp != curveref)
                            {
                                if ((refBeginPoint == curvecmp.EndPoint) || (refBeginPoint == curvecmp.HeadPoint))
                                {
                                    if (refBeginPoint == curvecmp.EndPoint)
                                    {
                                        refBeginPoint = curvecmp.HeadPoint;
                                        clockwise--;
                                    }
                                    else
                                    {
                                        refBeginPoint = curvecmp.EndPoint;
                                        clockwise++;
                                    }
                                    curveref = curvecmp;

                                    geometryBoundaryLoop.CurveList.Add(curvecmp);
                                    tempLoop.CurveList.Remove(curvecmp);
                                    curveAdded = true;
                                    break;
                                }
                            }
                        }
                    }
                    if (clockwise > 0)
                    {
                        //make counterclockwise
                        geometryBoundaryLoop.CurveList.Reverse();
                    }
                }
                catch (Exception e)
                {
                    Debug.Assert(false, "Could not change order boundary" + e.ToString());
                }
                stabilityModel.Geometry.Create(geometryBoundaryLoop);

                // step 5: Set geometry bounds
                if (stabilityModel.Geometry.Points.Count > 0)
                {
                    double left = Double.MaxValue;
                    double right = Double.MinValue;
                    double bottom = Double.MaxValue;

                    foreach (var loop in stabilityModel.Geometry.Loops)
                    {
                        foreach (var point in loop.Points)
                        {
                            if (point.X < left)
                            {
                                left = point.X;
                            }

                            if (point.X > right)
                            {
                                right = point.X;
                            }

                            if (point.Z < bottom)
                            {
                                bottom = point.Z;
                            }
                        }
                    }

                    stabilityModel.Geometry.Left = left;
                    stabilityModel.Geometry.Right = right;
                    stabilityModel.Geometry.Bottom = bottom;
                }
                else
                {
                    stabilityModel.Geometry.Left = GeometryConstants.DefaultLeftLimitGeometry;
                    stabilityModel.Geometry.Right = GeometryConstants.DefaultRightLimitGeometry;
                    stabilityModel.Geometry.Bottom = GeometryConstants.DefaultBottomLimitGeometry;
                }

                // step 6: Delete loose lines
                stabilityModel.Geometry.DeleteLooseCurves();
                stabilityModel.Geometry.DeleteLoosePoints();
            }
            catch (GeometryException e)
            {
                Debug.Assert(false, e.ToString());
            }
        }

        /// <summary>
        /// RegenerateCurveList
        /// After the CreateSplitCurves function, new curves are created. The whole geometry needs to be rebuild:
        /// Old curves pointers should be removed, old curve objects should be deleted and possible refId conflicts should be resolved 
        /// For the renumbering of curves, the direction of the polygon is important.
        /// </summary>
        /// <param name="aPublisher">Required for creating new Geometry Objects</param>
        /// <param name="aCurveDictionary">LookupTable<Old curveid, list of new curve ids></param>
        /// <param name="aSurfaceList">List of surfaces (containing a list of top boundary curves, bottomboundary curves</param>
        /// <param name="aMStabProject">New Geometry objects</param>
        private void RegenerateCurveList(Dictionary<GeometryCurve, List<GeometryCurve>> aCurveDictionary, ref List<List<List<GeometryCurve>>> aSurfaceList, ref StabilityModel aMStabProject)
        {
            foreach (var topAndBotCurves in aSurfaceList)
            {
                for (int topCurve = 0; topCurve < topAndBotCurves[0].Count; topCurve++)
                {
                    if (aCurveDictionary.ContainsKey(topAndBotCurves[0][topCurve]))
                    {
                        GeometryCurve curveId = topAndBotCurves[0][topCurve]; //Remember curve
                        topAndBotCurves[0].RemoveAt(topCurve);

                        //Replace curves
                        for (int i = 0; i < aCurveDictionary[curveId].Count; i++)
                        {
                            topAndBotCurves[0].Insert(topCurve + i, aCurveDictionary[curveId][i]);
                        }
                    }
                }

                for (int bottomCurve = 0; bottomCurve < topAndBotCurves[1].Count; bottomCurve++)
                {
                    if (aCurveDictionary.ContainsKey(topAndBotCurves[1][bottomCurve]))
                    {
                        GeometryCurve curveId = topAndBotCurves[1][bottomCurve]; //Remember curve
                        topAndBotCurves[1].RemoveAt(bottomCurve);

                        //Replace curves (reverse)
                        for (int i = 0; i < aCurveDictionary[curveId].Count; i++)
                        {
                            topAndBotCurves[1].Insert(bottomCurve + i, aCurveDictionary[curveId][aCurveDictionary[curveId].Count - 1 - i]);
                        }
                    }
                }
            }
            //Remove all unused curves
            foreach (var item in aCurveDictionary)
            {
                aMStabProject.Geometry.Curves.Remove(aMStabProject.Geometry.Curves.Find(x => x == item.Key));
            }
        }

        private void CreateSoilSurface(ClassicGeometryFunctionDefinitions.GeometryLayers classicGeometryLayers,
                                       ClassicGeometryFunctionDefinitions.GeometryBoundaries classicGeometryBoundaries,
                                       Dictionary<long, GeometryCurve> classic2NewCurves,
                                       ref StabilityModel project,
                                       ref Dictionary<GeometrySurface, int> surfaceDictionary, //Required for waternet
                                       ref Dictionary<long, List<GeometryCurve>> boundaryDictionary, //Required for waternet
                                       ref Dictionary<GeometryCurve, List<GeometryCurve>> curveDictionary)
        {
            if (surfaceDictionary == null)
            {
                surfaceDictionary = new Dictionary<GeometrySurface, int>();
            }

            if (boundaryDictionary == null)
            {
                boundaryDictionary = new Dictionary<long, List<GeometryCurve>>();
            }

            List<List<List<GeometryCurve>>> surfaceList = GetSurfaces(classicGeometryLayers, classicGeometryBoundaries, classic2NewCurves, boundaryDictionary, project);

            curveDictionary = CreateSplitCurves(surfaceList, ref project);

            RegenerateCurveList(curveDictionary, ref surfaceList, ref project);
            int classicLayerIndex = -1;
            foreach (var surface in surfaceList)
            {
                //Keep track of old layer index
                classicLayerIndex++;

                List<List<GeometryCurve>> multipleSurfaces = DetectMultipleSurfaces(surface[0], surface[1], project);
                for (int i = 0; i < multipleSurfaces.Count; i++)
                {
                    List<GeometryCurve> newSurface = multipleSurfaces[i];

                    RemoveOpenEndCurves(project, newSurface);

                    if (newSurface.Count > 0)
                    {
                        GeometryLoop geometryLoop = CreateGeometryLoop(project, newSurface);
                        if (geometryLoop.HasArea())
                        {
                            project.Geometry.Create(geometryLoop);
                            //create surface
                            var geoSurface = new GeometrySurface();
                            geoSurface.SetOuterLoop(geometryLoop);
                            project.Geometry.Create(geoSurface);

                            //Add soil material to surface
                            var soilsurface = new SoilLayer2D();
                            // For some reason, old sti files can handle soils regardless of their case so do that here too.
                            soilsurface.Soil = project.SoilModel.Soils.Find(x => x.Description.ToLower().Equals(classicGeometryLayers.layers[classicLayerIndex].materialOfLayer.ToLower()));
                            soilsurface.GeometrySurface = geoSurface;
                            project.Create(soilsurface);
                            //Translate new Surfaces into old layers
                            surfaceDictionary.Add(geoSurface, classicLayerIndex);
                        }
                    }
                }
            }
        }

        /// <summary>
        /// GetSurfaces
        /// -Retrieve the classic Layer layout
        /// -Return an index indepedent (refIds instead of indices) surfaces structure:
        /// -Each surface consists of two lists TopBoundary list and BottomBoundary list (no polygons formed yet)
        /// -The top boundary is ordered from left to right
        /// -the bottom boundary is ordered from right to left
        /// </summary>
        /// <param name="aClassicGeometryLayers">Classical Geometry layer layout with absolute indices</param>
        /// <param name="aClassicGeometryBoundaries">TopBoundaries and Bottomboundaries per layer</param>
        /// <param name="aMStabProject">New Geometry layout with surfaces</param>
        /// <returns>Returns a list of Botom/Top curvelists per Surface</returns>
        ///      
        private List<List<List<GeometryCurve>>> GetSurfaces(ClassicGeometryFunctionDefinitions.GeometryLayers aClassicGeometryLayers,
                                                            ClassicGeometryFunctionDefinitions.GeometryBoundaries aClassicGeometryBoundaries,
                                                            Dictionary<long, GeometryCurve> aClassic2NewCurves,
                                                            Dictionary<long, List<GeometryCurve>> aBoundaries,
                                                            StabilityModel aMStabProject)
        {
            var surfaces = new List<List<List<GeometryCurve>>>();

            int boundaryIndex = 0;
            for (int layerIndex = 0; layerIndex < aClassicGeometryLayers.numberOfLayers; layerIndex++)
            {
                var topAndBottomCurves = new List<List<GeometryCurve>>();

                //Top side  
                boundaryIndex = aClassicGeometryLayers.layers[layerIndex].boundaryTopLayer;
                int numberCurvesTop = aClassicGeometryBoundaries.boundaries[boundaryIndex].numberOfCurvesOnBoundary;
                var curveListTop = new List<GeometryCurve>();
                int absIndex = 0;
                for (int curvesTop = 0; curvesTop < numberCurvesTop; curvesTop++)
                {
                    absIndex = aClassicGeometryBoundaries.boundaries[boundaryIndex].curvesIndices[curvesTop] - 1;
                    curveListTop.Add(aClassic2NewCurves[absIndex]);
                }

                //Check directions: TopBoundary Left - Right
                if (aMStabProject.Geometry.Curves.Find(x => x == curveListTop[0]).HeadPoint.X >
                    aMStabProject.Geometry.Curves.Find(x => x == curveListTop[curveListTop.Count - 1]).EndPoint.X)
                {
                    curveListTop.Reverse();
                }

                aBoundaries[boundaryIndex] = curveListTop;

                topAndBottomCurves.Add(curveListTop);
                //Bottom side                    
                var curveListBottom = new List<GeometryCurve>();
                boundaryIndex = aClassicGeometryLayers.layers[layerIndex].boundaryBottomLayer;
                int numberCurvesBottom = aClassicGeometryBoundaries.boundaries[boundaryIndex].numberOfCurvesOnBoundary;
                for (int curvesBottom = 0; curvesBottom < numberCurvesBottom; curvesBottom++)
                {
                    absIndex = aClassicGeometryBoundaries.boundaries[boundaryIndex].curvesIndices[curvesBottom] - 1;
                    curveListBottom.Add(aClassic2NewCurves[absIndex]);
                }

                //Check directions: BottomBoundary Right - Left
                if (aMStabProject.Geometry.Curves.Find(x => x == curveListBottom[0]).HeadPoint.X <
                    aMStabProject.Geometry.Curves.Find(x => x == curveListBottom[curveListBottom.Count - 1]).EndPoint.X)
                {
                    curveListBottom.Reverse();
                }
                topAndBottomCurves.Add(curveListBottom);

                aBoundaries[boundaryIndex] = curveListBottom;

                surfaces.Add(topAndBottomCurves);
            }

            return surfaces;
        }

        /// <summary>
        /// CreateSplitCurves
        /// Detects intersections between TopBoundaries and bottomboundaries
        /// If an intersection occures the bottomboundary curves will be cut into smaller pieces.
        /// The refid of the bottomboundary acts as lookuptable, which points to the newly created curves.
        /// Newly created curves are added into the StabilityModel. Due to this, classic indexing is pretty much worthless now.             
        /// </summary>
        /// <param name="aPublisher">Publisher, required for the creation of new geometry objects</param>
        /// <param name="aSurfaceList">List of Botom/Top curvelists per Surface</param>
        /// <param name="aMStabProject">Contains new Geometry objects</param>
        /// <returns>Lookuptable<old ref id, list of new ref ids></returns>
        private Dictionary<GeometryCurve, List<GeometryCurve>> CreateSplitCurves(List<List<List<GeometryCurve>>> aSurfaceList, ref StabilityModel aMStabProject)
        {
            var curveDictionary = new Dictionary<GeometryCurve, List<GeometryCurve>>();
            var regenerateCurves = new Dictionary<GeometryCurve, List<GeometryPoint>>();
            try
            {
                for (int surface = 0; surface < aSurfaceList.Count; surface++)
                {
                    List<GeometryCurve> curveListTop = aSurfaceList[surface][0];
                    List<GeometryCurve> curveListBottom = aSurfaceList[surface][1];
                    for (int curveTop = 0; curveTop < curveListTop.Count; curveTop++)
                    {
                        var bPointTop = (GeometryPoint) aMStabProject.Geometry.Curves.Find(x => x == curveListTop[curveTop]).HeadPoint;
                        var ePointTop = (GeometryPoint) aMStabProject.Geometry.Curves.Find(x => x == curveListTop[curveTop]).EndPoint;

                        GeometryPoint dummyPoint = null;
                        if (bPointTop.X > ePointTop.X)
                        {
                            dummyPoint = bPointTop;
                            bPointTop = ePointTop;
                            ePointTop = dummyPoint;
                        }

                        for (int curveBottom = 0; curveBottom < curveListBottom.Count; curveBottom++)
                        {
                            var bPointBottom = (GeometryPoint) aMStabProject.Geometry.Curves.Find(x => x == curveListBottom[curveBottom]).HeadPoint;
                            var ePointBottom = (GeometryPoint) aMStabProject.Geometry.Curves.Find(x => x == curveListBottom[curveBottom]).EndPoint;

                            if (bPointBottom.X > ePointBottom.X)
                            {
                                dummyPoint = bPointBottom;
                                bPointBottom = ePointBottom;
                                ePointBottom = dummyPoint;
                            }
                            var intersectPoint = new GeometryPoint();
                            if (GeometryData.CheckIfIntersectStricktly(bPointTop, ePointTop, ePointBottom, bPointBottom, ref intersectPoint))
                            {
                                GeometryPoint point = null;
                                //Split top curve
                                double d1x = Math.Abs(bPointTop.X - intersectPoint.X);
                                double d1y = Math.Abs(bPointTop.Y - intersectPoint.Y);
                                double d1z = Math.Abs(bPointTop.Z - intersectPoint.Z);
                                double d2x = Math.Abs(ePointTop.X - intersectPoint.X);
                                double d2y = Math.Abs(ePointTop.Y - intersectPoint.Y);
                                double d2z = Math.Abs(ePointTop.Z - intersectPoint.Z);

                                double sum1 = d1x + d1y + d1z;
                                double sum2 = d2x + d2y + d2z;

                                if ((sum1 > 1e-4) && (sum2 > 1e-4))
                                {
                                    point = (GeometryPoint) aMStabProject.Geometry.Points.Find(x => (Math.Abs((x.X - intersectPoint.X)) < 1e-4) &&
                                                                                                    (Math.Abs((x.Y - intersectPoint.Y)) < 1e-4) &&
                                                                                                    (Math.Abs((x.Z - intersectPoint.Z)) < 1e-4));
                                    if (point == null)
                                    {
                                        aMStabProject.Geometry.Create(intersectPoint);
                                    }
                                    else
                                    {
                                        intersectPoint = point;
                                    }

                                    if (regenerateCurves.ContainsKey(curveListTop[curveTop]))
                                    {
                                        regenerateCurves[curveListTop[curveTop]].Add(intersectPoint);
                                    }
                                    else
                                    {
                                        var newpointList = new List<GeometryPoint>();
                                        newpointList.Add(intersectPoint);
                                        regenerateCurves.Add(curveListTop[curveTop], newpointList);
                                    }
                                }

                                //Split bottom curve
                                d1x = Math.Abs(bPointBottom.X - intersectPoint.X);
                                d1y = Math.Abs(bPointBottom.Y - intersectPoint.Y);
                                d1z = Math.Abs(bPointBottom.Z - intersectPoint.Z);
                                d2x = Math.Abs(ePointBottom.X - intersectPoint.X);
                                d2y = Math.Abs(ePointBottom.Y - intersectPoint.Y);
                                d2z = Math.Abs(ePointBottom.Z - intersectPoint.Z);

                                sum1 = d1x + d1y + d1z;
                                sum2 = d2x + d2y + d2z;

                                if ((sum1 > 1e-4) && (sum2 > 1e-4))
                                {
                                    point = (GeometryPoint) aMStabProject.Geometry.Points.Find(x => (Math.Abs((x.X - intersectPoint.X)) < 1e-4) &&
                                                                                                    (Math.Abs((x.Y - intersectPoint.Y)) < 1e-4) &&
                                                                                                    (Math.Abs((x.Z - intersectPoint.Z)) < 1e-4));
                                    if (point == null)
                                    {
                                        aMStabProject.Geometry.Create(intersectPoint);
                                    }
                                    else
                                    {
                                        intersectPoint = point;
                                    }

                                    if (regenerateCurves.ContainsKey(curveListBottom[curveBottom]))
                                    {
                                        regenerateCurves[curveListBottom[curveBottom]].Add(intersectPoint);
                                    }
                                    else
                                    {
                                        var newpointList = new List<GeometryPoint>();
                                        newpointList.Add(intersectPoint);
                                        regenerateCurves.Add(curveListBottom[curveBottom], newpointList);
                                    }
                                }
                            }
                        }
                    }
                }
                //Add Curves
                foreach (var regCurve in regenerateCurves)
                {
                    var newPointList = new List<GeometryPoint>();

                    newPointList.Add((GeometryPoint) aMStabProject.Geometry.Curves.Find(x => x == regCurve.Key).HeadPoint);
                    newPointList.Add((GeometryPoint) aMStabProject.Geometry.Curves.Find(x => x == regCurve.Key).EndPoint);
                    for (int newPoints = 0; newPoints < regCurve.Value.Count; newPoints++)
                    {
                        newPointList.Add(regCurve.Value[newPoints]);
                    }
                    newPointList.Sort(delegate(GeometryPoint p1, GeometryPoint p2) { return p1.X.CompareTo(p2.X); });
                    //Add new curves
                    for (int newCurve = 0; newCurve < newPointList.Count - 1; newCurve++)
                    {
                        //Check if curves already exist                    
                        var curve1 = (GeometryCurve) aMStabProject.Geometry.Curves.Find(x => (x.HeadPoint == newPointList[newCurve]) && (x.EndPoint == newPointList[newCurve + 1]));
                        var curve2 = (GeometryCurve) aMStabProject.Geometry.Curves.Find(x => (x.HeadPoint == newPointList[newCurve + 1]) && (x.HeadPoint == newPointList[newCurve]));

                        GeometryCurve addcurve = null;
                        if ((curve1 != null) || (curve2 != null))
                        {
                            if (curve1 != null)
                            {
                                addcurve = curve1;
                            }
                            else
                            {
                                addcurve = curve2;
                            }
                        }
                        else
                        {
                            addcurve = new GeometryCurve();
                            addcurve.HeadPoint = newPointList[newCurve];
                            addcurve.EndPoint = newPointList[newCurve + 1];
                            aMStabProject.Geometry.Create(addcurve);
                        }

                        if (curveDictionary.ContainsKey(regCurve.Key))
                        {
                            curveDictionary[regCurve.Key].Add(addcurve);
                        }
                        else
                        {
                            curveDictionary[regCurve.Key] = new List<GeometryCurve>();
                            curveDictionary[regCurve.Key].Add(addcurve);
                        }
                    }
                }

                //Check items in curveDictionary for curves that should be replaced!
                foreach (var keyvalues in curveDictionary)
                {
                    for (int indexI = 0; indexI < keyvalues.Value.Count; indexI++)
                    {
                        if (curveDictionary.ContainsKey(keyvalues.Value[indexI]))
                        {
                            GeometryCurve replacekey = keyvalues.Value[indexI];
                            curveDictionary[keyvalues.Key].RemoveAt(indexI);

                            for (int indexJ = 0; indexJ < curveDictionary[replacekey].Count; indexJ++)
                            {
                                curveDictionary[keyvalues.Key].Insert(indexI + indexJ, curveDictionary[replacekey][indexJ]);
                            }
                        }
                    }
                }
            }
            catch (Exception e)
            {
                Trace.Assert(false, e.ToString());
            }

            return curveDictionary;
        }

        /// <summary>
        /// DetectMultipleSurfaces
        /// In the new geometry, one Layer can be splitted in two or even more different surfaces
        /// These changes have a major effect on the classic geometry (waternets, consolidation matrices and so on)
        /// The place where the top boundary intersects a bottomboundary occure are detected and act as startingpoint
        /// for a new surface.
        /// </summary>
        /// <param name="topCurves">List of Top boundary curves</param>
        /// <param name="bottomCurves">List of bottom boundary curves</param>
        /// <param name="project">New geometry objects</param>
        /// <returns>Returns a list of closed polygons</returns>
        private List<List<GeometryCurve>> DetectMultipleSurfaces(List<GeometryCurve> topList, List<GeometryCurve> bottomList, StabilityModel project)
        {
            try
            {
                var topCurves = new List<GeometryCurve>(topList);
                var bottomCurves = new List<GeometryCurve>(bottomList);

                //Remove curves bottom side that occur in top side
                var deleteCurves = new List<GeometryCurve>();
                foreach (var element in bottomCurves)
                {
                    if (topCurves.Contains(element))
                    {
                        deleteCurves.Add(element);
                    }
                }

                foreach (var element in deleteCurves)
                {
                    topCurves.Remove(element);
                    bottomCurves.Remove(element);
                }

                var surfaces = new List<List<GeometryCurve>>();

                if (CurvesAreCoinciding(topCurves, bottomCurves))
                {
                    return surfaces;
                }

                var surface = new List<GeometryCurve>();

                // Add an additional curve along the left boundary
                if (!topCurves[0].AreConnected(bottomCurves[bottomCurves.Count - 1]))
                {
                    GeometryPoint topConnectionPoint = topCurves[0].HeadPoint.X < topCurves[0].EndPoint.X ? topCurves[0].HeadPoint : topCurves[0].EndPoint;
                    GeometryPoint bottomConnectionPoint = bottomCurves[bottomCurves.Count - 1].HeadPoint.X < bottomCurves[bottomCurves.Count - 1].EndPoint.X ? bottomCurves[bottomCurves.Count - 1].HeadPoint : bottomCurves[bottomCurves.Count - 1].EndPoint;

                    var connectionCurve = new GeometryCurve(bottomConnectionPoint, topConnectionPoint);
                    surface.Add(connectionCurve);
                    project.Geometry.Curves.Add(connectionCurve);
                }

                // Add all curves along the top boundary
                for (int i = 0; i < topCurves.Count; i++)
                {
                    var point1 = project.Geometry.Curves.Find(x => x == topCurves[i]).HeadPoint;
                    var point2 = project.Geometry.Curves.Find(x => x == topCurves[i]).EndPoint;
                    GeometryPoint dummyPoint;
                    if (point1.X > point2.X)
                    {
                        dummyPoint = point1;
                        point1 = point2;
                        point2 = dummyPoint;
                    }

                    bool intersect = false;

                    for (int j = 0; j < bottomCurves.Count; j++)
                    {
                        var pointc1 = project.Geometry.Curves.Find(x => x == bottomCurves[j]).HeadPoint;
                        var pointc2 = project.Geometry.Curves.Find(x => x == bottomCurves[j]).EndPoint;

                        if (pointc1.X < pointc2.X)
                        {
                            dummyPoint = pointc1;
                            pointc1 = pointc2;
                            pointc2 = dummyPoint;
                        }

                        if ((point2 == pointc1) || (point2 == pointc2))
                        {
                            //Intersection occured now check RC

                            double dz1 = (point2.Z - point1.Z);
                            double dx1 = (point2.X - point1.X);
                            double dz2 = (pointc2.Z - pointc1.Z);
                            double dx2 = (pointc2.X - pointc1.X);

                            double a1 = dz1/dx1;
                            double a2 = dz2/dx2;

                            if (a1 != a2)
                            {
                                intersect = true;

                                //Close surface
                                surface.Add(topCurves[i]);

                                int addBottomCurves = bottomCurves.Count - j;
                                var tel = j;
                                if (point2.X < pointc1.X)
                                {
                                    addBottomCurves = addBottomCurves - 1;
                                    tel = j + 1;
                                }
                                for (int k = tel; k < bottomCurves.Count; k++)
                                {
                                    surface.Add(bottomCurves[k]);
                                }

                                for (int k = 0; k < addBottomCurves; k++)
                                {
                                    bottomCurves.RemoveAt(bottomCurves.Count - 1);
                                }

                                surfaces.Add(surface);

                                //Create new surface
                                surface = new List<GeometryCurve>();
                                break;
                            }
                        }
                    }

                    if (!intersect == true)
                    {
                        surface.Add(topCurves[i]);
                    }
                }

                // Add an additional curve along the right boundary
                if (bottomCurves.Count > 0 && !topCurves[topCurves.Count - 1].AreConnected(bottomCurves[0]))
                {
                    GeometryPoint topConnectionPoint = topCurves[topCurves.Count - 1].HeadPoint.X > topCurves[topCurves.Count - 1].EndPoint.X ? topCurves[topCurves.Count - 1].HeadPoint : topCurves[topCurves.Count - 1].EndPoint;
                    GeometryPoint bottomConnectionPoint = bottomCurves[0].HeadPoint.X > bottomCurves[0].EndPoint.X ? bottomCurves[0].HeadPoint : bottomCurves[0].EndPoint;

                    var connectionCurve = new GeometryCurve(topConnectionPoint, bottomConnectionPoint);
                    surface.Add(connectionCurve);
                    project.Geometry.Curves.Add(connectionCurve);
                }

                // Add all curves along the bottom boundary
                surface.AddRange(bottomCurves);

                // Add the surface to the list of surfaces
                surfaces.Add(surface);

                return surfaces;
            }
            catch (Exception e)
            {
                Trace.Assert(false, String.Format("Error found: {0}", e.ToString()));
            }
            return null;
        }

        /// <summary>
        /// Detects whether the curve lists coincide
        /// </summary>
        /// <param name="topCurves">first curve list</param>
        /// <param name="bottomCurves">second curve list</param>
        /// <returns>Coincide</returns>
        private bool CurvesAreCoinciding(List<GeometryCurve> topCurves, List<GeometryCurve> bottomCurves)
        {
            var topPoints = new List<GeometryPoint>();
            var bottomPoints = new List<GeometryPoint>();

            foreach (var topCurve in topCurves)
            {
                if (!topPoints.Contains(topCurve.HeadPoint))
                {
                    topPoints.Add(topCurve.HeadPoint);
                }
                if (!topPoints.Contains(topCurve.EndPoint))
                {
                    topPoints.Add(topCurve.EndPoint);
                }
            }

            foreach (var bottomCurve in bottomCurves)
            {
                if (!bottomPoints.Contains(bottomCurve.HeadPoint))
                {
                    bottomPoints.Add(bottomCurve.HeadPoint);
                }
                if (!bottomPoints.Contains(bottomCurve.EndPoint))
                {
                    bottomPoints.Add(bottomCurve.EndPoint);
                }
            }

            foreach (var topPoint in topPoints.ToArray())
            {
                if (bottomPoints.Contains(topPoint))
                {
                    topPoints.Remove(topPoint);
                    bottomPoints.Remove(topPoint);
                }
            }

            if (topPoints.Count == 0 && bottomPoints.Count == 0)
            {
                return true;
            }

            foreach (var point in topPoints)
            {
                bool containsPoint = false;
                foreach (var curve in bottomCurves)
                {
                    containsPoint |= curve.ContainsPoint(point.GetPoint3D(), 0.0000001);
                }

                if (!containsPoint)
                {
                    return false;
                }
            }

            foreach (var point in bottomPoints)
            {
                bool containsPoint = false;
                foreach (var curve in topCurves)
                {
                    containsPoint |= curve.ContainsPoint(point.GetPoint3D(), 0.0000001);
                }

                if (!containsPoint)
                {
                    return false;
                }
            }

            return true;
        }

        /// <summary>
        /// RemoveOpenEndCurves
        /// Check if Polygon is valid (open end points should not occure here)
        /// End point 1 of curve 1 should belong to begin point 1 of curve 2
        /// A Surface should have an area
        /// </summary>
        /// <param name="curves">Surface with list of curves</param>
        /// <param name="project">Geometry objects</param>
        private void RemoveOpenEndCurves(StabilityModel project, List<GeometryCurve> curves)
        {
            bool stopLoop = false;
            List<GeometryCurve> surface = curves;
            while (true)
            {
                for (int i = 0; i < curves.Count; i++)
                {
                    int curvelistCount = curves.Count; //Meight change during processing 

                    var point1 = project.Geometry.Curves.Find(x => x == surface[i]).HeadPoint;
                    var point2 = project.Geometry.Curves.Find(x => x == surface[i]).EndPoint;
                    int counter1 = 0;
                    int counter2 = 0;
                    for (int j = 0; j < surface.Count; j++)
                    {
                        if (curves[j] != curves[i])
                        {
                            var pointc1 = project.Geometry.Curves.Find(x => x == surface[j]).HeadPoint;
                            var pointc2 = project.Geometry.Curves.Find(x => x == surface[j]).EndPoint;

                            if (Routines2D.AreEqual(point1.X, pointc1.X, Epsilon) || Routines2D.AreEqual(point1.X, pointc2.X, Epsilon))
                            {
                                counter1++;
                            }
                            if (Routines2D.AreEqual(point2.X, pointc1.X, Epsilon) || Routines2D.AreEqual(point2.X, pointc2.X, Epsilon))
                            {
                                counter2++;
                            }
                        }
                    }
                    if (counter1 == 0 || counter2 == 0)
                    {
                        GeometryCurve deleteCurve = surface[i];
                        curves.RemoveAll(x => x == deleteCurve);
                        break;
                    }
                    if (i == curvelistCount - 1)
                    {
                        stopLoop = true;
                    }
                }
                if ((stopLoop) || (curves.Count == 0))
                {
                    break;
                }
            }

            //Pass surface back to main function
            //Check angle of begin curve and end curve

            //var l1point1 = (GeometryPoint)project.Geometry.Curves.Find(x => x == surface[0]).HeadPoint;
            //var l1point2 = (GeometryPoint)project.Geometry.Curves.Find(x => x == surface[0]).EndPoint;
            //double avgpos1 = l1point1.Z + ((l1point2.Z - l1point1.Z) / 2.0);
            //var l2point1 = (GeometryPoint)project.Geometry.Curves.Find(x => x == surface[surface.Count - 1]).HeadPoint;
            //var l2point2 = (GeometryPoint)project.Geometry.Curves.Find(x => x == surface[surface.Count - 1]).EndPoint;
            //double avgpos2 = l2point1.Z + ((l2point2.Z - l2point1.Z) / 2.0);

            //if (avgpos1 == avgpos2)
            //    surface.Clear();

            curves = surface;
        }

        /// <summary>
        /// CreateGeometryLoop
        /// Adds an extra curve on each side (left and right) to make it compatible with the new Geometry
        /// Each curve has two points. The distance between and and begin curve should be zero. If this is not the case, the curve is open 
        /// and should be closed.
        /// Add curves to GeometryLoop.
        /// </summary>
        /// <param name="project"></param>
        /// <param name="curves"></param>
        private GeometryLoop CreateGeometryLoop(StabilityModel project, List<GeometryCurve> curves)
        {
            var geometryLoop = new GeometryLoop();

            //Define starting point                                     
            GeometryCurve curve = project.Geometry.Curves.Find(x => x == curves[0]);
            geometryLoop.CurveList.Add(curve);

            var l1p1 = project.Geometry.Curves.Find(x => x == curves[0]).HeadPoint;
            var l1p2 = project.Geometry.Curves.Find(x => x == curves[0]).EndPoint;

            GeometryPoint dummyPoint;
            if (l1p1.X > l1p2.X)
            {
                dummyPoint = l1p1;
                l1p1 = l1p2;
                l1p2 = dummyPoint;
            }
            else if (l1p1.X == l1p2.X)
            {
                //Assume end point heigher then begin point
                if (l1p1.Z > l1p2.Z)
                {
                    dummyPoint = l1p1;
                    l1p1 = l1p2;
                    l1p2 = dummyPoint;
                }
            }

            //Validate other Curves
            for (int curveIndex = 1; curveIndex < curves.Count + 1; curveIndex++)
            {
                GeometryPoint l2p1;
                GeometryPoint l2p2;
                if (curveIndex < curves.Count)
                {
                    l2p1 = project.Geometry.Curves.Find(x => x == curves[curveIndex]).HeadPoint;
                    l2p2 = project.Geometry.Curves.Find(x => x == curves[curveIndex]).EndPoint;
                }
                else
                {
                    l2p1 = project.Geometry.Curves.Find(x => x == curves[0]).HeadPoint;
                    l2p2 = project.Geometry.Curves.Find(x => x == curves[0]).EndPoint;
                }

                //Check endpoint l1p2 relative to other two points of curve 2
                double d1 = Math.Sqrt(((l1p2.X - l2p1.X)*(l1p2.X - l2p1.X)) + ((l1p2.Z - l2p1.Z)*(l1p2.Z - l2p1.Z)));
                double d2 = Math.Sqrt(((l1p2.X - l2p2.X)*(l1p2.X - l2p2.X)) + ((l1p2.Z - l2p2.Z)*(l1p2.Z - l2p2.Z)));
                double lowVal = 0;
                if (d2 < d1)
                {
                    //Swap points
                    dummyPoint = l2p1;
                    l2p1 = l2p2;
                    l2p2 = dummyPoint;
                    lowVal = d2;
                }
                else
                {
                    lowVal = d1;
                }

                if (lowVal != 0)
                {
                    //Add new curve
                    //Create new Curve
                    curve = new GeometryCurve();
                    curve.HeadPoint = l1p2;
                    curve.EndPoint = l2p1;
                    project.Geometry.Create(curve);
                    //Add curve in loop                            
                    geometryLoop.CurveList.Add(curve);
                    curves.Insert(curveIndex, curve);
                    //Now End point previous curve is begin point next curve
                    l1p1 = l1p2;
                    l1p2 = l2p1;
                }
                else
                {
                    //Now End point previous curve is begin point next curve
                    l1p1 = l2p1;
                    l1p2 = l2p2;
                    if (curveIndex < curves.Count)
                    {
                        curve = project.Geometry.Curves.Find(x => x == curves[curveIndex]);
                        geometryLoop.CurveList.Add(curve);
                    }
                }
            }
            return geometryLoop;
        }

        #region Nested type: PointType

        /// <summary>
        /// Create new Geometry from classic MStab inputfile
        /// </summary>
        private enum PointType
        {
            ptNotUsed = 0,
            ptInCurve,
            ptWaterPoint
        };

        #endregion
    }
}