//-----------------------------------------------------------------------
// <copyright file="PipingCalculator.cs" company="Deltares">
//  Copyright (c) 2009 Deltares. All rights reserved.
// </copyright>
// <author>B.S.T. The</author>
// <email>tom.the@deltares.nl</email>
// <date>25-8-2009</date>
// <summary>Base class for piping calculators Bligh and Sellmeijer</summary>
//-----------------------------------------------------------------------

using Deltares.Geometry;
using Deltares.Geotechnics.Soils;
using Deltares.Geotechnics.SurfaceLines;
using Deltares.Probabilistic;
using Deltares.Uplift;

namespace Deltares.Dam.Data
{
    using System;
    using System.Collections.Generic;

    using Deltares.Piping.Data;
    using Deltares.Standard;
    using Deltares.Geotechnics;

    public struct PipingCalculatorInput
    {
        public double Length;
        public double WhitesConstant;
        public double WaterViscosity;
        public double DInBetweenAquiferlayer;
        public double DBottomAquiferlayer;
        public double PermeabilityInBetweenAquiferlayer;
        public double PermeabilityBottomAquiferlayer;
        public double BeddingAngle;
        public double D70;
        public ProbabilisticStruct ProbPermeabilityInBetweenAquiferlayer;
        public ProbabilisticStruct ProbPermeabilityBottomAquiferlayer;
        public ProbabilisticStruct ProbWhitesConstant;
        public ProbabilisticStruct ProbBeddingAngle;
        public ProbabilisticStruct ProbD70;
    }

    public struct PipingProbabilisticParameters
    {
        public PipingProbabilisticParameters(DistributionType layerHeightDistribution, double layerHeightDeviation)
        {
            LayerHeightDistribution = layerHeightDistribution;
            LayerHeightDeviation = layerHeightDeviation;
        }
        public DistributionType LayerHeightDistribution;
        public double LayerHeightDeviation;
    }

    public class PipingCalculator
    {
        public const string PipingFilenameExtension = ".txt";
        public const string DamParametersSectionId = "[DamParameters]";
        public const string ResultsSectionId = "[Results]";
        public const double D50DividedByD70 = 0.81;
        protected const double CDefaultFluidisationGradient = 0.3;
        protected static double CGammaPDividedByGammaW = 16.5 / Physics.UnitWeightOfwater; // Gamma korrel is now assumed 16.5 kN/m3
        protected double CDefaultWaterViscosity = Physics.WaterViscosity;
        protected const double cDegreeToRadian = Math.PI / 180.0;
        protected ModelParametersForPLLines modelParametersForPLLines;
        protected double requiredSafetyFactor;
        protected const double cToleranceHead = 0.00000001;
        private IList<GaugePLLine> gaugePLLines;
        private IList<Gauge> gauges;
        public double HCritical { get; set; }
        public double? HeaveFactor { get; set; }
        public PipingProbabilisticParameters? PipingProbabilisticParameters { get; set; }
        protected double upliftCriterion;
        protected bool hasProbabilisticParameters = false;
        public string CalculationModelIdentifier = "Piping";
        private bool isHydraulicShortcut = false;
        public string PipingCalculationDirectory { get; set; }

        public PipingCalculator(ModelParametersForPLLines modelParametersForPLLines, double requiredSafetyFactor,
            IList<GaugePLLine> gaugePLLines, IList<Gauge> gauges, PipingProbabilisticParameters? pipingProbabilisticParameters,
            double upliftCriterion)
        {
            this.modelParametersForPLLines = modelParametersForPLLines;
            this.requiredSafetyFactor = requiredSafetyFactor;
            this.gaugePLLines = gaugePLLines;
            this.gauges = gauges;
            this.HCritical = 0.0;
            this.PipingProbabilisticParameters = pipingProbabilisticParameters;
            this.FilenameCalculation = "";
            this.PipingCalculationDirectory = "";
            this.upliftCriterion = upliftCriterion;
            hasProbabilisticParameters = (pipingProbabilisticParameters != null);
        }

        public const double cDefaultMaxReturnValue = 90.0;
        public const double cDefaultMinReturnValue = 0.0;

        protected PLLines PLLines { get; set; }
        protected SoilProfile1D SoilProfile { get; set; }
        public UpliftLocationAndResult UpliftLocationAndResult { get; set; }
        public string FilenameCalculation { get; set; }
        protected GeometryPoint entryPoint { get; set; }

        /// <summary>
        /// Determines the height of the cover layer (distance between surfacelevel and top of aquifer where uplift occurs).
        /// </summary>
        /// <param name="surfaceLevel">The surface level.</param>
        /// <returns></returns>
        protected double DetermineHeightCoverLayer(double surfaceLevel)
        {
            var topLevelAquifer = SoilProfile.GetLayerWithId(UpliftLocationAndResult.LayerWhereUpliftOccuresId).TopLevel;
            topLevelAquifer = Math.Min(topLevelAquifer, surfaceLevel);
            var d = surfaceLevel - topLevelAquifer;

            // if d negative is negative then top of aquifer is higher then surface layer. 
            // This means that the aquifer is exposed on the surface. 
            // In this case d = 0
            d = Math.Max(0, d);
            return d;
        }

        /// <summary>
        /// Determines the piping factor based on h actual and hCritical, setting upperlimit of cDefaultMaxReturnValue.
        /// </summary>
        /// <param name="hActual">The h actual.</param>
        /// <param name="hCritical">The h critical.</param>
        /// <returns></returns>
        protected double DeterminePipingFactorWithUpperLimit(double hActual, double hCritical)
        {
            var pipingFactor = cDefaultMaxReturnValue;
            if (hActual > cToleranceHead)
            {
                var hc = hCritical;
                pipingFactor = hc / hActual;
            }

            return (pipingFactor > cDefaultMaxReturnValue) ? cDefaultMaxReturnValue : pipingFactor;
            
        }

        /// <summary>
        /// Chek  whether there is hydraulic shortcut.
        /// </summary>
        /// <value>
        ///   <c>true</c> if [is hydraulic shortcut]; otherwise, <c>false</c>.
        /// </value>
        public bool IsHydraulicShortcut
        {
            get { return isHydraulicShortcut; }
            set { isHydraulicShortcut = value; }
        }

        /// <summary>
        /// Calculate desigh values for piping channel length and shoulderheight at a given point
        /// </summary>
        /// <param name="location"></param>
        /// <param name="surfaceLine"></param>
        /// <param name="soilProfile"></param>
        /// <param name="waterLevel"></param>
        /// <param name="point"></param>
        /// <returns></returns>
        public virtual PipingDesign CalculateDesignAtPoint(Location location, SurfaceLine2 surfaceLine, SoilProfile1D soilProfile, double waterLevel, GeometryPoint point)
        {
            SetupCalculationAtGivenPoint(location, surfaceLine, soilProfile, waterLevel, point);
            return new PipingDesign(0.0, 0.0);
        }

        /// <summary>
        /// Base for Calculate piping factor. This always return null as the overriding method(s) must provide the real value.
        /// </summary>
        /// <param name="soilProfile"></param>
        /// <param name="waterLevel"></param>
        /// <param name="location"></param>
        /// <param name="surfaceLine"></param>
        /// <returns>Piping factor.</returns>
        public virtual double? CalculatePipingFactor(Location location, SurfaceLine2 surfaceLine, SoilProfile1D soilProfile, double waterLevel)
        {
            this.HCritical = 0.0;
            SetupCalculation(location, surfaceLine, soilProfile, waterLevel);
            return (null);
        }

        /// <summary>
        /// Base for Calculate reliability index. This always return null as the overriding method(s) must provide the real value.
        /// </summary>
        /// <param name="location"></param>
        /// <param name="surfaceLine"></param>
        /// <param name="soilProfile"></param>
        /// <param name="waterLevel"></param>
        /// <returns></returns>
        /// <returns>Reliability index.</returns>
        public virtual double? CalculateReliabilityIndex(Location location, SurfaceLine2 surfaceLine, SoilProfile1D soilProfile, 
            ProbabilisticStruct waterLevel)
        {
            SetupCalculation(location, surfaceLine, soilProfile, waterLevel.Mean);
            return (null);
        }

        /// <summary>
        /// Base for Calculate probability of failure. This always return null as the overriding method(s) must provide the real value.
        /// </summary>
        /// <param name="location"></param>
        /// <param name="surfaceLine"></param>
        /// <param name="soilProfile"></param>
        /// <param name="waterLevel"></param>
        /// <returns></returns>
        /// <returns>Probability of failure.</returns>
        public virtual double? CalculatePipingFailureProbability(Location location, SurfaceLine2 surfaceLine, SoilProfile1D soilProfile, 
            ProbabilisticStruct waterLevel)
        {
            return (null);
        }

        /// <summary>
        /// Base for Calculates the new shoulder height when possible. This always return null as the overriding method(s) must provide the real value.
        /// </summary>
        /// <param name="location"></param>
        /// <param name="surfaceLine"></param>
        /// <param name="soilProfile"></param>
        /// <param name="waterLevel"></param>
        /// <returns></returns>
        public virtual double? CalculateDesignShoulderHeight(Location location, SurfaceLine2 surfaceLine, SoilProfile1D soilProfile, 
            double waterLevel)
        {
            SetupCalculation(location, surfaceLine, soilProfile, waterLevel);
            if (UpliftLocationAndResult != null)
            {
                return CalculateExtraShoulderHeight();
            }
            else
            {
                return null;
            }
        }

        /// <summary>
        /// Prepare for calculation
        /// </summary>
        protected virtual void SetupCalculation(Location location, SurfaceLine2 surfaceLine, SoilProfile1D soilProfile, 
            double waterLevel)
        {
            this.SoilProfile = soilProfile;
            ThrowHelper.ThrowWhenConditionIsTrue<DamFailureMechanismeCalculatorException>
                (
                    string.Format(ThrowHelper.GetResourceString(StringResourceNames.NoSoilProfileDefinedForLocation), ""),
                    () => (this.SoilProfile == null || this.SoilProfile.Layers.Count < 1)
                );

            ThrowHelper.ThrowWhenParameterIsMissing(this, "SurfaceLine", surfaceLine);

            entryPoint = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver);
            ThrowHelper.ThrowWhenParameterIsMissing(surfaceLine, "DikeToeAtRiver", entryPoint);

            ThrowHelper.ThrowIfArgumentNull(this.SoilProfile.BottomAquiferLayer, StringResourceNames.TwoSandlayersRequiredInSoilProfile);

            this.PLLines = CreatePLLines(location, surfaceLine, soilProfile, waterLevel);

            UpliftLocationDeterminator upliftLocationDeterminator = new UpliftLocationDeterminator
            {
                PLLines = this.PLLines,
                SoilProfile = this.SoilProfile,
                SurfaceLine = surfaceLine,
                DikeEmbankmentMaterial = location.GetDikeEmbankmentSoil(),
                XSoilGeometry2DOrigin = location.XSoilGeometry2DOrigin
            };
            UpliftLocationAndResult = upliftLocationDeterminator.GetLocationAndResult(upliftCriterion);

            // next lines are outcommented because they are not in the Release 15.1
            // but they are needed to get the same results as in the new 18.1 version
            // start block for debugging
//            const double upliftCriterionTolerance = 0.000000001;
//            double upliftCriterionTmp = upliftCriterion - upliftCriterionTolerance;
//            UpliftLocationAndResult = upliftLocationDeterminator.GetLocationAndResult(upliftCriterionTmp);
            // end block for debugging
        }

        /// <summary>
        /// Prepare for calculation
        /// </summary>
        protected virtual void SetupCalculationAtGivenPoint(Location location, SurfaceLine2 surfaceLine, SoilProfile1D soilProfile,
            double waterLevel, GeometryPoint givenPoint)
        {
            this.SoilProfile = soilProfile;
            ThrowHelper.ThrowWhenConditionIsTrue<DamFailureMechanismeCalculatorException>
                (
                    string.Format(ThrowHelper.GetResourceString(StringResourceNames.NoSoilProfileDefinedForLocation), ""),
                    () => (this.SoilProfile == null || this.SoilProfile.Layers.Count < 1)
                );

            ThrowHelper.ThrowWhenParameterIsMissing(this, "SurfaceLine", surfaceLine);

            entryPoint = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver);
            ThrowHelper.ThrowWhenParameterIsMissing(surfaceLine, "DikeToeAtRiver", entryPoint);

            ThrowHelper.ThrowIfArgumentNull(this.SoilProfile.BottomAquiferLayer, StringResourceNames.TwoSandlayersRequiredInSoilProfile);

            this.PLLines = CreatePLLines(location, surfaceLine, soilProfile, waterLevel);

            UpliftLocationDeterminator upliftLocationDeterminator = new UpliftLocationDeterminator
            {
                PLLines = this.PLLines,
                SoilProfile = this.SoilProfile,
                SurfaceLine = surfaceLine,
                DikeEmbankmentMaterial = location.GetDikeEmbankmentSoil(),
                XSoilGeometry2DOrigin = location.XSoilGeometry2DOrigin
            };
            UpliftLocationAndResult = upliftLocationDeterminator.GetUpliftFactorAtPoint(givenPoint); 
        }

        /// <summary>
        /// Create PLLines with waterlevel as input
        /// </summary>
        /// <param name="waterLevel"></param>
        /// <returns>created pllines</returns>
        protected PLLines CreatePLLines(Location location, SurfaceLine2 surfaceLine, SoilProfile1D soilProfile, double waterLevel)
        {
            var plLineCreator = new PLLinesCreator
            {
                WaterLevelRiverHigh = waterLevel,
                SurfaceLine = surfaceLine,
                WaterLevelPolder = location.PolderLevel,
                HeadInPLLine2 = location.HeadPL2,
                HeadInPLLine3 = location.HeadPl3,
                HeadInPLLine4 = location.HeadPl4,
                ModelParametersForPLLines = modelParametersForPLLines,
                SoilProfile = soilProfile,
                GaugePLLines = this.gaugePLLines,
                Gauges = this.gauges,
                IsAdjustPL3AndPL4SoNoUpliftWillOccurEnabled = false, // for piping this must be set to false
                PlLineOffsetBelowDikeTopAtRiver = location.PlLineOffsetBelowDikeTopAtRiver,
                PlLineOffsetBelowDikeTopAtPolder = location.PlLineOffsetBelowDikeTopAtPolder,
                DikeEmbankmentMaterial = location.GetDikeEmbankmentSoil(),
                IsHydraulicShortcut = this.IsHydraulicShortcut,
                XSoilGeometry2DOrigin = location.XSoilGeometry2DOrigin
            };

            return plLineCreator.CreateAllPLLines(location);
        }

        /// <summary>
        /// Throws exception D70 not in correct ranfe.
        /// </summary>
        /// <param name="d70">The D70.</param>
        protected void ThrowIfInvalidD70(double d70)
        {
            const double dMinValue = 1e-8;
            const double dMaxValue = 0.1;

            if ((d70 < dMinValue) || (d70 > dMaxValue))
            {
                throw new PipingCalculationException(String.Format("D70 must be between {0} and {1}, but is {2}", dMinValue, dMaxValue, d70));    
            }
        }

        /// <summary>
        /// Defines the piping input from the location parameters.
        /// </summary>
        /// <returns></returns>
        protected PipingCalculatorInput DefinePipingInputFromModel()
        {
            PipingCalculatorInput pipingCalculatorInput = new PipingCalculatorInput();

            Soil inBetweenAquiferlayerSoil = this.SoilProfile.BottomAquiferLayer.Soil;
            Soil bottomAquiferlayerSoil = this.SoilProfile.BottomAquiferLayer.Soil;
            double heightInBetweenAquiferlayer = 0.0;
            double heightBottomAquiferlayer = this.SoilProfile.GetLayerHeight(this.SoilProfile.BottomAquiferLayer);
            if (this.SoilProfile.InBetweenAquiferLayer != null)
            {
                inBetweenAquiferlayerSoil = this.SoilProfile.InBetweenAquiferLayer.Soil;
                heightInBetweenAquiferlayer = this.SoilProfile.GetLayerHeight(this.SoilProfile.InBetweenAquiferLayer);
            }
            else
            {
                heightBottomAquiferlayer = heightBottomAquiferlayer / 2;
                heightInBetweenAquiferlayer = heightBottomAquiferlayer;
            }

            ThrowIfInvalidD70(inBetweenAquiferlayerSoil.DiameterD70);

            pipingCalculatorInput.Length = UpliftLocationAndResult.X - entryPoint.X;
            double topLevel = SoilProfile.GetLayerWithId(UpliftLocationAndResult.LayerWhereUpliftOccuresId).TopLevel;
            pipingCalculatorInput.DInBetweenAquiferlayer = heightInBetweenAquiferlayer;
            pipingCalculatorInput.DBottomAquiferlayer = heightBottomAquiferlayer;
            pipingCalculatorInput.WhitesConstant = inBetweenAquiferlayerSoil.WhitesConstant;
            pipingCalculatorInput.WaterViscosity = CDefaultWaterViscosity;
            pipingCalculatorInput.BeddingAngle = inBetweenAquiferlayerSoil.BeddingAngle;
            pipingCalculatorInput.D70 = Physics.FactorMeterToMicroMeter * inBetweenAquiferlayerSoil.DiameterD70;
            pipingCalculatorInput.PermeabilityInBetweenAquiferlayer = inBetweenAquiferlayerSoil.PermeabKx;
            pipingCalculatorInput.PermeabilityBottomAquiferlayer = bottomAquiferlayerSoil.PermeabKx;

            if (hasProbabilisticParameters)
            {
                ThrowIfInvalidD70(inBetweenAquiferlayerSoil.DiameterD70Stochast.Deviation);
                pipingCalculatorInput.ProbPermeabilityInBetweenAquiferlayer = new ProbabilisticStruct(inBetweenAquiferlayerSoil.PermeabKx, inBetweenAquiferlayerSoil.PermeabKxStochast.Deviation, (int)inBetweenAquiferlayerSoil.PermeabKxStochast.DistributionType, true);
                pipingCalculatorInput.ProbPermeabilityBottomAquiferlayer = new ProbabilisticStruct(bottomAquiferlayerSoil.PermeabKx, bottomAquiferlayerSoil.PermeabKxStochast.Deviation, (int)bottomAquiferlayerSoil.PermeabKxStochast.DistributionType, true);
                pipingCalculatorInput.ProbWhitesConstant = new ProbabilisticStruct(inBetweenAquiferlayerSoil.WhitesConstant, inBetweenAquiferlayerSoil.WhitesConstantStochast.Deviation, (int)inBetweenAquiferlayerSoil.WhitesConstantStochast.DistributionType, true);
                pipingCalculatorInput.ProbBeddingAngle = new ProbabilisticStruct(inBetweenAquiferlayerSoil.BeddingAngle, inBetweenAquiferlayerSoil.BeddingAngleStochast.Deviation, (int)inBetweenAquiferlayerSoil.BeddingAngleStochast.DistributionType, true);
                pipingCalculatorInput.ProbD70 = new ProbabilisticStruct(
                    Physics.FactorMeterToMicroMeter * inBetweenAquiferlayerSoil.DiameterD70,
                    Physics.FactorMeterToMicroMeter * inBetweenAquiferlayerSoil.DiameterD70Stochast.Deviation, 
                    (int) inBetweenAquiferlayerSoil.DiameterD70Stochast.DistributionType, true);
            }
            return pipingCalculatorInput;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="permeabilityPipingLayer"></param>
        /// <returns></returns>
        static protected double KIntrinisc(double permeabilityPipingLayer, double waterViscosity)
        {
            double result = permeabilityPipingLayer * waterViscosity / Physics.GravityConstant;
            return result;
        }

        /// <summary>
        /// Check on precondition
        /// </summary>
        protected void ThrowIfNoPLLinesDefined()
        {
            if (PLLines == null)
                throw new UpliftLocationDeterminatorException("Required pllines not found");
        }

        /// <summary>
        /// Calculates the height of the extra shoulder.
        /// </summary>
        /// <returns></returns>
        protected double CalculateExtraShoulderHeight()
        {
            var calculator = new UpliftCalculator();
            calculator.SoilProfile = this.SoilProfile;
            calculator.SurfaceLevel = UpliftLocationAndResult.Z;
            calculator.PhreaticLevel = this.PLLines.Lines[PLLineType.PL1].ZFromX(UpliftLocationAndResult.X);
            calculator.TopOfLayerToBeEvaluated = SoilProfile.GetLayerWithId(UpliftLocationAndResult.LayerWhereUpliftOccuresId).TopLevel;
            PLLine plLine;
            if (UpliftLocationAndResult.LayerWhereUpliftOccuresId == SoilProfile.BottomAquiferLayer.Id)
            {
                plLine = this.PLLines.Lines[PLLineType.PL3];
            }
            else
            {
                plLine = this.PLLines.Lines[PLLineType.PL4];
            }
            return calculator.CalculateExtraHeight(plLine.ZFromX(UpliftLocationAndResult.X), upliftCriterion);
        }
    }
}