Index: DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs =================================================================== diff -u -r5284 -r5290 --- DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs (.../PlLinesToWaternetConverter.cs) (revision 5284) +++ DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs (.../PlLinesToWaternetConverter.cs) (revision 5290) @@ -20,6 +20,7 @@ // All rights reserved. using System; +using System.Collections; using System.Collections.Generic; using System.Data; using System.Linq; @@ -29,7 +30,6 @@ using Deltares.DamEngine.Data.Geometry; using Deltares.DamEngine.Data.Geotechnics; using Deltares.DamEngine.Data.Standard; -using Deltares.MacroStability.CSharpWrapper.Input; using CharacteristicPointType = Deltares.DamEngine.Data.Geotechnics.CharacteristicPointType; using SoilProfile = Deltares.DamEngine.Data.Geotechnics.SoilProfile; @@ -39,13 +39,18 @@ { public enum LayerType { - BottomAquifer, - TopLayerInBetweenAquiferCluster, - BottomLayerInBetweenAquiferCluster, - HighestAquifer, + BottomAquiferCluster, + InBetweenAquiferCluster, + HighestAquiferCluster, LowestLayer } + private enum BoundaryType + { + Top, + Bottom + } + private const double deviationX = 1e-06; private const double toleranceAlmostEqual = 1e-09; @@ -147,8 +152,8 @@ ThrowWhenPlLinesIsNull(plLines); ThrowWhenSoilProfileIsNull(soilProfile); - Point2D[] bottomAquiferCoordinates = DetermineLayerBoundaryPoints(LayerType.BottomAquifer, soilProfile).ToArray(); - if (!bottomAquiferCoordinates.Any()) + DetermineLayerBoundaryPoints(LayerType.BottomAquiferCluster, soilProfile, out Point2D[] bottomAquiferCoordinates, out _); + if (bottomAquiferCoordinates == null) { throw new NoNullAllowedException(string.Format(Resources.NoBottomAquiferLayer, soilProfile.Name)); } @@ -270,9 +275,8 @@ waternet.HeadLineList.Add(headLine); for (var i = 0; i < inBetweenAquiferCount; i++) { - Point2D[] inBetweenAquiferUpperCoordinates = DetermineLayerBoundaryPoints(LayerType.TopLayerInBetweenAquiferCluster, soilProfile, i).ToArray(); - Point2D[] inBetweenAquiferLowerCoordinates = DetermineLayerBoundaryPoints(LayerType.BottomLayerInBetweenAquiferCluster, soilProfile, i).ToArray(); - + DetermineLayerBoundaryPoints(LayerType.InBetweenAquiferCluster, soilProfile, out Point2D[] inBetweenAquiferUpperCoordinates, out Point2D[] inBetweenAquiferLowerCoordinates, i); + if (inBetweenAquiferUpperCoordinates.Any() && inBetweenAquiferLowerCoordinates.Any()) { WaternetLine waternetLineTop = CreateWaternetLine(inBetweenAquiferUpperCoordinates); @@ -348,20 +352,21 @@ /// /// Determines all the points of the layer boundary for the specified layer type: - /// - For LayerType = BottomAquifer, it will return the top boundary of the bottom aquifer (= waternet line for PL3) - /// - For TopLayerInBetweenAquiferCluster, it will return the top boundary of the in-between aquifer (= waternet line for PL4) + /// - For LayerType = BottomAquifer, layerBoundaryTop is the top boundary of the bottom aquifer (= waternet line for PL3) + /// - For LayerType = InBetweenAquiferCluster, it will return the top boundary of the in-between aquifer (= waternet line for PL4) /// - For BottomLayerInBetweenAquiferCluster, it will return the bottom boundary of the in-between aquifer (= waternet line for PL4) - /// - For HighestAquifer, it will return the top boundary of the highest aquifer (= waternet line for PL1 for "Hydrostatic" waternet) - /// - For LowestLayer, it will return the bottom boundary of the lowest layer (= waternet line for PL1 for "Full hydrostatic" waternet) + /// - For LayerType = HighestAquifer, it will return the top boundary of the highest aquifer (= waternet line for PL1 for "Hydrostatic" waternet) + /// - For LayerType = LowestLayer, it will return the bottom boundary of the lowest layer (= waternet line for PL1 for "Full hydrostatic" waternet) /// /// The type of layer boundary. /// The soil profile 2D. + /// All the points of the layer top boundary. + /// All the points of the layer bottom boundary. /// In case of several in-between aquifers, the index of the in-between aquifer must be specified. /// All the points of the layer boundary. - internal static IEnumerable DetermineLayerBoundaryPoints(LayerType layerType, SoilProfile2D soilProfile, int indexInBetweenAquifer = -1) + internal static void DetermineLayerBoundaryPoints(LayerType layerType, SoilProfile2D soilProfile, out Point2D[] layerBoundaryTop, out Point2D[] layerBoundaryBottom, int indexInBetweenAquifer = -1) { double[] xCoordinates = DetermineAllXCoordinatesOfSoilProfile(soilProfile); - SoilLayer1D previousAquiferLayer = null; var xCoordinatesAll = new List(); foreach (double xCoordinate in xCoordinates) { @@ -376,55 +381,79 @@ } } - var layerBoundaryPoints = new List(); + var layerBoundaryTopPoints = new List(); + var layerBoundaryBottomPoints = new List(); + double previousAquiferTop = double.NaN; + double previousAquiferBottom = double.NaN; foreach (double xCoordinate in xCoordinatesAll) { SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(xCoordinate); - // Determine if the layer is in range of the previous layer - // If not, return empty coordinates, because the layer is interrupted - SoilLayer1D currentAquifer = GetSoilLayer1D(layerType, crossSection, indexInBetweenAquifer); - if (previousAquiferLayer != null && currentAquifer != null - && !AreHorizontallyConnected(previousAquiferLayer, currentAquifer)) + // Determine if the cluster of layers is in range of the previous cluster of layers. + // If not, return empty coordinates, because the cluster of layer is interrupted. + double currentAquiferTop = GetLevel(layerType, BoundaryType.Top, crossSection, indexInBetweenAquifer); + double currentAquiferBottom = GetLevel(layerType, BoundaryType.Bottom, crossSection, indexInBetweenAquifer); + + if (!double.IsNaN(previousAquiferTop) && !double.IsNaN(previousAquiferBottom) + && !double.IsNaN(currentAquiferTop) + && !double.IsNaN(currentAquiferBottom) + && !AreHorizontallyConnected(previousAquiferTop, previousAquiferBottom, currentAquiferTop, currentAquiferBottom)) { - return Enumerable.Empty(); + layerBoundaryTop = null; + layerBoundaryBottom = null; + return; } - - if (currentAquifer != null) + + if (!double.IsNaN(currentAquiferTop) && !double.IsNaN(currentAquiferBottom)) { - previousAquiferLayer = currentAquifer; - layerBoundaryPoints.Add(layerType is LayerType.BottomLayerInBetweenAquiferCluster or LayerType.LowestLayer ? new Point2D(xCoordinate, currentAquifer.BottomLevel) : new Point2D(xCoordinate, currentAquifer.TopLevel)); + previousAquiferTop = currentAquiferTop; + previousAquiferBottom = currentAquiferBottom; + layerBoundaryTopPoints.Add(new Point2D(xCoordinate, currentAquiferTop)); + layerBoundaryBottomPoints.Add(new Point2D(xCoordinate, currentAquiferBottom)); } } // Perform a short validation that the coordinates are fully defined from the beginning to the end // of the profile. If not, the layer is not continuous. - if (!IsLayerBoundaryContinuous(layerBoundaryPoints, xCoordinates.First(), xCoordinates.Last())) + if (!IsLayerBoundaryContinuous(layerBoundaryTopPoints, xCoordinates.First(), xCoordinates.Last()) + || !IsLayerBoundaryContinuous(layerBoundaryBottomPoints, xCoordinates.First(), xCoordinates.Last())) { - return Enumerable.Empty(); + layerBoundaryTop = null; + layerBoundaryBottom = null; + return; } - - return layerBoundaryPoints; + layerBoundaryTop = layerBoundaryTopPoints.ToArray(); + layerBoundaryBottom = layerBoundaryBottomPoints.ToArray(); } - private static bool IsLayerBoundaryContinuous(List boundaryPoints, double leftGeometryBoundary, double rightGeometryBoundary) + private static bool IsLayerBoundaryContinuous(IEnumerable boundaryPoints, double leftGeometryBoundary, double rightGeometryBoundary) { - return boundaryPoints.Any() && (Math.Abs(boundaryPoints.First().X - leftGeometryBoundary) < toleranceAlmostEqual) - && Math.Abs(boundaryPoints.Last().X - rightGeometryBoundary) < toleranceAlmostEqual; + IEnumerable point2Ds = boundaryPoints.ToList(); + return point2Ds.Any() && (Math.Abs(point2Ds.First().X - leftGeometryBoundary) < toleranceAlmostEqual) + && Math.Abs(point2Ds.Last().X - rightGeometryBoundary) < toleranceAlmostEqual; } - private static SoilLayer1D GetSoilLayer1D(LayerType layerType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer) + private static double GetLevel(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer) { + SoilLayer1D layer = GetLayer(layerType, boundaryType, soilProfile1D, indexInBetweenAquifer); + if (layer != null) + { + return boundaryType == BoundaryType.Top ? layer.TopLevel : layer.BottomLevel; + } + + return double.NaN; + } + + private static SoilLayer1D GetLayer(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer) + { switch (layerType) { - case LayerType.BottomAquifer: - return soilProfile1D.BottomAquiferLayer; - case LayerType.TopLayerInBetweenAquiferCluster: - return soilProfile1D.GetInBetweenAquiferClusters[indexInBetweenAquifer].Item1; - case LayerType.BottomLayerInBetweenAquiferCluster: - return soilProfile1D.GetInBetweenAquiferClusters[indexInBetweenAquifer].Item2; - case LayerType.HighestAquifer: - return soilProfile1D.GetHighestAquifer(); + case LayerType.BottomAquiferCluster: + return boundaryType == BoundaryType.Top ? soilProfile1D.BottomAquiferLayer : soilProfile1D.DeepestAquiferLayer; + case LayerType.InBetweenAquiferCluster: + return boundaryType == BoundaryType.Top ? soilProfile1D.GetInBetweenAquiferClusters[indexInBetweenAquifer].Item1 : soilProfile1D.GetInBetweenAquiferClusters[indexInBetweenAquifer].Item2; + case LayerType.HighestAquiferCluster: + return boundaryType == BoundaryType.Top ? soilProfile1D.GetHighestAquifer() : soilProfile1D.GetLowestLayerOfHighestAquiferCluster(); case LayerType.LowestLayer: return soilProfile1D.Layers.Last(); } @@ -511,7 +540,7 @@ } case IntrusionVerticalWaterPressureType.HydroStatic: { - Point2D[] highestAquiferCoordinates = DetermineLayerBoundaryPoints(LayerType.HighestAquifer, soilProfile2D).ToArray(); + DetermineLayerBoundaryPoints(LayerType.HighestAquiferCluster, soilProfile2D, out Point2D[] highestAquiferCoordinates, out _); return CreateWaternetLine(highestAquiferCoordinates); } case IntrusionVerticalWaterPressureType.Linear: @@ -521,7 +550,7 @@ } case IntrusionVerticalWaterPressureType.FullHydroStatic: { - Point2D[] lowestBoundaryCoordinates = DetermineLayerBoundaryPoints(LayerType.LowestLayer, soilProfile2D).ToArray(); + DetermineLayerBoundaryPoints(LayerType.LowestLayer, soilProfile2D, out _, out Point2D[] lowestBoundaryCoordinates); return CreateWaternetLine(lowestBoundaryCoordinates); } case null: @@ -613,48 +642,39 @@ { SoilProfile1D crossSectionLeft = soilProfile.GetSoilProfile1D(xCoordinate - deviationX); SoilProfile1D crossSectionRight = soilProfile.GetSoilProfile1D(xCoordinate + deviationX); - if (crossSectionLeft != null && crossSectionRight != null) + if (crossSectionLeft == null || crossSectionRight == null) { - SoilLayer1D currentAquiferLeft = GetSoilLayer1D(layerType, crossSectionLeft, indexInBetweenAquifer); - SoilLayer1D currentAquiferRight = GetSoilLayer1D(layerType, crossSectionRight, indexInBetweenAquifer); - if (currentAquiferLeft != null && currentAquiferRight != null) - { - if (layerType is LayerType.BottomAquifer or LayerType.TopLayerInBetweenAquiferCluster) - { - return Math.Abs(currentAquiferLeft.TopLevel - currentAquiferRight.TopLevel) > GeometryConstants.Accuracy; - } + return false; + } - return Math.Abs(currentAquiferLeft.BottomLevel - currentAquiferRight.BottomLevel) > GeometryConstants.Accuracy; - } + double aquiferLeftTop = GetLevel(layerType, BoundaryType.Top, crossSectionLeft, indexInBetweenAquifer); + double aquiferLefBottom = GetLevel(layerType, BoundaryType.Bottom, crossSectionLeft, indexInBetweenAquifer); + double aquiferRightTop = GetLevel(layerType, BoundaryType.Top, crossSectionRight, indexInBetweenAquifer); + double aquiferRightBottom = GetLevel(layerType, BoundaryType.Bottom, crossSectionRight, indexInBetweenAquifer); + if (!double.IsNaN(aquiferLeftTop) && !double.IsNaN(aquiferLefBottom) && !double.IsNaN(aquiferRightTop) && !double.IsNaN(aquiferRightBottom)) + { + return Math.Abs(aquiferLeftTop - aquiferRightTop) > GeometryConstants.Accuracy || Math.Abs(aquiferLefBottom - aquiferRightBottom) > GeometryConstants.Accuracy; } return false; } - private static bool AreHorizontallyConnected(SoilLayer1D leftSoilLayer, SoilLayer1D rightSoilLayer) + private static bool AreHorizontallyConnected(double leftLayerTop, double leftLayerBottom, double rightLayerTop, double rightLayerBottom) { - // The layers are identical - if (ReferenceEquals(leftSoilLayer, rightSoilLayer)) - { - return true; - } - // Left soil layer envelopes whole right soil layer - if (leftSoilLayer.BottomLevel <= rightSoilLayer.BottomLevel - && leftSoilLayer.TopLevel >= rightSoilLayer.TopLevel) + if (leftLayerBottom <= rightLayerBottom && leftLayerTop >= rightLayerTop) { return true; } // Right soil layer envelopes whole left soil layer - if (rightSoilLayer.BottomLevel <= leftSoilLayer.BottomLevel - && rightSoilLayer.TopLevel >= leftSoilLayer.TopLevel) + if (rightLayerBottom <= leftLayerBottom && rightLayerTop >= leftLayerTop) { return true; } - return (rightSoilLayer.TopLevel <= leftSoilLayer.TopLevel && rightSoilLayer.TopLevel >= leftSoilLayer.BottomLevel) // Top level lies inbetween the left soil layer - || (rightSoilLayer.BottomLevel >= leftSoilLayer.BottomLevel && rightSoilLayer.BottomLevel <= leftSoilLayer.TopLevel); // Bottom level lies inbetween the left soil layer + return (rightLayerTop <= leftLayerTop && rightLayerTop >= leftLayerBottom) // Top level lies in-between the left soil layer + || (rightLayerBottom >= leftLayerBottom && rightLayerBottom <= leftLayerTop); // Bottom level lies in-between the left soil layer } private static bool IsBelowSoilProfile(SoilProfile1D soilProfile, PlLine line)