Index: DamEngine/trunk/src/Deltares.DamEngine.TestHelpers/Factories/FactoryForSoilProfiles.cs =================================================================== diff -u -r5241 -r5251 --- DamEngine/trunk/src/Deltares.DamEngine.TestHelpers/Factories/FactoryForSoilProfiles.cs (.../FactoryForSoilProfiles.cs) (revision 5241) +++ DamEngine/trunk/src/Deltares.DamEngine.TestHelpers/Factories/FactoryForSoilProfiles.cs (.../FactoryForSoilProfiles.cs) (revision 5251) @@ -1122,7 +1122,68 @@ return soilProfile2D; } - + + /// + /// -50 -20/-19.999/-19.998 -10 + /// |-----------------|-|-|---------------| Level 10 m + /// | surface 1 |2|3| surface 4 | + /// |-----------------|-|-|---------------| Level 0 m + /// | surface 5 |6| surface 7 | + /// | | | | + /// |-------------------|-|---------------| Level -15 m + /// +public static SoilProfile2D CreateSoilProfile2DWithFourSurfacesAndThreeVerticalSeparationDistantFrom1Mm(bool areConsecutiveSoilsIdentical) + { + var soilProfile2D = new SoilProfile2D + { + Geometry = new GeometryData + { + Left = -50, + Right = -10, + Bottom = -15 + } + }; + var soil1 = new Soil("Soil1", 11, 12); + Soil soil2 = areConsecutiveSoilsIdentical ? soil1 : new Soil("Soil2", 13, 14); + Soil soil3 = areConsecutiveSoilsIdentical ? soil1 : new Soil("Soil3", 13.1, 14.1); + Soil soil4 = areConsecutiveSoilsIdentical ? soil1 : new Soil("Soil4", 13.2, 14.2); + var soil5 = new Soil("Soil5", 15, 16); + Soil soil6 = areConsecutiveSoilsIdentical ? soil5 : new Soil("Soil6", 17, 18); + Soil soil7 = areConsecutiveSoilsIdentical ? soil5 : new Soil("Soil7", 17.1, 18.1); + + SoilLayer2D soilLayer1 = CreateRectangularSoilLayer2D(10, 0, -50, -20, soil1, soilProfile2D); + SoilLayer2D soilLayer2 = CreateRectangularSoilLayer2D(10, 0, -20, -19.999, soil2, soilProfile2D); + SoilLayer2D soilLayer3 = CreateRectangularSoilLayer2D(10, 0, -19.999, -19.998, soil3, soilProfile2D); + SoilLayer2D soilLayer4 = CreateRectangularSoilLayer2D(10, 0, -19.998, -10, soil4, soilProfile2D); + SoilLayer2D soilLayer5 = CreatePentagonSoilLayer2D(new Point2D(-50, 0), new Point2D(-20, 0), new Point2D(-19.999, 0), new Point2D(-19.999, -15), new Point2D(-50, -15), soil5, soilProfile2D); + SoilLayer2D soilLayer6 = CreateRectangularSoilLayer2D(0, -15, -19.999, -19.998, soil6, soilProfile2D); + SoilLayer2D soilLayer7 = CreateRectangularSoilLayer2D(0, -15, -19.998, -10, soil7, soilProfile2D); + + soilProfile2D.Surfaces.Add(soilLayer1); + soilProfile2D.Surfaces.Add(soilLayer2); + soilProfile2D.Surfaces.Add(soilLayer3); + soilProfile2D.Surfaces.Add(soilLayer4); + soilProfile2D.Surfaces.Add(soilLayer5); + soilProfile2D.Surfaces.Add(soilLayer6); + soilProfile2D.Surfaces.Add(soilLayer7); + + soilProfile2D.Geometry.Loops.Add(soilLayer1.GeometrySurface.OuterLoop); + soilProfile2D.Geometry.Surfaces.Add(soilLayer1.GeometrySurface); + soilProfile2D.Geometry.Loops.Add(soilLayer2.GeometrySurface.OuterLoop); + soilProfile2D.Geometry.Surfaces.Add(soilLayer2.GeometrySurface); + soilProfile2D.Geometry.Loops.Add(soilLayer3.GeometrySurface.OuterLoop); + soilProfile2D.Geometry.Surfaces.Add(soilLayer3.GeometrySurface); + soilProfile2D.Geometry.Loops.Add(soilLayer4.GeometrySurface.OuterLoop); + soilProfile2D.Geometry.Surfaces.Add(soilLayer4.GeometrySurface); + soilProfile2D.Geometry.Loops.Add(soilLayer5.GeometrySurface.OuterLoop); + soilProfile2D.Geometry.Surfaces.Add(soilLayer5.GeometrySurface); + soilProfile2D.Geometry.Loops.Add(soilLayer6.GeometrySurface.OuterLoop); + soilProfile2D.Geometry.Surfaces.Add(soilLayer6.GeometrySurface); + soilProfile2D.Geometry.Loops.Add(soilLayer7.GeometrySurface.OuterLoop); + soilProfile2D.Geometry.Surfaces.Add(soilLayer7.GeometrySurface); + + return soilProfile2D; +} private static string GetNewUniqueLayerId(SoilProfile1D soilProfile1D) { var num = 0; Index: DamEngine/trunk/src/Deltares.DamEngine.Data/Standard/ComputeDoubles.cs =================================================================== diff -u -r5044 -r5251 --- DamEngine/trunk/src/Deltares.DamEngine.Data/Standard/ComputeDoubles.cs (.../ComputeDoubles.cs) (revision 5044) +++ DamEngine/trunk/src/Deltares.DamEngine.Data/Standard/ComputeDoubles.cs (.../ComputeDoubles.cs) (revision 5251) @@ -101,7 +101,7 @@ return false; } - private static bool IsNearEqual(this double aValue1, double aValue2, double aTolerance) + public static bool IsNearEqual(this double aValue1, double aValue2, double aTolerance) { if (Math.Abs(aValue1 - aValue2) < aTolerance) { Index: DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/LayerDetector.cs =================================================================== diff -u -r4897 -r5251 --- DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/LayerDetector.cs (.../LayerDetector.cs) (revision 4897) +++ DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/LayerDetector.cs (.../LayerDetector.cs) (revision 5251) @@ -148,22 +148,21 @@ // finds intersections with regard to publisherEventArgs curve private void FindIntersections(Point2D aHeadpoint, Point2D aEndPoint, SoilLayer2D aSurfaceRefKey) { - var maxvalue = 999999999.0; - double minvalue = -999999999.0; + const double maxvalue = 999999999.0; + const double minvalue = -999999999.0; var verticalHeadPoint = new Point2D(verticalXCoordinate, maxvalue); var verticalEndPoint = new Point2D(verticalXCoordinate, minvalue); - Point2D intersectionPoint; LineIntersection intersectionResult = Routines2D.DetermineIf2DLinesIntersectStrickly(verticalHeadPoint, - verticalEndPoint, aHeadpoint, aEndPoint, out intersectionPoint); + verticalEndPoint, aHeadpoint, aEndPoint, out Point2D intersectionPoint); if (intersectionResult == LineIntersection.Intersects) { AddPointToList(new LayerIntersectionPoint(intersectionPoint, aSurfaceRefKey)); } else if (intersectionResult == LineIntersection.NoIntersection || intersectionResult == LineIntersection.Parallel) { - const double cEpsilon = 1.0e-3; + const double cEpsilon = 1.0e-7; if ((Routines2D.DoesPointExistInLine(verticalHeadPoint, verticalEndPoint, aHeadpoint, cEpsilon)) && (Routines2D.DoesPointExistInLine(verticalHeadPoint, verticalEndPoint, aEndPoint, cEpsilon))) Index: DamEngine/trunk/src/Deltares.DamEngine.Data.Tests/Geotechnics/SoilProfile2DTests.cs =================================================================== diff -u -r5205 -r5251 --- DamEngine/trunk/src/Deltares.DamEngine.Data.Tests/Geotechnics/SoilProfile2DTests.cs (.../SoilProfile2DTests.cs) (revision 5205) +++ DamEngine/trunk/src/Deltares.DamEngine.Data.Tests/Geotechnics/SoilProfile2DTests.cs (.../SoilProfile2DTests.cs) (revision 5251) @@ -19,7 +19,6 @@ // Stichting Deltares and remain full property of Stichting Deltares at all times. // All rights reserved. -using System.Collections.Generic; using System.Linq; using Deltares.DamEngine.Data.Geotechnics; using Deltares.DamEngine.TestHelpers.Factories; @@ -31,6 +30,8 @@ [TestFixture] public class SoilProfile2DTests { + private double precison10Decimals = 0.000000000051; + [Test] public void GivenProfile2DWhenCloningThenClonedProfile2DEqualsOriginal() { @@ -54,7 +55,44 @@ // Then CompareSoilProfile2D(soilProfile2D, cloneSoilProfile2D); } + + + [Test] + [TestCase( true, -20.000, -20.000)] + [TestCase( false, -20.000, -19.997)] + [TestCase( true, -50.000, -49.999)] + [TestCase( true, -10.000, -10.001)] + public void GivenProfile2DWithVerticalSeparationLayerWhenDeterminingSoilProfile1DAtSeparationThenReturnsExpectedSoilProfile1D(bool areConsecutiveSoilsIdentical, double givenX, double expectedX) + { + // Given + SoilProfile2D soilProfile2D = FactoryForSoilProfiles.CreateSoilProfile2DWithFourSurfacesAndThreeVerticalSeparationDistantFrom1Mm(areConsecutiveSoilsIdentical); + soilProfile2D.Name = "Soil Profile 2D Test"; + + // When + SoilProfile1D soilProfile1D = soilProfile2D.GetSoilProfile1D(givenX); + // Then + Assert.That(soilProfile1D.Layers, Has.Count.EqualTo(2)); + Assert.Multiple(() => + { + Assert.That(soilProfile1D.Layers[0].TopLevel, Is.EqualTo(10).Within(precison10Decimals)); + Assert.That(soilProfile1D.Layers[1].TopLevel, Is.EqualTo(0).Within(precison10Decimals)); + Assert.That(soilProfile1D.BottomLevel, Is.EqualTo(-15).Within(precison10Decimals)); + Assert.That(soilProfile1D.Name, Is.EqualTo("Generated at x = " + expectedX.ToString("F3") + " from Soil Profile 2D Test")); + + if (areConsecutiveSoilsIdentical) + { + Assert.That(soilProfile1D.Layers[0].Soil.Name, Is.EqualTo("Soil1")); + Assert.That(soilProfile1D.Layers[1].Soil.Name, Is.EqualTo("Soil5")); + } + else + { + Assert.That(soilProfile1D.Layers[0].Soil.Name, Is.EqualTo("Soil4")); + Assert.That(soilProfile1D.Layers[1].Soil.Name, Is.EqualTo("Soil7")); + } + }); + } + private static void CompareSoilProfile2D(SoilProfile2D sourceSoilProfile2D, SoilProfile2D targetSoilProfile2D) { Assert.That(targetSoilProfile2D, Is.Not.SameAs(sourceSoilProfile2D)); Index: DamEngine/trunk/src/Deltares.DamEngine.Calculators/PlLinesCreator/PlLinesCreator.cs =================================================================== diff -u -r5250 -r5251 --- DamEngine/trunk/src/Deltares.DamEngine.Calculators/PlLinesCreator/PlLinesCreator.cs (.../PlLinesCreator.cs) (revision 5250) +++ DamEngine/trunk/src/Deltares.DamEngine.Calculators/PlLinesCreator/PlLinesCreator.cs (.../PlLinesCreator.cs) (revision 5251) @@ -897,26 +897,22 @@ /// The surface point. private void ConfigureUpliftCalculator(UpliftCalculator upliftCalculator, GeometryPoint surfacePoint) { - // Offset to solve issue MWDAM-764 - const double offset = 0.01; upliftCalculator.SurfaceLevel = surfacePoint.Z; if (currentPl1Line != null) { upliftCalculator.PhreaticLevel = currentPl1Line.ZFromX(surfacePoint.X); - // set phreatic level to calculate upliftfactor + // set phreatic level to calculate uplift factor } else { upliftCalculator.PhreaticLevel = WaterLevelPolder; // if not PL1 created then assume phreatic level is same as waterlevel polder } - SoilProfile1D actualSoilProfile = GetSoilProfileBelowPoint(surfacePoint.X + offset); - // At the end of a geometry, there are no layers to be found beyond that end. In that case - // get the layers just before the end of the geometry. + SoilProfile1D actualSoilProfile = GetSoilProfileBelowPoint(surfacePoint.X); if (actualSoilProfile.LayerCount == 0) { - actualSoilProfile = GetSoilProfileBelowPoint(surfacePoint.X - offset); + actualSoilProfile = GetSoilProfileBelowPoint(surfacePoint.X); } AdaptSoilProfileForSurfaceLevelAccuracy(actualSoilProfile, upliftCalculator.SurfaceLevel); @@ -1782,20 +1778,17 @@ private SoilProfile1D GetRelevantSoilProfileForAquiferLayersSearch() { - // Offset to solve issue MWDAM-557 - const double offset = 0.01; - GeometryPoint relevantPoint = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder); if (relevantPoint == null) { throw new PlLinesCreatorException("Characteristic point \"dike toe at polder\" is not defined."); } - return GetSoilProfileBelowPoint(relevantPoint.X + offset); + return GetSoilProfileBelowPoint(relevantPoint.X); } /// - /// Throws the when water level above dike. + /// Throws the exception when water level above dike. /// /// The water level river. /// The surface line. Index: DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/SoilProfile2D.cs =================================================================== diff -u -r5234 -r5251 --- DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/SoilProfile2D.cs (.../SoilProfile2D.cs) (revision 5234) +++ DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/SoilProfile2D.cs (.../SoilProfile2D.cs) (revision 5251) @@ -21,8 +21,9 @@ using System; using System.Collections.Generic; -using System.Runtime.CompilerServices; +using System.Linq; using Deltares.DamEngine.Data.Geometry; +using Deltares.DamEngine.Data.Standard; using Deltares.DamEngine.Data.Standard.Language; using Deltares.DamEngine.Data.Standard.Validation; @@ -35,6 +36,7 @@ { private readonly List surfaces = new List(); private const double deviation = 1E-05; + private const double toleranceAlmostEqual = 1e-07; /// /// Initializes a new instance of the class. @@ -80,38 +82,34 @@ Geometry.Right = Geometry.MaxGeometryPointsX; Geometry.Left = Geometry.MinGeometryPointsX; } - - var soilProfile = new SoilProfile1D + + // At the end of a geometry, there are no layers to be found beyond that end. In that case get the layers just before + // the end of the geometry. + if (x.IsGreaterThanOrEqualTo(Geometry.Right, toleranceAlmostEqual)) { - Name = "Generated at x = " + x + " from " + Name - }; - var detector = new LayerDetector(Surfaces); - if (x > Geometry.Right) - { x = Geometry.Right - diff; } - if (x < Geometry.Left) + if (Geometry.Left.IsGreaterThanOrEqualTo(x, toleranceAlmostEqual)) { x = Geometry.Left + diff; } - detector.DetermineMaterials(x); - - if (detector.LayerList.Count > 0) + // If the X coincides with a vertical layer separation and both profiles at left and at right of the separation are different, + // then the X is shifted by 1 mm to the right. + bool isXShiftedToRight; + double xRelevant = x; + do { - soilProfile.BottomLevel = detector.LayerList[detector.LayerList.Count - 1].EndPoint.Z; - for (var i = 0; i < detector.LayerList.Count; i++) + isXShiftedToRight = HasLayerAVerticalPartAtGivenX(xRelevant) && !AreSoilProfilesIdentical(xRelevant - diff / 2, xRelevant + diff / 2); + if (isXShiftedToRight) { - var layer = new SoilLayer1D(detector.LayerList[i].Soil, detector.LayerList[i].StartPoint.Z) - { - IsAquifer = detector.LayerList[i].IsAquifer - }; - soilProfile.Layers.Add(layer); + xRelevant += diff; } - } - return soilProfile; + } while (isXShiftedToRight && xRelevant < Geometry.Right); + + return DetermineSoilProfile1DAtX(xRelevant); } /// @@ -372,4 +370,73 @@ } return isPointInSurface; } + + private SoilProfile1D DetermineSoilProfile1DAtX(double x) + { + var soilProfile = new SoilProfile1D + { + Name = "Generated at x = " + x.ToString("F3") + " from " + Name + }; + + var detector = new LayerDetector(Surfaces); + detector.DetermineMaterials(x); + + if (detector.LayerList.Count > 0) + { + soilProfile.BottomLevel = detector.LayerList[detector.LayerList.Count - 1].EndPoint.Z; + for (var i = 0; i < detector.LayerList.Count; i++) + { + var layer = new SoilLayer1D(detector.LayerList[i].Soil, detector.LayerList[i].StartPoint.Z) + { + IsAquifer = detector.LayerList[i].IsAquifer, + WaterpressureInterpolationModel = detector.LayerList[i].WaterpressureInterpolationModel, + Name = detector.LayerList[i].Name, + SoilName = detector.LayerList[i].SoilName + }; + soilProfile.Layers.Add(layer); + } + } + + return soilProfile; + } + + private bool HasLayerAVerticalPartAtGivenX(double xCoordinate) + { + foreach (SoilLayer2D surface in Surfaces) + { + foreach (GeometryCurve curve in surface.GeometrySurface.OuterLoop.CurveList) + { + if (curve.HeadPoint.X.IsNearEqual(xCoordinate, toleranceAlmostEqual) && + curve.EndPoint.X.IsNearEqual(xCoordinate, toleranceAlmostEqual)) + { + return true; + } + } + } + + return false; + } + + private bool AreSoilProfilesIdentical(double x1, double x2) + { + SoilProfile1D soilProfile1 = DetermineSoilProfile1DAtX(x1); + SoilProfile1D soilProfile2 = DetermineSoilProfile1DAtX(x2); + + bool isIdentical = soilProfile1.Layers.Count == soilProfile2.Layers.Count; + if (!isIdentical) + { + return false; + } + + for (var i = 0; i < soilProfile1.Layers.Count; i++) + { + isIdentical = isIdentical && soilProfile1.Layers[i].TopLevel.IsNearEqual(soilProfile2.Layers[i].TopLevel, toleranceAlmostEqual); + isIdentical = isIdentical && soilProfile1.Layers[i].Soil == soilProfile2.Layers[i].Soil; + isIdentical = isIdentical && soilProfile1.Layers[i].SoilName == soilProfile2.Layers[i].SoilName; + isIdentical = isIdentical && soilProfile1.Layers[i].IsAquifer == soilProfile2.Layers[i].IsAquifer; + isIdentical = isIdentical && soilProfile1.Layers[i].WaterpressureInterpolationModel == soilProfile2.Layers[i].WaterpressureInterpolationModel; + } + + return isIdentical; + } } \ No newline at end of file