Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/ExitPointDeterminator.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/ExitPointDeterminator.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/ExitPointDeterminator.cs (revision 253) @@ -0,0 +1,630 @@ +using System; +using System.Collections.Generic; +using System.Linq; + +namespace Deltares.DamPiping +{ + /// + /// Class to determine the exit point + /// + public class ExitPointDeterminator + { + private double leakageLength; + private double phiPolder; + + private double phiToe; + private double rExit = PipingConstants.RExitDefault; + private double rToe = PipingConstants.RToeDefault; + + //TODO: this is not officially documented, functionality required by DAM + private double upliftFactor; + +// //TODO: this is not officially documented, functionality required by Ringtoets +// private readonly List upliftLocationAndResultList = new List(); + private double volumetricWeightOfWater = 9.81; +// private LanguageType language = LanguageType.Dutch; +// +// /// +// /// Language for (validation) messages +// /// +// public LanguageType Language +// { +// get { return language; } +// set { language = value; } +// } + + /// + /// Piezometric head at polder side as input + /// + public double PhiPolder + { + get + { + return phiPolder; + } + set + { + phiPolder = value; + } + } + + /// + /// Damping factor at Toe as input (default = 1). + /// + public double RToe + { + get + { + return rToe; + } + set + { + rToe = value; + } + } + + /// + /// River water level as input + /// + public double HRiver { get; set; } + + /// + /// Leakage length as input + /// + public double LeakageLength + { + get + { + return leakageLength; + } + set + { + leakageLength = value; + } + } + + /// + /// 2D surface line with all relevant points labelled (toe, ditch) as input. + /// + public PipingSurfaceLine SurfaceLine { get; set; } + + /// + /// 1D soil profile with all layers and neccessary material parameters as input + /// + public PipingProfile SoilProfile { get; set; } + + /// + /// Required factor of safety for Uplift as input + /// + public double RequiredFactorOfSafety { get; set; } + + /// + /// Exit point as result + /// + public PipingPoint ExitPoint { get; set; } + + /// + /// Volumetric weight of water as input + /// + public double VolumetricWeightOfWater + { + get + { + return volumetricWeightOfWater; + } + set + { + volumetricWeightOfWater = value; + } + } + + /// + /// Model factor Uplift as input + /// + public double ModelFactorUplift { get; set; } + + /// + /// Determines whether the dry oven weight should be used instead of the weight above phreatic level as input + /// + public bool IsUseOvenDryUnitWeight { get; set; } + + /// + /// Minimum thickness of the cover layer as input + /// + public double MinimumThicknessCoverLayer { get; set; } + + /// + /// damping factor at the exit point (default = 1) as input + /// + public double RExit + { + get + { + return rExit; + } + set + { + rExit = value; + } + } + + /// + /// Gets or sets the uplift factor. + /// + /// + /// The uplift factor. + /// + //TODO: this is not officially documented, functionality required by DAM + public double UpliftFactor + { + get + { + return upliftFactor; + } + set + { + upliftFactor = value; + } + } + + /// + /// Gets or sets the layer where uplift occures identifier. + /// + /// + /// The layer where uplift occures identifier. + /// + //TODO: this is not officially documented, functionality required by DAM + public string LayerWhereUpliftOccuresId { get; set; } + +// /// +// /// List of all uplift locations with upliftfactor as result +// /// TODO: this is not officially documented, functionality required by Ringtoets +// /// +// /// +// public List GetUpliftLocationAndResultList() +// { +// return upliftLocationAndResultList; +// } + + /// + /// Start the calculation of all output parameters. + /// + public virtual void Calculate() + { +// var restoreLanguage = LocalizationManager.CurrentLanguage; +// LocalizationManager.CurrentLanguage = Language; +// var dumper = new XmlDumper(); + try + { +// DumpThis(dumper); + + DeterminePhiToe(); + var relevantSurfacePointsList = DetermineRelevantSurfacePointList(); + var surfacePointsList = relevantSurfacePointsList as PipingPoint[] ?? + relevantSurfacePointsList.ToArray(); +// FillUpliftLocationAndResultList(surfacePointsList); + foreach (var surfacePoint in surfacePointsList) + { + var upliftLocationAndResult = GetUpliftFactorAtPoint(surfacePoint); + if (upliftLocationAndResult != null) + { + if (upliftLocationAndResult.UpliftFactor < RequiredFactorOfSafety && + upliftLocationAndResult.Zu < 0) + { + ExitPoint = surfacePoint; + upliftFactor = upliftLocationAndResult.UpliftFactor; + LayerWhereUpliftOccuresId = upliftLocationAndResult.LayerWhereUpliftOccuresId; + break; + } + } + } + } + finally + { +// DumpThis(dumper); +// LocalizationManager.CurrentLanguage = restoreLanguage; + } + } + + private void LogParameterIsNaN(IList list, string paramName) + { + var msg = string.Format(Helper.Translate("NaNParameterError"), paramName); + list.Add(msg); + } + + /// + /// Validates the input + /// + /// a filled list when errors are found else an empty list + public virtual List Validate() + { +// var restoreLanguage = LocalizationManager.CurrentLanguage; +// LocalizationManager.CurrentLanguage = Language; + var errors = new List(); +// var dumper = new XmlDumper(); + try + { +// DumpThis(dumper); + + if (double.IsNaN(PhiPolder)) + LogParameterIsNaN(errors, "PhiPolder"); + if (double.IsNaN(RToe)) + LogParameterIsNaN(errors, "RToe"); + if (double.IsNaN(HRiver)) + LogParameterIsNaN(errors, "HRiver"); + if (double.IsNaN(LeakageLength)) + LogParameterIsNaN(errors, "LeakageLength"); + if (double.IsNaN(RequiredFactorOfSafety)) + LogParameterIsNaN(errors, "RequiredFactorOfSafety"); + if (double.IsNaN(VolumetricWeightOfWater)) + LogParameterIsNaN(errors, "VolumetricWeightOfWater"); + if (double.IsNaN(ModelFactorUplift)) + LogParameterIsNaN(errors, "ModelFactorUplift"); + if (double.IsNaN(MinimumThicknessCoverLayer)) + LogParameterIsNaN(errors, "MinimumThicknessCoverLayer"); + if (double.IsNaN(RExit)) + LogParameterIsNaN(errors, "RExit"); + if (double.IsNaN(UpliftFactor)) + LogParameterIsNaN(errors, "UpliftFactor"); + + if (SoilProfile == null) + { + errors.Add(Helper.Translate("SoilprofileNotDefined")); + return errors; + } + if (SurfaceLine == null) + { + errors.Add(Helper.Translate("SurfaceLineNotDefined")); + return errors; + } + if (SurfaceLine.Points.Count < 2) + { + errors.Add(Helper.Translate("SurfaceLineAtLeastTwoPoints")); + return errors; + } + + try + { + SurfaceLine.Validate(); + } + catch (Exception e) + { + errors.Add(e.Message); + } + + string layerName; + var point = GetLocalDikeToe() ?? SurfaceLine.Points[0]; + DeterminePhiToe(); + var phi = DeterminePhiAtX(point.X); + var upliftCalculator = PrepareUpliftCalculator(point, phi, out layerName); + var upliftErrors = upliftCalculator.Validate(); + errors.AddRange(upliftErrors); + try + { + SoilProfile.Validate(); + } + catch (Exception e) + { + errors.Add(e.Message); + } + double? topOfLayerToBeEvaluated = null; + if (SoilProfile.BottomAquiferLayer != null) + { + topOfLayerToBeEvaluated = SoilProfile.BottomAquiferLayer.TopLevel; + } + if (SoilProfile.TopAquiferLayer != null) + { + topOfLayerToBeEvaluated = SoilProfile.TopAquiferLayer.TopLevel; + } + if (topOfLayerToBeEvaluated == null) + { + errors.Add(Helper.Translate("SoilProfileHasNoAquifer")); + return errors; + } + // EffectiveThicknessCalcualtor should be validated here, but has no Validate. So we use validator of SoilVolumicMassCalculator + var svmCalculator = PrepareSoilVolumicMassCalculator(point, topOfLayerToBeEvaluated.Value); + var massErrors = svmCalculator.Validate(); + errors.AddRange(massErrors); + } + finally + { +// DumpThis(dumper); +// LocalizationManager.CurrentLanguage = restoreLanguage; + } + + return errors; + } + +// /// +// /// TODO: this is not officially documented, functionality required by Ringtoets +// /// fill list with all relevant surface point locations with uplift factor +// /// +// /// +// private void FillUpliftLocationAndResultList(IEnumerable relevantSurfacePointsList) +// { +// foreach (var surfacePoint in relevantSurfacePointsList) +// { +// var upliftLocationAndResult = GetUpliftFactorAtPoint(surfacePoint); +// if (upliftLocationAndResult != null) +// { +// upliftLocationAndResultList.Add(upliftLocationAndResult); +// } +// } +// } + + /// + /// Determine the (initial) piezometric head at the toe. + /// + /// + private void DeterminePhiToe() + { + phiToe = PhiPolder + (RToe*(HRiver - PhiPolder)); + } + + /// + /// Returns the dike toe as defined for this kernel. If the surface line has no defined toe than its + /// assumed that the toe is the first point of the surfaceline. + /// + /// + private PipingPoint GetLocalDikeToe() + { + var point = SurfaceLine.GetDikeToeInward(); + return point; + } + + /// + /// Determines the relevant points for the calculation of the possible exit points + /// + /// + private IEnumerable DetermineRelevantSurfacePointList() + { + PipingPoint startSurfacePoint = GetLocalDikeToe(); + + IEnumerable relevantSurfacePointsList = from PipingPoint point in SurfaceLine.Points + where point.X >= startSurfacePoint.X + orderby point.X ascending + select point; + relevantSurfacePointsList = GetCheckedSurfaceLine(relevantSurfacePointsList); + return relevantSurfacePointsList; + } + + /// + /// Ensures that the points on the surface line are never more than cDiff (0.5) apart. + /// + /// + /// + private IEnumerable GetCheckedSurfaceLine(IEnumerable originalLine) + { + const double cDiff = 0.5; + var newLine = new List(); + if (originalLine == null) + { + return newLine; + } + var pipingPoints = originalLine as PipingPoint[] ?? originalLine.ToArray(); + var x = pipingPoints.First().X; + foreach (var point in pipingPoints) + { + while (point.X > x + cDiff) + { + var newPoint = new PipingPoint(point) + { + X = x + cDiff + }; + if (newPoint.X > newLine.Last().X) + { + newPoint.Z = newLine.Last().Z + ((newPoint.X - newLine.Last().X)/(point.X - newLine.Last().X))* + (point.Z - newLine.Last().Z); + newLine.Add(newPoint); + } + x = newPoint.X; + } + newLine.Add(point); + } + return newLine; + } + + /// + /// Calculate upliftfactor for given point + /// + /// + /// location and upliftfactor + private UpliftLocationAndResult GetUpliftFactorAtPoint(PipingPoint point) + { + string layerName; + var phi = DeterminePhiAtX(point.X); + var upliftCalculator = PrepareUpliftCalculator(point, phi, out layerName); + + upliftCalculator.Calculate(); + var res = new UpliftLocationAndResult(point) + { + UpliftFactor = upliftCalculator.FoSu, + LayerWhereUpliftOccuresId = layerName, + Zu = upliftCalculator.Zu, + Phicu = upliftCalculator.DeltaPhiCu + upliftCalculator.HExit, + Phi = phi + }; + return res; + } + + /// + /// Prepares the WTIUpliftCalculator for calculation + /// + /// + /// + /// + /// A filled calculator ready for calculation/validation + private UpliftCalculator PrepareUpliftCalculator(PipingPoint point, double phi, out string layerName) + { + var upliftCalculator = new UpliftCalculator + { + HExit = DeterminePhreaticLevelAtPoint(point), + HRiver = HRiver, + RExit = RExit, + PhiPolder = PhiPolder, + PhiExit = phi, + ModelFactorUplift = ModelFactorUplift, + VolumetricWeightOfWater = VolumetricWeightOfWater + }; + + var bottomStress = 0.0; + layerName = ""; + if (SoilProfile.BottomAquiferLayer != null) + { + bottomStress = DetermineEffectiveStressAtPoint(point); + upliftCalculator.EffectiveStress = bottomStress; + if (SoilProfile.BottomAquiferLayer != null) + { + layerName = SoilProfile.BottomAquiferLayer.Name; + } + } + + double? inbetweenStress = null; + if (SoilProfile.TopAquiferLayer != null) + { + inbetweenStress = DetermineEffectiveStressAtPoint(point); + } + if (inbetweenStress != null) + { + if (inbetweenStress.Value <= bottomStress) + { + upliftCalculator.EffectiveStress = inbetweenStress.Value; + if (SoilProfile.TopAquiferLayer != null) + { + layerName = SoilProfile.TopAquiferLayer.Name; + } + } + } + else + { + upliftCalculator.EffectiveStress = bottomStress; + } + return upliftCalculator; + } + + /// + /// Determines the effective stress at the given point + /// + /// + /// + private double DetermineEffectiveStressAtPoint(PipingPoint point) + { + var effectiveThicknessCalculator = new EffectiveThicknessCalculator + { + ExitPointXCoordinate = point.X, + SoilProfile = SoilProfile, + SurfaceLine = SurfaceLine, + PhreaticLevel = DeterminePhreaticLevelAtPoint(point), + VolumicWeightOfWater = VolumetricWeightOfWater + }; + effectiveThicknessCalculator.Calculate(); + return effectiveThicknessCalculator.EffectiveStress; + } + + /// + /// Prepares the SoilVolumicMassCalculator for calculation; pure for the validation part + /// + /// + /// + /// A filled calculator ready for calculation/validation + private SoilVolumicMassCalculator PrepareSoilVolumicMassCalculator(PipingPoint point, double top) + { + var svmCalculator = new SoilVolumicMassCalculator + { + SoilProfile = SoilProfile, + PhreaticLevel = DeterminePhreaticLevelAtPoint(point), + VolumicWeightOfWater = VolumetricWeightOfWater, + SurfaceLevel = point.Z, + TopOfLayerToBeEvaluated = top, + MinimumThicknessCoverLayer = MinimumThicknessCoverLayer, + IsUseOvenDryUnitWeight = IsUseOvenDryUnitWeight + }; + return svmCalculator; + } + + /// + /// Determines the phreatic level at the given point + /// Phreaticlevel is at least polder level when calculation a ditch point and else the phreatic level is the surface level + /// + /// + /// + private double DeterminePhreaticLevelAtPoint(PipingPoint point) + { + return IsPointPartOfDitch(point) ? Math.Max(PhiPolder, point.Z) : point.Z; + } + + /// + /// Determines whether the point is situated within the ditch + /// + /// + /// + private bool IsPointPartOfDitch(PipingPoint point) + { + return SurfaceLine.HasDitch() && SurfaceLine.DitchDikeSide.X <= point.X && + SurfaceLine.DitchPolderSide.X >= point.X; + } + + /// + /// Determination of the piezometric head at a given X location + /// + /// + /// + private double DeterminePhiAtX(double x) + { + var xToe = GetLocalDikeToe().X; + if (leakageLength <= 0 || Double.IsNaN(leakageLength)) + { + // if leekagelength is unknown, return phi polder + return phiPolder; + } + var exp = Math.Exp((xToe - x)/leakageLength); + if (Double.IsNaN(exp)) + { + throw new PipingException(Helper.Translate("NaNResultError")); + } + return phiPolder + ((phiToe - phiPolder)*exp); + } + + // #region DebugOut + // private bool isDebugOutputEnabled = false; + // + // protected void DumpThis(XmlDumper dumper) + // { + // if (isDebugOutputEnabled) + // { + // dumper.DumpObject(this); + // } + // } + // + // /// + // /// Enables debug-dumping to XML + // /// + // public void EnableDebugOutput() + // { + // isDebugOutputEnabled = true; + // } + // + // /// + // /// Fills the data from debug file. + // /// + // public ExitPointDeterminator FillFromDebugFile(string content) + // { + // var serializer = new Deltares.Standard.IO.Xml.XmlDeserializer(); + // if (!string.IsNullOrEmpty(content)) + // { + // var res = (ExitPointDeterminator)serializer.XmlDeserializeFromString(content); + // return res; + // } + // return null; + // } + // + // /// + // /// Gets the main result + // /// + // public string GetMainResult() + // { + // return "ExitPoint.X = " + ExitPoint.X.ToString() + " ExitPoint.Z = " + ExitPoint.Z.ToString(); + // } + // #endregion + } + } \ No newline at end of file Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/UpliftCalculator.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/UpliftCalculator.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/UpliftCalculator.cs (revision 253) @@ -0,0 +1,305 @@ +using System; +using System.Collections.Generic; + +namespace Deltares.DamPiping +{ + /// + /// Exception class for UpliftCalculator + /// + public class UpliftCalculatorException : Exception + { + public UpliftCalculatorException(string message) : base(message) {} + } + + /// + /// Class for calculating the Uplift output parameters + /// + public class UpliftCalculator + { + // Input parameters + + private double deltaPhiCu = Double.NaN; + private double foSu = Double.NaN; + private double hcu = Double.NaN; + private double volumetricWeightOfWater = PipingConstants.UnitWeightOfWater; + private double zu = Double.NaN; +// private LanguageType language = LanguageType.Dutch; + +// /// +// /// Language for (validation) messages +// /// +// public LanguageType Language +// { +// get { return language; } +// set { language = value; } +// } + + /// + /// Gets or sets the volumic weight of water. + /// + /// + /// The volumic weight of water. + /// + public double VolumetricWeightOfWater + { + get + { + return volumetricWeightOfWater; + } + set + { + volumetricWeightOfWater = value; + } + } + + /// + /// Gets or sets the model factor uplift. + /// + /// + /// The model factor uplift. + /// + public double ModelFactorUplift { get; set; } + + /// + /// Gets or sets the effective stress. + /// + /// + /// The effective stress. + /// + public double EffectiveStress { get; set; } + + /// + /// Gets or sets the h-river. + /// + /// + /// The h river. + /// + public double HRiver { get; set; } + + /// + /// Gets or sets the phi-exit. + /// + /// + /// The phi exit. + /// + public double PhiExit { get; set; } + + /// + /// Gets or sets the R-Exit. + /// + /// + /// The r exit. + /// + public double RExit { get; set; } + + /// + /// Gets or sets the H-exit. + /// + /// + /// The h exit. + /// + public double HExit { get; set; } + + public double PhiPolder { get; set; } + + // Output parameters + + /// + /// Gets the Zu. + /// + /// + /// The Zu. + /// + public double Zu + { + get + { + return zu; + } + set {} + } + + /// + /// Gets the Fo_su. + /// + /// + /// The fo su. + /// + public double FoSu + { + get + { + return foSu; + } + set {} + } + + /// + /// Gets the delta phi cu. + /// + /// + /// The delta phi cu. + /// + public double DeltaPhiCu + { + get + { + return deltaPhiCu; + } + set {} + } + + /// + /// Gets the hcu. + /// + /// + /// The hcu. + /// + public double Hcu + { + get + { + return hcu; + } + set {} + } + + /// + /// Calculates the output parameters from the input. + /// + public void Calculate() + { +// var restoreLanguage = LocalizationManager.CurrentLanguage; +// LocalizationManager.CurrentLanguage = Language; +// var dumper = new XmlDumper(); + try + { +// DumpThis(dumper); + + if (EffectiveStress < 0) + { + var message = Helper.Translate("EffectiveStressNegative"); + throw new UpliftCalculatorException(message); + } + if (Math.Abs(RExit) < PipingConstants.Epsilon) + { + var message = Helper.Translate("RexitIsZero"); + throw new PipingException(message); + } + deltaPhiCu = EffectiveStress/VolumetricWeightOfWater; + zu = ModelFactorUplift*deltaPhiCu - (PhiExit - HExit); + DetermineFactorOfSafety(); + hcu = (deltaPhiCu + HExit - PhiPolder) / RExit + PhiPolder; + } + finally + { +// DumpThis(dumper); +// LocalizationManager.CurrentLanguage = restoreLanguage; + } + } + + private void DetermineFactorOfSafety() + { + foSu = 0; + + var reducedDeltaPhiCu = ModelFactorUplift*deltaPhiCu; + if (Math.Abs(reducedDeltaPhiCu) > PipingConstants.Epsilon) + { + foSu = reducedDeltaPhiCu / Math.Max(0, PhiExit - HExit); + } + } + + private void LogParameterIsNaN(IList list, string paramName) + { + var msg = string.Format(Helper.Translate("NaNParameterError"), paramName); + list.Add(msg); + } + + /// + /// Validates the input + /// + /// a filled list when errors are found else an empty list + public List Validate() + { +// var restoreLanguage = LocalizationManager.CurrentLanguage; +// LocalizationManager.CurrentLanguage = Language; + var errors = new List(); +// var dumper = new XmlDumper(); + try + { +// DumpThis(dumper); + + if (double.IsNaN(VolumetricWeightOfWater)) + LogParameterIsNaN(errors, "VolumetricWeightOfWater"); + if (double.IsNaN(ModelFactorUplift)) + LogParameterIsNaN(errors, "ModelFactorUplift"); + if (double.IsNaN(EffectiveStress)) + LogParameterIsNaN(errors, "EffectiveStress"); + if (double.IsNaN(HRiver)) + LogParameterIsNaN(errors, "HRiver"); + if (double.IsNaN(PhiExit)) + LogParameterIsNaN(errors, "PhiExit"); + if (double.IsNaN(RExit)) + LogParameterIsNaN(errors, "RExit"); + if (double.IsNaN(HExit)) + LogParameterIsNaN(errors, "HExit"); + if (double.IsNaN(PhiPolder)) + LogParameterIsNaN(errors, "PhiPolder"); + + if (Math.Abs(RExit) < PipingConstants.Epsilon) + { + errors.Add(Helper.Translate("RexitIsZero")); + } + } + finally + { +// DumpThis(dumper); +// LocalizationManager.CurrentLanguage = restoreLanguage; + } + return errors; + } + +// #region DebugOut +// private bool isDebugOutputEnabled = false; +// +// private void DumpThis(XmlDumper dumper) +// { +// if (isDebugOutputEnabled) +// { +// dumper.DumpObject(this); +// } +// } +// +// /// +// /// Enables debug-dumping to XML +// /// +// public void EnableDebugOutput() +// { +// isDebugOutputEnabled = true; +// } +// +// /// +// /// Fills the data from debug file. +// /// +// public WTIUpliftCalculator FillFromDebugFile(string content) +// { +// var serializer = new Deltares.Standard.IO.Xml.XmlDeserializer(); +// if (!string.IsNullOrEmpty(content)) +// { +// var res = (WTIUpliftCalculator) serializer.XmlDeserializeFromString(content); +// return res; +// } +// return null; +// } +// +// /// +// /// Gets the main result +// /// +// public string GetMainResult() +// { +// return "FoS = " + foSu.ToString(); +// } +// +// #endregion + } +} \ No newline at end of file Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/UpliftLocationAndResult.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/UpliftLocationAndResult.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/UpliftLocationAndResult.cs (revision 253) @@ -0,0 +1,46 @@ +namespace Deltares.DamPiping +{ + /// + /// Class for the storage of the uplift results per location. + /// + public class UpliftLocationAndResult : PipingPoint + { + public UpliftLocationAndResult() + { + UpliftFactor = 0; + LayerWhereUpliftOccuresId = ""; + } + + public UpliftLocationAndResult(PipingPoint point) + { + X = point.X; + Y = point.Y; + Z = point.Z; + } + + /// + /// The uplift factor + /// + public double UpliftFactor { get; set; } + + /// + /// The limit state value + /// + public double Zu { get; set; } + + /// + /// The name (ID) of the layer where uplift occurs + /// + public string LayerWhereUpliftOccuresId { get; set; } + + /// + /// Critical head + /// + public double Phicu { get; set; } + + /// + /// Calculated head + /// + public double Phi { get; set; } + } +} \ No newline at end of file Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Helper.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Helper.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Helper.cs (revision 253) @@ -0,0 +1,29 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace Deltares.DamPiping +{ + public class Helper + { + public static string Translate(string textId, string id = "") + { + return textId + ": " + id; + } + + public static bool AlmostEquals(double double1, double double2, double precision) + { + if (Double.IsNaN(double1) && Double.IsNaN(double2)) + return true; + if (Double.IsNaN(double1) || Double.IsNaN(double2)) + return false; + return Math.Abs(double1 - double2) <= precision; + } + public static bool AlmostEquals(double double1, double double2) + { + return AlmostEquals(double1, double2, 1E-07); + } + + } +} Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingSurfaceLine.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingSurfaceLine.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingSurfaceLine.cs (revision 253) @@ -0,0 +1,362 @@ +using System; +using System.Collections.Generic; +using System.Linq; +//using Deltares.Mathematics; + +namespace Deltares.DamPiping +{ + /// + /// PipingSurfaceLine Exception class + /// + public class PipingSurfaceLineException : Exception + { + /// + /// Initializes a new instance of the class. + /// + /// The message that describes the error. + public PipingSurfaceLineException(string message) : base(message) { } + } + + /// + /// Class to hold the surface line for the Piping kernel + /// + public class PipingSurfaceLine + { + private List points = new List(); + + /// + /// Gets or sets the name. + /// + /// + /// The name. + /// + public string Name {get; set;} + + /// + /// Gets or sets the points. + /// + /// + /// The points. + /// + public List Points + { + get + { + return points; + } + set + { + points = value; + } + } + + /// + /// Determines whether this instance has ditch. + /// + /// + public bool HasDitch() + { + var i = 0; + var foundDitch = false; + while (!foundDitch && i < Points.Count) + { + foundDitch = points[i].Type == PipingCharacteristicPointType.BottomDitchDikeSide; + i++; + } + return foundDitch; + } + + /// + /// Gets the ditch dike side. + /// + /// + /// The ditch dike side. + /// + public PipingPoint DitchDikeSide + { + get + { + var point = points.Where(x => x.Type == PipingCharacteristicPointType.DitchDikeSide); + return point.FirstOrDefault(); + } + } + + /// + /// Gets the bottom ditch dike side. + /// + /// + /// The bottom ditch dike side. + /// + public PipingPoint BottomDitchDikeSide + { + get + { + var point = points.Where(x => x.Type == PipingCharacteristicPointType.BottomDitchDikeSide); + return point.FirstOrDefault(); + } + } + + /// + /// Gets the bottom ditch polder side. + /// + /// + /// The bottom ditch polder side. + /// + public PipingPoint BottomDitchPolderSide + { + get + { + var point = points.Where(x => x.Type == PipingCharacteristicPointType.BottomDitchPolderSide); + return point.FirstOrDefault(); + } + } + + /// + /// Gets the ditch polder side. + /// + /// + /// The ditch polder side. + /// + public PipingPoint DitchPolderSide + { + get + { + var point = points.Where(x => x.Type == PipingCharacteristicPointType.DitchPolderSide); + return point.FirstOrDefault(); + } + } + + /// + /// Gets the shoulder base inside. + /// + /// + /// The shoulder base inside. + /// + public PipingPoint ShoulderBaseInside + { + get + { + var point = points.Where(x => x.Type == PipingCharacteristicPointType.ShoulderBaseInside); + return point.FirstOrDefault(); + } + } + + /// + /// Gets the dike toe at polder. + /// + /// + /// The dike toe at polder. + /// + public PipingPoint DikeToeAtPolder + { + get + { + var point = points.Where(x => x.Type == PipingCharacteristicPointType.DikeToeAtPolder); + return point.FirstOrDefault(); + } + } + + /// + /// If shoulder is present then the toe of the shoulder (DikeToeAtPolder) is returned + /// Else if present the toe of the dike (ShoulderBaseInside) is returned + /// Else the first point of the surfaceline is returned if present. + /// + /// toe of the dike + public virtual PipingPoint GetDikeToeInward() + { + if (ShoulderBaseInside != null) + { + return ShoulderBaseInside; + } + if (DikeToeAtPolder != null) + { + return DikeToeAtPolder; + } + return Points.FirstOrDefault(); + } + + /// + /// Gets the z at x. + /// + /// The x. + /// + public virtual double GetZatX(double x) + { + for (int i = 0; i < Points.Count - 1; i++) + { + var current = Points[i]; + var next = Points[i + 1]; + + var leftOffset = x - current.X; + var rightOffset = next.X - x; + + if (Math.Abs(leftOffset) < PipingConstants.Accuracy) + { + return current.Z; + } + if (Math.Abs(rightOffset) < PipingConstants.Accuracy) + { + return next.Z; + } + if (leftOffset >= 0 && rightOffset >= 0) + { + var fraction = leftOffset / (leftOffset + rightOffset); + + return (1.0 - fraction) * current.Z + fraction * next.Z; + } + } + + return Double.NaN; + } + + /// + /// Intersections the points xz with line xz. + /// + /// The begin. + /// The end. + /// + public IList IntersectionPointsXzWithLineXz(PipingPoint begin, PipingPoint end) + { + var intersectionPointsWithLine = new List(); + + for (int pointIndex = 0; pointIndex < Points.Count - 1; pointIndex++) + { + DoIntersectAndAddToCollection(begin, end, Points[pointIndex], Points[pointIndex + 1], intersectionPointsWithLine); + } + return intersectionPointsWithLine; + } + + /// + /// Validates this instance. + /// + /// + internal void Validate() + { + if (!AreAllCharacteristicPointsOrdered()) + { + var format = Helper.Translate("CharacteristicPointXNotAscending"); + var message = String.Format(format, Name); + throw new PipingSurfaceLineException(message); + } + if (!IsDitchCorrect()) + { + var format = Helper.Translate("DitchIncorrect"); + var message = String.Format(format, Name); + throw new PipingSurfaceLineException(message); + } + if (!ArePointsAscending()) + { + var format = Helper.Translate("PointXNotAscending"); + var message = String.Format(format, Name); + throw new PipingSurfaceLineException(message); + } + } + + /// + /// Checks if points are ascending. + /// + /// true if points are ascending, otherwise false + private bool ArePointsAscending() + { + for (int i = 1; i < points.Count; i++) + { + if (points[i].X < points[i - 1].X) + { + return false; + } + } + return true; + } + + + /// + /// Does the intersect and add to collection. + /// + /// a point1. + /// a point2. + /// The begin. + /// The end. + /// The intersection points with line. + private static void DoIntersectAndAddToCollection(PipingPoint aPoint1, PipingPoint aPoint2, PipingPoint begin, PipingPoint end, ICollection intersectionPointsWithLine) + { + Point2D intersectionPoint2D = null; + Routines2D.DetermineIf2DLinesIntersectStrickly(new Point2D(aPoint1.X, aPoint1.Z), new Point2D(aPoint2.X, aPoint2.Z), new Point2D(begin.X, begin.Z), new Point2D(end.X, end.Z), ref intersectionPoint2D); + if (intersectionPoint2D != null) + { + var intersectionPoint = new PipingPoint(intersectionPoint2D.X, 0, intersectionPoint2D.Y); + if (NoPointSameXzLocation(intersectionPointsWithLine, intersectionPoint)) + { + intersectionPointsWithLine.Add(intersectionPoint); + } + } + } + + /// + /// Determines if there is a point with same xz location. + /// + /// The collection. + /// The point. + /// + private static bool NoPointSameXzLocation(IEnumerable collection, PipingPoint point) + { + return !collection.Any( + p => Math.Abs(p.X - point.X) < PipingConstants.Accuracy && + Math.Abs(p.Z - point.Z) < PipingConstants.Accuracy); + } + + /// + /// Determines whether the ditch (when defined) is defined correct. It is correct if all + /// four points are there and they are in the proper left to right order. + /// + /// + /// true if its is else false + /// + private bool IsDitchCorrect() + { + if (DitchDikeSide != null || BottomDitchDikeSide != null || BottomDitchPolderSide != null || DitchPolderSide != null) + { + // At least on ditch point given, so they must now all exist + if (DitchDikeSide == null || BottomDitchDikeSide == null || BottomDitchPolderSide == null || DitchPolderSide == null) + { + return false; + } + // they must be left to right + if (BottomDitchDikeSide.X < DitchDikeSide.X || BottomDitchPolderSide.X < BottomDitchDikeSide.X || + DitchPolderSide.X < BottomDitchPolderSide.X) + { + return false; + } + } + // no ditch or ditch is ok + return true; + } + + /// + /// Determines whether the characteristic points (when defined) are ordered correct. + /// Note that the ditch is checked separatly. + /// + /// + private bool AreAllCharacteristicPointsOrdered() + { + if (ShoulderBaseInside != null) + { + if (DikeToeAtPolder != null && DikeToeAtPolder.X < ShoulderBaseInside.X) + { + return false; + } + if (DitchDikeSide != null && DitchDikeSide.X < ShoulderBaseInside.X) + { + return false; + } + } + if (DikeToeAtPolder != null) + { + if (DitchDikeSide != null && DitchDikeSide.X < DikeToeAtPolder.X) + { + return false; + } + } + return true; + } + } +} Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingPoint.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingPoint.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingPoint.cs (revision 253) @@ -0,0 +1,71 @@ +namespace Deltares.DamPiping +{ + /// + /// Class to hold points (3D) for the Piping kernel + /// + public class PipingPoint + { + /// + /// Initializes a new instance of the class. + /// + public PipingPoint() + { + Type = PipingCharacteristicPointType.None; + } + + /// + /// Initializes a new instance of the class. + /// + /// The point. + /// + public PipingPoint(PipingPoint point, PipingCharacteristicPointType type = PipingCharacteristicPointType.None) + { + X = point.X; + Y = point.Y; + Z = point.Z; + Type = type; + } + + /// + /// Initializes a new instance of the class. + /// + /// The x. + /// The y. + /// The z. + /// + public PipingPoint(double x, double y, double z, PipingCharacteristicPointType type = PipingCharacteristicPointType.None) + { + X = x; + Y = y; + Z = z; + Type = type; + } + + /// + /// Gets or sets the x. + /// + /// + /// The x. + /// + public double X { get; set; } + + /// + /// Gets or sets the y. + /// + /// + /// The y. + /// + public double Y { get; set; } + + /// + /// Gets or sets the z. + /// + /// + /// The z. + /// + public double Z { get; set; } + + public PipingCharacteristicPointType Type { get; set; } + + } +} Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Vector3D.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Vector3D.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Vector3D.cs (revision 253) @@ -0,0 +1,247 @@ +using System; +using System.ComponentModel; +using System.Globalization; + + +namespace Deltares.DamPiping +{ + [TypeConverter(typeof(Vector3D))] + public class Vector3D : ExpandableObjectConverter + { + private double x; + private double y; + private double z; + + public double X + { + get + { + return this.x; + } + set + { + this.x = value; + } + } + + public double Y + { + get + { + return this.y; + } + set + { + this.y = value; + } + } + + public double Z + { + get + { + return this.z; + } + set + { + this.z = value; + } + } + + public Vector3D() + { + this.x = this.y = this.z = 0.0; + } + + public Vector3D(double aX, double aY, double aZ) + { + this.x = aX; + this.y = aY; + this.z = aZ; + } + + public Vector3D(Point3D aPoint) + { + this.x = aPoint.X; + this.y = aPoint.Y; + this.z = aPoint.Z; + } + + public Vector3D(double aValue1, Vector3D aVector1, double aValue2, Vector3D aVector2) + { + this.x = aValue1 * aVector1.x + aValue2 * aVector2.x; + this.y = aValue1 * aVector1.y + aValue2 * aVector2.y; + this.z = aValue1 * aVector1.z + aValue2 * aVector2.z; + } + + public Vector3D(double aValue1, Vector3D aVector1, double aValue2, Vector3D aVector2, double aValue3, Vector3D aVector3, double aValue4, Vector3D aVector4) + { + this.x = aValue1 * aVector1.X + aValue2 * aVector2.X + aValue3 * aVector3.X + aValue4 * aVector4.X; + this.y = aValue1 * aVector1.Y + aValue2 * aVector2.Y + aValue3 * aVector3.Y + aValue4 * aVector4.Y; + this.z = aValue1 * aVector1.Z + aValue2 * aVector2.Z + aValue3 * aVector3.Z + aValue4 * aVector4.Z; + } + + public Vector3D(Vector3D aVector) + { + this.x = aVector.x; + this.y = aVector.y; + this.z = aVector.z; + } + + public static Vector3D operator -(Vector3D aVector) + { + return new Vector3D(-aVector.x, -aVector.y, -aVector.z); + } + + public static Vector3D operator +(Vector3D aVector2, Vector3D aVector1) + { + return new Vector3D(aVector2.x + aVector1.x, aVector2.y + aVector1.y, aVector2.z + aVector1.z); + } + + public static Vector3D operator -(Vector3D aVector2, Vector3D aVector1) + { + return new Vector3D(aVector2.x - aVector1.x, aVector2.y - aVector1.y, aVector2.z - aVector1.z); + } + + public static Vector3D operator *(double aValue, Vector3D aVector) + { + return new Vector3D(aValue * aVector.x, aValue * aVector.y, aValue * aVector.z); + } + + public static Vector3D operator *(Vector3D aVector, double aValue) + { + return new Vector3D(aValue * aVector.x, aValue * aVector.y, aValue * aVector.z); + } + + public static double operator ~(Vector3D aVector) + { + return Math.Sqrt(aVector.x * aVector.x + aVector.y * aVector.y + aVector.z * aVector.z); + } + + public static double operator |(Vector3D aVector2, Vector3D aVector1) + { + return aVector2.x * aVector1.x + aVector2.y * aVector1.y + aVector2.z * aVector1.z; + } + + public static Vector3D operator *(Vector3D aVector1, Vector3D aVector2) + { + return new Vector3D(aVector1.y * aVector2.z - aVector1.z * aVector2.y, aVector1.z * aVector2.x - aVector1.x * aVector2.z, aVector1.x * aVector2.y - aVector1.y * aVector2.x); + } + + public static Point3D operator *(Point3D aPoint, Vector3D aVector) + { + return new Point3D(aPoint.Y * aVector.z - aPoint.Z * aVector.y, aPoint.Z * aVector.x - aPoint.X * aVector.z, aPoint.X * aVector.y - aPoint.Y * aVector.x); + } + + public static Point3D operator +(Point3D aPoint, Vector3D aVector) + { + return new Point3D(aPoint.X + aVector.X, aPoint.Y + aVector.Y, aPoint.Z + aVector.Z); + } + + public static Point3D operator +(Vector3D aVector, Point3D aPoint) + { + return new Point3D(aVector.X + aPoint.X, aVector.Y + aPoint.Y, aVector.Z + aPoint.Z); + } + + public static Point3D operator -(Point3D aPoint, Vector3D aVector) + { + return new Point3D(aPoint.X - aVector.X, aPoint.Y - aVector.Y, aPoint.Z - aVector.Z); + } + + public static Vector3D operator !(Vector3D aVector) + { + Vector3D vector3D = new Vector3D(aVector); + vector3D.Normalize(); + return vector3D; + } + + public void Init(double aX, double aY, double aZ) + { + this.x = aX; + this.y = aY; + this.z = aZ; + } + + public void Init(Vector3D aVector) + { + this.x = aVector.x; + this.y = aVector.y; + this.z = aVector.z; + } + + public void Normalize() + { + double num = Math.Sqrt(this.x * this.x + this.y * this.y + this.z * this.z); + if (num == 0.0) + return; + this.x = this.x / num; + this.y = this.y / num; + this.z = this.z / num; + } + + public Vector3D Add(Vector3D aVector) + { + this.x += aVector.x; + this.y += aVector.y; + this.z += aVector.z; + return this; + } + + public bool IsZero() + { + if (Math.Abs(this.x - 0.0) < 0.001 && Math.Abs(this.y - 0.0) < 0.001) + return Math.Abs(this.z - 0.0) < 0.001; + return false; + } + + public Vector3D RotateVectorAroundAxis(Vector3D aAxis, double aFi) + { + Vector3D vector3D1 = (this | aAxis) * aAxis; + Vector3D vector3D2 = this - vector3D1; + Vector3D vector3D3 = aAxis * vector3D2; + double num1 = Math.Cos(aFi); + double num2 = Math.Sin(aFi); + this.Init(vector3D2 * num1 + vector3D3 * num2 + vector3D1); + return this; + } + + public bool FromString(string aString, params char[] aSeparator) + { + string[] strArray = aString.Split(aSeparator); + if (strArray.Length != 3) + return false; + this.x = Convert.ToDouble(strArray[0]); + this.y = Convert.ToDouble(strArray[1]); + this.z = Convert.ToDouble(strArray[2]); + return true; + } + + public string ToStringSep(char aSeparator) + { + return this.x.ToString() + (object)aSeparator + this.y.ToString() + (object)aSeparator + this.z.ToString(); + } + + public string ToStringSep(IFormatProvider aFormatProvider, char aSeparator) + { + return this.x.ToString(aFormatProvider) + (object)aSeparator + this.y.ToString(aFormatProvider) + (object)aSeparator + this.z.ToString(aFormatProvider); + } + + public override bool Equals(object aObject) + { + Vector3D vector3D = aObject as Vector3D; + if (Math.Abs(this.x - vector3D.x) < 0.001 && Math.Abs(this.y - vector3D.y) < 0.001) + return Math.Abs(this.z - vector3D.z) < 0.001; + return false; + } + + public override object ConvertTo(ITypeDescriptorContext aContext, CultureInfo aCulture, object aValue, Type aDestinationType) + { + return (object)""; + } + + public override int GetHashCode() + { + return base.GetHashCode(); + } + } +} Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Point3D.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Point3D.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Point3D.cs (revision 253) @@ -0,0 +1,222 @@ +using System; +using System.ComponentModel; +using System.Globalization; + +namespace Deltares.DamPiping +{ + [TypeConverter(typeof(Point3D))] + public class Point3D : ExpandableObjectConverter, IComparable + { + private double x; + private double y; + private double z; + + [Description("Set X coordinate")] + public double X + { + get + { + return this.x; + } + set + { + this.x = value; + } + } + + [Description("Set Y coordinate")] + public double Y + { + get + { + return this.y; + } + set + { + this.y = value; + } + } + + [Description("Set Z coordinate")] + public double Z + { + get + { + return this.z; + } + set + { + this.z = value; + } + } + + public Point3D() + { + this.x = 0.0; + this.y = 0.0; + this.z = 0.0; + } + + public Point3D(Point3D point) + { + this.x = point.X; + this.y = point.Y; + this.z = point.Z; + } + + public Point3D(double aX, double aY, double aZ) + { + this.x = aX; + this.y = aY; + this.z = aZ; + } + + public Point3D(double aValue1, Point3D aPoint1, double aValue2, Point3D aPoint2) + { + this.x = aValue1 * aPoint1.X + aValue2 * aPoint2.X; + this.y = aValue1 * aPoint1.Y + aValue2 * aPoint2.Y; + this.z = aValue1 * aPoint1.Z + aValue2 * aPoint2.Z; + } + + public Point3D(double aValue1, Point3D aPoint1, double aValue2, Point3D aPoint2, double aValue3, Point3D aPoint3) + { + this.x = aValue1 * aPoint1.X + aValue2 * aPoint2.X + aValue3 * aPoint3.X; + this.y = aValue1 * aPoint1.Y + aValue2 * aPoint2.Y + aValue3 * aPoint3.Y; + this.z = aValue1 * aPoint1.Z + aValue2 * aPoint2.Z + aValue3 * aPoint3.Z; + } + + public Point3D(double aValue1, Point3D aPoint1, double aValue2, Point3D aPoint2, double aValue3, Point3D aPoint3, double aValue4, Point3D aPoint4) + { + this.x = aValue1 * aPoint1.X + aValue2 * aPoint2.X + aValue3 * aPoint3.X + aValue4 * aPoint4.X; + this.y = aValue1 * aPoint1.Y + aValue2 * aPoint2.Y + aValue3 * aPoint3.Y + aValue4 * aPoint4.Y; + this.z = aValue1 * aPoint1.Z + aValue2 * aPoint2.Z + aValue3 * aPoint3.Z + aValue4 * aPoint4.Z; + } + + public static Point3D operator +(Point3D aPoint1, Point3D aPoint2) + { + return new Point3D(aPoint1.X + aPoint2.x, aPoint1.Y + aPoint2.y, aPoint1.Z + aPoint2.Z); + } + + public static Point3D operator -(Point3D aPoint1, Point3D aPoint2) + { + return new Point3D(aPoint1.X - aPoint2.x, aPoint1.Y - aPoint2.y, aPoint1.Z - aPoint2.Z); + } + + public static Point3D operator *(double scale, Point3D point) + { + return point * scale; + } + + public static Point3D operator *(Point3D point, double scale) + { + return new Point3D(point.x * scale, point.y * scale, point.z * scale); + } + + public static Point3D operator /(double scale, Point3D point) + { + return point / scale; + } + + public static Point3D operator /(Point3D point, double scale) + { + return point * (1.0 / scale); + } + + public bool Compare(Point3D aPoint, double aDistance) + { + return Math.Pow(this.x - aPoint.X, 2.0) + Math.Pow(this.y - aPoint.Y, 2.0) + Math.Pow(this.z - aPoint.Z, 2.0) <= Math.Pow(aDistance, 2.0); + } + + public void Init(double aX, double aY, double aZ) + { + this.x = aX; + this.y = aY; + this.z = aZ; + } + + public void Init(Point3D aPoint) + { + this.x = aPoint.X; + this.y = aPoint.Y; + this.z = aPoint.Z; + } + + public void Add(Point3D aPoint) + { + this.x += aPoint.X; + this.y += aPoint.Y; + this.z += aPoint.Z; + } + + public void Mul(double aValue) + { + this.x *= aValue; + this.y *= aValue; + this.z *= aValue; + } + + public bool FromString(string aString, params char[] aSeperator) + { + string[] strArray = aString.Split(aSeperator); + if (strArray.Length != 3) + return false; + this.x = Convert.ToDouble(strArray[0]); + this.y = Convert.ToDouble(strArray[1]); + this.z = Convert.ToDouble(strArray[2]); + return true; + } + + public Point2D GetPointXZ() + { + return new Point2D(this.X, this.Z); + } + + public double LengthSquared() + { + return this.DotProduct(this); + } + + public double Length() + { + return Math.Sqrt(this.LengthSquared()); + } + + public double DotProduct(Point3D p) + { + return this.X * p.X + this.Y * p.Y + this.Z * p.Z; + } + + public Point3D CrossProduct(Point3D p) + { + return new Point3D(this.Y * p.Z - this.Z * p.Y, this.Z * p.X - this.X * p.Z, this.X * p.Y - this.Y * p.X); + } + + public override bool Equals(object aObject) + { + Point3D point3D = aObject as Point3D; + if (point3D == null || Math.Abs(this.x - point3D.x) >= 0.001 || Math.Abs(this.y - point3D.y) >= 0.001) + return false; + return Math.Abs(this.z - point3D.z) < 0.001; + } + + public override object ConvertTo(ITypeDescriptorContext aContext, CultureInfo aCulture, object aValue, Type aDestinationType) + { + return (object)""; + } + + public override string ToString() + { + return this.X.ToString("F2") + " " + this.Y.ToString("F2") + " " + this.Z.ToString("F2"); + } + + public override int GetHashCode() + { + return base.GetHashCode(); + } + + public int CompareTo(object obj) + { + return this.X.CompareTo(((Point3D)obj).X); + } + } +} Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Deltares.DamPiping.csproj =================================================================== diff -u -r251 -r253 --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Deltares.DamPiping.csproj (.../Deltares.DamPiping.csproj) (revision 251) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Deltares.DamPiping.csproj (.../Deltares.DamPiping.csproj) (revision 253) @@ -41,9 +41,27 @@ + + + + + + + + + + + + + + + + + + Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingProfile.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingProfile.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingProfile.cs (revision 253) @@ -0,0 +1,217 @@ +using System; +using System.Collections.Generic; +using System.Linq; + +namespace Deltares.DamPiping +{ + /// + /// PipingProfile Exception class + /// + public class PipingProfileException : Exception + { + /// + /// Initializes a new instance of the class. + /// + /// The message that describes the error. + public PipingProfileException(string message) : base(message) { } + } + + /// + /// Class to hold profiles (1D) for the Piping kernel + /// + public class PipingProfile + { + private List layers = new List(); + + /// + /// Gets or sets the name. + /// + /// + /// The name. + /// + public string Name { get; set; } + + /// + /// Gets or sets the layers. + /// + /// + /// The layers. + /// + public List Layers + { + get + { + return layers; + } + set + { + layers = value; + } + } + + /// + /// Assigns the specified profile. + /// + /// The profile. + public void Assign(PipingProfile profile) + { + BottomLevel = profile.BottomLevel; + foreach (var pipingLayer in profile.Layers) + { + var layer = new PipingLayer(); + Layers.Add(layer); + layer.Assign(pipingLayer); + } + } + + /// + /// Gets the top level of the profile. + /// + public virtual double TopLevel + { + get + { + return layers.Any() ? layers.First().TopLevel : BottomLevel; + } + } + + /// + /// Gets or sets the bottom level of the profile. + /// + public double BottomLevel { get; set; } + + /// + /// Gets the top aquifer layer. + /// + /// + /// The top aquifer layer. + /// + public PipingLayer TopAquiferLayer + { + get + { + IList sortedLayers = Layers.OrderByDescending(l => l.TopLevel).ToList(); + return sortedLayers.FirstOrDefault(layer => layer.IsAquifer); + } + } + + /// + /// Gets the bottom aquifer layer. + /// + /// + /// The bottom aquifer layer. + /// + public PipingLayer BottomAquiferLayer + { + get + { + IList sortedLayers = Layers.OrderBy(l => l.TopLevel).ToList(); + PipingLayer aquiferLayer = null; + var aquiferIndex = -1; + + // Search deepest aquifer layer + for (var layerIndex = 0; layerIndex < sortedLayers.Count; layerIndex++) + { + var layer = sortedLayers[layerIndex]; + if (layer.IsAquifer) + { + aquiferIndex = layerIndex; + aquiferLayer = layer; + break; + } + } + + // aquifer may consists of more then 1 connected (aquifer) layers + // Search all layers above the first found aquifer to find top aquifer layer + if (aquiferIndex >= 0) + { + for (var layerIndex = aquiferIndex + 1; layerIndex < sortedLayers.Count; layerIndex++) + { + var layer = sortedLayers[layerIndex]; + if (layer.IsAquifer) + { + aquiferLayer = layer; + } + else + { + break; + } + } + } + return aquiferLayer; + } + } + + /// + /// Gets (calculates) the height for a given layer in the profile + /// + /// The layer to process + /// The height + public double GetLayerHeight(PipingLayer soilLayer) + { + return GetLayerHeight(layers.IndexOf(soilLayer)); + } + + internal void Validate() + { + if (layers.Count == 0) + { + var message = Helper.Translate("ProfileHasNoLayers"); + throw new PipingProfileException(message); + } + if (!AreLayersOrderedTopToBottom()) + { + var message = Helper.Translate("LayersNotOrderedFromTopToBottom"); + throw new PipingProfileException(message); + } + if (!IsBottomLevelDeepEnough()) + { + var format = Helper.Translate("BottomLevelNotDeepEnough"); + var topLevel = layers[layers.Count - 1].TopLevel; + var message = String.Format(format, BottomLevel, PipingConstants.MinimumLayerThickness, topLevel); + throw new PipingProfileException(message); + } + } + + /// + /// Gets the height of the layer. + /// + /// Index of the layer. + /// The height + private double GetLayerHeight(int layerIndex) + { + var soilLayer = Layers[layerIndex]; + var soilLayerBelow = (layerIndex < layers.Count - 1) ? layers[layerIndex + 1] : null; + var levelBelow = (soilLayerBelow != null) ? soilLayerBelow.TopLevel : BottomLevel; + return soilLayer.TopLevel - levelBelow; + } + + /// + /// Determines whether the bottom level is deep. + /// + /// + /// true when it is deep enough else false. + /// + private bool IsBottomLevelDeepEnough() + { + return layers[layers.Count - 1].TopLevel - BottomLevel >= PipingConstants.MinimumLayerThickness; + } + + /// + /// Determines whether the layers are ordered top to bottom. + /// + /// + /// true if they are else false. + /// + private bool AreLayersOrderedTopToBottom() + { + for (int i = 1; i < layers.Count; i++) + { + if (layers[i-1].TopLevel - layers[i].TopLevel < PipingConstants.MinimumLayerThickness) + return false; + } + return true; + } + + } +} Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/EffectiveThicknessCalculator.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/EffectiveThicknessCalculator.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/EffectiveThicknessCalculator.cs (revision 253) @@ -0,0 +1,515 @@ +using System; +using System.Collections.Generic; + +namespace Deltares.DamPiping +{ + /// + /// Exception class for SoilVolumicMassCalculator + /// + public class EffectiveThicknessCalculatorException : Exception + { + public EffectiveThicknessCalculatorException(string message) : base(message) {} + } + + /// + /// Class for calculation of effective height and effective stress + /// + public class EffectiveThicknessCalculator + { + // Input parameters + private double bottomAquiferLayerEffectiveHeight = Double.NaN; + private double bottomAquiferLayerEffectiveStress = Double.NaN; + private double inbetweenAquiferLayerEffectiveHeight = Double.NaN; + private double inbetweenAquiferLayerEffectiveStress = Double.NaN; +// private LanguageType language = LanguageType.Dutch; + +// /// +// /// Language for (validation) messages +// /// +// public LanguageType Language +// { +// get { return language; } +// set { language = value; } +// } + + /// + /// Initializes a new instance of the class. + /// + public EffectiveThicknessCalculator() + { + VolumicWeightOfWater = PipingConstants.UnitWeightOfWater; + } + + /// + /// Gets or sets the volumic weight of water. + /// + /// + /// The volumic weight of water. + /// + public double VolumicWeightOfWater { get; set; } + + /// + /// Gets or sets the soil profile. + /// + /// + /// The soil profile. + /// + public PipingProfile SoilProfile { get; set; } + + /// + /// Gets or sets the surface line. + /// + /// + /// The surface line. + /// + public PipingSurfaceLine SurfaceLine { get; set; } + + /// + /// Gets or sets the exit point x coordinate. + /// + /// + /// The exit point x coordinate. + /// + public double ExitPointXCoordinate { get; set; } + + /// + /// Gets or sets the phreatic level. + /// + /// + /// The phreatic level. + /// + public double PhreaticLevel { get; set; } + + // In the docs this is called the PolderLevel + + // Output parameters + /// + /// Gets the effective stress. + /// + /// + /// The effective stress. + /// + public double EffectiveStress + { + // For now we use the highest aquifer for the determination of EffectiveStress + get + { + return SoilProfile.TopAquiferLayer != null ? inbetweenAquiferLayerEffectiveStress : bottomAquiferLayerEffectiveStress; + } + } + + /// + /// Gets the height of the effective. + /// + /// + /// The height of the effective. + /// + public double EffectiveHeight + { + // For now we use the highest aquifer for the determination of EffectiveHeight + get + { + return SoilProfile.TopAquiferLayer != null ? inbetweenAquiferLayerEffectiveHeight : bottomAquiferLayerEffectiveHeight; + } + } + + /// + /// Calculates the output parameters from the input. + /// + public void Calculate() + { +// var restoreLanguage = LocalizationManager.CurrentLanguage; +// LocalizationManager.CurrentLanguage = Language; +// var dumper = new XmlDumper(); + try + { +// DumpThis(dumper); + + ThrowWhenSoilProfileIsNull(); + ThrowWhenSurfaceLineIsNull(); + inbetweenAquiferLayerEffectiveStress = Double.NaN; + inbetweenAquiferLayerEffectiveHeight = Double.NaN; + if (SoilProfile.TopAquiferLayer != null) + { + // Calculate for the inbetween aquifer + CalculateEffectiveHeightAndStress(SoilProfile.TopAquiferLayer.TopLevel, + out inbetweenAquiferLayerEffectiveHeight, out inbetweenAquiferLayerEffectiveStress); + } + // Calculate for the bottom aquifer + CalculateEffectiveHeightAndStress(SoilProfile.BottomAquiferLayer.TopLevel, + out bottomAquiferLayerEffectiveHeight, out bottomAquiferLayerEffectiveStress); + } + finally + { +// DumpThis(dumper); +// LocalizationManager.CurrentLanguage = restoreLanguage; + } + } + + /// + /// Validates the input + /// + /// a filled list when errors are found else an empty list + public virtual List Validate() + { +// var restoreLanguage = LocalizationManager.CurrentLanguage; +// LocalizationManager.CurrentLanguage = Language; + var errors = new List(); +// var dumper = new XmlDumper(); + try + { +// DumpThis(dumper); + + if (double.IsNaN(VolumicWeightOfWater)) + LogParameterIsNaN(errors, "VolumicWeightOfWater"); + if (double.IsNaN(PhreaticLevel)) + LogParameterIsNaN(errors, "PhreaticLevel"); + if (double.IsNaN(ExitPointXCoordinate)) + LogParameterIsNaN(errors, "ExitPointXCoordinate"); + + + if (SoilProfile == null) + { + errors.Add(Helper.Translate("SoilprofileNotDefined")); + return errors; + } + if (SurfaceLine == null) + { + errors.Add(Helper.Translate("SurfaceLineNotDefined")); + return errors; + } + if (SurfaceLine.Points.Count < 2) + { + errors.Add(Helper.Translate("SurfaceLineAtLeastTwoPoints")); + return errors; + } + + try + { + SurfaceLine.Validate(); + } + catch (Exception e) + { + errors.Add(e.Message); + } + + try + { + SoilProfile.Validate(); + } + catch (Exception e) + { + errors.Add(e.Message); + } + double? topOfLayerToBeEvaluated = null; + if (SoilProfile.BottomAquiferLayer != null) + { + topOfLayerToBeEvaluated = SoilProfile.BottomAquiferLayer.TopLevel; + } + if (SoilProfile.TopAquiferLayer != null) + { + topOfLayerToBeEvaluated = SoilProfile.TopAquiferLayer.TopLevel; + } + if (topOfLayerToBeEvaluated == null) + { + errors.Add(Helper.Translate("SoilProfileHasNoAquifer")); + return errors; + } + } + finally + { +// DumpThis(dumper); +// LocalizationManager.CurrentLanguage = restoreLanguage; + } + + return errors; + } + + /// + /// Gets the effective width ditch. + /// Note: Assumes a valid ditch (i.e. Xcoors for diketop, dikebottom, polderbottom, poldertop are all ascending + /// Made public pure to be able to unit test. + /// + /// the effective width ditch + public double GetEffectiveWidthDitch() + { + if (SurfaceLine.DitchPolderSide.Z > + SurfaceLine.DitchDikeSide.Z) + { + // Ditch at polder side higher than dike side so intersect horizontaly from dike to polder side + // p1,p2 = line at polderside + var p1 = new Point2D(SurfaceLine.DitchPolderSide.X, + SurfaceLine.DitchPolderSide.Z); + var p2 = new Point2D(SurfaceLine.BottomDitchPolderSide.X, + SurfaceLine.BottomDitchPolderSide.Z); + // p3, p4 = line from dike side level horizontal to polderside + var p3 = new Point2D(SurfaceLine.DitchDikeSide.X, + SurfaceLine.DitchDikeSide.Z); + var p4 = new Point2D(SurfaceLine.DitchPolderSide.X, + SurfaceLine.DitchDikeSide.Z); + Point2D ip = null; + Routines2D.DetermineIf2DLinesIntersectStrickly(p1, p2, p3, p4, ref ip); + return ip.X - SurfaceLine.DitchDikeSide.X; + } + else + { + if (SurfaceLine.DitchPolderSide.Z < + SurfaceLine.DitchDikeSide.Z) + { + // Ditch at dike side higher than polder side so intersect horizontaly from polder to dike side + // p1,p2 = line at dikeside + var p1 = new Point2D(SurfaceLine.DitchDikeSide.X, + SurfaceLine.DitchDikeSide.Z); + var p2 = new Point2D(SurfaceLine.BottomDitchDikeSide.X, + SurfaceLine.BottomDitchDikeSide.Z); + // p3, p4 = line from polder side level horizontal to dikeside + var p3 = new Point2D(SurfaceLine.DitchPolderSide.X, + SurfaceLine.DitchPolderSide.Z); + var p4 = new Point2D(SurfaceLine.DitchDikeSide.X, + SurfaceLine.DitchPolderSide.Z); + Point2D ip = null; + Routines2D.DetermineIf2DLinesIntersectStrickly(p1, p2, p3, p4, ref ip); + return SurfaceLine.DitchPolderSide.X - ip.X; + } + } + return SurfaceLine.DitchPolderSide.X - SurfaceLine.DitchDikeSide.X; + } + + + private void LogParameterIsNaN(IList list, string paramName) + { + var msg = string.Format(Helper.Translate("NaNParameterError"), paramName); + list.Add(msg); + } + + /// + /// Calculates the effective height and stress at a given topLevel. + /// + /// The top level. + /// Height of the effective. + /// The effective stress. + private void CalculateEffectiveHeightAndStress(double topLevel, out double effectiveHeight, out double effectiveStress) + { + var exitPoint = new PipingPoint(ExitPointXCoordinate, 0, 0); // Y and Z are not important, because only X will be used + var surfaceLevel = SurfaceLine.GetZatX(ExitPointXCoordinate); + // The minimum thickness of the cover layer is only to be used for point that are actual situated in the ditch. + // Points besides the ditch should be calculated as is. + if (IsPointPartOfDitch(exitPoint)) + { + effectiveHeight = DetermineThicknessCoverLayerAtDitch(topLevel); + } + else + { + effectiveHeight = surfaceLevel - topLevel; + } + + var soilVolumeMassCalculator = new SoilVolumicMassCalculator + { + MinimumThicknessCoverLayer = effectiveHeight, + PhreaticLevel = PhreaticLevel, + SurfaceLevel = topLevel + effectiveHeight, + TopOfLayerToBeEvaluated = topLevel, + SoilProfile = SoilProfile, + VolumicWeightOfWater = VolumicWeightOfWater + }; + effectiveStress = soilVolumeMassCalculator.CalculateEffectiveStress(); + } + + /// + /// Determines whether the point is situated within the ditch + /// + /// + /// + private bool IsPointPartOfDitch(PipingPoint point) + { + if (SurfaceLine.HasDitch() && SurfaceLine.DitchDikeSide.X <= point.X && + SurfaceLine.DitchPolderSide.X >= point.X) + { + return true; + } + return false; + } + + /// + /// Determines the value for the thickness of the coverLayer at a ditch. + /// This is described in paragraph 3.3.3 of the Requirements and Functional Design + /// + /// + /// the thickness of the coverLayer at a ditch + private double DetermineThicknessCoverLayerAtDitch(double topAquifer) + { + double minimumThickness = 0.0; + if (SurfaceLine.HasDitch()) + { + // find the middle of the bottom of the ditch to determine the start point of the check (startPointUnderDitch). + var ditchBottomCenterX = (SurfaceLine.BottomDitchDikeSide.X + + SurfaceLine.BottomDitchPolderSide.X)/2.0; + // Determine total cover layer thickness outside the ditch + double hOutside = Math.Min(SurfaceLine.DitchDikeSide.Z, + SurfaceLine.DitchPolderSide.Z) - topAquifer; + // Determine width of ditch at surface level + double bSurfaceLevel = GetEffectiveWidthDitch(); + // Determine width of ditch at bottom level + double bBottom = SurfaceLine.BottomDitchPolderSide.X - + SurfaceLine.BottomDitchDikeSide.X; + // Determine cover layer thickness at centre of ditch + double hCentre = SurfaceLine.GetZatX(ditchBottomCenterX) - topAquifer; + + if (bSurfaceLevel < hOutside) + { + // case 1 in the description + minimumThickness = hOutside; + } + else + { + if (bBottom > hCentre) + { + // case 2 in the description + minimumThickness = hCentre; + } + else + { + // case 3 in the description + minimumThickness = DetermineMinimumThicknessInBetween(topAquifer, ditchBottomCenterX); + } + } + } + return minimumThickness; + } + + /// + /// Determines the value for the minimum thickness of the coverLayer. Instead of using just the actual thickness of the coverLayer (which can be determined + /// from the profile and surfaceline), the thickness of the coverlayer has a prescribed minimum value which can exceed the actual thickness. The procedure + /// is to set up lines using an gradient of 2:1 from the top of the aquifer just under the center of the ditch to the left (dike side) and right (polderside). + /// Then determine where these lines intersects the surface lines. If this is to the left and/or to the right of the ditch, keep the top of the ditch at that + /// side as the result for that side. If it intersects the ditch, keep the lowest height of intersection points (there can be more than one) as result for + /// that side. The overall result is the lowest of the results left and right. + /// This is described in paragraph 3.3.3 of the Requirements and Functional Design + /// + /// + /// + /// the value for the minimum thickness of the coverLayer + private double DetermineMinimumThicknessInBetween(double topAquifer, double ditchBottomCenterX) + { + var startPointUnderDitch = new PipingPoint(ditchBottomCenterX, 0, topAquifer); + // set a fictive point left (dikeside) of the startpoint with a gradient of 2:1 and create the line from this point to the startpoint + var pointLeft = new PipingPoint(ditchBottomCenterX - 1000, 0, topAquifer + 2000); + // Determine the intersection(s) of the line with the surfaceline and from that the intersection point with the lowest Z value + // If the intersection point lies left of ditch, take the top of ditch at dike side as result. + var pointsLeft = SurfaceLine.IntersectionPointsXzWithLineXz(pointLeft, startPointUnderDitch); + double minHeightLeft; + var lowestPoint = new PipingPoint(0, 0, double.MaxValue); + foreach (var pipingPoint in pointsLeft) + { + if (pipingPoint.Z <= lowestPoint.Z) + { + lowestPoint = pipingPoint; + } + } + if (lowestPoint.X < SurfaceLine.DitchDikeSide.X) + { + minHeightLeft = SurfaceLine.DitchDikeSide.Z; + } + else + { + minHeightLeft = lowestPoint.Z; + } + + // Do the same at the right side of the ditch + var pointRight = new PipingPoint(ditchBottomCenterX + 1000, 0, topAquifer + 2000); + var pointsRight = SurfaceLine.IntersectionPointsXzWithLineXz(pointRight, startPointUnderDitch); + double minHeightRight; + lowestPoint = new PipingPoint(0, 0, double.MaxValue); + foreach (var pipingPoint in pointsRight) + { + if (pipingPoint.Z <= lowestPoint.Z) + { + lowestPoint = pipingPoint; + } + } + + if (lowestPoint.X > SurfaceLine.DitchPolderSide.X) + { + minHeightRight = SurfaceLine.DitchPolderSide.Z; + } + else + { + minHeightRight = lowestPoint.Z; + } + + // Take the lowest of right and left + double minimumThickness = Math.Min(minHeightLeft, minHeightRight); + minimumThickness = minimumThickness - topAquifer; + // Never lower than 0 + minimumThickness = Math.Max(minimumThickness, 0.0); + return minimumThickness; + } + + /// + /// Check precondition: Throws the when soil profile is null. + /// + /// The soilprofile is not defined + private void ThrowWhenSoilProfileIsNull() + { + if (SoilProfile == null) + { + throw new EffectiveThicknessCalculatorException(Helper.Translate("SoilprofileNotDefined")); + } + } + + /// + /// Check precondition: Throws the when surfaceline is null. + /// + /// The soilprofile is not defined + private void ThrowWhenSurfaceLineIsNull() + { + if (SurfaceLine == null) + { + throw new EffectiveThicknessCalculatorException(Helper.Translate("SurfaceLineNotDefined")); + } + } + +// #region DebugOut +// private bool isDebugOutputEnabled = false; +// +// protected void DumpThis(XmlDumper dumper) +// { +// if (isDebugOutputEnabled) +// { +// dumper.DumpObject(this); +// } +// } +// +// /// +// /// Enables debug-dumping to XML +// /// +// public void EnableDebugOutput() +// { +// isDebugOutputEnabled = true; +// } +// +// /// +// /// Fills the data from debug file. +// /// +// public EffectiveThicknessCalculator FillFromDebugFile(string content) +// { +// var serializer = new Deltares.Standard.IO.Xml.XmlDeserializer(); +// if (!string.IsNullOrEmpty(content)) +// { +// var res = (EffectiveThicknessCalculator)serializer.XmlDeserializeFromString(content); +// return res; +// } +// return null; +// } +// +// /// +// /// Gets the main result +// /// +// public string GetMainResult() +// { +// return "Eff height = " + EffectiveHeight.ToString(); +// } +// +// #endregion + } +} \ No newline at end of file Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingConstants.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingConstants.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingConstants.cs (revision 253) @@ -0,0 +1,22 @@ +namespace Deltares.DamPiping +{ + /// + /// Class holding all constant default values specific to the WTI Piping kernel + /// + public class PipingConstants + { + public const double UnitWeightOfWater = 9.81; + + public const double GammaSubParticles = 16.5; + public const double BeddingAngleSellmeijerOriginal = 41.0; + public const double BeddingAngleSellmeijerRevised = 37.0; + public const double WhitesDragCoefficient = 0.25; + public const double D70Mean = 2.08e-4; + public const double RExitDefault = 1.00; + public const double RToeDefault = 1.00; + public const double RcDefault = 0.3; + public const double Accuracy = 0.0005; + public const double MinimumLayerThickness = 0.001; + public const double Epsilon = 1e-8; + } +} \ No newline at end of file Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingException.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingException.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingException.cs (revision 253) @@ -0,0 +1,30 @@ +using System; + +namespace Deltares.DamPiping +{ + /// + /// Generic Piping Exception class, based on BaseException + /// + /// + [Serializable] + public class PipingException : BaseException + { + /// + /// PipingException, based on BaseException + /// + public PipingException() {} + + /// + /// PipingException, based on BaseException + /// + /// + public PipingException(string message) : base(message) {} + + /// + /// PipingException, based on BaseException + /// + /// + /// + public PipingException(string message, Exception exception) : base(message, exception) {} + } +} \ No newline at end of file Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingCharacteristicPointType.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingCharacteristicPointType.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingCharacteristicPointType.cs (revision 253) @@ -0,0 +1,13 @@ +namespace Deltares.DamPiping +{ + public enum PipingCharacteristicPointType + { + None, + ShoulderBaseInside, // Insteek binnenberm + DikeToeAtPolder, // Teen dijk binnenwaarts + DitchDikeSide, // Insteek sloot dijkzijde + BottomDitchDikeSide, // Slootbodem dijkzijde + BottomDitchPolderSide, // Slootbodem polderzijde + DitchPolderSide, // Insteek sloot polderzijde + } +} Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Point2D.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Point2D.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Point2D.cs (revision 253) @@ -0,0 +1,118 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace Deltares.DamPiping +{ + public class Point2D + { + public const double Epsilon = 1E-10; + + public double X { get; set; } + + public double Y { get; set; } + + public Point2D() + { + this.X = 0.0; + this.Y = 0.0; + } + + public Point2D(double aX, double aY) + { + this.X = aX; + this.Y = aY; + } + + public Point2D(double aValue1, Point2D aPoint1, double aValue2, Point2D aPoint2) + { + this.X = aValue1 * aPoint1.X + aValue2 * aPoint2.X; + this.Y = aValue1 * aPoint1.Y + aValue2 * aPoint2.Y; + } + + public Point2D(double aValue1, Point2D aPoint1, double aValue2, Point2D aPoint2, double aValue3, Point2D aPoint3, double aValue4, Point2D aPoint4) + { + this.X = aValue1 * aPoint1.X + aValue2 * aPoint2.X + aValue3 * aPoint3.X + aValue4 * aPoint4.X; + this.Y = aValue1 * aPoint1.Y + aValue2 * aPoint2.Y + aValue3 * aPoint3.Y + aValue4 * aPoint4.Y; + } + + public Point2D(Point2D aPoint) + { + this.X = aPoint.X; + this.Y = aPoint.Y; + } + + public static Point2D operator -(Point2D aPoint) + { + return new Point2D(-aPoint.X, -aPoint.Y); + } + + public static Point2D operator +(Point2D aPoint2, Point2D aPoint1) + { + return new Point2D(aPoint2.X + aPoint1.X, aPoint2.Y + aPoint1.Y); + } + + public static Point2D operator -(Point2D aPoint2, Point2D aPoint1) + { + return new Point2D(aPoint2.X - aPoint1.X, aPoint2.Y - aPoint1.Y); + } + + public static Point2D operator *(double aValue1, Point2D aPoint1) + { + return new Point2D(aValue1 * aPoint1.X, aValue1 * aPoint1.Y); + } + + public static Point2D operator *(Point2D aPoint1, double aValue1) + { + return new Point2D(aValue1 * aPoint1.X, aValue1 * aPoint1.Y); + } + + public void Init(double aX, double aY) + { + this.X = aX; + this.Y = aY; + } + + public void Init(Point2D aPoint) + { + this.X = aPoint.X; + this.Y = aPoint.Y; + } + + public Point2D Add(Point2D aPoint) + { + this.X += aPoint.X; + this.Y += aPoint.Y; + return this; + } + + public bool FromString(string aString, params char[] aSeparator) + { + string[] strArray = aString.Split(aSeparator); + if (strArray.Length != 2) + return false; + this.X = Convert.ToDouble(strArray[0]); + this.Y = Convert.ToDouble(strArray[1]); + return true; + } + + public override bool Equals(object obj) + { + Point2D point2D = obj as Point2D; + if (point2D != null && Helper.AlmostEquals(X, point2D.X, 1E-10)) + return Helper.AlmostEquals(Y, point2D.Y, 1E-10); + return false; + } + + public override int GetHashCode() + { + return this.GetType().GetHashCode() + 11 * (int)Math.Round(this.X * 1000000000) + 29 * (int)Math.Round(this.Y * 1000000000); + } + + public override string ToString() + { + return this.X.ToString("F3") + ", " + this.Y.ToString("F3"); + } + } +} Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/BaseException.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/BaseException.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/BaseException.cs (revision 253) @@ -0,0 +1,38 @@ +using System; + +namespace Deltares.DamPiping +{ + /// + /// Generic BaseException class, based on Exception + /// + /// + [Serializable] + public class BaseException : Exception + { + /* + For guidelines and best practices regarding the creation of new exception types, see + http://msdn.microsoft.com/en-us/library/seyhszts.aspx + http://msdn.microsoft.com/en-us/library/87cdya3t.aspx + http://msdn.microsoft.com/en-us/library/vstudio/ms229064(v=vs.100).aspx + http://msdn.microsoft.com/en-us/library/vstudio/ms229007(v=vs.100).aspx + */ + + /// + /// Custom Exception, based on Exception + /// + public BaseException() {} + + /// + /// Custom Exception, based on Exception + /// + /// + public BaseException(string message) : base(message + Environment.NewLine + string.Format("Exception of type {0}", typeof(T).Name)) {} + + /// + /// Custom Exception, based on Exception + /// + /// + /// + public BaseException(string message, Exception inner) : base(message + Environment.NewLine + string.Format("Exception of type {0}", typeof(T).Name), inner) {} + } +} \ No newline at end of file Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/SoilVolumicMassCalculator.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/SoilVolumicMassCalculator.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/SoilVolumicMassCalculator.cs (revision 253) @@ -0,0 +1,369 @@ +using System; +using System.Collections.Generic; +using System.Linq; + +namespace Deltares.DamPiping +{ + /// + /// Exception class for SoilVolumicMassCalculator + /// + public class SoilVolumicMassCalculatorException : Exception + { + public SoilVolumicMassCalculatorException(string message) : base(message) {} + } + + /// + /// Calculation class to determine weight of the soil between TopOfLayerToBeEvaluated and SurfaceLevel + /// + public class SoilVolumicMassCalculator + { + private const double cTolerance = 0.00000001; + + /// + /// Initializes a new instance of the class. + /// + public SoilVolumicMassCalculator() + { + VolumicWeightOfWater = PipingConstants.UnitWeightOfWater; + IsUseOvenDryUnitWeight = false; + } + + public bool IsUseOvenDryUnitWeight { get; set; } + public double VolumicWeightOfWater { get; set; } + public double PhreaticLevel { get; set; } + public double SurfaceLevel { get; set; } + public double TopOfLayerToBeEvaluated { get; set; } + public PipingProfile SoilProfile { get; set; } + public double MinimumThicknessCoverLayer { get; set; } + + /// + /// Calculates the total weight of the soil between SurfaceLevel and TopOfLayerToBeEvaluated, + /// taking into account whether the soil is submerged or not + /// + /// total mass + public double CalculateTotalMass() + { + var weight = double.NaN; +// var dumper = new XmlDumper(); + try + { +// DumpThis(dumper); + + ThrowWhenSoilProfileIsNull(); + ThrowWhenSurfaceLevelNotInsideProfile(); + ThrowWhenTopLayerToBeEvaluatedNotInsideProfile(); + + weight = SoilProfile.Layers + .Where(layer => layer.TopLevel > TopOfLayerToBeEvaluated) + .Sum(layer => GetWeight(layer)); + + // Is there water above the soil layer? If so add the volumic weight + // of water to the soil mass + if (PhreaticLevel > SurfaceLevel) + { + var height = PhreaticLevel - SurfaceLevel; + weight += height*VolumicWeightOfWater; + } + } + finally + { +// DumpThis(dumper); + } + + return weight; + } + + /// + /// Calculates the effective stress of the soil between TopOfLayerToBeEvaluated and SurfaceLevel + /// + /// + public double CalculateEffectiveStress() + { + var weight = double.NaN; +// var dumper = new XmlDumper(); + try + { +// DumpThis(dumper); + + ThrowWhenSoilProfileIsNull(); + ThrowWhenSoilProfileHasNoLayers(); + ThrowWhenSurfaceLevelNotInsideProfile(); + ThrowWhenTopLayerToBeEvaluatedNotInsideProfile(); + + var involvedLayers = SoilProfile.Layers.Where(layer => layer.TopLevel > TopOfLayerToBeEvaluated); + weight = involvedLayers.Sum(layer => GetEffectiveStress(layer)); + // Check to see if thickness of all involved layers is is smaller than the given minimum thickness. If so, + // then add the weight of the "missing" soil(s) by multiplying the difference by the weight of the toplevel. + // Using the weight of the toplevel might be theoreticaly incorrect but this is a good and simple approach + // which has been accorded by Erik Vastenburg. + var thicknessInvolvedLayers = involvedLayers.Sum(layer => GetLayerHeight(layer)); + if (thicknessInvolvedLayers < MinimumThicknessCoverLayer) + { + var layer = involvedLayers.First(); + double bottomLevel = Math.Min(involvedLayers.First().TopLevel, SurfaceLevel); + double height = MinimumThicknessCoverLayer - thicknessInvolvedLayers; + double topLevel = bottomLevel + height; + double factorWet; + double factorDry; + DetermineHeightAndDryAndWetFraction(topLevel, bottomLevel, out factorWet, out factorDry); + weight += GetSoilUnitWeightDry(layer)*factorDry*height + + (layer.BelowPhreaticLevel - VolumicWeightOfWater)*factorWet*height; + } + } + finally + { +// DumpThis(dumper); + } + + return weight; + } + + private void LogParameterIsNaN(IList list, string paramName) + { + var msg = string.Format(Helper.Translate("NaNParameterError"), paramName); + list.Add(msg); + } + + /// + /// Validates the input + /// + /// a filled list when errors are found else an empty list + public List Validate() + { + var errors = new List(); +// var dumper = new XmlDumper(); + try + { +// DumpThis(dumper); + + if (double.IsNaN(VolumicWeightOfWater)) + LogParameterIsNaN(errors, "VolumicWeightOfWater"); + if (double.IsNaN(PhreaticLevel)) + LogParameterIsNaN(errors, "PhreaticLevel"); + if (double.IsNaN(SurfaceLevel)) + LogParameterIsNaN(errors, "SurfaceLevel"); + if (double.IsNaN(TopOfLayerToBeEvaluated)) + LogParameterIsNaN(errors, "TopOfLayerToBeEvaluated"); + if (double.IsNaN(MinimumThicknessCoverLayer)) + LogParameterIsNaN(errors, "MinimumThicknessCoverLayer"); + + if (SoilProfile == null) + { + errors.Add(Helper.Translate("SoilprofileNotDefined")); + return errors; + } + if ((SurfaceLevel - SoilProfile.TopLevel) > cTolerance || + (SurfaceLevel - SoilProfile.BottomLevel) < -cTolerance) + { + errors.Add(Helper.Translate("SurfacelevelNotInsideSoilProfile")); + } + if ((TopOfLayerToBeEvaluated - SoilProfile.TopLevel) > cTolerance || + (TopOfLayerToBeEvaluated - SoilProfile.BottomLevel) < -cTolerance) + { + errors.Add(Helper.Translate("TopLayerNotInsideSoilProfile")); + } + } + finally + { +// DumpThis(dumper); + } + + return errors; + } + + /// + /// Gets the contribution of the weight of a layer, taking in account the phreatic level + /// + /// The layer. + /// + private double GetWeight(PipingLayer layer) + { + double topLevel, bottomLevel; + double height = GetLayerHeight(layer, out topLevel, out bottomLevel); + double factorWet; + double factorDry; + DetermineHeightAndDryAndWetFraction(topLevel, bottomLevel, out factorWet, out factorDry); + + return GetSoilUnitWeightDry(layer)*factorDry*height + layer.BelowPhreaticLevel*factorWet*height; + } + + /// + /// Gets the contribution to the effective stress of a layer, taking in account the phreatic level + /// + /// The layer. + /// + private double GetEffectiveStress(PipingLayer layer) + { + double topLevel, bottomLevel; + double height = GetLayerHeight(layer, out topLevel, out bottomLevel); + double factorWet; + double factorDry; + DetermineHeightAndDryAndWetFraction(topLevel, bottomLevel, out factorWet, out factorDry); + + // If above phreatic line use the dry weight of the soil + // if below phreatic line use the effective submerged weight (gamma_sat - gamma_water) + return GetSoilUnitWeightDry(layer)*factorDry*height + (layer.BelowPhreaticLevel - VolumicWeightOfWater)*factorWet*height; + } + + /// + /// Determines the height and dry and wet fraction of the layer. + /// + /// + /// + /// The factor wet. + /// The factor dry. + private void DetermineHeightAndDryAndWetFraction(double topLevel, double bottomLevel, out double factorWet, out double factorDry) + { + double height = topLevel - bottomLevel; + factorWet = 0.0; + factorDry = 0.0; + + if (topLevel < PhreaticLevel || Helper.AlmostEquals(topLevel, PhreaticLevel)) + { + factorWet = 1.0; + } + else if (bottomLevel > PhreaticLevel || Helper.AlmostEquals(bottomLevel, PhreaticLevel)) + { + factorDry = 1.0; + } + else + { + var pLevel = (topLevel - PhreaticLevel); + factorDry = Helper.AlmostEquals(pLevel, 0.0) ? 0.0 : pLevel/height; + factorWet = 1.0 - factorDry; + } + } + + /// + /// Calculates the height of a layer taking into account the surfacelevel and the TopOfLayerToBeEvaluated + /// + /// + /// the height of the layer + private double GetLayerHeight(PipingLayer layer) + { + var layerHeight = SoilProfile.GetLayerHeight(layer); + double topLevel; + if (layer.TopLevel > SurfaceLevel) + { + topLevel = SurfaceLevel; + layerHeight = Math.Max(layerHeight - (layer.TopLevel - SurfaceLevel), 0); + } + else + { + topLevel = layer.TopLevel; + } + + double bottomLevel = topLevel - layerHeight; + bottomLevel = bottomLevel > TopOfLayerToBeEvaluated ? bottomLevel : TopOfLayerToBeEvaluated; + return topLevel - bottomLevel; + } + + /// + /// Gets the dry unit weight for the layer. + /// + /// The layer. + /// the dry unit weight + private double GetSoilUnitWeightDry(PipingLayer layer) + { + if (IsUseOvenDryUnitWeight) + { + return layer.DryUnitWeight; + } + return layer.AbovePhreaticLevel; + } + + /// + /// Gets the height of the layer. + /// + /// The layer. + /// The top level. + /// The bottom level. + /// + private double GetLayerHeight(PipingLayer layer, out double topLevel, out double bottomLevel) + { + var layerHeight = SoilProfile.GetLayerHeight(layer); + if (layer.TopLevel > SurfaceLevel) + { + topLevel = SurfaceLevel; + layerHeight = Math.Max(layerHeight - (layer.TopLevel - SurfaceLevel), 0); + } + else + { + topLevel = layer.TopLevel; + } + + bottomLevel = topLevel - layerHeight; + bottomLevel = bottomLevel > TopOfLayerToBeEvaluated ? bottomLevel : TopOfLayerToBeEvaluated; + return topLevel - bottomLevel; + } + + /// + /// Check precondition: Throws the when soil profile is null. + /// + /// The soilprofile is not defined + private void ThrowWhenSoilProfileIsNull() + { + if (SoilProfile == null) + { + throw new SoilVolumicMassCalculatorException(Helper.Translate("SoilprofileNotDefined")); + } + } + + /// + /// Throws the when soil profile has no layers. + /// + /// The soilprofile is not defined + private void ThrowWhenSoilProfileHasNoLayers() + { + if (SoilProfile.Layers.Count == 0) + { + throw new SoilVolumicMassCalculatorException(Helper.Translate("ProfileHasNoLayers")); + } + } + + /// + /// Check precondition: Throws the when surface level is not inside profile. + /// + /// Surfaceline is not inside soil profile + private void ThrowWhenSurfaceLevelNotInsideProfile() + { + if ((SurfaceLevel - SoilProfile.TopLevel) > cTolerance || + (SurfaceLevel - SoilProfile.BottomLevel) < -cTolerance) + { + throw new SoilVolumicMassCalculatorException(Helper.Translate("SurfacelevelNotInsideSoilProfile")); + } + } + + /// + /// Check precondition: Throws the when top layer to be evaluated is not inside profile. + /// + /// The top layer to be evaluated is not inside the soil profile + private void ThrowWhenTopLayerToBeEvaluatedNotInsideProfile() + { + if ((TopOfLayerToBeEvaluated - SoilProfile.TopLevel) > cTolerance || (TopOfLayerToBeEvaluated - SoilProfile.BottomLevel) < -cTolerance) + { + throw new SoilVolumicMassCalculatorException(Helper.Translate("TopLayerNotInsideSoilProfile")); + } + } + +// #region DebugOut +// private bool isDebugOutputEnabled = false; +// +// private void DumpThis(XmlDumper dumper) +// { +// if (isDebugOutputEnabled) +// { +// dumper.DumpObject(this); +// } +// } +// +// /// +// /// Enables debug-dumping to XML +// /// +// public void EnableDebugOutput() +// { +// isDebugOutputEnabled = true; +// } +// #endregion + } +} \ No newline at end of file Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Routines2D.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Routines2D.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/Routines2D.cs (revision 253) @@ -0,0 +1,516 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; + +namespace Deltares.DamPiping +{ + public enum LineIntersection + { + NoIntersection, + Intersects, + Parallel, + } + public static class Routines2D + { +// private const float cCoincideTolerance = 0.0001f; +// private const double cEpsilon = 0.0001; +// + public static double Compute2DDistance(double aX1, double aY1, double aX2, double aY2) + { + double num1 = aX1 - aX2; + double num2 = aY1 - aY2; + return Math.Sqrt(num1 * num1 + num2 * num2); + } + +// public static PointLocation RelativePointPositionToLine(double aLine1X, double aLine1Y, double aLine2X, double aLine2Y, double aPointX, double aPointY) +// { +// double num = Routines2D.CrossProduct(aLine2X - aLine1X, aLine2Y - aLine1Y, aPointX - aLine1X, aPointY - aLine1Y); +// if (Math.Abs(num) < 0.0001) +// return PointLocation.On; +// return num < 0.0 ? PointLocation.Right : PointLocation.Left; +// } +// + public static LineIntersection DetermineIf2DLinesIntersectStrickly(Point2D point1, Point2D point2, Point2D point3, Point2D point4, ref Point2D intersectionPoint) + { + return Routines2D.DetermineIf2DLinesIntersectStrickly(point1, point2, point3, point4, out intersectionPoint, 0.0001); + } + + public static LineIntersection DetermineIf2DLinesIntersectStrickly(Point2D point1, Point2D point2, Point2D point3, Point2D point4, out Point2D intersectionPoint, double tolerance) + { + LineIntersection lineIntersection = Routines2D.DetermineIf2DLinesIntersectWithExtrapolation(Routines2D.CalculateNormalLineConstants(point1, point2), Routines2D.CalculateNormalLineConstants(point3, point4), out intersectionPoint); + switch (lineIntersection) + { + case LineIntersection.Intersects: + if (!Routines2D.DoesPointExistInLine(point1, point2, intersectionPoint, tolerance) || !Routines2D.DoesPointExistInLine(point3, point4, intersectionPoint, tolerance)) + { + intersectionPoint = (Point2D)null; + lineIntersection = LineIntersection.NoIntersection; + break; + } + break; + case LineIntersection.Parallel: + if (!Routines2D.DoLinesAtLeastPartialyOverlap(point1, point2, point3, point4, tolerance)) + { + if (Routines2D.DetermineIfPointsCoincide(point1, point3, tolerance) || Routines2D.DetermineIfPointsCoincide(point1, point4, tolerance)) + { + intersectionPoint = point1; + lineIntersection = LineIntersection.Intersects; + } + if (Routines2D.DetermineIfPointsCoincide(point2, point3, tolerance) || Routines2D.DetermineIfPointsCoincide(point2, point4, tolerance)) + { + intersectionPoint = point2; + lineIntersection = LineIntersection.Intersects; + break; + } + break; + } + break; + } + return lineIntersection; + } + + public static bool DoLinesAtLeastPartialyOverlap(Point2D point1, Point2D point2, Point2D point3, Point2D point4, double tolerance) + { + bool flag = Routines2D.AreLinesParallel(point1, point2, point3, point4, tolerance); + if (!flag) + return false; + if (Math.Abs(point1.X - point2.X) < double.Epsilon) + { + double num1 = Math.Max(point1.Y, point2.Y); + double num2 = Math.Min(point1.Y, point2.Y); + double num3 = Math.Max(point3.Y, point4.Y); + double num4 = Math.Min(point3.Y, point4.Y); + if (num1 <= num4 || num2 >= num3) + flag = false; + } + else + { + double num1 = Math.Max(point1.X, point2.X); + double num2 = Math.Min(point1.X, point2.X); + double num3 = Math.Max(point3.X, point4.X); + double num4 = Math.Min(point3.X, point4.X); + if (num1 <= num4 || num2 >= num3) + flag = false; + } + return flag; + } + + public static bool AreLinesParallel(Point2D point1, Point2D point2, Point2D point3, Point2D point4) + { + return Routines2D.AreLinesParallel(point1, point2, point3, point4, 0.0001); + } + + public static bool AreLinesParallel(Point2D point1, Point2D point2, Point2D point3, Point2D point4, double tolerance) + { + Vector3D vector3D1 = new Vector3D(point2.X, point2.Y, 0.0) - new Vector3D(point1.X, point1.Y, 0.0); + Vector3D vector3D2 = new Vector3D(point4.X, point4.Y, 0.0) - new Vector3D(point3.X, point3.Y, 0.0); + vector3D1.Normalize(); + vector3D2.Normalize(); + return Routines2D.AreLinesParallel(new Point2D(vector3D1.X, vector3D1.Y), new Point2D(vector3D2.X, vector3D2.Y), tolerance); + } + + public static LineIntersection DetermineIf2DLinesIntersectWithExtrapolation(Point2D aPoint1, Point2D aPoint2, Point2D aPoint3, Point2D aPoint4, ref Point2D aIntersectionPoint) + { + return Routines2D.DetermineIf2DLinesIntersectWithExtrapolation(Routines2D.CalculateNormalLineConstants(aPoint1, aPoint2), Routines2D.CalculateNormalLineConstants(aPoint3, aPoint4), out aIntersectionPoint); + } + + public static bool DetermineIfPointsCoincide(Point2D point1, Point2D point2, double tolerance) + { + return Math.Abs(point1.X - point2.X) < tolerance && Math.Abs(point1.Y - point2.Y) < tolerance; + } + + public static bool DoesPointExistInLine(Point2D linePoint1, Point2D linePoint2, Point2D point, double tolerance) + { + double x1 = linePoint1.X; + double y1 = linePoint1.Y; + double x2 = linePoint2.X; + double y2 = linePoint2.Y; + if (linePoint1.X == linePoint2.X && linePoint1.Y == linePoint2.Y) + return Routines2D.Compute2DDistance(point.X, point.Y, x1, y1) < tolerance; + double num1 = Math.Max(x1, x2); + double num2 = Math.Min(x1, x2); + double num3 = Math.Max(y1, y2); + double num4 = Math.Min(y1, y2); + double num5 = num2 - tolerance; + double num6 = num4 - tolerance; + double num7 = num1 + tolerance; + double num8 = num3 + tolerance; + double x3 = point.X; + double y3 = point.Y; + if (x3 <= num5 || x3 >= num7 || (y3 <= num6 || y3 >= num8)) + return false; + double num9 = y1 - y2; + double num10 = x2 - x1; + double num11 = -(num9 * x1 + num10 * y1); + return Math.Abs((num9 * x3 + num10 * y3 + num11) / Routines2D.Compute2DDistance(x1, y1, x2, y2)) < tolerance; + } + +// public static double LinInpolY(double aXStart, double aYStart, double aXEnd, double aYEnd, double aXInterpolation) +// { +// if (Math.Abs(aXEnd - aXStart) < 1E-15) +// return (aYEnd + aYStart) * 0.5; +// return GeneralMathRoutines.LinearInterpolate(aYStart, aYEnd, (aXInterpolation - aXStart) / (aXEnd - aXStart)); +// } +// +// public static double CalculateSquaredDistance(Point2D point, Point2D segmentPoint1, Point2D segmentPoint2) +// { +// Point2D point2D1 = segmentPoint2 - segmentPoint1; +// Point2D point2D2 = point - segmentPoint1; +// Func func = (Func)((first, second) => first.X * second.X + first.Y * second.Y); +// double num1 = func(point2D2, point2D1); +// double num2 = func(point2D1, point2D1); +// Point2D point2D3; +// if (num1 <= 0.0) +// point2D3 = point - segmentPoint1; +// else if (num2 <= num1) +// { +// point2D3 = point - segmentPoint2; +// } +// else +// { +// double num3 = num1 / num2; +// Point2D point2D4 = segmentPoint1 + num3 * point2D1; +// point2D3 = point - point2D4; +// } +// return func(point2D3, point2D3); +// } +// +// public static double CalculateSquaredDistance(Point2D segment1Point1, Point2D segment1Point2, Point2D segment2Point1, Point2D segment2Point2) +// { +// return ((IEnumerable)new double[4] +// { +// Routines2D.CalculateSquaredDistance(segment1Point1, segment2Point1, segment2Point2), +// Routines2D.CalculateSquaredDistance(segment1Point2, segment2Point1, segment2Point2), +// Routines2D.CalculateSquaredDistance(segment2Point1, segment1Point1, segment1Point2), +// Routines2D.CalculateSquaredDistance(segment2Point2, segment1Point1, segment1Point2) +// }).Min(); +// } +// +// public static double CalculateDistanceToLine(double pointX, double pointY, double segmentStartX, double segmentStartY, double segmentEndX, double segmentEndY) +// { +// return Math.Sqrt(Routines2D.CalculateSquaredDistance(new Point2D(pointX, pointY), new Point2D(segmentStartX, segmentStartY), new Point2D(segmentEndX, segmentEndY))); +// } +// +// public static void FindParallelLine(double x1, double y1, double x2, double y2, double distance, out double resultX1, out double resultY1, out double resultX2, out double resultY2) +// { +// if (Math.Abs(distance) > 0.0) +// { +// double num1 = x2 - x1; +// double num2 = y2 - y1; +// double num3 = Routines2D.Compute2DDistance(x1, y1, x2, y2); +// double num4 = distance / num3; +// double num5 = num4 * num1; +// double num6 = num4 * num2; +// x1 -= num6; +// x2 -= num6; +// y1 += num5; +// y2 += num5; +// } +// resultX1 = x1; +// resultX2 = x2; +// resultY1 = y1; +// resultY2 = y2; +// } +// +// public static PointInPolygon CheckIfPointIsInPolygon(List polygon, double x, double y) +// { +// PointInPolygon pointInPolygon = PointInPolygon.OutsidePolygon; +// if (polygon.Count > 0) +// polygon.Add(polygon[0]); +// double count = (double)polygon.Count; +// if (count > 2.0) +// { +// double x1 = polygon[0].X - x; +// double y1 = polygon[0].Y - y; +// if (Math.Abs(x1) < 1E-10 && Math.Abs(y1) < 1E-10) +// return PointInPolygon.OnPolygonEdge; +// double num1 = Math.Atan2(y1, x1); +// double num2 = 0.0; +// for (int index = 1; (double)index < count; ++index) +// { +// double x2 = polygon[index].X - x; +// double y2 = polygon[index].Y - y; +// if (Math.Abs(x2) < 1E-10 && Math.Abs(y2) < 1E-10) +// return PointInPolygon.OnPolygonEdge; +// double num3 = Math.Atan2(y2, x2); +// double num4 = num3 - num1; +// if (num4 < -1.0 * Math.PI) +// num4 += 2.0 * Math.PI; +// if (num4 > Math.PI) +// num4 -= 2.0 * Math.PI; +// if (Math.PI - num4 < 1E-10 || Math.PI + num4 < 1E-10) +// return PointInPolygon.OnPolygonEdge; +// num2 += num4; +// num1 = num3; +// } +// pointInPolygon = num2 > 19.0 * Math.PI / 10.0 || num2 < -19.0 * Math.PI / 10.0 ? PointInPolygon.InsidePolygon : PointInPolygon.OutsidePolygon; +// } +// return pointInPolygon; +// } +// +// public static void GetPointOnLineClosestTo(double aPointX, double aPointY, double aLine1X, double aLine1Y, double aLine2X, double aLine2Y, out double aResultX, out double aResultY) +// { +// double x = Routines2D.Compute2DDistance(aLine1X, aLine1Y, aLine2X, aLine2Y); +// double num = (Math.Pow(Routines2D.Compute2DDistance(aLine1X, aLine1Y, aPointX, aPointY), 2.0) - Math.Pow(Routines2D.Compute2DDistance(aLine2X, aLine2Y, aPointX, aPointY), 2.0) + Math.Pow(x, 2.0)) / (2.0 * x); +// if (num <= 0.0) +// { +// aResultX = aLine1X; +// aResultY = aLine1Y; +// } +// else if (num >= x) +// { +// aResultX = aLine2X; +// aResultY = aLine2Y; +// } +// else +// { +// aResultX = aLine1X + num / x * (aLine2X - aLine1X); +// aResultY = aLine1Y + num / x * (aLine2Y - aLine1Y); +// } +// } +// +// public static List IntersectCircleline(double aX, double aY, double aR, double aX1, double aX2, double aY1, double aY2) +// { +// double num1 = aX2 - aX1; +// double x1 = aX1 - aX; +// double num2 = aY2 - aY1; +// double x2 = aY1 - aY; +// List point2DList = new List(); +// if (Math.Abs(num1) > 1E-08 || Math.Abs(num2) > 1E-08) +// { +// double num3 = num1 * num1 + num2 * num2; +// double num4 = 2.0 * (num1 * x1 + num2 * x2); +// double num5 = x1 * x1 + x2 * x2 - aR * aR; +// double d = num4 * num4 - 4.0 * num3 * num5; +// if (d > 1E-08) +// { +// double num6 = (-num4 + Math.Sqrt(d)) / (2.0 * num3); +// if (num6 >= -1E-08 && num6 <= 1.00000001) +// point2DList.Add(new Point2D(aX1 + num6 * num1, aY1 + num6 * num2)); +// double num7 = (-num4 - Math.Sqrt(d)) / (2.0 * num3); +// if (num7 >= -1E-08 && num7 <= 1.00000001) +// point2DList.Add(new Point2D(aX1 + num7 * num1, aY1 + num7 * num2)); +// } +// else if (Math.Abs(d) <= 1E-08) +// { +// double num6 = -num4 / (2.0 * num3); +// if (num6 >= -1E-08 && num6 <= 1.00000001) +// point2DList.Add(new Point2D(aX1 + num6 * num1, aY1 + num6 * num2)); +// } +// } +// else if (Math.Abs(Math.Pow(x1, 2.0) + Math.Pow(x2, 2.0) - Math.Pow(aR, 2.0)) < 1E-08) +// point2DList.Add(new Point2D(aX1, aY1)); +// return point2DList; +// } +// +// public static Clockwise IsClockWise(IEnumerable aPolygon) +// { +// Point2D[] array = aPolygon.Distinct().ToArray(); +// if (array.Length < 3) +// return Clockwise.NotEnoughUniquePoints; +// double num1 = 0.0; +// for (int index = 0; index < array.Length - 1; ++index) +// { +// double num2 = (array[index + 1].X - array[index].X) * (array[index + 1].Y + array[index].Y); +// num1 += num2; +// } +// double num3 = (array[0].X - array[array.Length - 1].X) * (array[0].Y + array[array.Length - 1].Y); +// int num4 = Math.Sign(num1 + num3); +// if (num4 == 0) +// return Clockwise.PointsOnLine; +// return (double)num4 <= 0.0 ? Clockwise.AntiClockwise : Clockwise.IsClockwise; +// } +// +// public static double ComputeTriangleArea(double aX1, double aY1, double aX2, double aY2, double aX3, double aY3) +// { +// return 0.5 * ((aX3 - aX1) * (aY2 - aY1) - (aY3 - aY1) * (aX2 - aX1)); +// } +// +// public static bool AreEqual(double x1, double x2, double tolerance) +// { +// return Math.Abs(x1 - x2) < tolerance; +// } +// +// public static bool AreEqual(Point2D p1, Point2D p2, int tolerance) +// { +// if (Routines2D.AreEqual(p1.X, p2.X, (double)tolerance)) +// return Routines2D.AreEqual(p1.Y, p2.Y, (double)tolerance); +// return false; +// } +// +// public static bool AreEqual(int x1, int x2, int tolerance) +// { +// return Math.Abs(x1 - x2) < tolerance; +// } +// +// public static bool DetermineIfPointsCoincide(double aX1, double aY1, double aX2, double aY2, double aTolerance) +// { +// if (Math.Abs(aX2 - aX1) < aTolerance) +// return Math.Abs(aY2 - aY1) < aTolerance; +// return false; +// } +// +// public static double Determine2DPolygonArea(List closedPolygon) +// { +// int count = closedPolygon.Count; +// double num = 0.0; +// if (!Routines2D.DetermineIfPointsCoincide(closedPolygon[0], closedPolygon[count - 1], 9.99999974737875E-05)) +// { +// closedPolygon.Add(closedPolygon[0]); +// count = closedPolygon.Count; +// } +// for (int index = 0; index < count - 1; ++index) +// num += Routines2D.CrossProduct(closedPolygon[index].X, closedPolygon[index].Y, closedPolygon[index + 1].X, closedPolygon[index + 1].Y); +// return num / 2.0; +// } +// +// public static Point2D DeterminePolygonCentroid(List closedPolygon) +// { +// int count = closedPolygon.Count; +// if (!Routines2D.DetermineIfPointsCoincide(closedPolygon[0], closedPolygon[count - 1], 9.99999974737875E-05)) +// closedPolygon.Add(closedPolygon[0]); +// else +// --count; +// double num1 = Routines2D.Determine2DPolygonArea(closedPolygon); +// if (num1 > 0.0) +// { +// double num2 = 0.0; +// double num3 = 0.0; +// for (int index = 0; index < count; ++index) +// { +// double num4 = Routines2D.CrossProduct(closedPolygon[index].X, closedPolygon[index].Y, closedPolygon[index + 1].X, closedPolygon[index + 1].Y); +// num2 += (closedPolygon[index].X + closedPolygon[index + 1].X) * num4; +// num3 += (closedPolygon[index].Y + closedPolygon[index + 1].Y) * num4; +// } +// return new Point2D(num2 / (6.0 * num1), num3 / (6.0 * num1)); +// } +// double num5 = 0.0; +// double num6 = 0.0; +// for (int index = 0; index < count; ++index) +// { +// num5 += closedPolygon[index].X; +// num6 += closedPolygon[index].Y; +// } +// return new Point2D(num5 / (double)count, num6 / (double)count); +// } +// +// public static MinMax GetMinMax(List aPointList, ref double aMin, ref double aMax, MinMaxCordinate aMinMaxCoordinate) +// { +// if (aPointList.Count < 2) +// return MinMax.Error; +// aMin = aPointList.Min((Func)(x => +// { +// if (aMinMaxCoordinate != MinMaxCordinate.X) +// return x.Y; +// return x.X; +// })); +// aMax = aPointList.Max((Func)(x => +// { +// if (aMinMaxCoordinate != MinMaxCordinate.X) +// return x.Y; +// return x.X; +// })); +// return MinMax.Calculated; +// } +// +// public static double GetAngle(double x1, double y1, double x2, double y2) +// { +// double aValue1_1 = y2 - y1; +// double aValue1_2 = x2 - x1; +// if (aValue1_2.IsNearEqual(1E-08)) +// return 90.0; +// if (aValue1_1.IsNearEqual(1E-08)) +// return 0.0; +// return Math.Atan(aValue1_1 / aValue1_2) * 180.0 / Math.PI; +// } +// +// public static bool IsPointInRectangularBounds(double aXLeft, double aXRight, double aYTop, double aYBottom, Point2D aPoint) +// { +// bool flag1 = aPoint.X.IsGreaterThanOrEqualTo(aXLeft) && aPoint.X.IsLessThanOrEqualTo(aXRight); +// bool flag2 = aPoint.Y.IsGreaterThanOrEqualTo(aYBottom) && aPoint.Y.IsLessThanOrEqualTo(aYTop); +// if (flag1) +// return flag2; +// return false; +// } +// +// public static double GetMinimumDistanceAmongPoints(List pointList) +// { +// if (pointList.Count < 2) +// return 0.0; +// List doubleList = new List(); +// for (int index1 = 0; index1 < pointList.Count; ++index1) +// { +// for (int index2 = index1 + 1; index2 < pointList.Count; ++index2) +// doubleList.Add(Routines2D.Compute2DDistance(pointList[index1].X, pointList[index1].Y, pointList[index2].X, pointList[index2].Y)); +// } +// doubleList.Sort(); +// return doubleList[0]; +// } +// +// public static Point2D RotateAroundPoint(Point2D rotatePoint, Point2D rotateCenter, double theta) +// { +// double num1 = rotatePoint.X - rotateCenter.X; +// double num2 = rotatePoint.Y - rotateCenter.Y; +// return new Point2D() +// { +// X = Math.Cos(theta) * num1 - Math.Sin(theta) * num2 + rotateCenter.X, +// Y = Math.Sin(theta) * num1 + Math.Cos(theta) * num2 + rotateCenter.Y +// }; +// } +// +// public static bool IsBetween(double x, double x1, double x2) +// { +// if (x1 <= x && x2 >= x) +// return true; +// if (x1 >= x) +// return x2 <= x; +// return false; +// } +// + private static double CrossProduct(double pointAx, double pointAy, double pointBx, double pointBy) + { + return pointAx * pointBy - pointBx * pointAy; + } + + private static Vector3D CalculateNormalLineConstants(Point2D aPoint1, Point2D aPoint2) + { + return new Vector3D() + { + X = aPoint2.Y - aPoint1.Y, + Y = -(aPoint2.X - aPoint1.X), + Z = (aPoint2.Y - aPoint1.Y) * aPoint1.X - (aPoint2.X - aPoint1.X) * aPoint1.Y + }; + } + + private static bool AreLinesParallel(Point2D aLine1Constant, Point2D aLine2Constant, double tolerance) + { + double x1 = aLine1Constant.X; + double y1 = aLine1Constant.Y; + double x2 = aLine2Constant.X; + double y2 = aLine2Constant.Y; + return Math.Abs(Routines2D.CrossProduct(x1, x2, y1, y2)) < tolerance; + } + + private static LineIntersection DetermineIf2DLinesIntersectWithExtrapolation(Vector3D aLine1Constant, Vector3D aLine2Constant, out Point2D aIntersectionPoint) + { + aIntersectionPoint = new Point2D(0.0, 0.0); + double x1 = aLine1Constant.X; + double y1 = aLine1Constant.Y; + double z1 = aLine1Constant.Z; + double x2 = aLine2Constant.X; + double y2 = aLine2Constant.Y; + double z2 = aLine2Constant.Z; + if (Routines2D.AreLinesParallel(new Point2D(x1, x2), new Point2D(y1, y2), 0.0001)) + { + aIntersectionPoint = (Point2D)null; + return LineIntersection.Parallel; + } + double num1 = (y2 * z1 - y1 * z2) / (x1 * y2 - x2 * y1); + double num2 = (z1 * x2 - z2 * x1) / (x2 * y1 - x1 * y2); + aIntersectionPoint.X = num1; + aIntersectionPoint.Y = num2; + return LineIntersection.Intersects; + } + } +} Index: dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingLayer.cs =================================================================== diff -u --- dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingLayer.cs (revision 0) +++ dam failuremechanisms/damPiping/trunk/src/Deltares.DamPiping/PipingLayer.cs (revision 253) @@ -0,0 +1,71 @@ +namespace Deltares.DamPiping +{ + /// + /// Class to hold layers (1D) for the Piping kernel + /// + public class PipingLayer + { + /// + /// Assigns the specified layer. + /// + /// The layer. + public void Assign(PipingLayer layer) + { + AbovePhreaticLevel = layer.AbovePhreaticLevel; + BelowPhreaticLevel = layer.BelowPhreaticLevel; + DryUnitWeight = layer.DryUnitWeight; + TopLevel = layer.TopLevel; + IsAquifer = layer.IsAquifer; + Name = layer.Name; + } + + /// + /// Gets or sets the unit weight above phreatic level. + /// + /// + /// The unit weight above phreatic level. + /// + public double AbovePhreaticLevel { get; set; } + + /// + /// Gets or sets the unit weight below phreatic level. + /// + /// + /// The unit weight below phreatic level. + /// + public double BelowPhreaticLevel { get; set; } + + /// + /// Gets or sets the dry unit weight. + /// + /// + /// The dry unit weight. + /// + public double DryUnitWeight { get; set; } + + /// + /// Gets or sets the top level. + /// + /// + /// The top level. + /// + public double TopLevel { get; set; } + + /// + /// Gets or sets a value indicating whether this instance is aquifer. + /// + /// + /// true if this instance is aquifer; otherwise, false. + /// + public bool IsAquifer { get; set; } + + /// + /// Gets or sets the name. + /// + /// + /// The name. + /// + public string Name { get; set; } + + } +}