Index: DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs =================================================================== diff -u -r5567 -r5923 --- DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs (.../PlLinesToWaternetConverter.cs) (revision 5567) +++ DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs (.../PlLinesToWaternetConverter.cs) (revision 5923) @@ -235,34 +235,31 @@ return; } - if (!SoilProfile2DHelper.AreInBetweenAquiferClustersPresent(soilProfile, out int inBetweenAquiferCount)) + SoilProfile2DHelper.DetermineInBetweenAquifersLayerBoundaryPoints(soilProfile, out List inBetweenAquiferUpperCoordinates, + out List inBetweenAquiferLowerCoordinates); + + if (!inBetweenAquiferUpperCoordinates.Any() || !inBetweenAquiferLowerCoordinates.Any()) { return; } var headLine = CreateLine(plLine4, headLine4Name); waternet.HeadLineList.Add(headLine); - for (var i = 0; i < inBetweenAquiferCount; i++) + for (var i = 0; i < inBetweenAquiferUpperCoordinates.Count; i++) { - SoilProfile2DHelper.DetermineAquiferLayerBoundaryPoints(SoilProfile2DHelper.LayerType.InBetweenAquiferCluster, soilProfile, out Point2D[] inBetweenAquiferUpperCoordinates, out Point2D[] inBetweenAquiferLowerCoordinates, i); - - if (inBetweenAquiferUpperCoordinates.Any() && inBetweenAquiferLowerCoordinates.Any()) - { - WaternetLine waternetLineTop = CreateWaternetLine(inBetweenAquiferUpperCoordinates); - waternetLineTop.HeadLine = headLine; - waternetLineTop.Name = waternetLine4TopName + " (" + (i + 1) + ")"; - waternet.WaternetLineList.Add(waternetLineTop); - - WaternetLine waternetLineBottom = CreateWaternetLine(inBetweenAquiferLowerCoordinates); - waternetLineBottom.HeadLine = headLine; - waternetLineBottom.Name = waternetLine4BottomName + " (" + (i + 1) + ")"; - waternet.WaternetLineList.Add(waternetLineBottom); - } - else - { - throw new NoNullAllowedException(string.Format(Resources.PlLinesToWaternetConverter_ErrorDuringDeterminationCoordinatesInBetweenAquifer, soilProfile.Name)); - } + WaternetLine waternetLineTop = CreateWaternetLine(inBetweenAquiferUpperCoordinates[i].CalcPoints); + waternetLineTop.HeadLine = headLine; + waternetLineTop.Name = waternetLine4TopName + " (" + (i + 1) + ")"; + waternet.WaternetLineList.Add(waternetLineTop); } + + for (var i = 0; i < inBetweenAquiferLowerCoordinates.Count; i++) + { + WaternetLine waternetLineBottom = CreateWaternetLine(inBetweenAquiferLowerCoordinates[i].CalcPoints); + waternetLineBottom.HeadLine = headLine; + waternetLineBottom.Name = waternetLine4BottomName + " (" + (i + 1) + ")"; + waternet.WaternetLineList.Add(waternetLineBottom); + } } private static TLineType CreateLine(PlLine plLine, string name) Index: DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/GeometrySurface.cs =================================================================== diff -u -r5634 -r5923 --- DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/GeometrySurface.cs (.../GeometrySurface.cs) (revision 5634) +++ DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/GeometrySurface.cs (.../GeometrySurface.cs) (revision 5923) @@ -183,7 +183,8 @@ int currentIndex = OuterLoop.CalcPoints.IndexOf(startTopPoint); FillTopOrBottomPoints(topPoints, endTopPoint, currentIndex); - + topPoints.SyncPoints(); + return topPoints; } @@ -209,9 +210,10 @@ FillTopOrBottomPoints(bottomPoints, startBottomPoint, currentIndex); - // As the curves are sorted clockwise, the bottomcurves are always orientated from right to left + // As the curves are sorted clockwise, the bottom curves are always orientated from right to left // while this result of this method must be from left to right so reverse the curves bottomPoints.CalcPoints.Reverse(); + bottomPoints.SyncPoints(); return bottomPoints; } Index: DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/Common/SoilProfile2DHelper.cs =================================================================== diff -u -r5919 -r5923 --- DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/Common/SoilProfile2DHelper.cs (.../SoilProfile2DHelper.cs) (revision 5919) +++ DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/Common/SoilProfile2DHelper.cs (.../SoilProfile2DHelper.cs) (revision 5923) @@ -32,7 +32,6 @@ public enum LayerType { BottomAquiferCluster, - InBetweenAquiferCluster, HighestAquiferCluster, LowestLayer } @@ -45,29 +44,7 @@ private const double deviationX = 1e-06; private const double toleranceAlmostEqual = 1e-09; - - internal static bool AreInBetweenAquiferClustersPresent(SoilProfile2D soilProfile, out int count) - { - double[] xCoordinates = DetermineAllXCoordinatesOfSoilProfile(soilProfile); - var currentCount = 0; - var previousCount = 0; - for (var i = 0; i < xCoordinates.Length; i++) - { - SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(xCoordinates[i]); - currentCount = crossSection.GetInBetweenAquiferClusters?.Count ?? 0; - if (i > 0 && currentCount != previousCount) - { - count = 0; - return false; - } - previousCount = currentCount; - } - - count = currentCount; - return count > 0; - } - /// /// Check if at least one isolated in-between aquifer is present. /// @@ -79,7 +56,7 @@ { return true; } - + // The inner loops are removed otherwise in case of "Aquitard inner loop in aquifer" the method // IsAtLeastOneIsolatedInBetweenAquiferPresent will return True because the number of aquifers in the first crossSection // cutting the inner loop is higher than in the previousCrossSection outside the inner loop. @@ -91,7 +68,7 @@ { SoilProfile1D crossSection = soilProfileWithoutInnerLoops.GetSoilProfile1D(xCoordinates[i]); int currentCount = crossSection.GetInBetweenAquiferClusters?.Count ?? 0; - + if (i > 0 && (currentCount != previousCount) && IsAtLeastOneIsolatedInBetweenAquiferPresent(currentCount, previousCount, crossSection, previousCrossSection)) { return true; @@ -171,24 +148,21 @@ /// /// Determines all the points of the layer boundary for the specified layer type: /// - For LayerType = BottomAquiferCluster, layerBoundaryTop is the top boundary of the bottom aquifer (= waternet line for PL3) - /// - For LayerType = InBetweenAquiferCluster, layerBoundaryTop and layerBoundaryBottom are both boundaries of the in-between aquifer (= waternet line for PL4) /// - For LayerType = HighestAquiferCluster, layerBoundaryTop is the top of the highest aquifer (= waternet line for PL1 for "Hydrostatic" waternet) /// - For LayerType = LowestLayer, layerBoundaryBottom is the bottom boundary of the lowest layer (= waternet line for PL1 for "Full hydrostatic" waternet) - /// When the aquifer is not continuous, null top and bottom boundaries are returned. /// /// 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 void DetermineAquiferLayerBoundaryPoints(LayerType layerType, SoilProfile2D soilProfile, out Point2D[] layerBoundaryTop, out Point2D[] layerBoundaryBottom, int indexInBetweenAquifer = -1) + internal static void DetermineAquiferLayerBoundaryPoints(LayerType layerType, SoilProfile2D soilProfile, out Point2D[] layerBoundaryTop, out Point2D[] layerBoundaryBottom) { double[] xCoordinates = DetermineAllXCoordinatesOfSoilProfile(soilProfile); var xCoordinatesAll = new List(); foreach (double xCoordinate in xCoordinates) { - if (HasAquiferAVerticalPartAtGivenX(layerType, xCoordinate, soilProfile, indexInBetweenAquifer)) + if (HasAquiferAVerticalPartAtGivenX(layerType, xCoordinate, soilProfile)) { xCoordinatesAll.Add(xCoordinate - deviationX); xCoordinatesAll.Add(xCoordinate + deviationX); @@ -201,36 +175,20 @@ var layerBoundaryTopPoints = new List(); var layerBoundaryBottomPoints = new List(); - double previousAquiferTop = double.NaN; - double previousAquiferBottom = double.NaN; - foreach (double xCoordinate in xCoordinatesAll) + foreach (double t in xCoordinatesAll) { - SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(xCoordinate); + SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(t); + double currentAquiferTop = GetLevel(layerType, BoundaryType.Top, crossSection); + double currentAquiferBottom = GetLevel(layerType, BoundaryType.Bottom, crossSection); - // Determine if the cluster of layers is in range of the previous cluster of layers. - // If not, return empty coordinates, because the cluster of layers 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)) - { - layerBoundaryTop = null; - layerBoundaryBottom = null; - return; - } - if (!double.IsNaN(currentAquiferTop) && !double.IsNaN(currentAquiferBottom)) { - previousAquiferTop = currentAquiferTop; - previousAquiferBottom = currentAquiferBottom; - layerBoundaryTopPoints.Add(new Point2D(xCoordinate, currentAquiferTop)); - layerBoundaryBottomPoints.Add(new Point2D(xCoordinate, currentAquiferBottom)); + layerBoundaryTopPoints.Add(new Point2D(t, currentAquiferTop)); + layerBoundaryBottomPoints.Add(new Point2D(t, 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(layerBoundaryTopPoints, xCoordinates.First(), xCoordinates.Last()) @@ -245,7 +203,115 @@ layerBoundaryBottom = layerBoundaryBottomPoints.ToArray(); } + private static void AddAquiferContainingPointIfNotYetAdded(SoilProfile2D soilProfile, double xCoord, double zCoord, ref List relevantAquifers) + { + foreach (GeometrySurface surface in soilProfile.Geometry.Surfaces) + { + var polygon = new GeometryLoop(); + foreach (Point2D point in surface.OuterLoop.CalcPoints) + { + polygon.CalcPoints.Add(point); + } + + if (Routines2D.CheckIfPointIsInPolygon(polygon, xCoord, zCoord) == PointInPolygon.OutsidePolygon) + { + continue; + } + + if (!relevantAquifers.Contains(surface)) + { + relevantAquifers.Add(surface); + } + + break; + } + } + /// + /// Determines all the points of the layer boundaries (top and bottom) of all the in-between aquifers. + /// + /// The soil profile 2D. + /// All the points of the layer top boundary of all the in-between aquifers. + /// All the points of the layer bottom boundary of all the in-between aquifers. + /// All the points of the layer boundary of all the in-between aquifers. + internal static void DetermineInBetweenAquifersLayerBoundaryPoints(SoilProfile2D soilProfile, out List aquifersBoundaryTop, out List aquifersBoundaryBottom) + { + // The inner loops are removed because they can't be in-between aquifers + SoilProfile2D soilProfileWithoutInnerLoops = CreateSoilProfileWithoutInnerLoops(soilProfile); + double[] xCoordinates = DetermineAllXCoordinatesOfSoilProfile(soilProfileWithoutInnerLoops); + + var relevantAquifers = new List(); + foreach (double xCoordinate in xCoordinates) + { + SoilProfile1D crossSection = soilProfileWithoutInnerLoops.GetSoilProfile1D(xCoordinate); + if (crossSection.GetInBetweenAquiferClusters.Count > 0) + { + for (var i = 0; i < crossSection.GetInBetweenAquiferClusters.Count; i++) + { + double zLevel = 0.5 * (crossSection.GetInBetweenAquiferClusters[i].Item1.TopLevel + crossSection.GetInBetweenAquiferClusters[i].Item2.BottomLevel); + AddAquiferContainingPointIfNotYetAdded(soilProfileWithoutInnerLoops, xCoordinate, zLevel, ref relevantAquifers); + } + } + } + + List isolatedLayerBoundaryTop = []; + List isolatedLayerBoundaryBottom = []; + aquifersBoundaryTop = []; + aquifersBoundaryBottom = []; + if (relevantAquifers != null) + { + foreach (GeometrySurface surface in relevantAquifers) + { + isolatedLayerBoundaryTop.Add(surface.DetermineTopGeometrySurface()); + isolatedLayerBoundaryBottom.Add(surface.DetermineBottomGeometrySurface()); + } + + aquifersBoundaryTop = ConnectPolylines(isolatedLayerBoundaryTop); + aquifersBoundaryBottom = ConnectPolylines(isolatedLayerBoundaryBottom); + } + } + + private static List ConnectPolylines(List polylines) + { + var connectedPolylines = new List(); + var usedPolylines = new HashSet(); + + foreach (GeometryPointString polyline in polylines) + { + if (usedPolylines.Contains(polyline)) + continue; + + var currentPolyline = new GeometryPointString(); + currentPolyline.CalcPoints.AddRange(polyline.CalcPoints); + usedPolylines.Add(polyline); + + bool foundConnection; + do + { + foundConnection = false; + foreach (GeometryPointString nextPolyline in polylines) + { + if (usedPolylines.Contains(nextPolyline)) + continue; + + if (currentPolyline.CalcPoints.Last().LocationEquals(nextPolyline.CalcPoints.First())) + { + currentPolyline.CalcPoints.AddRange(nextPolyline.CalcPoints.Skip(1)); + usedPolylines.Add(nextPolyline); + foundConnection = true; + break; + } + } + } while (foundConnection); + + currentPolyline.SyncPoints(); + connectedPolylines.Add(currentPolyline); + } + + return connectedPolylines; + } + + /// /// Determine all the xCoordinates to make cross-sections. /// /// The soil profile 2D. @@ -256,17 +322,17 @@ double[] xCoordinates = points.Select(point => point.X).OrderBy(x => x).Distinct().ToArray(); return xCoordinates; } - + private static bool IsLayerBoundaryContinuous(IEnumerable boundaryPoints, double leftGeometryBoundary, double rightGeometryBoundary) { IEnumerable point2Ds = boundaryPoints.ToList(); return point2Ds.Any() && (Math.Abs(point2Ds.First().X - leftGeometryBoundary) < toleranceAlmostEqual) && Math.Abs(point2Ds.Last().X - rightGeometryBoundary) < toleranceAlmostEqual; } - - private static double GetLevel(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer) + + private static double GetLevel(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D) { - SoilLayer1D layer = GetLayer(layerType, boundaryType, soilProfile1D, indexInBetweenAquifer); + SoilLayer1D layer = GetLayer(layerType, boundaryType, soilProfile1D); if (layer != null) { return boundaryType == BoundaryType.Top ? layer.TopLevel : layer.BottomLevel; @@ -275,14 +341,12 @@ return double.NaN; } - private static SoilLayer1D GetLayer(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer) + private static SoilLayer1D GetLayer(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D) { switch (layerType) { 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: @@ -292,7 +356,7 @@ return null; } - private static bool HasAquiferAVerticalPartAtGivenX(LayerType layerType, double xCoordinate, SoilProfile2D soilProfile, int indexInBetweenAquifer) + private static bool HasAquiferAVerticalPartAtGivenX(LayerType layerType, double xCoordinate, SoilProfile2D soilProfile) { SoilProfile1D crossSectionLeft = soilProfile.GetSoilProfile1D(xCoordinate - deviationX); SoilProfile1D crossSectionRight = soilProfile.GetSoilProfile1D(xCoordinate + deviationX); @@ -301,33 +365,15 @@ return false; } - 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); + double aquiferLeftTop = GetLevel(layerType, BoundaryType.Top, crossSectionLeft); + double aquiferLefBottom = GetLevel(layerType, BoundaryType.Bottom, crossSectionLeft); + double aquiferRightTop = GetLevel(layerType, BoundaryType.Top, crossSectionRight); + double aquiferRightBottom = GetLevel(layerType, BoundaryType.Bottom, crossSectionRight); 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(double leftLayerTop, double leftLayerBottom, double rightLayerTop, double rightLayerBottom) - { - // Left soil layer envelopes whole right soil layer - if (leftLayerBottom <= rightLayerBottom && leftLayerTop >= rightLayerTop) - { - return true; - } - - // Right soil layer envelopes whole left soil layer - if (rightLayerBottom <= leftLayerBottom && rightLayerTop >= leftLayerTop) - { - return true; - } - - 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 - } } \ No newline at end of file Index: DamEngine/trunk/src/Deltares.DamEngine.Calculators.Tests/KernelWrappers/Common/SoilProfile2DHelperTests.cs =================================================================== diff -u -r5382 -r5923 --- DamEngine/trunk/src/Deltares.DamEngine.Calculators.Tests/KernelWrappers/Common/SoilProfile2DHelperTests.cs (.../SoilProfile2DHelperTests.cs) (revision 5382) +++ DamEngine/trunk/src/Deltares.DamEngine.Calculators.Tests/KernelWrappers/Common/SoilProfile2DHelperTests.cs (.../SoilProfile2DHelperTests.cs) (revision 5923) @@ -33,12 +33,11 @@ public class SoilProfile2DHelperTests { private const double precision5Decimals = 0.0000051; - - private static readonly Random random = new(21); - private static readonly double leftCoordinate = random.NextDouble(); - private static readonly double rightCoordinate = leftCoordinate + 1; - private static readonly double middleXCoordinate = 0.5 * (leftCoordinate + rightCoordinate); + private const double leftCoordinate = -50; + private const double rightCoordinate = 50; + private const double middleXCoordinate = 0.5 * (leftCoordinate + rightCoordinate); + /// /// Different soil profiles 2D are tested, see the pictures in method GetSoilProfilesWithContinuousBottomAquiferLayer, /// with a top clay layer and bottom aquifer layers always forming a continuous bottom layer. @@ -185,22 +184,25 @@ soilProfile.Surfaces.Add(soilLayerAquiferInBetweenRight); soilProfile.Surfaces.Add(soilLayerInBetween); soilProfile.Surfaces.Add(soilLayerAquiferBottom); - soilProfile.Geometry.Surfaces.Add(new GeometrySurface()); + soilProfile.Geometry.Surfaces.Add(soilLayer.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferInBetweenLeft.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferInBetweenRight.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerInBetween.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferBottom.GeometrySurface); // Call - bool isInBetweenLayerPresent = SoilProfile2DHelper.AreInBetweenAquiferClustersPresent(soilProfile, out int count); - SoilProfile2DHelper.DetermineAquiferLayerBoundaryPoints(SoilProfile2DHelper.LayerType.InBetweenAquiferCluster, soilProfile, out Point2D[] topLevelInBetweenAquifer, out Point2D[] bottomLevelInBetweenAquifer, 0); + SoilProfile2DHelper.DetermineInBetweenAquifersLayerBoundaryPoints(soilProfile, out List topLevelInBetweenAquifer, out List bottomLevelInBetweenAquifer); // Assert Assert.Multiple(() => { - Assert.That(isInBetweenLayerPresent, Is.True); - Assert.That(count, Is.EqualTo(1)); + Assert.That(topLevelInBetweenAquifer, Has.Count.EqualTo(1)); + Assert.That(bottomLevelInBetweenAquifer, Has.Count.EqualTo(1)); }); GeometryPoint[] expectedTop = [ new(leftCoordinate, -10), - new(middleXCoordinate, -10), + new(middleXCoordinate, -10), new(rightCoordinate, -10) ]; GeometryPoint[] expectedBottom = @@ -209,8 +211,8 @@ new(middleXCoordinate, -20), new(rightCoordinate, -20) ]; - AssertGeometry(topLevelInBetweenAquifer, expectedTop); - AssertGeometry(bottomLevelInBetweenAquifer, expectedBottom); + AssertGeometry(topLevelInBetweenAquifer[0].Points, expectedTop); + AssertGeometry(bottomLevelInBetweenAquifer[0].Points, expectedBottom); } /// @@ -246,17 +248,19 @@ soilProfile.Surfaces.Add(soilLayerAquiferInBetween); soilProfile.Surfaces.Add(soilLayerInBetween); soilProfile.Surfaces.Add(soilLayerAquiferBottom); - soilProfile.Geometry.Surfaces.Add(new GeometrySurface()); + soilProfile.Geometry.Surfaces.Add(soilLayer.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferInBetween.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerInBetween.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferBottom.GeometrySurface); // Call - bool isInBetweenLayerPresent = SoilProfile2DHelper.AreInBetweenAquiferClustersPresent(soilProfile, out int count); - SoilProfile2DHelper.DetermineAquiferLayerBoundaryPoints(SoilProfile2DHelper.LayerType.InBetweenAquiferCluster, soilProfile, out Point2D[] topLevelInBetweenAquifer, out Point2D[] bottomLevelInBetweenAquifer, 0); + SoilProfile2DHelper.DetermineInBetweenAquifersLayerBoundaryPoints(soilProfile, out List topLevelInBetweenAquifer, out List bottomLevelInBetweenAquifer); // Assert Assert.Multiple(() => { - Assert.That(isInBetweenLayerPresent, Is.True); - Assert.That(count, Is.EqualTo(1)); + Assert.That(topLevelInBetweenAquifer, Has.Count.EqualTo(1)); + Assert.That(bottomLevelInBetweenAquifer, Has.Count.EqualTo(1)); }); GeometryPoint[] expectedTop = [ @@ -268,8 +272,8 @@ new(leftCoordinate, -20), new(rightCoordinate, -20) ]; - AssertGeometry(topLevelInBetweenAquifer, expectedTop); - AssertGeometry(bottomLevelInBetweenAquifer, expectedBottom); + AssertGeometry(topLevelInBetweenAquifer[0].Points, expectedTop); + AssertGeometry(bottomLevelInBetweenAquifer[0].Points, expectedBottom); } /// @@ -338,34 +342,46 @@ soilProfile.Surfaces.Add(soilLayerAquiferInBetweenRight); soilProfile.Surfaces.Add(soilLayerInBetween); soilProfile.Surfaces.Add(soilLayerAquiferBottom); - soilProfile.Geometry.Surfaces.Add(new GeometrySurface()); + soilProfile.Geometry.Surfaces.Add(soilLayer.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferInBetweenLeft.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerGap.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferInBetweenRight.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerInBetween.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferBottom.GeometrySurface); // Call - bool isInBetweenLayerPresent = SoilProfile2DHelper.AreInBetweenAquiferClustersPresent(soilProfile, out int count); - SoilProfile2DHelper.DetermineAquiferLayerBoundaryPoints(SoilProfile2DHelper.LayerType.InBetweenAquiferCluster, soilProfile, out Point2D[] topLevelInBetweenAquifer, out Point2D[] bottomLevelInBetweenAquifer, 0); + SoilProfile2DHelper.DetermineInBetweenAquifersLayerBoundaryPoints(soilProfile, out List topLevelInBetweenAquifer, out List bottomLevelInBetweenAquifer); // Assert Assert.Multiple(() => { - Assert.That(isInBetweenLayerPresent, Is.True); - Assert.That(count, Is.EqualTo(1)); + Assert.That(topLevelInBetweenAquifer, Has.Count.EqualTo(2)); + Assert.That(bottomLevelInBetweenAquifer, Has.Count.EqualTo(2)); }); - GeometryPoint[] expectedTop = + GeometryPoint[] expectedTop1 = [ new(leftCoordinate, -10), - new(middleXCoordinate, -10), + new(middleXCoordinate, -10) + ]; + GeometryPoint[] expectedTop2 = + [ new(middleXCoordinate, -15), new(rightCoordinate, -15) ]; - GeometryPoint[] expectedBottom = + GeometryPoint[] expectedBottom1 = [ new(leftCoordinate, -20), - new(middleXCoordinate, -20), + new(middleXCoordinate, -20) + ]; + GeometryPoint[] expectedBottom2 = + [ new(middleXCoordinate, -21), new(rightCoordinate, -21) ]; - AssertGeometry(topLevelInBetweenAquifer, expectedTop); - AssertGeometry(bottomLevelInBetweenAquifer, expectedBottom); + AssertGeometry(topLevelInBetweenAquifer[0].Points, expectedTop1); + AssertGeometry(topLevelInBetweenAquifer[1].Points, expectedTop2); + AssertGeometry(bottomLevelInBetweenAquifer[0].Points, expectedBottom1); + AssertGeometry(bottomLevelInBetweenAquifer[1].Points, expectedBottom2); } /// @@ -430,34 +446,40 @@ soilProfile.Surfaces.Add(soilLayerAquiferInBetweenRight); soilProfile.Surfaces.Add(soilLayerInBetween); soilProfile.Surfaces.Add(soilLayerAquiferBottom); - soilProfile.Geometry.Surfaces.Add(new GeometrySurface()); + soilProfile.Geometry.Surfaces.Add(soilLayer.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferInBetweenLeft.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferInBetweenRight.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerInBetween.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferBottom.GeometrySurface); // Call - bool isInBetweenLayerPresent = SoilProfile2DHelper.AreInBetweenAquiferClustersPresent(soilProfile, out int count); - SoilProfile2DHelper.DetermineAquiferLayerBoundaryPoints(SoilProfile2DHelper.LayerType.InBetweenAquiferCluster, soilProfile, out Point2D[] topLevelInBetweenAquifer, out Point2D[] bottomLevelInBetweenAquifer, 0); + SoilProfile2DHelper.DetermineInBetweenAquifersLayerBoundaryPoints(soilProfile, out List topLevelInBetweenAquifer, out List bottomLevelInBetweenAquifer); // Assert Assert.Multiple(() => { - Assert.That(isInBetweenLayerPresent, Is.True); - Assert.That(count, Is.EqualTo(1)); + Assert.That(topLevelInBetweenAquifer, Has.Count.EqualTo(2)); + Assert.That(bottomLevelInBetweenAquifer, Has.Count.EqualTo(1)); }); - GeometryPoint[] expectedTop = + GeometryPoint[] expectedTop1 = [ new(leftCoordinate, -10), - new(middleXCoordinate, -10), + new(middleXCoordinate, -10) + ]; + GeometryPoint[] expectedTop2 = + [ new(middleXCoordinate, -5), new(rightCoordinate, -5) ]; GeometryPoint[] expectedBottom = [ new(leftCoordinate, -20), new(middleXCoordinate, -20), - new(middleXCoordinate, -20), new(rightCoordinate, -20) ]; - AssertGeometry(topLevelInBetweenAquifer, expectedTop); - AssertGeometry(bottomLevelInBetweenAquifer, expectedBottom); + AssertGeometry(topLevelInBetweenAquifer[0].Points, expectedTop1); + AssertGeometry(topLevelInBetweenAquifer[1].Points, expectedTop2); + AssertGeometry(bottomLevelInBetweenAquifer[0].Points, expectedBottom); } /// @@ -548,34 +570,45 @@ soilProfile.Surfaces.Add(soilLayerAquiferInBetweenRight); soilProfile.Surfaces.Add(soilLayerInBetween); soilProfile.Surfaces.Add(soilLayerAquiferBottom); - soilProfile.Geometry.Surfaces.Add(new GeometrySurface()); + soilProfile.Geometry.Surfaces.Add(soilLayer.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferInBetweenLeft.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferInBetweenRight.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerInBetween.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferBottom.GeometrySurface); // Call - bool isInBetweenLayerPresent = SoilProfile2DHelper.AreInBetweenAquiferClustersPresent(soilProfile, out int count); - SoilProfile2DHelper.DetermineAquiferLayerBoundaryPoints(SoilProfile2DHelper.LayerType.InBetweenAquiferCluster, soilProfile, out Point2D[] topLevelInBetweenAquifer, out Point2D[] bottomLevelInBetweenAquifer, 0); + SoilProfile2DHelper.DetermineInBetweenAquifersLayerBoundaryPoints(soilProfile, out List topLevelInBetweenAquifer, out List bottomLevelInBetweenAquifer); // Assert Assert.Multiple(() => { - Assert.That(isInBetweenLayerPresent, Is.True); - Assert.That(count, Is.EqualTo(1)); + Assert.That(topLevelInBetweenAquifer, Has.Count.EqualTo(2)); + Assert.That(bottomLevelInBetweenAquifer, Has.Count.EqualTo(2)); }); - GeometryPoint[] expectedTop = + GeometryPoint[] expectedTop1 = [ new(leftCoordinate, -10), - new(middleXCoordinate, -10), + new(middleXCoordinate, -10) + ]; + GeometryPoint[] expectedTop2 = + [ new(middleXCoordinate, -15), new(rightCoordinate, -15) ]; - GeometryPoint[] expectedBottom = + GeometryPoint[] expectedBottom1 = [ new(leftCoordinate, -20), - new(middleXCoordinate, -20), + new(middleXCoordinate, -20) + ]; + GeometryPoint[] expectedBottom2 = + [ new(middleXCoordinate, -17), new(rightCoordinate, -17) ]; - AssertGeometry(topLevelInBetweenAquifer, expectedTop); - AssertGeometry(bottomLevelInBetweenAquifer, expectedBottom); + AssertGeometry(topLevelInBetweenAquifer[0].Points, expectedTop1); + AssertGeometry(topLevelInBetweenAquifer[1].Points, expectedTop2); + AssertGeometry(bottomLevelInBetweenAquifer[0].Points, expectedBottom1); + AssertGeometry(bottomLevelInBetweenAquifer[1].Points, expectedBottom2); } /// @@ -666,84 +699,125 @@ soilProfile.Surfaces.Add(soilLayerAquiferInBetweenRight); soilProfile.Surfaces.Add(soilLayerInBetween); soilProfile.Surfaces.Add(soilLayerAquiferBottom); - soilProfile.Geometry.Surfaces.Add(new GeometrySurface()); + soilProfile.Geometry.Surfaces.Add(soilLayer.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferInBetweenLeft.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferInBetweenRight.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerInBetween.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferBottom.GeometrySurface); // Call - bool isInBetweenLayerPresent = SoilProfile2DHelper.AreInBetweenAquiferClustersPresent(soilProfile, out int count); - SoilProfile2DHelper.DetermineAquiferLayerBoundaryPoints(SoilProfile2DHelper.LayerType.InBetweenAquiferCluster, soilProfile, out Point2D[] topLevelInBetweenAquifer, out Point2D[] bottomLevelInBetweenAquifer, 0); + SoilProfile2DHelper.DetermineInBetweenAquifersLayerBoundaryPoints(soilProfile, out List topLevelInBetweenAquifer, out List bottomLevelInBetweenAquifer); // Assert Assert.Multiple(() => { - Assert.That(isInBetweenLayerPresent, Is.True); - Assert.That(count, Is.EqualTo(1)); + Assert.That(topLevelInBetweenAquifer, Has.Count.EqualTo(2)); + Assert.That(bottomLevelInBetweenAquifer, Has.Count.EqualTo(2)); }); - GeometryPoint[] expectedTop = + GeometryPoint[] expectedTop1 = [ new(leftCoordinate, -10), - new(middleXCoordinate, -10), + new(middleXCoordinate, -10) + ]; + GeometryPoint[] expectedTop2 = + [ new(middleXCoordinate, -5), new(rightCoordinate, -5) ]; - GeometryPoint[] expectedBottom = + GeometryPoint[] expectedBottom1 = [ new(leftCoordinate, -20), - new(middleXCoordinate, -20), + new(middleXCoordinate, -20) + ]; + GeometryPoint[] expectedBottom2 = + [ new(middleXCoordinate, -21), new(rightCoordinate, -21) ]; - AssertGeometry(topLevelInBetweenAquifer, expectedTop); - AssertGeometry(bottomLevelInBetweenAquifer, expectedBottom); + AssertGeometry(topLevelInBetweenAquifer[0].Points, expectedTop1); + AssertGeometry(topLevelInBetweenAquifer[1].Points, expectedTop2); + AssertGeometry(bottomLevelInBetweenAquifer[0].Points, expectedBottom1); + AssertGeometry(bottomLevelInBetweenAquifer[1].Points, expectedBottom2); } - + /// - /// --------------------------------------------------- Level 0 m - /// top layer - /// -------------------------|------------------------- Level -10 m - /// in-between layer left | in-between layer right - /// -------------------------|------------------------- Level -20 m - /// in-between layer - /// --------------------------------------------------- Level -25 m - /// bottom aquifer layer - /// --------------------------------------------------- Level -30 m + /// ----------------------------------------------\ E Level 0 m + /// D /-------------------\ F Level -2 m + /// Top aquitard / /------------------\ G Level -4 m + /// B / / C | + /// ---------------------/ / | Level -10 m + /// in-between aquifer / Aquitard | + /// ---------------------/ A | Level -13 m + /// ------------------------------------------------| Level -15 m + /// bottom aquifer | + /// ------------------------------------------------| Level -20 m /// [Test] - [TestCase(true, false)] - [TestCase(false, true)] - public void Given2DSoilProfileWithDiscontinuousInBetweenAquifer_WhenDeterminingInBetweenAquiferCount_ThenReturns0( - bool isLeftInBetweenLayerAquifer, bool isRightInBetweenLayerAquifer) + public void Given2DSoilProfileNotVerticalAtRightBoundaryAndWithInBetweenAquiferWithJump_WhenDeterminingLayerBoundaryPointsOfInBetweenAquifer_ThenExpectedCoordinatesReturned() { + var pointA = new Point2D(middleXCoordinate - 10, -13); + var pointB = new Point2D(middleXCoordinate - 10, -10); + var pointC = new Point2D(middleXCoordinate + 10, -4); + var pointD = new Point2D(middleXCoordinate + 10, -2); + var pointE = new Point2D(rightCoordinate - 2, 0); + var pointF = new Point2D(rightCoordinate - 1, -2); + var pointG = new Point2D(rightCoordinate, -4); + // Setup - SoilLayer2D soilLayer = CreateRectangularSoilLayer2D(0, -10, leftCoordinate, rightCoordinate, false); - SoilLayer2D soilLayerAquiferInBetweenLeft = CreateRectangularSoilLayer2D(-10, -20, leftCoordinate, middleXCoordinate, isLeftInBetweenLayerAquifer); - SoilLayer2D soilLayerAquiferInBetweenRight = CreateRectangularSoilLayer2D(-10, -20, middleXCoordinate, rightCoordinate, isRightInBetweenLayerAquifer); - SoilLayer2D soilLayerInBetween = CreateRectangularSoilLayer2D(-20, -25, leftCoordinate, rightCoordinate, false); - SoilLayer2D soilLayerAquiferBottom = CreateRectangularSoilLayer2D(-25, -30, leftCoordinate, rightCoordinate, true); + SoilLayer2D soilLayerTop = FactoryForSoilProfiles.CreateHexagonSoilLayer2D(new Point2D(leftCoordinate, 0), pointE, pointF, pointD, pointB, new Point2D(leftCoordinate, -10)); + SoilLayer2D soilLayerInBetweenAquifer = FactoryForSoilProfiles.CreatePolygoneSoilLayer2D([ + ..new[] + { + new Point2D(leftCoordinate, -10), pointB, pointD, pointF, pointG, pointC, pointA, new Point2D(leftCoordinate, -13) + } + ], null, true); + SoilLayer2D soilLayerAquitard = FactoryForSoilProfiles.CreateHexagonSoilLayer2D(new Point2D(leftCoordinate, -13), pointA, pointC, pointG, new Point2D(rightCoordinate, -15), new Point2D(leftCoordinate, -15)); + SoilLayer2D soilLayerAquiferBottom = CreateRectangularSoilLayer2D(-15, -20, leftCoordinate, rightCoordinate, true); var soilProfile = new SoilProfile2D { Geometry = new GeometryData { Left = leftCoordinate, Right = rightCoordinate, - Bottom = -30 + Bottom = -20 } }; - soilProfile.Surfaces.Add(soilLayer); - soilProfile.Surfaces.Add(soilLayerAquiferInBetweenLeft); - soilProfile.Surfaces.Add(soilLayerAquiferInBetweenRight); - soilProfile.Surfaces.Add(soilLayerInBetween); + soilProfile.Surfaces.Add(soilLayerTop); + soilProfile.Surfaces.Add(soilLayerInBetweenAquifer); + soilProfile.Surfaces.Add(soilLayerAquitard); soilProfile.Surfaces.Add(soilLayerAquiferBottom); + soilProfile.Geometry.Surfaces.Add(soilLayerTop.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerInBetweenAquifer.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquitard.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerAquiferBottom.GeometrySurface); // Call - bool isInBetweenLayerPresent = SoilProfile2DHelper.AreInBetweenAquiferClustersPresent(soilProfile, out int count); + SoilProfile2DHelper.DetermineInBetweenAquifersLayerBoundaryPoints(soilProfile, out List topLevelInBetweenAquifer, out List bottomLevelInBetweenAquifer); // Assert Assert.Multiple(() => { - Assert.That(isInBetweenLayerPresent, Is.False); - Assert.That(count, Is.EqualTo(0)); + Assert.That(topLevelInBetweenAquifer, Has.Count.EqualTo(1)); + Assert.That(bottomLevelInBetweenAquifer, Has.Count.EqualTo(1)); }); + GeometryPoint[] expectedTop = + [ + new(leftCoordinate, -10), + new(pointB.X, pointB.Z), + new(pointD.X, pointD.Z), + new(pointF.X, pointF.Z), + new(pointG.X, pointG.Z) + ]; + GeometryPoint[] expectedBottom = + [ + new(leftCoordinate, -13), + new(pointA.X, pointA.Z), + new(pointC.X, pointC.Z), + new(pointG.X, pointG.Z) + ]; + AssertGeometry(topLevelInBetweenAquifer[0].Points, expectedTop); + AssertGeometry(bottomLevelInBetweenAquifer[0].Points, expectedBottom); } /// Index: DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/Routines2D.cs =================================================================== diff -u -r5741 -r5923 --- DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/Routines2D.cs (.../Routines2D.cs) (revision 5741) +++ DamEngine/trunk/src/Deltares.DamEngine.Data/Geometry/Routines2D.cs (.../Routines2D.cs) (revision 5923) @@ -231,8 +231,8 @@ { // This is done by calculating and adding the angles from the point to all points // of the polygon (using the ArcTan2 function). If the sum of these angles is - // aproximate zero the point lies outside the polygon. - // If the sum in aproximate + or - 2Pi the point lies in the polygon; + // approximate zero the point lies outside the polygon. + // If the sum in approximate + or - 2Pi the point lies in the polygon; // Performance tests have proven that this algorithm performs better than the // polynomial winding number algorithm described at: // http://geomalgorithms.com/a03-_inclusion.html Index: DamEngine/trunk/src/Deltares.DamEngine.Calculators.Tests/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverterTests.cs =================================================================== diff -u -r5562 -r5923 --- DamEngine/trunk/src/Deltares.DamEngine.Calculators.Tests/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverterTests.cs (.../PlLinesToWaternetConverterTests.cs) (revision 5562) +++ DamEngine/trunk/src/Deltares.DamEngine.Calculators.Tests/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverterTests.cs (.../PlLinesToWaternetConverterTests.cs) (revision 5923) @@ -565,7 +565,22 @@ soilProfile.Surfaces.Add(aquifer7A); soilProfile.Surfaces.Add(aquifer7B); soilProfile.Surfaces.Add(aquifer7C); - soilProfile.Geometry.Surfaces.Add(new GeometrySurface()); + soilProfile.Geometry.Surfaces.Add(aquitard1.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquifer2A.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquifer2B.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquifer2C.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquitard3.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquifer4A.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquifer4B.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquifer4C.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquitard5.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquifer6A.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquifer6B.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquifer6C.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquifer7A.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquifer7B.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(aquifer7C.GeometrySurface); + soilProfile.Geometry.SurfaceLine.CalcPoints.Add(new Point2D(leftCoordinate, 20)); soilProfile.Geometry.SurfaceLine.CalcPoints.Add(new Point2D(rightCoordinate, 20)); soilProfile.Geometry.SurfaceLine.SyncPoints();