Index: DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/Common/SoilProfile2DHelper.cs =================================================================== diff -u -r5923 -r5928 --- DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/Common/SoilProfile2DHelper.cs (.../SoilProfile2DHelper.cs) (revision 5923) +++ DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/Common/SoilProfile2DHelper.cs (.../SoilProfile2DHelper.cs) (revision 5928) @@ -32,6 +32,7 @@ public enum LayerType { BottomAquiferCluster, + InBetweenAquiferCluster, HighestAquiferCluster, LowestLayer } @@ -158,34 +159,19 @@ /// All the points of the layer boundary. 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)) - { - xCoordinatesAll.Add(xCoordinate - deviationX); - xCoordinatesAll.Add(xCoordinate + deviationX); - } - else - { - xCoordinatesAll.Add(xCoordinate); - } - } - + double[] xCoordinates = DetermineXCoordinatesOfSoilProfileIncludingVerticalParts(layerType, soilProfile); var layerBoundaryTopPoints = new List(); var layerBoundaryBottomPoints = new List(); - foreach (double t in xCoordinatesAll) + foreach (double xCoordinate in xCoordinates) { - SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(t); - double currentAquiferTop = GetLevel(layerType, BoundaryType.Top, crossSection); - double currentAquiferBottom = GetLevel(layerType, BoundaryType.Bottom, crossSection); + SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(xCoordinate); + double currentAquiferTop = GetLevel(layerType, BoundaryType.Top, crossSection, -1); + double currentAquiferBottom = GetLevel(layerType, BoundaryType.Bottom, crossSection, -1); - if (!double.IsNaN(currentAquiferTop) && !double.IsNaN(currentAquiferBottom)) { - layerBoundaryTopPoints.Add(new Point2D(t, currentAquiferTop)); - layerBoundaryBottomPoints.Add(new Point2D(t, currentAquiferBottom)); + layerBoundaryTopPoints.Add(new Point2D(xCoordinate, currentAquiferTop)); + layerBoundaryBottomPoints.Add(new Point2D(xCoordinate, currentAquiferBottom)); } } @@ -238,18 +224,21 @@ { // The inner loops are removed because they can't be in-between aquifers SoilProfile2D soilProfileWithoutInnerLoops = CreateSoilProfileWithoutInnerLoops(soilProfile); - double[] xCoordinates = DetermineAllXCoordinatesOfSoilProfile(soilProfileWithoutInnerLoops); + double[] xCoordinates = DetermineXCoordinatesOfSoilProfileIncludingVerticalParts(LayerType.InBetweenAquiferCluster, soilProfileWithoutInnerLoops); - var relevantAquifers = new List(); + var relevantTopAquifers = new List(); + var relevantBottomAquifers = 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); + double zLevelTop = crossSection.GetInBetweenAquiferClusters[i].Item1.TopLevel - GeometryConstants.Accuracy; + double zLevelBottom = crossSection.GetInBetweenAquiferClusters[i].Item2.BottomLevel + GeometryConstants.Accuracy; + AddAquiferContainingPointIfNotYetAdded(soilProfileWithoutInnerLoops, xCoordinate, zLevelTop, ref relevantTopAquifers); + AddAquiferContainingPointIfNotYetAdded(soilProfileWithoutInnerLoops, xCoordinate, zLevelBottom, ref relevantBottomAquifers); } } } @@ -258,57 +247,85 @@ List isolatedLayerBoundaryBottom = []; aquifersBoundaryTop = []; aquifersBoundaryBottom = []; - if (relevantAquifers != null) + if (relevantTopAquifers != null && relevantBottomAquifers != null) { - foreach (GeometrySurface surface in relevantAquifers) + isolatedLayerBoundaryTop.AddRange(relevantTopAquifers.Select(surface => surface.DetermineTopGeometrySurface())); + isolatedLayerBoundaryBottom.AddRange(relevantBottomAquifers.Select(surface => surface.DetermineBottomGeometrySurface())); + aquifersBoundaryTop = ConnectPolyLines(isolatedLayerBoundaryTop); + aquifersBoundaryBottom = ConnectPolyLines(isolatedLayerBoundaryBottom); + } + } + + private static double[] DetermineXCoordinatesOfSoilProfileIncludingVerticalParts(LayerType layerType, SoilProfile2D soilProfile) + { + double[] xCoordinates = DetermineAllXCoordinatesOfSoilProfile(soilProfile); + var xCoordinatesAll = new List(); + foreach (double xCoordinate in xCoordinates) + { + SoilProfile1D crossSection = soilProfile.GetSoilProfile1D(xCoordinate); + int count = layerType == LayerType.InBetweenAquiferCluster ? crossSection.GetInBetweenAquiferClusters.Count : 1; { - isolatedLayerBoundaryTop.Add(surface.DetermineTopGeometrySurface()); - isolatedLayerBoundaryBottom.Add(surface.DetermineBottomGeometrySurface()); + for (var i = 0; i < count; i++) + { + if (HasAquiferAVerticalPartAtGivenX(layerType, xCoordinate, soilProfile, i)) + { + xCoordinatesAll.Add(xCoordinate - deviationX); + xCoordinatesAll.Add(xCoordinate + deviationX); + } + else + { + xCoordinatesAll.Add(xCoordinate); + } + } } - - aquifersBoundaryTop = ConnectPolylines(isolatedLayerBoundaryTop); - aquifersBoundaryBottom = ConnectPolylines(isolatedLayerBoundaryBottom); } + return xCoordinatesAll.ToArray(); } - private static List ConnectPolylines(List polylines) + /// + /// Connect the poly-lines in the list when possible (i.e. when the end point of a poly-line coincides with the start point + /// of another poly-line. + /// + /// + /// A list of connected poly-lines. + internal static List ConnectPolyLines(List polyLines) { - var connectedPolylines = new List(); - var usedPolylines = new HashSet(); + var connectedPolyLines = new List(); + var usedPolyLines = new HashSet(); - foreach (GeometryPointString polyline in polylines) + foreach (GeometryPointString polyLine in polyLines) { - if (usedPolylines.Contains(polyline)) + if (usedPolyLines.Contains(polyLine)) continue; - var currentPolyline = new GeometryPointString(); - currentPolyline.CalcPoints.AddRange(polyline.CalcPoints); - usedPolylines.Add(polyline); + var currentPolyLine = new GeometryPointString(); + currentPolyLine.CalcPoints.AddRange(polyLine.CalcPoints); + usedPolyLines.Add(polyLine); - bool foundConnection; + bool isConnectionFound; do { - foundConnection = false; - foreach (GeometryPointString nextPolyline in polylines) + isConnectionFound = false; + foreach (GeometryPointString nextPolyLine in polyLines) { - if (usedPolylines.Contains(nextPolyline)) + if (usedPolyLines.Contains(nextPolyLine)) continue; - if (currentPolyline.CalcPoints.Last().LocationEquals(nextPolyline.CalcPoints.First())) + if (currentPolyLine.CalcPoints.Last().LocationEquals(nextPolyLine.CalcPoints.First())) { - currentPolyline.CalcPoints.AddRange(nextPolyline.CalcPoints.Skip(1)); - usedPolylines.Add(nextPolyline); - foundConnection = true; + currentPolyLine.CalcPoints.AddRange(nextPolyLine.CalcPoints.Skip(1)); + usedPolyLines.Add(nextPolyLine); + isConnectionFound = true; break; } } - } while (foundConnection); + } while (isConnectionFound); - currentPolyline.SyncPoints(); - connectedPolylines.Add(currentPolyline); + currentPolyLine.SyncPoints(); + connectedPolyLines.Add(currentPolyLine); } - return connectedPolylines; + return connectedPolyLines; } /// @@ -330,9 +347,9 @@ && Math.Abs(point2Ds.Last().X - rightGeometryBoundary) < toleranceAlmostEqual; } - private static double GetLevel(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D) + private static double GetLevel(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer) { - SoilLayer1D layer = GetLayer(layerType, boundaryType, soilProfile1D); + SoilLayer1D layer = GetLayer(layerType, boundaryType, soilProfile1D, indexInBetweenAquifer); if (layer != null) { return boundaryType == BoundaryType.Top ? layer.TopLevel : layer.BottomLevel; @@ -341,12 +358,14 @@ return double.NaN; } - private static SoilLayer1D GetLayer(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D) + private static SoilLayer1D GetLayer(LayerType layerType, BoundaryType boundaryType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer) { 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: @@ -356,7 +375,7 @@ return null; } - private static bool HasAquiferAVerticalPartAtGivenX(LayerType layerType, double xCoordinate, SoilProfile2D soilProfile) + private static bool HasAquiferAVerticalPartAtGivenX(LayerType layerType, double xCoordinate, SoilProfile2D soilProfile, int indexInBetweenAquifer) { SoilProfile1D crossSectionLeft = soilProfile.GetSoilProfile1D(xCoordinate - deviationX); SoilProfile1D crossSectionRight = soilProfile.GetSoilProfile1D(xCoordinate + deviationX); @@ -365,10 +384,10 @@ return false; } - 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); + 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; Index: DamEngine/trunk/src/Deltares.DamEngine.Calculators.Tests/KernelWrappers/Common/SoilProfile2DHelperTests.cs =================================================================== diff -u -r5923 -r5928 --- DamEngine/trunk/src/Deltares.DamEngine.Calculators.Tests/KernelWrappers/Common/SoilProfile2DHelperTests.cs (.../SoilProfile2DHelperTests.cs) (revision 5923) +++ DamEngine/trunk/src/Deltares.DamEngine.Calculators.Tests/KernelWrappers/Common/SoilProfile2DHelperTests.cs (.../SoilProfile2DHelperTests.cs) (revision 5928) @@ -19,7 +19,6 @@ // Stichting Deltares and remain full property of Stichting Deltares at all times. // All rights reserved. -using System; using System.Collections.Generic; using System.Linq; using Deltares.DamEngine.Calculators.KernelWrappers.Common; @@ -151,23 +150,26 @@ /// /// --------------------------------------------------------------- Level 0 m - /// top layer + /// top aquitard /// --------------------------------------------------------------- Level -10 m /// in-between aquifer layer left | in-between aquifer layer right /// --------------------------------------------------------------- Level -20 m - /// in-between layer + /// in-between aquifer layer + /// --------------------------------------------------------------- Level -22 m + /// in-between aquitard /// --------------------------------------------------------------- Level -25 m /// bottom aquifer layer /// --------------------------------------------------------------- Level -30 m /// [Test] - public void Given2DSoilProfileWithContinuousInBetweenAquiferLayerConsistingOfTwoParts_WhenDeterminingLayerBoundaryPointsOfInBetweenAquifer_ThenExpectedCoordinatesReturned() + public void Given2DSoilProfileWithContinuousInBetweenAquiferClusterConsistingOfThreeParts_WhenDeterminingLayerBoundaryPointsOfInBetweenAquifer_ThenExpectedCoordinatesReturned() { // Setup SoilLayer2D soilLayer = CreateRectangularSoilLayer2D(0, -10, leftCoordinate, rightCoordinate, false); SoilLayer2D soilLayerAquiferInBetweenLeft = CreateRectangularSoilLayer2D(-10, -20, leftCoordinate, middleXCoordinate, true); SoilLayer2D soilLayerAquiferInBetweenRight = CreateRectangularSoilLayer2D(-10, -20, middleXCoordinate, rightCoordinate, true); - SoilLayer2D soilLayerInBetween = CreateRectangularSoilLayer2D(-20, -25, leftCoordinate, rightCoordinate, false); + SoilLayer2D soilLayerInBetweenAquifer = CreateRectangularSoilLayer2D(-20, -22, leftCoordinate, rightCoordinate, true); + SoilLayer2D soilLayerInBetweenAquitard = CreateRectangularSoilLayer2D(-22, -25, leftCoordinate, rightCoordinate, false); SoilLayer2D soilLayerAquiferBottom = CreateRectangularSoilLayer2D(-25, -30, leftCoordinate, rightCoordinate, true); var soilProfile = new SoilProfile2D @@ -182,12 +184,14 @@ soilProfile.Surfaces.Add(soilLayer); soilProfile.Surfaces.Add(soilLayerAquiferInBetweenLeft); soilProfile.Surfaces.Add(soilLayerAquiferInBetweenRight); - soilProfile.Surfaces.Add(soilLayerInBetween); + soilProfile.Surfaces.Add(soilLayerInBetweenAquifer); + soilProfile.Surfaces.Add(soilLayerInBetweenAquitard); soilProfile.Surfaces.Add(soilLayerAquiferBottom); 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(soilLayerInBetweenAquifer.GeometrySurface); + soilProfile.Geometry.Surfaces.Add(soilLayerInBetweenAquitard.GeometrySurface); soilProfile.Geometry.Surfaces.Add(soilLayerAquiferBottom.GeometrySurface); // Call @@ -202,14 +206,13 @@ GeometryPoint[] expectedTop = [ new(leftCoordinate, -10), - new(middleXCoordinate, -10), + new(middleXCoordinate, -10), new(rightCoordinate, -10) ]; GeometryPoint[] expectedBottom = [ - new(leftCoordinate, -20), - new(middleXCoordinate, -20), - new(rightCoordinate, -20) + new(leftCoordinate, -22), + new(rightCoordinate, -22) ]; AssertGeometry(topLevelInBetweenAquifer[0].Points, expectedTop); AssertGeometry(bottomLevelInBetweenAquifer[0].Points, expectedBottom); @@ -279,7 +282,7 @@ /// /// ------------------------------------------------------------- Level 0 m /// top layer - /// ------------------------------------------------------------- Level -10 m + /// -----------------------------|------------------------------- Level -10 m /// | Gap layer /// in-between aquifer |------------------------------- Level -15 m /// layer left | @@ -771,8 +774,8 @@ 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); + 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 { @@ -1010,6 +1013,54 @@ }).SetName("Right aquifer fully enveloped by left aquifer"); } + [Test] + public void GivenAListOfPolyLines_WhenConnectingThem_ThenExpectedConnectedPolyLinesReturned() + { + // Setup + var polyLine1 = new GeometryPointString(); + polyLine1.CalcPoints.Add(new Point2D(-1, 2)); + polyLine1.CalcPoints.Add(new Point2D(3, 4)); + polyLine1.CalcPoints.Add(new Point2D(5, 6)); + + var polyLine2 = new GeometryPointString(); + polyLine2.CalcPoints.Add(new Point2D(5, 6)); + polyLine2.CalcPoints.Add(new Point2D(9, 1)); + polyLine2.CalcPoints.Add(new Point2D(15, 26)); + + var polyLine3 = new GeometryPointString(); + polyLine3.CalcPoints.Add(new Point2D(3, 4)); + polyLine3.CalcPoints.Add(new Point2D(9, 5)); + polyLine3.CalcPoints.Add(new Point2D(6, 8)); + + var polyLines = new List + { + polyLine1, + polyLine2, + polyLine3 + }; + + // Call + List connectedPolyLines = SoilProfile2DHelper.ConnectPolyLines(polyLines); + + // Assert - Expected is that poly-lines 1 and 2 are connected but poly-line 3 not + Assert.Multiple(() => + { + Assert.That(connectedPolyLines, Has.Count.EqualTo(2)); + Assert.That(connectedPolyLines[0].CalcPoints, Has.Count.EqualTo(5)); + Assert.That(connectedPolyLines[1].CalcPoints, Has.Count.EqualTo(3)); + Assert.That(connectedPolyLines[0].CalcPoints[0].X, Is.EqualTo(-1)); + Assert.That(connectedPolyLines[0].CalcPoints[0].Z, Is.EqualTo(2)); + Assert.That(connectedPolyLines[0].CalcPoints[1].X, Is.EqualTo(3)); + Assert.That(connectedPolyLines[0].CalcPoints[1].Z, Is.EqualTo(4)); + Assert.That(connectedPolyLines[0].CalcPoints[2].X, Is.EqualTo(5)); + Assert.That(connectedPolyLines[0].CalcPoints[2].Z, Is.EqualTo(6)); + Assert.That(connectedPolyLines[0].CalcPoints[3].X, Is.EqualTo(9)); + Assert.That(connectedPolyLines[0].CalcPoints[3].Z, Is.EqualTo(1)); + Assert.That(connectedPolyLines[0].CalcPoints[4].X, Is.EqualTo(15)); + Assert.That(connectedPolyLines[0].CalcPoints[4].Z, Is.EqualTo(26)); + }); + } + private static SoilLayer2D CreateSoilLayer2D(Point2D topLeftCoordinate, Point2D topRightCoordinate, Point2D bottomRightCoordinate, Point2D bottomLeftCoordinate, bool isAquifer) {